這篇隨筆原本是我上實驗課時候的筆記，2023 年 7 月曾經在 CSDN 平臺上 釋出過。
今天恰好有朋友跟我問起 MATLAB 自帶的求解器輸入很不直觀的問題，我開啟這個文章發給他的時候發現自己一年前寫的 LaTeX 公式依託答辯，所以重打了一遍。再加上由於 CSDN 平臺的持續擺爛，終於是用不下去了。原文裡的公式使用的諸如
array一類的環境本來就不是很魯棒，被 CSDN 這個大概是沒最佳化過（應該是沒最佳化過，因為在部落格園就沒什麼問題）的渲染一折騰，文章裡的數學公式現在極其混亂，徹底沒法看了。考慮到這個破文裡的內容可能確實還有點用，這裡搬運到部落格園再發一遍。由於文章內容已經是一年多前的東西了，那時候寫的東西可能很多方面沒有估計到，不太完善，還請讀者見諒。
- 前言
- linprog() 的用法以及 MATLAB 的標準型
- 函式功能的實現
- 運輸問題
- 運輸問題的程式碼實現
前言
前段時間我們學校開始上 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() 或者非標準的運輸問題。