本推送總結
數學建模
中常用的
數學規劃
模型，並附詳細的MATLAB求解案例。分為四個模組：
求解數學模型的一般步驟如下：
讀題+理解模型；
設計
決策變數
；
列出
目標函式
及
約束條件
；
表示（畫）出約束條件代表的可行域；
在可行域內求目標函式的最優解及最優值；
如果你還不太理解，跟著例子往下看...
1 線性規劃問題（LP）
求
線性目標函式
在
線性約束條件
下的最大值或最小值的問題，統稱為
線性規劃
問題。
舉一個簡單案例：
例：
某機床廠生產甲、乙兩種機床，每臺銷售後的利潤分別為4000元與3000元。生產甲機床需用A、B機器加工，加工時間分別為每臺2小時和1小時；生產乙機床需用A、B、C三種機器加工，加工時間為每臺各一小時。若每天可用於加工的機器時數分別為A機器10小時、B機器8小時和C機器7小時，問該廠應生產甲、乙機床各幾臺，才能使總利潤最大？
根據上述描述建立數學模型：
設該廠生產
x1
臺甲機床和
x2
乙機床時總利潤最大，則
x1
，
x2
應滿足：
上式中：
x1
，
x2
稱之為決策變數； -（1）式被稱為問題的目標函式； -（2）中的幾個不等式是問題的約束條件，記為
s.t.
（即 subject to）；
上面的目標函式及約束條件均為線性函式，故被稱為線性規劃問題。
在解決實際問題時，把問題歸結成一個線性規劃數學模型非常重要，但往往也最困難，模型建立得是否恰當，直接影響到求解。
選適當的決策變數，是建立有效模型的關鍵之一。
線性規劃的目標函式可以是求最大值，也可以是求最小值，約束條件的不等號可以是小於號也可以是大於號。為了避免這種形式多樣性帶來的不便，
MATLAB中規定線性規劃的標準形式為
：
其中
c
和
x
為 n 維列向量，
A
、
Aeq
為適當維數的矩陣，
b
、
beq
為適當維數的列向量；lb與ub分別為x的下界和上界，維度與x對應。
模型已經表示為MATLAB的標準形式，下面就可以採用MATLAB進行求解了，此處採用
linprog
函式。
風格1
function lp_
% MATLAB求解線性規劃問題，linprog
%
c = -［4000；3000］； % 目標函式係數
A = ［2 1；1 1；0 1］；
b = ［10；8；7］；
Aeq = ［］；
beq = ［］；
lb = ［0 0］； % 變數下界
x = linprog（c，A，b，Aeq，beq，lb）
fval = -c‘*x
需要注意的就是將約束條件轉換為合適維度的矩陣！
如果你不太想進行矩陣（係數）轉換，也可以採用另一種程式設計風格：
風格2
% 返回的fval是負的，
% 即使問題是正的，
% 在matlab內部，prob2struct將最大化問題轉化為目標函式負的最小化問題
x1 = optimvar（’x1‘，’LowerBound‘，0，’UpperBound‘，inf）；
x2 = optimvar（’x2‘，’LowerBound‘，0，’UpperBound‘，7）；
prob = optimproblem（’Objective‘，4000*x1 + 3000*x2，’ObjectiveSense‘，’max‘）；
prob。Constraints。c1 = 2*x1 + x2 <= 10；
prob。Constraints。c2 = x1 + x2 <= 8；
problem = prob2struct（prob）；
［sol，fval，~，output］ = linprog（problem）
相信透過上述的簡單例子，你已經對線性問題的建模和求解有了一定的瞭解，下面我們看更復雜的模型。
2 非線性規劃
若目標函式·或·約束條件中包含非線性函式，就稱這種規劃問題為非線性規劃問題。
非線性規劃問題有無約束和有約束兩種，它們的求解思路不太一致，具體的方法可以參考“運籌學”或者“最最佳化”的書籍；其中一類比較特殊的問題是
二次規劃問題
。
一般說來，解非線性規劃要比解線性規劃問題困難得多。而且，也不像線性規劃有
單純形法
這一通用方法，非線性規劃目前還沒有適於各種問題的一般演算法，各方法都有特定的適用範圍。
線性規劃與非線性規劃的區別在於
：如果線性規劃的最優解存在，則其最優解只能在其可行域的邊界上達到（特別是可行域的頂點上達到），而非線性規劃的最優解（如果存在）可能在其可行域的任意一點達到。下方線性問題的圖解法很容易反映解在頂點這一性質：
非線性規劃問題在MATLAB裡的標準形式為：
在非線性求解時常用的函式為
fmincon
，舉一個示例模型來看看它的用法：
下方為
fmincon
的呼叫格式：
不包含非線性約束
［X，FVAL，EXITFLAG，OUTPUT，LAMBDA］ = fmincon（FUN，X0，A，B，Aeq，Beq，LB，UB）
包含非線性約束
［X，FVAL，EXITFLAG，OUTPUT，LAMBDA］ = fmincon（FUN，X0，A，B，Aeq，Beq，LB，UB，NONLCON）
自己定義目標函式和非線性約束
［X，FVAL，EXITFLAG，OUTPUT，LAMBDA］ = fmincon（@（X）MYOBJ（X），X0，A，B，Aeq，Beq，LB，UB，@（X）MYCON（X））性約束函式
% 目標函式
function F = MYOBJ（X）
F = ……
% 非線性約束函式
function ［G，Heq］ = MYCON（X）
G = 。。。。。 % 非線性不等式約束條件，G<=0
Heq = 。。。。。 % 非線性等式約束條件，Heq == 0
（該模型的程式碼在文末下載）
採用fmicon的求解結果為：
除了採用預設的
內點法
外，也可以根據需要為fmincon函式設定不同的尋優演算法：
接下來是動態規劃問題。
3 動態規劃
動態規劃基本思路
將多階段決策過程劃分階段，恰當地選擇狀態變數、決策變數以定義最優指標函式，從而把問題化成一族同類型的子問題，然後逐個求解；
求解時從邊界條件開始，逆序過程行進，逐段遞推尋優．在每一個子問題求解時，都要使用它前面已求出的子問題的最優結果．最後一個子問題的最優解，就是整個問題的最優解；
動態規劃方法是既將當前一段與未來各段分開，又把當前效益和未來效益結合起來考慮的一種最最佳化方法，因此每段的最優決策選取是從全域性考慮的，與該段的最優選擇一般是不同的；
學習動態規劃，我們首先要了解多階段決策問題。
比如：
生產決策問題
：企業在生產過程中，由於需求是隨時間變化的，因此企業為了獲得全年的最佳生產效益，就要在整個生產過程中逐月或逐季度地根據庫存和需求決定生產計劃。
機器負荷分配問題
：某種機器可以在高低兩種不同的負荷下進行生產。要求制定一個五年計劃，在每年開始時，決定如何重新分配完好的機器在兩種不同的負荷下生產的數量，使在五年內產品的總產量達到最高。
太空梭飛行控制問題
：由於太空梭的運動的環境是不斷變化的，因此就要根據太空梭飛行在不同環境中的情況，不斷地決定太空梭的飛行方向和速度（狀態），使之能最省燃料和完成飛行任務（如軟著陸）。
針對多階段決策過程的最最佳化問題，美國數學家Bellman等人在20世紀50年代初提出了著名的最最佳化原理，把多階段決策問題轉化為一系列單階段最最佳化問題，從而逐個求解，創立了解決這類過程最佳化問題的新方法：
動態規劃
。對最佳路徑（最佳決策過程）所經過的各個階段，其中每個階段始點到全過程終點的路徑，必定是該階段始點到全過程終點的一切可能路徑中的最佳路徑（最優決策），這就是Bellman提出的著名的
最最佳化原理
。即：
一個最優策略的子策略必然也是最優的。
我們還是直接看一個求解案例：
例
：設有m=4個目標，目標價值（重要性和危害性）各不相同，用數值Ak（k=1，2，…，m）表示，計劃用n=5枚導彈突襲，導彈擊毀目標地機率
中a_k是常數，取決於導彈地特性與目標地性質；u_k為向目標發射地導彈數，如何做出方案使預期地突擊效果最大？
首先需要
建立問題的模型
：
根據模型和動態規劃演算法建立函式檔案：
主函式
clear
n = 5；
x1 = ［n；nan*ones（n，1）］；
x2 = ［0：n］’；
x = ［x1 x2 x2 x2］
［p f］ = dynprog（x，‘DecisF1’，‘SubObjF1’，‘TransF1’，‘ObjF1’）
disp（‘p的第三列是決策變數值’）
決策函式
function u = DecisF1（k，x）
% 決策函式
if k ==4
u =x；
else
u = 0：x；
end
end
狀態轉移方程
function y = TransF1（k，x，u）
% 狀態轉移方程
y = x-u；
end
階段k的指標函式
function v = SubObjF1（k，x，u）
% 階段k的指標函式
a = ［0。2 0。3 0。5 0。9］；
A = ［8 7 6 3］；
v = A（k）*（1-exp（-a（k）*u））；
v = -v；
end
基本方程中的函式
function y = ObjF1（v，f）
% 基本方程中的函式
y = v+f；
% y = -y；
end
動態規劃dynprog函式
function ［p_opt， fval］ = dynprog（x，DecisFun， SubObjFun， TransFun， ObjFun）
% 動態規劃逆序演算法
% x是狀態變數，一列代表一個階段狀態；
% DecisFun（k，x）由階段k的狀態變數x求出相應的允許決策變數；
% SubObjFun（k，x，u）是階段指標函式；
% TransFun（k，x，u）是狀態轉移函式，其中x是階段k的某狀態變數，u是相應的決策變數；
% ObjFun（v，f）是第k階段至最後階段指標函式，當ObjFun（v，f） = v+f時，輸入的ObjFun可以省略；
% fval是一個列向量，各元素分別表示p_opt各最優策略組對應始端狀態x的最優函式值
完整的程式碼在文末
在動態規劃問題中，
路徑規劃
也非常重要，下方列出了關於
A星演算法
和
機率路線圖（基於dijkstra演算法）
的求解案例：
A星演算法
function path=AStar（obstacle，map）
path=［］； %用於儲存路徑
open=［］； %OpenList，考慮的節點
close=［］； %CloseList，不考慮的節點
findFlag=false； %findFlag用於判斷while迴圈是否結束
%% 1。將起始點放在Openlist中
% open：［節點座標，代價值F=G+H，代價值G，父節點座標］
open =［map。start（1）， map。start（2） ，。。。
0+h（map。start，map。goal） ，。。。
0 ， 。。。
map。start（1） ， map。start（2）］；
%與當前節點的座標差值（前兩列）
%與當前節點的距離值（最後一列）
NEXT = ［-1，1，14；。。。
0，1，10；。。。
1，1，14；。。。
-1，0，10；。。。
1，0，10；。。。
-1，-1，14；。。。
0，-1，10；。。。
1，-1，14］；
%% 2。 迴圈以下過程
while ~findFlag
%——————————首先判斷是否達到目標點，或無路徑——-
if isempty（open（：，1））
disp（‘沒有目標路徑’）；
return；
end
%判斷目標點是否出現在open列表中
［isopenFlag，Id］=isopen（map。goal，open）；
if isopenFlag
disp（‘找到目標’）；
close = ［open（Id，：）；close］；
findFlag=true；
break；
end
%按照Openlist中的第三列（代價函式F）進行排序，查詢F值最小的節點
［Y，I］ = sort（open（：，3））； %對OpenList中第三列排序
open=open（I，：）；%open中第一行節點是F值最小的
%將F值最小的節點（即open中第一行節點），放到close第一行（close是不斷積壓的），作為當前節點
close = ［open（1，：）；close］；
current = open（1，：）；
open（1，：）=［］；%因為已經從open中移除了，所以第一列需要為空
%對當前節點周圍的8個相鄰節點進行計算，演算法的主體：————————————
for in=1：length（NEXT（：，1））
m=［current（1，1）+NEXT（in，1） ， current（1，2）+NEXT（in，2） ， 0 ， 0 ， 0 ，0］；
m（4）=current（1，4）+NEXT（in，3）； % m（4） 相鄰節點G值
m（3）=m（4）+h（m（1：2），map。goal）；% m（3） 相鄰節點F值
if isObstacle（m，obstacle）
continue；
end
［flag，targetInd］=findIndex（m，open，close）；
if flag==1
continue；
elseif flag==2
m（5：6）=［current（1，1），current（1，2）］；
open = ［open；m］；
else
if m（3） < open（targetInd，3）
m（5：6）=［current（1，1），current（1，2）］；
open（targetInd，：） = m；
end
end
end
%繪製
pause（0。01）；
fillPlot（close，［204，159，42］/255）；
hold on；
fillPlot（open，［51，102，153］/255）；
end
%% 3。 追溯路徑
path=GetPath（close，map。start）；
end
基於dijkstra的機率路線圖
% 載入地圖
% —————————————————————————————————————
% PIC = imread（uigetfile（‘*。png’））； % 載入地圖
PIC = imread（‘map2。png’）； % 載入地圖
imshow（PIC）； % 顯示地圖
hold on；
% 關鍵引數
% —————————————————————————————————————
feasibleR = 51； % 可行區域的紅色值【0~255】
NUM = 200； % 離散點規模
% 一些樣式設定，可以註釋掉
% —————————————————————————————————————
ax = gca；
ax。Parent。Units = ‘normalized’；
ax。Parent。Color = ‘k’；
ax。Parent。Position = ［0。1 0。1 0。8 0。7］；
ax。Units = ‘normalized’；
ax。Position = ［0。1 0。001 0。8 0。998］；
ax。Color = ［51 51 51］/255；
% 選擇起/始點
% —————————————————————————————————————
StartPoint = ginput（1）；
EndPoint = ginput（1）；
scatter（［StartPoint（1） EndPoint（1）］，［StartPoint（2） EndPoint（2）］，。。。
‘filled’，‘SizeData’，150，‘MarkerFaceColor’，［199 237 53］/255，‘MarkerEdgeColor’，‘w’）；
% 繪製出起始點
pause（1）
% 生成並繪製可行點
%__________________________________________________________________________
Height = size（PIC，1）；
Width = size（PIC，2）；
randPoints = floor（［rand（NUM，1）*Height，。。。
rand（NUM，1）*Width］）+1； % 空間隨機撒點
allPoints = ［StartPoint；EndPoint］； % 初始化可行點列表
for i=1：NUM
if（PIC（randPoints（i，2），randPoints（i，1），1） == feasibleR ）
allPoints = ［allPoints；double（randPoints（i，：））］；
end
end% 得到所有可行點列表
for i = 1：size（allPoints，1）
scatter（allPoints（i，1），allPoints（i，2），‘filled’，‘SizeData’，20，‘MarkerFaceColor’，［199 237 53］/255）；
pause（。0005）
end
% 初始化無向圖
鄰接矩陣
%__________________________________________________________________________
AdjacencyMatrix = genAM（PIC，allPoints）；
% 對矩陣AdjacencyMatrix使用【Dijkstra最短路徑演算法】；
%__________________________________________________________________________
4 多目標規劃
上面的問題無論約束條件有多複雜，目標都是隻有一個。但是在現實中很多問題都能最終抽象為
多目標最佳化問題
。這是因為現實問題往往比較複雜，一般在達成目標時會權衡很多方面。這裡簡單介紹主流的幾種方法。
包含多個可能有衝突的目標函式
希望找到能夠很好平衡全部最佳化目標的解集
理解多目標最佳化，需要先理解
帕累托最優
（Pareto Optimal）
帕累托最優
帕雷託最優是指資源分配的一種理想狀態。給定固有的一群人和可分配的資源，如果從一種分配狀態到另一種狀態的變化中，在
沒有使任何人境況變壞的前提下，使得至少一個人變得更好
，這就是
帕雷託
改善。
帕雷託最優的狀態就是不可能再有更多的帕雷託改善的狀態；換句話說，不可能在不使任何其他人受損的情況下再改善某些人的境況。
支配（Dominace）
當x1和x2滿足如下條件時稱x1支配x2：
對於所有目標函式x1不比x2差； 至少在一個目標函式上，x1嚴格比x2要好
；
（上圖中點1支配點2；點5支配點1；點1和點4互不支配）
不可支配解集
（Non-dominated solution set） 當一個解集中任何一個解都不能被該集合中其他解支配，那麼就稱該解集為不可支配解集。
帕累托最優解集
（Pareto-optimal set） 所有可行解的不可支配解集被成為帕累托最優解集。
帕累托最優前沿面
（Pareto-optimal front） 帕累拖最優解集的邊界（boundary）被成為帕累托最優前沿面。
多目標最佳化問題的經典方法有哪些？
線性加權法
其中權重代表了每個目標函式的重要程度。
優點
是簡單；
缺點
是很難設定一個權重向量能夠去獲得
帕累托最優解
；在一些非凸情況不能夠保證獲得帕累托最優解。
約束轉化法
（注：以上關於多目標最佳化的部分配圖來源於：
https：//
hpzhao。github。io
）
多目標遺傳演算法
基於遺傳演算法的多目標最佳化就是利用遺傳演算法的原理來搜尋帕累托最優前沿面。關於遺傳演算法的介紹在往期推送已有，直接搜尋歷史推送即可。
其基本流程如下：
隨機產生初始群體；
計算各點的目標函式值和約束函式值；
根據目標函式值對群體分級；
根據約束函式值和分級結果計算各點的約束罰項、劣解罰項及總罰項；
根據各點的總罰項計算
適應度
；
根據各點的適應度，進行選擇、交叉和變異操作，生成新群體；
將總罰項為0的點放入非劣解集侯選表，對侯選表進行檢查，保留第1級非劣點 ， 刪除其它點；
檢查是否收斂，如沒有，轉到步驟2；
刪除侯選表中與其它點距離太近的點；
輸出侯選表中的Pareto最優解集及對應的目標函式值；
決策人根據個人偏好從Pareto最優解集中挑選出最適合該問題的解；
遺傳演算法相比與傳統演算法的優點是能夠得到一個最優解集，而不是單單一個最優解，這樣給我們更多的選擇。但計算複雜度可能稍高，而且裡面涉及的一些函式需要精心設計。
多目標最佳化是一個非常複雜的問題，這次肝不動了，之後的推送再介紹（100個“在看”），有一箇中文期刊裡的案例可以參考《飛機方案多目標最佳化的Pareto遺傳演算法》
交流群：1129425848
