1 整數規劃的MATLAB求解方法（一） 用MATLAB求解一般混合整數規劃問題由于MATLAB優化工具箱中并未提供求解純整數規劃和混合整數規劃的函數，因而需要自行根據需要和設定相關的算法來實現。現在有許多用戶發布的工具箱可以解決該類問題。這里我們給出開羅大學的Sherif和Tawfik在MATLAB Central上發布的一個用于求解一般混合整數規劃的程序，在此命名為intprog，在原程序的基礎上做了簡單的修改，將其選擇分枝變量的算法由自然序改造成分枝變量選擇原則中的一種，即：選擇與整數值相差最大的非整數變量首先進行分枝。intprog函數的調用格式如下： x,fval,exitflag=intprog(c,A,b,Aeq,beq,lb,ub,M,TolXInteger) 該函數解決的整數規劃問題為：在上述標準問題中，假設為維設計變量，且問題具有不等式約束個，等式約束個，那么：、均為維列向量，為維列向量，為維列向量，為維矩陣，為維矩陣。在該函數中，輸入參數有c,A,b,Aeq,beq,lb,ub,M和TolXInteger。其中c為目標函數所對應設計變量的系數，A為不等式約束條件方程組構成的系數矩陣，b為不等式約束條件方程組右邊的值構成的向量。Aeq為等式約束方程組構成的系數矩陣，beq為等式約束條件方程組右邊的值構成的向量。lb和ub為設計變量對應的上界和下界。M為具有整數約束條件限制的設計變量的序號，例如問題中設計變量為，要求和為整數，則M=2;3;6；若要求全為整數，則M=1:6，或者M=1;2;3;4;5;6。TolXInteger為判定整數的誤差限，即若某數x和最鄰近整數相差小于該誤差限，則認為x即為該整數。在該函數中，輸出參數有x, fval和exitflag。其中x為整數規劃問題的最優解向量，fval為整數規劃問題的目標函數在最優解向量x處的函數值，exitflag為函數計算終止時的狀態指示變量。例1 求解整數規劃問題：算法：c=-1;-1;A=-4 2;4 2;0 -2;b=-1;11;-1;lb=0;0;M=1;2; %均要求為整數變量Tol=1e-8; %判斷是否整數的誤差限x,fval=linprog(c,A,b,lb,) %求解原問題松弛線性規劃x1,fval1=intprog(c,A,b,lb,M,Tol) %求最優解整數解結果：x =%松弛線性規劃問題的最優解 1.5000 2.5000fval = -4.0000x1 =%整數規劃的最優解 2 1fval2 = -3.0000（二） 用MATLAB求解0-1規劃問題在MATLAB優化工具箱中，提供了專門用于求解0-1規劃問題的函數bintprog，其算法基礎即為分枝界定法，在MATLAB中調用bintprog函數求解0-1規劃時，需要遵循MATLAB中對0-1規劃標準性的要求。 0-1規劃問題的MATLAB標準型在上述模型中，目標函數f 需要極小化，以及需要滿足的約束條件，不等式約束一定要化為形式為“”。假設為維設計變量，且問題具有不等式約束個，等式約束個，那么：、均為維列向量，為維列向量，為維列向量，為維矩陣，為維矩陣。如果不滿足標準型的要求，則需要對原問題進行轉化，化為標準型之后才能使用相關函數，標準化的方法和線性規劃中的類似。0-1規劃問題的MATLAB求解函數 MATLAB優化工具箱中求解0-1規劃問題的命令為bintprog bintprog的調用格式x = bintprog(f)x = bintprog(f,A,b)x = bintprog(f,A,b,Aeq,beq)x = bintprog(f,A,b,Aeq,beq,x0)x = bintprog(f,A,b,Aeq,Beq,x0,options)x,fval = bintprog(.)x,fval,exitflag = bintprog(.)x,fval,exitflag,output = bintprog(.)命令詳解1）x = bintprog(f) 該函數調用格式求解如下形式的0-1規劃問題 2）x = bintprog(c,A,b) 該函數調用格式求解如下形式的0-1規劃問題 3）x = bintprog (c,A,b,Aeq,beq)該函數調用格式求解如下形式的0-1規劃問題：4）x = bintprog (c,A,b,Aeq,beq,x0) 該函數調用格式求解如下形式的0-1規劃問題在前一個調用格式的基礎上同時設置求解算法的初始解為x0，如果初始解x0不在0-1規劃問題的可行域中，算法將采用默認的初始解5）x = bintprog (c,A,b,Aeq,beq,x0,options)用options指定的優化參數進行最小化。可以使用optimset來設置這些參數 上面的函數調用格式僅設置了最優解這一輸出參數，如果需要更多的輸出參數，則可以參照下面的調用格式：x,fval = bintprog(.) 在優化計算結束之時返回整數規劃問題在解x處的目標函數值fval x,fval,exitflag = bintprog(.) 在優化計算結束之時返回exitflag值，描述函數計算的退出條件。 x,fval,exitflag,output = bintprog(.)在優化計算結束之時返回結構變量output 例2：求解0-1規劃問題為了程序的可讀性，我們用一維下標來表示設計變量，即用表示，用表示，用表示，用表示，于是約束條件和目標函數分別為： 算法：c=20;12;33;26;22;15;29;23;21;13;31;24;22;16;32;23;Aeq=1 1 1 1 0 0 0 0 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 0 0 0 0 1 1 1 1 0 0 0 0; 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 1 0 0 0; 0 1 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 1 0; 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1; ;beq=ones(1,8);x,fval=bintprog(c,Aeq,beq);B=reshape(x,4,4); %由于x是一列元素，為了使結果更加直觀，故排成與效率矩陣E相對應的形式Bfval結果：Optimization terminated.ans = 0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1fval = 85整數規劃的應用組件配套問題某機械產品需要由三個工廠開工一起生產后組裝完成。每件機械需要4個組件1和3個組件2。生產這兩種組件需要消耗兩種原材料A和B。已知這兩種原材料的供應量分別為400kg和600kg。由于三個工廠的生產條件和擁有設備工藝條件不同，每個工廠生產組件的能力和原材料的消耗也不盡相同，且每個工廠開工一次都是配套生產一定數量的組件1和組件2，其具體的數據如表所示。表11-2 各工廠生產能力和消耗原材料的數據表每個工廠消耗原材料的數量（公斤）每個工廠各組件的生產能力（件數）A材料 B材料組件1組件2工廠工廠工廠9647109879695現在的最優化問題是，這三個工廠應當如何安排生產，才能使該產品的配套數達到最大？（）組件配套問題的建模設和是三個工廠分別開工的次數，將其作為該問題的設計變量。由于每個工廠開工一次都是配套生產，故每次開工消耗的原材料一定，且生產的組件數也是一定的。在這個假設的前提之下，我們可以得出該問題的目標函數和約束條件。因為原材料的總量是有限的，根據工廠的開工次數，可得工廠將消耗材料，工廠將消耗材料，工廠將消耗材料，故有約束條件：同理，對于材料的總量約束條件可以表達為：然后再來分析三個工廠零件生產的情況，三個工廠生產的組件1的數量分別為和，故組件1的總數為：同理，組件2的總數為：下一步是分析目標函數，該問題要求的不是生產的各種組件總數最多，而是要求產品的配套數量最大，即能組成的機械數目最多。問題中已經給出了該種機械中兩種組件的配比為4:3，故組件1所能成套的數目滿足約束條件：同理，組件2所能成套的數目滿足約束條件：因而，所能組成的成品機械的數目應該為和中的較小值，即：那么，我們求解的目標即是使得能夠盡可能大，可以看出，這種形式并不是有關設計變量的線性函數，我們需要對其進行轉化，此時我們可以令一個人工設計變量等于目標函數值，則有：在該假設下，一定滿足關系式：且結合約束關系，對的約束可以表示為：同理，對的約束可以表示為：對的上述關系進行整理，可以得到關系：同理對也可以得到不等式關系為：上面兩個式子即為由組件的配比數得到的關于組件數目的約束條件，此時問題的目標就是如何使得取到最大值，由于開工的次數一定是一個非負整數，故也是一個約束條件，決定了該問題是一個純整數規劃問題。結合前面給出的原材料約束，可以得到如下的數學模型：（）組件配套問題的求解利用8.節中給出的函數對此問題求解，代碼和運行結果如下：算法：%目標函數所對應的設計變量的系數，為轉化為極小，故取原目標函數的相反數f=0;0;0;-1;%不等式約束A= 9 6 4 0; 7 10 9 0; -8 -7 -9 4; -6 -9 -5 3;b=400;600;0;0;%邊界約束，由于無上界，故設置ub=Inf;Inf;Inf;Inf;lb=0;0;0;0;ub=Inf;Inf;Inf;Inf;%所有變量均為整數變量，故將所有序號組成向量MM=1;2;3;4;%判定為整數的誤差限Tol=1e-8;%求最優解x和目標函數值fval，并返回狀態指示x,fval,exitflag=intprog(f,A,b,lb,ub,M,Tol)結果：x =最優解向量x 18 15 36 141fval = 在最優解向量x處，原目標函數值的相反數 -141.000exitflag= 算法終止指示變量，說明問題收斂到了最優解x 1由運行結果可知，工廠、和需要分別開工18、15和36次，這樣所能生產出來的成套的機械為141件。2 多目標規劃的MATLAB求解方法（一） 多目標規劃的MATLAB求解由于多目標規劃中的求解涉及到的方法非常多，故在MATLAB中可以利用不同的函數進行求解，例如在評價函數法中我們所得最后的評價函數為一線性函數，且約束條件也為線性函數，則我們可以利用MATLAB優化工具箱中提供的linprog函數進行求解，如果我們得到的評價函數為非線性函數，則可以利用MATLAB優化工具箱中提供的fmincon函數進行求解，如果我們采用最大最小法進行求解，則可以利用MATLAB優化工具箱中提供的fminimax函數進行求解。下面我們就結合理論求解的幾種方法，講解一下典型多目標規劃問題的MATLAB求解方法。例1 利用理想點法求解我們首先根據評價函數法中的理想點法的理論基礎，按照理想點法的求解思路分別對兩個單目標規劃問題進行求解：求解的MATLAB的代碼和相應的運行結果為：算法：c=2;-3;A=3 2;1 1;b=12;8lb=0;0x,fval=linprog(c,A,b,lb,)結果：x = 0.0000 6.0000fval = -18.0000于是可知。當時，單目標線性規劃的最優函數值為。求解的MATLAB的代碼和相應的運行結果為：算法：c=-5;-3;A=3 2;1 1;b=12;8lb=0;0x,fval=linprog(c,A,b,lb,)結果：Optimization terminated.x = 4.0000 0.0000fval =-20.0000于是可知，當時，單目標線性規劃的最優函數值為。由上述兩個單目標線性規劃的求解結果可知，因而是一個不可能達到的理想點，因而我們構造如下評價函數：構造描述該評價函數的M-函數文件objfun.m如下：function f=objfun(x)f=sqrt(2*x(1)-3*x(2)+18)2+(5*x(1)+3*x(2)-20)2);然后用非線性規劃的方式求解上述問題：算法：A=3 2;1 1;b=12;8;lb=0;0;x0=0;0;x,fval,exitflag=fmincon(objfun,x0,A,b,lb,)結果： Active inequalities (to within options.TolCon = 1e-006): lower upper ineqlin ineqnonlin 1 x = 0.0235 5.9647fval = 1.9941exitflag = 5由運行結果可知在該評價函數標準之下，問題的最優解為：此時，各目標函數的取值為：。它與理想點在評價函數標準下的最小距離為1.9941。例2 利用評價函數中的線性加權和法求解如下多目標規劃問題：其中權系數為。建立線性加權和法的評價函數為：將相應的權系數代入上式即整理出目標函數為：于是建立目標函數的M-函數文件objfun.m：function f=objfun(x)f=x(1)2+1.2*x(2)2+1.4*x(3)2;由于目標函數非線性函數且具有線性等式約束和邊界約束，因而我們調用MATLAB中求解非線性規劃的fmincon函數對此問題進行求解，同時注意如果只考慮第一個目標函數，由這種特殊形式，即在設計變量的和為一定值，需要求其平方和的最小值時，最優解必然是當這幾個設計變量的值相等時取得，于是我們可以將這個解設為問題的初始點，開始迭代。算法：Aeq=1 1 1;beq=3;lb=0;0;0;x0=1;1;1;x,fval=fmincon(objfun,x0,Aeq,beq,lb,)結果：No active inequalities.x = 1.1776 0.9812 0.8412fval = 3.5327結果說明，問題的最優解為：我們在求解多目標規劃問題時，可以采用評價函數法中的最大最小法，而MATLAB也為這種方法提供了專門的求解函數fminimax，在講解這方面的例題之前，我們首先介紹一下MATLAB優化工具箱中所提供的最大最小法的求解函數fminimax。最大最小法問題的MATLAB標準形式為：函數fminimax的調用方式和其他的最優化函數類似，其中所涉及的輸入參數和輸出參數的含義與非線性規劃的求解函數fmincon類似，使用方法也基本相同，細節問題讀者可以參考MATLAB的幫助文件。例3 求解最大最小問題：首先建立描述目標函數的M-函數文件objfun.m，注意到一共有五個目標函數，所求目標為這五個函數最大值中的最小值，代碼如下：function f = objfun(x)f(1)= 2*x(1)2+x(2)2-48*x(1)-40*x(2)+304; f(2)= -x(1)2 - 3*x(2)2;f(3)= x(1) + 3*x(2) -18;f(4)= -x(1)- x(2);f(5)= x(1) + x(2) - 8;然后設置求解的初始點為x0=0;0，調用fminimax求解該問題。算法：x0 = 0; 0;x,fval,maxfval = fminimax(objfun,x0)結果：Local minimum possible. Constraints satisfied.fminimax stopped because the predicted change in the objective functionis less than the default value of the function tolerance and constraints were satisfied to within the default value of the constraint tolerance.Active inequalities (to within options.TolCon = 1e-006): lower upper ineqlin ineqnonlin 1 5x = 4.0000 4.0000fval = 0.0000 -64.0006 -1.9999 -8.0000 -0.0000maxfval = 2.6897e-008上述結果說明當時，目標函數的最大值達到最小，這一組的函數值為0.0000，-64.0006，-1.9999，-8.0000，-0.0000，于是最大值為0。多目標規劃的應用投資問題（全國大學生數學建模競賽試題）假設市場上有種資產，比如股票、債券等可以供投資者選擇，某公司有數額為的一筆相當大的資金可用作一個時間的投資。通過財務人員對種資產進行評估，大概可以估算出在這一時期內購買資產的平均收益率為，并預測出購買的損失率為。考慮到投資越分散，總的風險越小，公司決定當用這筆資金購買若干種資產時，總體風險可用所投資的中的最大一個風險來度量。購買要付交易費，費率為，并且當購買額不超過給定值時，交易費按購買計算（不買當然無須付費）。另外，假定同期銀行存款利率是，且既無交易費又無風險(5)。已知時的相關數據如下表所示：表1 投資各種資產的參數值（元）282.51103211.52198235.54.552252.66.540試給該公司設計一種投資組合方案，即用給定的資金，有選擇地購買若干種資產或存銀行生息，使凈收益盡可能大，而總體風險盡可能小。（）投資問題的建模 為了建立數學模型，首先對模型進行一些必要的假設及符號說明：(1) 模型的假設 在一個時期內所給出的，保持不變。在一個時間內所購買的各種資產(如股票、證券等)不進行買賣交易，即在買入后不再賣出。 每種投資是否收益是相互獨立的。 在投資過程中，無論盈利與否必須先付交易費。(2)符號說明 (元)：公司現有投資總金額；：欲購買的第種資產種類（其中表示存入銀行）：公司購買的金額；：公司購買的平均收益率；：公司購買的平均損失率；：公司購買超過時所付交易費率。下面來建立模型。設購買的金額為，所付的交易費，則由于投資額相當大，所以總可以假定對每個的投資。這時上面的大括號公式可簡化為：對投資的凈收益為：對的風險為：對投資所需資金為投資金額與所需的手續費之和，即：當購買的金額為，投資組合的凈收益總額為：整體風險為：資金約束為：根據題目要求，以凈收益總額最大，而整體風險最小為目標建立模型如下：很顯然，這是一個多目標規劃問題。（）投資問題的求解在此我們采用主要目標法對該問題進行求解，即根據問題的實際情況，確定一個目標為主要目標，而把其余目標作為次要目標，并且根據經驗，選取一定的界限值。這樣就可以把次要目標作為約束來處理，于是就將原來的多目標問題轉化為一個在新的約束下的單目標最優化問題。在上述例子中如果以收益為主要目標，則可以固定風險水平，給定風險一個界限，講問題轉化稱為求最大風險不超過時的最大收益，即下面的線性規劃模型： （1）若投資者希望總盈利至少達到水平以上，則可以在風險最小的情況下尋找相應的投資組合，從而將原模型轉化成為下列的線性規劃模型進行求解： （2）根據上面的分析，我們利用主