助力國賽 | 第3彈 規劃問題（MATLAB版）

助力國賽 | 第3彈 規劃問題（MATLAB版）

來自專欄 MATLAB學習15 人贊了文章

前言

上次介紹了回歸模型，用到的方法主要是擬合和插值，如果把回歸問題看成是預測問題的一類，其實我們還可以使用灰色預測和人工神經網路（ANN）等更高級的方法來進行求解。

本次主要介紹另一種比較常用的模型——規劃模型。規劃類問題是常見的數學建模問題，把一些優化問題進行離散化處理後就可以使用規劃模型來進行求解，再加上這類問題佔據了數學建模的「大片江山」，那麼快速求解規劃類問題就成為了參賽隊員的基本素養。MATLAB作為一款強大的通用型數學軟體，在處理這類問題上，也封裝了一些標準函數，可以簡單快速的得到所要的結果，但有一個缺點就是不直觀、甚至有些難以理解，另一款專門用於求解規劃問題的軟體——Lingo就在方面做的好多了，這將會是下一次主要介紹的對象。所以，對於一個規劃問題，首先選擇的是Lingo，如果涉及到演算法改進等方面才會去考慮用MATLAB，不過掌握基本的操作還是很必要的

這次主要介紹常見規劃模型的MATLAB求解，包括線性規劃、整數規劃和非線性規劃三個部分。一般說來，如果能熟練掌握這幾部分的操作，可以解決絕大部分的規劃模型的求解問題。對於像規劃、優化這一類問題，除了有這些基本操作，還存在更「高級」的方法，如智能演算法和計算機模擬等，這些會在後續介紹到，那就敬請期待吧。

• 線性規劃
• 非線性規劃
• 整數規劃

線性規劃

數學規劃是研究如何在有限資源的情況下取到最大效益的問題，是運籌學的一個重要分支，而線性規劃(Linear Programming 簡記LP)則是數學規劃的一個重要分支。自從1947 年G. B. Dantzig 提出求解線性規劃的單純形方法以來，線性規劃在理論上趨向成熟，在實用中日益廣泛與深入。特別是在計算機能處理成千上萬個約束條件和決策變數的線性規劃問題之後，線性規劃的適用領域更為廣泛了，已成為現代管理中經常採用的基本方法之一。

事實上，線性規劃對於我們來說已經是多年的「老朋友」了，在初高中就已經與它相遇了。那時，老師告訴我們，首先在紙上畫一「十字架」，然後把各函數圖像畫上去，用筆或者尺子比劃一下就可以得出答案了。承認這種方法簡單方便，但是我們都已經不再是當年的「少年」了，我們有更好的工具了，追求的是自動化和系統化了，只希望幾行代碼就解決事。MATLAB在這方面做的還不錯。由於MATLAB已經把規劃模型給標準化了，我們就必須了解它是如何「理解」線性規劃問題的。

線性規劃的 MATLAB 標準形式

規劃問題無非就兩類問題，目標函數的最大值和最小值問題。MATLAB很「聰明」，把兩類問題合成了一類問題，也就是說在MATLAB中求解出來的是目標函數的最小值問題，那如果是最大值問題呢？目標函數前面添負號。

MATLAB中規定線性規劃的標準形式為：

基本函數形式為linprog(c,A,b)，它的返回值是向量x的值。還有其它的一些函數調用形式（在Matlab 指令窗運行help linprog 可以看到所有的函數調用形式），如：

[x,fval]=linprog(c,A,b,Aeq,beq,LB,UB,X0,OPTIONS)

這裡fval 返回目標函數的值；Aeq和beq對應等式Ax = b；LB 和UB 分別是變數x的下界和上界；x0是x的初始值；OPTIONS 是控制參數。

如果是求解目標函數最大值問題

可以改寫成這樣的MATLAB標準型

應用

看如下一個例子：

求解線性規劃問題：

code:

c = [2;3;1];a = [1,4,2;3,2,0];b = [8,6];[x,y] = linprog(c,-a,-b,[],[],zeros(3,1))

輸出結果為：

再來看一個例子：

先將其化成求解最小值問題

code:

format longf = [-170.8582 17.7254 -41.2582 -2.2182 -131.8182 500000];A = [1 -0.17037 -0.5324 0 1 0;0 0.17037 0.5324 0 0 0;1 0.32 1 0 0 0;0 1 0 0 0 0;0 0 1 1 0 0;0 0 0 -1 -1 0];b = [0;888115;166805;521265.625;683400;-660000];Aeq = [0 0 0 0 0 1];beq = [1];lb = [0;0;0;0;0;0];[x,fval,exitflag,output,lambda] = linprog(f,A,b,Aeq,beq,lb,[])

結果為：

即：

非線性規劃

如果目標函數或約束條件中包含非線性函數，就稱這種規劃問題為非線性規劃問題，簡稱NP。一般說來，解非線性規劃要比解線性規劃問題困難得多。而且，也不像線性規劃有單純形法這一通用方法，非線性規劃目前還沒有適於各種問題的一般演算法，各個方法都有自己特定的適用範圍。

MATLAB 中非線性規劃的數學模型寫成以下形式：

其中f (x)是標量函數，A,B,Aeq,Beq是相應維數的矩陣和向量，C(x),Ceq(x)是非線性向量函數。

MATLAB中的命令是

X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)

它的返回值是向量x，其中FUN 是用M 文件定義的函數f(x)；X0 是x

的初始值；A,B,Aeq,Beq 定義了線性約束A* X ≤ B, Aeq * X = Beq ，如果沒有線性約束，則A=[],B=[],Aeq=[],Beq=[]；LB 和UB 是變數x的下界和上界，如果上界和下界沒有約束，則LB=[]，UB=[]，如果x無下界，則LB 的各分量都為-inf，如果x無上界，則UB的各分量都為inf；NONLCON 是用M 文件定義的非線性向量函數C(x),Ceq(x)；OPTIONS定義了優化參數，可以使用Matlab 預設的參數設置。

看如下一個例子：

編寫M文件，定義目標函數：

function f = fun00(x)f = -(sqrt(x(1)) + sqrt(x(2))+ sqrt(x(3))+ sqrt(x(4)));end

編寫M文件，定義約束條件：

function[g,ceq] = mycon0(x)g(1) = x(1)-400;g(2) = 1.1*x(1) + x(2) - 440;g(3) = 1.21*x(1) +1.1* x(2) + x(3) - 480;g(4) = 1.331*x(1) +1.21*x(2)+1.1*x(3)+x(4) - 532.4;ceq = 0;end

編寫主程序：

x0 = [1;1;1;1];lb = [0;0;0;0];ub =[];A = [];b = [];Aeq = [];beq =[];[x,fval] = fmincon(fun00,x0,A,b,Aeq,beq,lb,ub,mycon0)

結果為：

二次規劃

若某非線性規劃的目標函數為自變數x 的二次函數，約束條件又全是線性的，就稱這種規劃為二次規劃，這一類規劃問題比較特殊，MATLAB專門有函數求解，簡單介紹一下。

MATLAB 中二次規劃的數學模型可表述如下:

這裡H是實對稱矩陣，f,b是列向量，A是相應維數的矩陣。

Matlab中求解二次規劃的命令是

[X,FVAL]= QUADPROG(H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS)

返回值X 是決策向量x的值，返回值FVAL 是目標函數在x處的值。（具體細節可以參看在Matlab 指令中運行help quadprog 後的幫助）。

看如下一個例子：

求解代碼如下：

h = [4,-4;-4,8];f = [-6;-3];a = [1,1;4,1];b = [3;9];[x,value] = quadprog(h,f,a,b,[],[],zeros(2,1))

結果：

罰函數法

利用罰函數法，可將非線性規劃問題的求解，轉化為求解一系列無約束極值問題，

因而也稱這種方法為序列無約束最小化技術，簡記為SUMT (Sequential Unconstrained Minization Technique)。

罰函數法求解非線性規劃問題的思想是，利用問題中的約束函數作出適當的罰函

數，由此構造出帶參數的增廣目標函數，把問題轉化為無約束非線性規劃問題。主要有

兩種形式，一種叫外罰函數法，另一種叫內罰函數法，下面介紹外罰函數法。

考慮問題：

取一個充分大的數M > 0，構造函數

Matlab中可以直接利用max、min和sum函數。則以增廣目標函數P(x, M )為目標函數的無約束極值問題

的最優解x也是原問題的最優解。

看如下一個例子：

編寫M文件test0.m

function g =test0(x)M = 50000;f = x(1)^2 + x(2)^2 +8;g = f - M*min(x(1),0) - M*min(x(2),0) - M*min(x(1)^2 - x(2),0 + M*abs(-x(1) - x(2)^2 + 2));end

在MATLAB命令行窗口輸入：

[x,y] = fminunc(test0,rand(2,1))

即可求出答案。

整數規劃

規劃中的變數（部分或全部）限制為整數時，稱為整數規劃。若在線性規劃模型中，變數限制為整數，則稱為整數線性規劃。目前所流行的求解整數規劃的方法，往往只適用於整數線性規劃。目前還沒有一種方法能有效地求解一切整數規劃。因此下面簡單介紹幾種整數規劃類型和求解方法。

0-1整數規劃

0 ?1型整數規劃是整數規劃中的特殊情形，它的變數xj僅取值0 或1。這時xj稱為0 ?1變數，或稱二進位變數。xj僅取值0 或1 這個條件可由下述約束條件：0 ≤ xj≤ 1，整數所代替，是和一般整數規劃的約束條件形式一致的。在實際問題中，如果引入0 ?1變數，就可以把有各種情況需要分別討論的線性規劃問題統一在一個問題中討論了。

看如下一個例子並給出了思路：

求解思路：

• 先試探性求一個可行解，易看出(x1, x2, x3)= (1 0 , 0 , )滿足約束條件，故為一個可行解，且z = 3。
• 因為是求極大值問題，故求最優解時，凡是目標值z < 3的解不必檢驗是否滿足約束條件即可刪除，因它肯定不是最優解，於是應增加一個約束條件（目標值下界）：
• 改進過濾條件。
• 由於對每個組合首先計算目標值以驗證過濾條件，故應優先計算目標值z大的組合，這樣可提前抬高過濾門檻，以減少計算量。

按上述思路，最後求解得：

隨機取樣法

非線性規劃的求解本身就是一個難題，更何況非線性整數規劃呢.目前對於非線性整數規劃尚未有一種成熟而準確的求解方法，但是考慮到整數解是有限個，因而可以考慮枚舉法。當然，當自變數維數很大和取值範圍很寬情況下，企圖用顯枚舉法（即窮舉法）計算出最優值是不現實的，但是應用概率理論可以證明，在一定的計算量的情況下，完全可以得出一個滿意解。

考慮如下一個問題：

首先編寫M 文件mente.m定義目標函數f 和約束向量函數g，程序如下：

function [f,g]=mente(x)f=x(1)^2+x(2)^2+3*x(3)^2+4*x(4)^2+2*x(5)-8*x(1)-2*x(2)-3*x(3)-...x(4)-2*x(5);g=[sum(x)-400x(1)+2*x(2)+2*x(3)+x(4)+6*x(5)-8002*x(1)+x(2)+6*x(3)-200x(3)+x(4)+5*x(5)-200];end

編寫M文件ainint.m如下求問題的解：

rand(state,sum(clock)); p0=0; tic for i=1:10^6 x=99*rand(5,1); x1=floor(x);x2=ceil(x); [f,g]=mente(x1); if sum(g<=0)==4 if p0<=f x0=x1;p0=f; end end [f,g]=mente(x2); if sum(g<=0)==4 if p0<=f x0=x2;p0=f; end endendx0,p0toc

由於是隨機取樣，因而每次運行的結果可以不一樣，但差別不大。

編輯不易，歡迎推廣

推薦閱讀：