【旧文重发】MATLAB 通过函数封装一劳永逸地解决线性规划与运输问题的linprog的标准化操作（附MATLAB代码）
这篇随笔原本是我上实验课时候的笔记，2023 年 7 月曾经在 CSDN 平台上 发布过。
今天恰好有朋友跟我问起 MATLAB 自带的求解器输入很不直观的问题，我打开这个文章发给他的时候发现自己一年前写的 LaTeX 公式依托答辩，所以重打了一遍。再加上由于 CSDN 平台的持续摆烂，终于是用不下去了。原文里的公式使用的诸如
array一类的环境本来就不是很鲁棒，被 CSDN 这个大概是没优化过（应该是没优化过，因为在博客园就没什么问题）的渲染一折腾，文章里的数学公式现在极其混乱，彻底没法看了。考虑到这个破文里的内容可能确实还有点用，这里搬运到博客园再发一遍。由于文章内容已经是一年多前的东西了，那时候写的东西可能很多方面没有估计到，不太完善，还请读者见谅。
前言
前段时间我们学校开始上 MATLAB 的实验课了，就开始学习做这个线性规划，还要写报告什么的。在 MATLAB 里面做线性规划就要用到 linprog() 这个函数。
众所周知，每一次使用 linprog() 之前都要先把线性规划问题化为 MATLAB 标准型。
这就很烦，看到线性规划问题之后不能直接扔给 MATLAB，非得先掏出草稿纸一通哐哐化简然后才能敲进脚本。然后去机房还要带一堆草稿纸，更何况我还要写作业报告，每一次都要在报告里面写一大堆公式解释标准化的过程，我就觉得很烦。遇到运输问题形式的就更烦。系数矩阵全是 0 和 1，看得我眼睛都花了。
所以我们干嘛不写个函数脚本让 MATLAB 自己自动地完成标准化呢？
这个应该很多人都写过。但是我在网上没搜到，所以我就自己写了。写的比较草率，但是能用。（学校的作业写完顺便来水一篇博客。）
linprog() 的用法以及 MATLAB 的标准型
linprog() 用法如下：
x = linprog(f,A,b)
x = linprog(f,A,b,Aeq,beq)
x = linprog(f,A,b,Aeq,beq,lb,ub)
x = linprog(f,A,b,Aeq,beq,lb,ub,options)
x = linprog(problem)
[x,fval] = linprog(___)
[x,fval,exitflag,output] = linprog(___)
[x,fval,exitflag,output,lambda] = linprog(___)
不多赘述。可以看见这里传入的参数中，f,A,b,Aeq,beq,lb,ub 这几个参数是从原始线性规划里面得出来的，所以我们关键就是要让函数脚本自己生成这几个参数。
MATLAB 中的标准型如下：
我们的思路是这样的：把我们要求解的线性规划问题写进一个矩阵，然后让 MATLAB 自己去切割矩阵，加正负号调整 max 和 min 什么的，返回 f,A,b,Aeq,beq,lb,ub 这几个参数。
打个比方说，假如现在有如下线性规划：
用矩阵 opt_mat 表示线性规划的形式，矩阵的第一行表示约束条件，按照小于、等于、大于的顺序从上往下写。第二行到倒数第三行表示约束条件， le、eq、ge 分别表示小于等于约束、等式约束和大于等于约束的个数。最后两行是上界和下界，其余部分补 0。
那么这个线性规划问题就可以被简化为：
函数功能的实现
可以用如下的代码实现标准化：
function [f, A, b, Aeq, beq, lb, ub] = ...
OptStandardization( opt_mat, le, eq, ge, max_or_min )
%OPTSTANDARDIZATION 用于线性规划问题 linprog() 的标准化
% 矩阵 opt_mat 是表示线性规划问题形式的矩阵
%
% 从上往下写
% 先写目标函数，把数字排成一排，最后一个数字补上 0
% 再写小于约束
% 然后是等于约束
% 然后是大于约束
% le, eq, ge,
% 分别表示小于、大于、等于约束的个数
% 然后下面一行写下界
% 再下面一行写上界
% max_or_min 传入 "min" 或者 "max"
% 表示约束的类型，默认是 "min"
[ row, col ] = size(opt_mat);
if max_or_min == "min"
f = ( opt_mat( 1, ( 1:(col-1) ) ) )';
elseif max_or_min == "max"
f = ( -opt_mat( 1, ( 1:(col-1) ) ) )';
else
f = ( opt_mat( 1, ( 1:(col-1) ) ) )'; %默认是 "min"
end
if le ~= 0
A = opt_mat( ( 2:(le+1) ), (1:( col-1 ) ) );
b = ( opt_mat( ( 2:(le+1) ), col ) )';
else
A = [];
b = [];
end
if eq ~= 0
Aeq = opt_mat( ( (le+2):(le+eq+1) ), (1:col-1) );
beq = ( opt_mat( ( (le+2):(le+eq+1) ), col ) )';
else
Aeq = [];
beq = [];
end
if ge ~= 0
A = [ A ; ...
(-opt_mat( ( (le + eq + 2): ...
(le + eq + ge + 1) ), (1:col-1) ) ) ];
b = [ b ; ...
( -opt_mat( ( (le + eq + 2): ...
(le + eq + ge + 1) ), col ) )' ];
end
lb = ( opt_mat( row-1, ( 1:( col-1 ) ) ) )';
ub = ( opt_mat( row, ( 1:( col-1 ) ) ) )';
end
接下来的事情就变得很简单了。当我们要求解这个线性规划问题的时候，只需要把问题写作：
opt_mat = [ 2 3 -5 0 ;
1 3 1 12 ;
1 1 1 7 ;
2 -5 1 10 ;
0 0 0 0 ;
inf inf inf 0 ];
然后先调用 OptStandardization()，再把返回的参数传递给 linprog()，就结束了。
可以参考下面的完整代码：
%% 线性规划模型计算
%% 清理工作区
clc; clear; close all;
%% 定义传入的矩阵
opt_mat = [ 2 3 -5 0 ;
1 3 1 12 ;
1 1 1 7 ;
2 -5 1 10 ;
0 0 0 0 ;
inf inf inf 0 ];
%% 设定小于等于、等于和大于三种约束的数量
num_of_le = 1;
num_of_eq = 1;
num_of_ge = 1;
%% 标准化
[f, A, b, Aeq, beq, lb, ub] = ...
OptStandardization( opt_mat, num_of_le, num_of_eq, num_of_ge, "max" );
%% 线性规划求解
[x,fval] = linprog(f, A, b, Aeq, beq, lb, ub);
max_fval = -fval;
运行 MATLAB 之后，解得：
完美。
运输问题
对于运输问题的情形，则要更加复杂一些，不过也好搞的。我们这里暂时只讨论产销平衡的标准运输问题。
原理的解释同样是举个例子，看下面这个问题：
已知某企业有甲、乙、丙三个分厂生产一种产品，其产量分别为 7、9、7 个单位，需运往 A、B、C、D 四个门市部，各，门市部需求量分别为 3、5、7、8 个单位已知单位运价如下表，试确定运输计划使总运费最少
| A | B | C | D | 产量 | |
|---|---|---|---|---|---|
| 甲 | 12 | 13 | 10 | 11 | 7 |
| 乙 | 10 | 12 | 14 | 10 | 9 |
| 丙 | 14 | 11 | 15 | 12 | 7 |
| 需求量 | 3 | 5 | 7 | 8 | 23 |
| : 运输问题中的常数 |
我们知道，这是一个 标准运输问题，符合 产销平衡的特点，可以转化为线性规划问题再用 MATLAB 中的 linprog() 求解。
众所周知，运输问题转化为标准线性规划问题的建模过程是非常辛苦且复杂的，主要体现为算式很长、变量很多的特点。建模过程如下：
我们设 \(x_{1,1}\) 到 \(x_{3,4}\) 总共 12 个变量，用 \(x_{i,j}\) 表示每个产地送往每个销地的商品量。用 \(X\) 表示所有变量组成的矩阵 \((x_{i,j})\)。同时，现在有一个运价矩阵 \(P\)：
那么目标函数就可以写作：
然后呢，对于产销平衡的情形，这里的约束条件有两类：
- 每个工厂运往每个销地的商品总量要等于该工厂生产的产品总量
- 各个销地收到的商品总量要等于销地的需求总量
如果我们用 \(a\) 和 \(b\) 两个向量来分别表示各个工厂生产的产品总量和各个销地需要的商品的需求量，那么问题的约束就可以表示为下面的约束：
把上面的建模过程展开之后，就是下面的标准线性规划的式子：
可以看见这个式子非常的复杂，光是把这个 \(L^AT_E\chi\) 公式写出来就花了我十几分钟的时间。
但是实际上，我们可以对这个式子进行一个简化，划归成 MATLAB 标准型。
可以看见在这个模型里面的 约束都是等式，而约束的 系数矩阵 可以表示为：
这个矩阵明显是有规律的。首先，矩阵整体上可以分为上下两部分：上半部分包含 3 个长度相当于原始运输表列数 4 的、呈现阶梯状的全为 1 的行；下半部分呈现波涛状对角线排列，可以分为 \(3\) 个 \(4 \times 4\) 的单位矩阵。
在一般情况下，标准运输问题矩阵可以表示为下面的形式：
假如我们设原始的运输表存在 \(m\) 个产地和 \(n\) 个销地，那么就可以发现：在这个矩阵中：
- 上半部分矩阵的每一个分块都具有 \(m\) 行 \(n\) 列，而且对于第 \(i\) 个块，从上往下数第 \(i\) 行元素全为 \(1\)。
- 下半部分矩阵的每一个分块都是 \(n \times n\) 的尺寸，而且对角线元素全为 1。
这不就方便了吗？
也就是说，实际上只要我们能够知道原始的运输矩阵的行和列，就能够直接写出系数矩阵 \(A_{eq}\)。
运输问题的代码实现
那么，在 MATLAB 中，我们就可以定义一个函数 STDBoolTransMat(m, n) 来生成系数矩阵 Aeq。代码如下：
function [bool_mat] = STDBoolTransMat(m, n)
%STDBOOLTRANSMAT 返回运输问题的标准系数矩阵
% function [bool_mat] = STDBoolTransMat(m, n)
% 便于把运输问题转化为标准的线性规划问题来求解
%
% 输入：
% m : 运价表的行数
% n : 运价表的列数
%
% 输出：
% bool_mat : 运输问题的标准系数矩阵，由 0 和 1 组成
%
% 注意事项：输入的是运价表的行数和列数
% 而不是完整的运输表，也就是说输入的是产地和销地的数量
% 不包含产地行和销地行
bool_mat = zeros((m+n), (m*n));
for i = 1:m
bool_mat( i:i, (n*i-(n-1)):(n*i) ) = ones(1, n);
end
for i = 1:n
for j = 1:m
bool_mat( (i+m):(i+m), (n*j-(n-i)):(n*j-(n-i)) ) = 1;
end
end
end
该函数只要输入参数 \(m\) 和 \(n\) 就能返回运输问题的系数矩阵。在命令行中输入命令测试这一功能：
>> STDBoolTransMat(3, 4)
ans =
1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0
0 0 0 0 0 0 0 0 1 1 1 1
1 0 0 0 1 0 0 0 1 0 0 0
0 1 0 0 0 1 0 0 0 1 0 0
0 0 1 0 0 0 1 0 0 0 1 0
0 0 0 1 0 0 0 1 0 0 0 1
可以看见命令行返回的输出结果是符合预期的。
我们为了实现把标准运输问题自动转化为标准线性规划问题，定义函数 STDTransMat2OptMat()：
function [opt_mat] = STDTransMat2OptMat(trans_mat)
%STDTRANSMAT2OPTMAT 用于把标准运输问题的运输表转化为线性规划的标准表
% 配合 OptStandardization.m 可以实现自动解运输问题
%
% 最后修改日期：2023/7/5
%% 获得原始运输表矩阵的行和列
[row, col] = size(trans_mat);
% 把 rol - 1 和 col - 1 分别标注为 rolt 和 colt
rowt = row - 1;
colt = col - 1;
% 代表除去需求量和产量行和列之后的运价表的尺寸
% 这是为了在下面切割矩阵的时候方便记忆
% 如果一直记着 rol - 1 和 col - 1 很容易写错代码
% 也可以说是为了方便后期的代码维护吧
%% 通过矩阵运算切割矩阵，把运输表转化为线性规划标准格式表
% 初始化运输矩阵为一个全是 0 的矩阵
% 行数为运输表的行列总数 + 3，
% 3 代表最上面的一行目标，最下面的两行上界和下界
rowo = (3 + rowt + colt);
colo = (rowt*colt + 1);
% 保存这两个变量也是为了下面方便记忆
opt_mat = zeros( rowo, colo );
%% 填写线性规划标准格式表
% 把标准格式表的第一行设置为运输问题的目标
opt_mat( 1, 1:( rowt*colt ) ) ...
= ( ( reshape( (trans_mat(1:rowt, 1:colt))' , [], 1)) )';
% 这里的逻辑是先把只包含运价的运价表割下来
% 然后用 reshape() 拉成一根直杆子
% 然后转置成行向量，赋值给运输表第一行
opt_mat(2:(rowo-2), 1:(colo-1)) = STDBoolTransMat(rowt, colt);
% 运输问题的上半部分是阶梯状的行向量排列
for i = 1:rowt
opt_mat( (i+1):(i+1), colo:colo ) ...
= trans_mat(i:i, col:col); % 在最后一列添加右端常数项
end
% 运输问题的下半部分是波涛状的对角线排列
for i = 1:4
% 填入矩阵右端的常数项
opt_mat( (row+1):(rowo-2), colo:colo ) = (trans_mat( row:row, 1:colt ))';
% row+1 在这里表示撇开运输表上半部分的行开始往下填
% 其实是 rowt + 1 + 1
% 为了方便计算机处理所以这样写
% rolo-2 在这里表示在最后两行之前截尾，最后两行是上下界
end
% 标准运输问题下界都是 0，不用管
% 上界是 inf
opt_mat(rowo:rowo, 1:(colo - 1)) = inf(1, (colo - 1));
end
用下面的命令测试一下这个函数：
>> trans_mat = [ 12, 13, 10, 11, 7 ;
10, 12, 14, 10, 9 ;
14, 11, 15, 12, 7 ;
3, 5, 7, 8, 23 ];
>> STDTransMat2OptMat(trans_mat)
ans =
12 13 10 11 10 12 14 10 14 11 15 12 0
1 1 1 1 0 0 0 0 0 0 0 0 7
0 0 0 0 1 1 1 1 0 0 0 0 9
0 0 0 0 0 0 0 0 1 1 1 1 7
1 0 0 0 1 0 0 0 1 0 0 0 3
0 1 0 0 0 1 0 0 0 1 0 0 5
0 0 1 0 0 0 1 0 0 0 1 0 7
0 0 0 1 0 0 0 1 0 0 0 1 8
0 0 0 0 0 0 0 0 0 0 0 0 0
Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf 0
可以看见函数的输出结果还是符合我们的预期。非常的完美。
这样一来，下一次要求解运输问题的时候，只要把运输表当作参数传入矩阵，就能直接求出用于 linprog() 的参数了。
唯一有点美中不足的就是这个函数里面要求传入的矩阵包含运输表的产量列和销量行。 当时写的时候没考虑到可以采取其他的更合理的方式，比如用额外的参数传入之类的。
尝试一下解决刚才那道题：我们在文件里面调用已经写好的 OptStandardization()，把返回参数 [f, A, b, Aeq, beq, lb, ub] 传入 linprog()，然后可以求出问题的线性规划结果：
clc; clear; close all;
%% 原始的运输矩阵
trans_mat = [ 12, 13, 10, 11, 7 ;
10, 12, 14, 10, 9 ;
14, 11, 15, 12, 7 ;
3, 5, 7, 8, 23 ];
%% 调用 STDTransMat2OptMat() 自动转化为标准线性规划问题的矩阵
[opt_mat] = STDTransMat2OptMat(trans_mat);
%% 设定小于等于、等于和大于三种约束的数量
num_of_le = 0;
num_of_eq = 7;
num_of_ge = 0;
%% 标准化
[f, A, b, Aeq, beq, lb, ub] =...
OptStandardization( opt_mat, num_of_le, num_of_eq, num_of_ge, "min" );
%% 线性规划求解
[x,fval] = linprog(f, A, b, Aeq, beq, lb, ub);
解得运输问题的最优运输分配方案为：
这个结果跟手敲矩阵求解的结果基本一致。下回把具体的操作步骤改一改还可以解决 intlinprog() 或者非标准的运输问题。
