求助：使用Matlab的linprog函数求解数学优化模型

问题概述 #

需调用Matlab的linprog(f,A,b,Aeq,beq,lb,ub)函数求解污泥运输与处理中心选址的线性规划模型，需将目标函数与约束条件转换为该函数可识别的格式。

参数定义 #

n = 3; 
m = 4; 
J = 5;
Pi = [1000; 1500; 2000]; % 污水处理厂i的污泥产量
Lk = [2000; 2500; 3000; 3500]; % 处理中心k的最大处理能力
Cik = [3 4 5 6; 2 3 4 5; 4 5 6 7]; % 从污水处理厂i到处理中心k的单位运输成本
Ckj = [5 6 7 8; 4 5 6 7; 6 7 8 9; 3 4 5 6; 2 3 4 5]; % 从处理中心k到农场j的单位运输成本
CAk = [1000; 1500; 2000; 2500]; % 处理中心k的年运营成本
CIk = [5000; 6000; 7000; 8000]; % 处理中心k的建设成本
Nj = [1500; 2000; 2500; 3000; 3500]; % 农场j的生物肥料需求量
alpha = 0.8; % 污泥转化为生物肥料的转化率

变量定义 #

• X_k：0-1决策变量，选择位置k建立处理中心则为1，否则为0
• Y_ik：连续变量，从污水处理厂i运往处理中心k的污泥量
• Z_kj：连续变量，从处理中心k运往农场j的生物肥料量

目标函数 #

原目标为总运输成本、处理中心建设与运营成本之和：
$$
\min \sum_{i=1}^n \sum_{k=1}^m C_{ik}Y_{ik}X_k + \sum_{k=1}^m (CI_k + CA_k)X_k + \sum_{k=1}^m \sum_{j=1}^J C_{kj}Z_{kj}X_k
$$
需转换为f'*x的线性形式（x为所有变量的向量）。

约束条件 #

1. 污泥产量约束：每个污水处理厂的污泥全部运出
   $$
   \sum_{k=1}^m Y_{ik} = P_i \quad (i=1,2,...,n)
   $$
2. 处理能力约束：处理中心k接收的污泥量不超过其最大处理能力（仅当X_k=1时生效）
   $$
   \sum_{i=1}^n Y_{ik} \leq L_k X_k \quad (k=1,2,...,m)
   $$
3. 肥料转化约束：处理中心k产出的肥料量等于接收污泥量乘以转化率（仅当X_k=1时生效）
   $$
   \sum_{j=1}^J Z_{kj} = \alpha \sum_{i=1}^n Y_{ik} \quad (k=1,2,...,m)
   $$
4. 肥料需求约束：农场j接收的肥料量不超过其需求量
   $$
   \sum_{k=1}^m Z_{kj} \leq N_j \quad (j=1,2,...,J)
   $$
5. 非负约束：
   $$
   Y_{ik} \geq 0, \quad Z_{kj} \geq 0, \quad X_k \in {0,1}
   $$

转换为linprog格式的实现步骤 #

linprog仅支持连续变量，因此先将0-1变量X_k的约束设为$0 \leq X_k \leq 1$，后续可对结果取整；若需严格0-1解，建议使用intlinprog。

1. 变量向量化 #

将所有变量按顺序排成向量x：

• 前m个元素：$X_1,X_2,...,X_m$
• 接下来n*m个元素：$Y_{11},Y_{12},...,Y_{1m},Y_{21},...,Y_{nm}$
• 最后mJ个元素：$Z_{11},Z_{12},...,Z_{1J},Z_{21},...,Z_{mJ}$
  总变量数：$m + nm + m*J = 36$

2. 构造目标函数系数向量f #

total_vars = m + n*m + m*J;
f = zeros(total_vars, 1);

% X_k的系数：建设+运营成本
f(1:m) = CIk + CAk;

% Y_ik的系数：单位运输成本
y_start = m + 1;
for i = 1:n
    for k = 1:m
        f(y_start + (i-1)*m + (k-1)) = Cik(i,k);
    end
end

% Z_kj的系数：单位运输成本
z_start = m + n*m + 1;
for k = 1:m
    for j = 1:J
        f(z_start + (k-1)*J + (j-1)) = Ckj(j,k);
    end
end

3. 构造约束矩阵 #

等式约束（Aeq, beq）

Aeq = [];
beq = [];

% 污泥产量约束
for i = 1:n
    row = zeros(1, total_vars);
    y_indices = m + (i-1)*m + 1 : m + i*m;
    row(y_indices) = 1;
    Aeq = [Aeq; row];
    beq = [beq; Pi(i)];
end

% 肥料转化约束
for k = 1:m
    row = zeros(1, total_vars);
    z_indices = m + n*m + (k-1)*J + 1 : m + n*m + k*J;
    row(z_indices) = 1;
    y_indices = m + 1 + (k-1) : m + n*m : m + n*m + (k-1);
    row(y_indices) = -alpha;
    Aeq = [Aeq; row];
    beq = [beq; 0];
end

不等式约束（A, b）

A = [];
b = [];

% 处理能力约束
for k = 1:m
    row = zeros(1, total_vars);
    row(k) = -Lk(k);
    y_indices = m + 1 + (k-1) : m + n*m : m + n*m + (k-1);
    row(y_indices) = 1;
    A = [A; row];
    b = [b; 0];
end

% 肥料需求约束
for j = 1:J
    row = zeros(1, total_vars);
    z_indices = m + n*m + j : m*J : m + n*m + m*J;
    row(z_indices) = 1;
    A = [A; row];
    b = [b; Nj(j)];
end

4. 设置变量上下界（lb, ub） #

lb = zeros(total_vars, 1);
ub = ones(total_vars, 1);

% 调整Y、Z的上界为合理大数
y_ub = 1e6;
z_ub = 1e6;
ub(1:m) = 1;
ub(m+1:m+n*m) = y_ub;
ub(m+n*m+1:end) = z_ub;

5. 调用linprog求解 #

options = optimoptions('linprog','Display','iter');
[x,fval,exitflag,output] = linprog(f,A,b,Aeq,beq,lb,ub,options);

内容的提问来源于stack exchange，提问作者Carlo Soares