1. 程式人生 > >Matlab數學建模(五):優化模型之標準模型

Matlab數學建模(五):優化模型之標準模型

一、學習目標

(1)瞭解最優化模型。

(2)掌握線性規劃的優化求解。

(3)掌握整數規劃的優化求解。

(4)瞭解Matlab的圖形化應用。

二、例項演練

     1、談談你對最優化模型的瞭解。

        最優化模型是數學建模大賽中最常見的問題型別之一。一般說來,凡是尋求最大、最小、最遠、最近、最經濟、最豐富、最高效、最耗時的目標,都可以劃入優化問題的範疇。MATLAB 優化工具箱和全域性優化工具箱對多個優化問題提供了完整的解決方案,前者涵蓋了線性規劃、混合整型線性規劃、二次規劃、非線性優化、非線性最小二乘的求解器,後者囊括了全域性搜尋、多初始點、模式搜尋、遺傳演算法等求解演算法。

       最優化即在一定的條件下,尋求使目標最小(大)的設計引數或決策。在優化問題中有兩個關鍵物件:目標函式約束條件。哈哈哈,這讓筆者想起來高中學到的線性規劃,其實,本節談到的最優化模型跟高中的線性規劃還真有點像。不過,高中的線性規劃問題約束條件很少,一般是通過作圖法來求解。本節的最優化模型約束條件一般比高中的線性規劃多很多,並且我們要通過寫程式碼去求解,而不是作圖。常規優化問題,其數學表達可以描述為:

其中x 為長度n的決策變數向量,f(x) 為目標函式,G(x) 為約束函式。

以上數學表示式看不明白也沒關係,因為下面我們會通過兩個具體的例子去講解分析。

求解目標函式的最小(大)值,一個高效而精確的解決方案不僅取決於約束條件和變數數量,更取決於目標函式和約束函式的特性。明確優化型別是確認優化方案的前提,讓我們看一下這些特性如何劃分:

常見的目標函式有:

線性規劃:被廣泛的應用於變數之間可線性表示的財務、能源、運營研究等現代管理領域中。

混合整數線性規劃:擴充套件了線性規劃問題,增加了最優解中部分或全部變數必須是整數的約束。例如,如果一個變數代表要認購的股票數量,則只應取整數值。同樣,如果一個變數代表發電機的開/關狀態,則只應取二進位制值(0 或 1)。

二次規劃:目標函式或約束函式為多元二次函式。此優化應用於財務金融中投資組合優化、發電廠發電優化、工程中設計優化等領域。

最小二乘:分為線性和非線性,通過最小化誤差的平方和尋找變數的最優函式匹配。非線性最小二乘優化還可用於曲線擬合。

對 MATLAB 提供的各類優化問題的演算法,我們稱之為求解器(Solver)。根據其求解目標,被分為四大組:

  • 極小值優化組:找到目標函數出發點 x0 附近的區域性極小值

  • 多目標優化組:找到最小化一組函式的最大值或指定的值

  • 方程求解組:找到非線性方程 f(x) = 0 出發點 x0 附近的解

  • 最小二乘法(曲線擬合)組:最小化平方和

僅優化工具箱就提供了近 20 種求解器,面對如此繁多的選項,使用者往往一頭霧水。幸好,MATLAB 提供了簡單明瞭的參考工具 —— 優化決策表。可謂一表在手,優化不愁:

上表中*表示演算法由全域性工具箱提供。

我的天呀,你是不是已經被上面的內容嚇傻了。看不太懂?沒關係。我們現在只需知道有那麼一回事即可,不必苦惱。等遇到相應的問題時我們再去結合具體的例子去深入學習和理解。

    2、已知目標函式和約束條件如下,試求解目標函式的最小值。

                             

       初一看,HPS、PP、EP、P1、I1,……等等,這些是什麼東東?別緊張,它們只是變數,把它們當成x,y,z這種常見的變數去看待即可。我數了一下,約束函式中共有19條數學表示式,涉及16個變數。是挺難搞的,用高中的作圖法怕是解決不了。別怕,我們有Matlab,通過程式碼去搞定它。

       好,現在我們開始一步步去求解它。

a. 首先,根據題目確認這是一個線性規劃問題。而線性規劃的通用數學表示式和MATLAB標準形式為:

這個標準形式很重要,上面的A,b,Aeq,beq,lb,ub我們後面要用到。

b. 對於線性規劃的優化求解步驟(也適用於其他優化方案),建議如下:

    1 ) 選擇優化求解器

    2 ) 將所有變數合併為一個向量

    3 ) 建立邊界約束(lb,ub)

    4 ) 建立線性不等式約束(A,b)

    5 ) 建立線性等式約束(Aeq,beq)

    6 ) 建立目標函式

    7 ) 優化問題求解

    8 ) 結果檢驗

上面的求解步驟非常重要,我們的程式碼整體框架跟求解步驟差不多。

現在我們按照上面的求解步驟去求解:

(1)選擇優化求解器。

這道題目是這是一個線性規劃問題。求解線性規劃問題,我們一般選用linprog。linprog在寫程式碼將要寫完才需用到。在第一步我們只需要知道一個線性規劃問題,程式碼按照求解線性規劃問題去寫即可。

(2)將所有變數合併為一個向量。

目標函式和約束條件中共有HPS、PP、EP、P1、I1等16個變數,我們需要將其合併為一個向量。除了合併為一個向量外,我們還將每個變數依次分別賦值1,2,3,4,……16。

%% 將所有變數合併為一個向量,共16個變數
variables = {'I1','I2','HE1','HE2','LE1','LE2','C','BF1','BF2','HPS','MPS','LPS','P1','P2','PP','EP'}
N = length(variables)
for v = 1:N
    eval([variables{v},'=',num2str(v),';'])
end

(3)建立邊界約束(lb,ub)

lb是指low boundary,即最低邊界。ub是指up boundary,即最高邊界。在這一步中,我們需要把約束條件裡的不等式(該不等式只含單個變數,如2500<=P1<=6250,3000<=P2<=9000)找出來,並根據它們的上下邊界寫程式碼。

%% 設定上下限約束(lb<=x<=ub)
% 設定下限約束,即lb<=x
lb = zeros(size(variables)); % 1x16的矩陣,每個數都是0
lb([P1,P2,MPS,LPS]) = [2500,3000,271536,100623];
% 設定上限約束,即x<=ub
ub = Inf(size(variables)); % 1想6的矩陣,每個數都是無窮大
ub([P1,P2,I1,I2,C,LE2])= [6250,9000,192000,244000,62000,142000];

(4)建立線性不等式約束(A,b)。

在這一步需要將約束條件裡的不等式(該不等式含兩個或兩個以上變數)找出來,並根據它們的上下邊界寫程式碼。

%% 建立線性不等式約束(A*x<=b)
A = zeros(3,16); % 3x16的矩陣,每個數均為0,為什麼是3x16,因為約束條件有3個不等式
% 由不等式I1-HE1<=132000得到下面一行程式碼
A(1,I1)=1; A(1,HE1)=-1; b(1) = 132000;
% 由不等式-EP-PP<=12000得到下面一行程式碼
A(2,EP)=-1; A(2,PP)=-1; b(2) = -12000;
% 由不等式-P1-P2-PP<=-24550得到下面一行程式碼
A(3,[P1,P2,PP])=[-1,-1,-1];b(3)=-24550;

(5)建立線性等式約束(Aeq,beq)。

在這一步需要把約束條件中的等式找出來,並通過移項,讓等式的右邊為0。

%% 建立線性等式約束(Aeq*x=beq)
Aeq=zeros(8,16);beq=zeros(8,1) % 約束條件中共有8個等式
% 把等式I2=LE2+HE2轉化為LE2+HE2-I2=0後,得到下面一行程式碼
Aeq(1,[LE2,HE2,I2])=[1,1,-1];
Aeq(2,[LE1,LE2,BF2,LPS])=[1,1,1,-1];
Aeq(3,[I1,I2,BF1,HPS])=[1,1,1,-1];
Aeq(4,[C,MPS,LPS,HPS])=[1,1,1,-1];
Aeq(5,[LE1,HE1,C,I1])=[1,1,1,-1];
Aeq(6,[HE1,HE2,BF1,BF2,MPS])=[1,1,1,-1,-1];
Aeq(7,[HE1,LE1,C,P1,I1])=[1267.8,1251.4,192,3413,-1359.8];
Aeq(8,[HE2,LE2,P2,I2])=[1267.8,1251.4,3413,-1359.8];

(6)建立目標函式。

%% 建立目標函式
f = zeros(size(variables));
% 由目標函式0.002614HPS+0.0239PP+0.009825EP
f([HPS,PP,EP]) = [0.002614,0.0239,0.009825];

(7) 求解問題

%% 由linprog實現線性規劃問題求解
options = optimoptions('linprog','Algorithm','dual-simplex');
% 將前面已經確定的各個引數傳入linprog()中
[x, fval] = linprog(f,A,b,Aeq,beq,lb,ub,options);
for d=1:N
    fprintf('%12.2f\t%s\n',x(d),variables{d})
end

fprintf函式是將求解後每個變數的打印出來。

求解結果如下:

下面把完整的原始碼貼上:

clc,clear,close all
%% 選擇優化求解器,線性規劃求解可由linprog實現

%% 將所有變數合併為一個向量,共16個變數
variables = {'I1','I2','HE1','HE2','LE1','LE2','C','BF1','BF2','HPS','MPS','LPS','P1','P2','PP','EP'}
N = length(variables)
for v = 1:N
    eval([variables{v},'=',num2str(v),';'])
end

%% 設定上下限約束(lb<=x<=ub)
% 設定下限約束,即lb<=x
lb = zeros(size(variables)); % 1x16的矩陣,每個數都是0
lb([P1,P2,MPS,LPS]) = [2500,3000,271536,100623];
% 設定上限約束,即x<=ub
ub = Inf(size(variables)); % 1想6的矩陣,每個數都是無窮大
ub([P1,P2,I1,I2,C,LE2])= [6250,9000,192000,244000,62000,142000];

%% 建立線性不等式約束(A*x<=b)
A = zeros(3,16); % 3x16的矩陣,每個數均為0,為什麼是3x16,因為約束條件有3個不等式
% 由不等式I1-HE1<=132000得到下面一行程式碼
A(1,I1)=1; A(1,HE1)=-1; b(1) = 132000;
% 由不等式-EP-PP<=12000得到下面一行程式碼
A(2,EP)=-1; A(2,PP)=-1; b(2) = -12000;
% 由不等式-P1-P2-PP<=-24550得到下面一行程式碼
A(3,[P1,P2,PP])=[-1,-1,-1];b(3)=-24550;

%% 建立線性等式約束(Aeq*x=beq)
Aeq=zeros(8,16);beq=zeros(8,1) % 約束條件中共有8個等式
% 把等式I2=LE2+HE2轉化為LE2+HE2-I2=0後,得到下面一行程式碼
Aeq(1,[LE2,HE2,I2])=[1,1,-1];
Aeq(2,[LE1,LE2,BF2,LPS])=[1,1,1,-1];
Aeq(3,[I1,I2,BF1,HPS])=[1,1,1,-1];
Aeq(4,[C,MPS,LPS,HPS])=[1,1,1,-1];
Aeq(5,[LE1,HE1,C,I1])=[1,1,1,-1];
Aeq(6,[HE1,HE2,BF1,BF2,MPS])=[1,1,1,-1,-1];
Aeq(7,[HE1,LE1,C,P1,I1])=[1267.8,1251.4,192,3413,-1359.8];
Aeq(8,[HE2,LE2,P2,I2])=[1267.8,1251.4,3413,-1359.8];

%% 建立目標函式
f = zeros(size(variables));
% 由目標函式0.002614HPS+0.0239PP+0.009825EP
f([HPS,PP,EP]) = [0.002614,0.0239,0.009825];

%% 由linprog實現線性規劃問題求解
options = optimoptions('linprog','Algorithm','dual-simplex');
% 將前面已經確定的各個引數傳入linprog()中
[x, fval] = linprog(f,A,b,Aeq,beq,lb,ub,options);
for d=1:N
    fprintf('%12.2f\t%s\n',x(d),variables{d})
end

  3、下面是一個整數規劃問題,已知目標函式和約束條件如下,求解目標函式的最大值。

                                      

     求解最大值問題和求解最小值問題本質上是一致的,求解最大值也可以轉換為求解最小值。

例如:本題要求解z=3*x1-2*x2+5*x3的最大值,也就是要求解y=-3*x1+2*x2-5*x3的最小值。求解線性規劃最小值問題我們在上面已經學過。本題還有一個比較特殊的問題是,約束條件中的三個變數均為整數,而且是0或1,這也是所謂的0-1規劃問題。

求解整值問題要用到專門的求解器 intlinprog。

clc,clear,close all
% 求z=3*x1 - 2*x2 + 5*x3的最大值轉化為求y=-3*x1 + 2*x2 - 5*x3的最小值。
f = [-3;2;-5]; % 建立目標函式
intcon=[1,2,3];
A=[1 2 -1; 1 4 1; 1 1 0; 0 4 1]; % 四個不等式中的變數係數
b=[2;4;3;6]; % 約束條件中不等式右邊的常數
lb=[0,0,0]; % x1,x2,x3=0
ub=[1,1,1]; % x1,x2,x3=1
Aeq=[0,0,0];
beq=0;
x=intlinprog(f,intcon,A,b,Aeq,beq,lb,ub)

我想提一下intcon=[1,2,3]這句程式碼是怎麼回事。

下面我舉一個簡單的例子來說明intcon的用法。

X=[x1,x2,x3,x4,x5,x6],其中x2, x3, x6只能取整數
intcon = [2,3,6]
如果所有變數都只能取整數,則:intcon = [1,2,3,4,5,6]; 比較方便的寫法是:intcon = 1:6
如果只有x4取整數,則:intcon = 4;就是約束整形變數

在本題中,除了x1,x2,x3=0或1外,沒有其它的等式了。故Aeq=[0,0,0];beq=0;

 4、圖形化應用

在數學建模競賽中,我們一般通過程式碼來進行求解問題,圖形化應用可以用來檢驗結果是否正確。

MATLAB 在資料分析領域如此受歡迎,除了其提供豐富的內建演算法集,還有各類友好的應用介面。在優化工具箱中,也有這麼一個強大的工具—— Optimization App,可以在 MATLAB Apps 視窗或者執行 optmitool 命令開啟。它是一個互動式的圖形化應用工具,無需手寫程式碼,直接在圖形介面中設定各類求解器、配置目標函式、約束條件,即可執行優化演算法並使中間結果和最終結果視覺化。

在 Optimization App 中,只需點選選單欄中的 File > Generate Code,即可將 App 中的各項設定自動生成 MATLAB 程式碼,使用者可實現演算法的複用和二次開發。