第九章最優(yōu)化方法

1、第九章 最優(yōu)化方法 本章主要介紹線性規(guī)劃、0-1規(guī)劃、非線性規(guī)劃等問(wèn)題的MATLAB求解。 9.1 線性規(guī)劃（Linear Programming，簡(jiǎn)寫(xiě)為L(zhǎng)P）問(wèn)題線性規(guī)劃問(wèn)題就是求多變量線性函數(shù)在線性約束條件下的最優(yōu)值。滿足約束條件的解稱為可行解，所有可行解構(gòu)成的集合稱為可行域，滿足目標(biāo)式的可行解稱為最優(yōu)解。MATLAB解決的線性規(guī)劃問(wèn)題的標(biāo)準(zhǔn)形式為：min 其中為列向量，為矩陣。其它形式的線性規(guī)劃問(wèn)題都可經(jīng)過(guò)適當(dāng)變換化為此標(biāo)準(zhǔn)形式。在MATLAB中求解線性規(guī)劃問(wèn)題函數(shù)為linprog，其使用格式為：x, fval, exitflag, output, lambda = linprog(f,

2、 A, b, Aeq, beq, lb, ub)輸入部分：其中各符號(hào)對(duì)應(yīng)線性規(guī)劃問(wèn)題標(biāo)準(zhǔn)形式中的向量和矩陣，如果約束條件中有缺少，則其相應(yīng)位置用空矩陣代替。輸出部分：其中x為最優(yōu)解，用列向量表示；fval為最優(yōu)值；exitflag為退出標(biāo)志，若exitflag=1表示函數(shù)有最優(yōu)解，若exitflag=0表示超過(guò)設(shè)定的迭代最大次數(shù)，若exitflag=-2，表示約束區(qū)域不可行，若exitflag=-3，表示問(wèn)題無(wú)解，若exitflag=-4，表示執(zhí)行迭代算法時(shí)遇到NaN，若exitflag=-5，表示原問(wèn)題和對(duì)偶問(wèn)題均不可行，若exitflag=-7，表示搜索方向太小，不能繼續(xù)前進(jìn)；output

3、表明算法和迭代情況；lambda表示存儲(chǔ)情況。例1 用linprog函數(shù)求下面的線性規(guī)劃問(wèn)題輸入如下：>> f = -5, -4,-6;A = 1 -1 1; 3 2 4; 3 2 0;b = 20; 42; 30;lb = zeros(3,1);x,fval,exitflag,output,lambda = linprog(f,A,b,lb)注意：由于該問(wèn)題沒(méi)有等式約束，所以輸入格式中相應(yīng)的位置用代替，變量沒(méi)有上限約束，所以u(píng)b也用代替，但由于其在最后，可以不寫(xiě)。輸出結(jié)果如下：Optimization terminated.x = % 最優(yōu)解 0.0000 15.0000 3.0

4、000fval = %最優(yōu)值 -78.0000exitflag = %函數(shù)收斂于解 1output = iterations: 6 algorithm: 'large-scale: interior point' cgiterations: 0 message: 'Optimization terminated.'lambda = ineqlin: 3x1 double eqlin: 0x1 double upper: 3x1 double lower: 3x1 double例2 一家家具公司生產(chǎn)桌子和椅子，用于生產(chǎn)的全部勞動(dòng)力共計(jì)450個(gè)工時(shí)，原料是400個(gè)單

5、位的木材。每張桌子使用15個(gè)工時(shí)的勞動(dòng)力，20個(gè)單位的木材，售價(jià)為80元。每張椅子使用10個(gè)工時(shí)，用料5個(gè)單位，售價(jià)45元。問(wèn)為達(dá)到最大效益，應(yīng)如何安排生產(chǎn)？解 設(shè)生產(chǎn)桌子個(gè)，椅子個(gè)，建立如下模型：輸入如下：>> f = -80,-45;A = 20, 5; 15, 10;b = 400, 450;lb = zeros(2,1);x, fval, exitflag = linprog(f,A,b,lb)結(jié)果如下：Optimization terminated.x = 14.0000 24.0000fval = -2.2000e+003exitflag = 1注意：由于linprog

6、是求目標(biāo)函數(shù)的最小值，如求目標(biāo)函數(shù)的最大值，可先求出的最小值fval，則-fval就是的最大值。本例只輸出最優(yōu)解、最優(yōu)值和退出標(biāo)志，可見(jiàn)生產(chǎn)14個(gè)桌子，24個(gè)椅子，可獲得最大利潤(rùn)2200元。9.2 0-1規(guī)劃0-1規(guī)劃是一種特殊形式的整數(shù)規(guī)劃，它的決策變量?jī)H取值0或1.一般用0表示放棄，1表示選取，故0-1規(guī)劃可以用來(lái)處理選址問(wèn)題、指派問(wèn)題、裝箱問(wèn)題、項(xiàng)目評(píng)價(jià)、資金分配、生產(chǎn)計(jì)劃安排等問(wèn)題。在MATALB中求解0-1規(guī)劃問(wèn)題函數(shù)為bintprog，其針對(duì)下述0-1規(guī)劃：min 其中為列向量，為矩陣。使用格式為：x, fval, exitflag, output = bintprog(f, A,

7、 b, Aeq, beq)輸入部分：其中各符號(hào)對(duì)應(yīng)0-1規(guī)劃問(wèn)題標(biāo)準(zhǔn)形式中的向量和矩陣，如果約束條件中有缺少，則其相應(yīng)位置用空矩陣代替。輸出部分：其中x為最優(yōu)解，用列向量表示；fval為最優(yōu)值；exitflag為退出標(biāo)志，若exitflag=1表示函數(shù)有最優(yōu)解，若exitflag=-1表示超過(guò)設(shè)定的迭代最大次數(shù)，若exitflag=-2，表示問(wèn)題不可行，若exitflag=-4，表示搜索節(jié)點(diǎn)數(shù)超過(guò)了設(shè)定的搜索節(jié)點(diǎn)最大個(gè)數(shù)，若exitflag=-5，表示搜索時(shí)間超過(guò)了設(shè)定的指令運(yùn)行的最大秒數(shù)，若exitflag=-6，表示LP求解器在某節(jié)點(diǎn)處求解LP松弛問(wèn)題時(shí)的迭代次數(shù)超過(guò)了設(shè)定的迭代次數(shù)；ou

8、tput包含使用算法、迭代次數(shù)、搜索過(guò)的節(jié)點(diǎn)數(shù)、算法執(zhí)行時(shí)間、算法終止原因。例3 求解下述0-1規(guī)劃問(wèn)題。利用bintprog函數(shù)求解如下：>> f=-1;2;2;-6;-4;A=3,2,-1,4,2; 2,4, -2,-1,-2;b=5;5 ;x,fval,exit,output=bintprog(f,A,b)輸出結(jié)果：Optimization terminated.x = 1 1 1 0 0fval = -5exit = 1output = iterations: 4 nodes: 1 time: 0.0781 algorithm: 'LP-based branch-a

9、nd-bound' branchStrategy: 'maximum integer infeasibility' nodeSrchStrategy: 'best node search' message: 'Optimization terminated.這表明，當(dāng)時(shí)，目標(biāo)函數(shù)取最大值5。例4 資金分配問(wèn)題，某企業(yè)在今后3年內(nèi)有5個(gè)工程考慮開(kāi)工，每項(xiàng)工程的期望收入和年度費(fèi)用如表所示。假定每一項(xiàng)已經(jīng)批準(zhǔn)的工程要在整3年內(nèi)完成。企業(yè)應(yīng)如何選擇工程，使企業(yè)總收入最大？工程費(fèi)用（千元）收入（千元）第一年第二年第三年1518202471040339220

10、4741155861030最大可用資金數(shù)252525解 設(shè)決策變量為：其中，表示放棄第項(xiàng)工程，表示選擇第項(xiàng)工程。該問(wèn)題的數(shù)學(xué)模型為：利用bintprog函數(shù)求解如下：f=-20;40;20;15;30;A=5,4,3,7,8;1,7,9,4,6;8,10,2,1,10;b=25;25;25;x,fv,exit=bintprog(f,A,b)輸出結(jié)果如下：Optimization terminated.x = 1 1 1 1 0fv = -95exit = 1上述結(jié)果表明，該企業(yè)選擇第1，第2，第3，第4項(xiàng)工程，可以獲得最大利潤(rùn)95千元。指派問(wèn)題：設(shè)有項(xiàng)工作分配給個(gè)人去做，每人做一項(xiàng)，由于個(gè)人的

11、工作效率不同，因而完成同一工作所需時(shí)間也不同，設(shè)人員完成工作所需時(shí)間為（稱為效率矩陣），問(wèn)如何分配工作，使得完成所有工作所用的總時(shí)間最少？這類問(wèn)題稱為指派問(wèn)題（Assignment Problem），也稱最優(yōu)匹配問(wèn)題，它是一類典型的0-1規(guī)劃問(wèn)題。例5 有4個(gè)工人被分派做4項(xiàng)工作，每項(xiàng)工作只能一人做，每人只能做一項(xiàng)工作。現(xiàn)設(shè)每個(gè)工人做每項(xiàng)工作的耗時(shí)如表所示，求總耗時(shí)最少的分配方案。表1時(shí)間表 時(shí)間單位：小時(shí) 工作耗時(shí)工人1234115182124219232218326171619419212317解 設(shè)決策變量為，只取0或1。表示工人完成工作；表示工人不做工作。于是，上述問(wèn)題的數(shù)學(xué)模型如下：

12、下面給出Matlab計(jì)算程序。>> f=15;18;21;24;19;23;22;18;26;17;16;19;19;21;23;17;o4=ones(1,4);z4=zeros(1,4);z8=zeros(1,8);z12=zeros(1,12);Aeq1=o4,z12;z4,o4,z8;z8,o4,z4;z12,o4;e4=eye(4);Aeq2=e4,e4,e4,e4;Aeq = Aeq1;Aeq2;beq=ones(8,1);x,fv,exit=bintprog(f,Aeq,beq);x=transpose(reshape(x,4,4)計(jì)算接可以得出：fv = 70exit

13、 = 1指派方陣x = 1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0它表示的最優(yōu)指派方案為：工人1做工作1；工人2做工作4；工人3做工作3；工人4做工作2。該指派方案耗時(shí)最少，為70小時(shí)。9.3 非線性規(guī)劃9.3.1 有約束的一元函數(shù)的極（最）小值在MATLAB中使用fminbnd函數(shù)求一元函數(shù)在區(qū)間上的極（最）小值。其調(diào)用格式有：x,fval,exitflag,output = fminbnd(fun,x1,x2,options)其中輸入?yún)?shù)fun為目標(biāo)函數(shù)的表達(dá)式字符串或MATLAB定義的M函數(shù)，x1和x2為自變量x的取值范圍，options為指定優(yōu)化參數(shù)選項(xiàng)；輸出參數(shù)x

14、為函數(shù)fun在區(qū)間上取最小值時(shí)的x值，fval為目標(biāo)函數(shù)的最小值，exitflag為終止迭代的條件，output為優(yōu)化信息，參見(jiàn)optimset函數(shù)。若參數(shù)exitflag>0，表示函數(shù)收斂于x，若exitflag=0，表示超過(guò)函數(shù)估計(jì)值或迭代的最大數(shù)字，exitflag<0表示函數(shù)不收斂于x；若參數(shù)output=iterations表示迭代次數(shù)，output=funccount表示函數(shù)賦值次數(shù)，output=algorithm表示所使用的算法。例6 計(jì)算函數(shù)在區(qū)間(0,1)內(nèi)的最小值。在命令窗口輸入：>> x,fval,exitflag,output=fminbnd(

15、'(x3+cos(x)+x*log(x)/exp(x)',0,1)輸出結(jié)果為：x = 0.5223fval = 0.3974exitflag = 1output = iterations: 8 funcCount: 9 algorithm: 'golden section search, parabolic interpolation' message: 1x112 char表明該函數(shù)在x =0.5223時(shí)取最小值0.3974。ezplot('(x3+cos(x)+x*log(x)/exp(x)',0,1)為簡(jiǎn)化輸出結(jié)果，可輸入：>>

16、 x,fval=fminbnd('(x3+cos(x)+x*log(x)/exp(x)',0,1)9.3.2 無(wú)約束多元函數(shù)極（最）小值在MATLAB中使用fminsearch函數(shù)和fminunc函數(shù)求無(wú)約束多元函數(shù)的極（最）小值。（1）fminsearch函數(shù)該函數(shù)的調(diào)用格式為：x,fval,exitflag,output= fminsearch(fun,x0,options)其中輸入?yún)?shù)fun為目標(biāo)函數(shù)的表達(dá)式字符串或MATLAB自定義的M函數(shù)，x0為初始點(diǎn)，options為指定優(yōu)化參數(shù)選項(xiàng)；輸出參數(shù)x為多元函數(shù)的最小值點(diǎn)，fval為最小值，exitflag和output與

17、單變量情形一致。例7 求的最小值點(diǎn)。在命令窗口輸入：>> x,fval=fminsearch('2*x(1)3+4*x(1)*x(2)3-10*x(1)*x(2)+x(2)2', 0,0)輸出結(jié)果為：x = 1.0016 0.8335fval = -3.3241表明當(dāng)x1=1.0016，x2=0.8335時(shí)函數(shù)的最小值為-3.3241。（2）fminunc函數(shù)該函數(shù)的調(diào)用格式為：x,fval,exitflag,output,grad,hessian = fminunc (fun,x0,options)其中輸入?yún)?shù)fun為目標(biāo)函數(shù)的表達(dá)式字符串或MATLAB自定義的M

18、函數(shù)，x0為初始點(diǎn)，options為指定優(yōu)化參數(shù)選項(xiàng)；輸出參數(shù)x為多元函數(shù)的最小值點(diǎn)，fval為最小值，exitflag和output與單變量情形一致，grad為函數(shù)在解x處的梯度值，hessian為目標(biāo)函數(shù)在解x處的海賽（Hessian）值。fminsearch利用了單純形法的原理，fminunc利用了擬牛頓法的原理。當(dāng)函數(shù)的階數(shù)大于2時(shí)，使用fminunc比f(wàn)minsearch更有效，但當(dāng)所選函數(shù)高度不連續(xù)時(shí)，使用fminsearch效果較好，這兩個(gè)函數(shù)都容易陷入局部?jī)?yōu)化，并且結(jié)果的正確與否還要取決于初值點(diǎn)x0的選取。利用fiminunc求例7的最小值。在名利窗口輸入：>>x,

19、fval=fminunc('2*x(1)3+4*x(1)*x(2)3-10*x(1)*x(2)+x(2)2',0.5,0.5)輸出結(jié)果為：x = 1.0016 0.8335fval = -3.3241此結(jié)果和fminsearch函數(shù)的求解結(jié)果一樣。9.3.3有約束的多元函數(shù)最小值非線性有約束的多元函數(shù)的標(biāo)準(zhǔn)形式為：其中：x,b,beq,lb,ub是向量，A,Aeq為矩陣，C(x)、Ceq(x)是返回向量的函數(shù)，f(x)為目標(biāo)函數(shù)，f(x)、C(x)、Ceq(x)可以是非線性函數(shù)。在MATLAB用fmincon函數(shù)，求有約束的多元函數(shù)的最小值，其調(diào)用格式有：x = fmincon

20、(fun,x0,A,b)x = fmincon(fun,x0,A,b,Aeq,beq)x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub)x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)x,fval = fmincon()x,fval,exitflag = fmincon()x,fval,exitflag,output = fmincon()x,fval,exitflag,output,lambda = fmincon()x,

21、fval,exitflag,output,lambda,grad = fmincon()x,fval,exitflag,output,lambda,grad,hessian = fmincon()其中fun為目標(biāo)函數(shù)，它可用前面的方法定義；x0為初始值；A、b滿足線性不等式約束，若沒(méi)有不等式約束，則取A=，b=；Aeq、beq滿足等式約束，若沒(méi)有，則取Aeq=，beq=；lb、ub滿足，若沒(méi)有界，可設(shè)lb=，ub=；nonlcon的作用是通過(guò)接受的向量x來(lái)計(jì)算非線性不等約束和等式約束分別在x處的估計(jì)C和Ceq，通過(guò)M文件建立，如：>>x = fmincon(myfun,x0,A,b

22、,Aeq,beq,lb,ub,mycon)，先建立非線性約束函數(shù)，并保存為mycon.m：function C,Ceq = mycon(x)C = % 計(jì)算x處的非線性不等約束的函數(shù)值。Ceq = % 計(jì)算x處的非線性等式約束的函數(shù)值。輸出參數(shù)x為多元函數(shù)的最小值點(diǎn)，fval為最小值，exitflag和output與單變量情形一致，lambda是Lagrange乘子，它體現(xiàn)哪一個(gè)約束有效，grad為函數(shù)在解x處的梯度值，hessian為目標(biāo)函數(shù)在解x處的海賽（Hessian）值。例8 求下面問(wèn)題在初始點(diǎn)x=(10, 10, 10)處的最優(yōu)解。min sub.to 解：約束條件的標(biāo)準(zhǔn)形式為sub

23、.to 在命令窗口輸入：>> fun= '-x(1)*x(2)*x(3)'x0=10,10,10;A=-1 -2 -2;1 2 2;b=0;72;x,fval=fmincon(fun,x0,A,b)輸出結(jié)果為：Optimization terminated: magnitude of directional derivative in search direction less than 2*options.TolFun and maximum constraint violation is less than options.TolCon.Active inequ

24、alities (to within options.TolCon = 1e-006): lower upper ineqlin ineqnonlin 2 x = 24.0000 12.0000 12.0000fval = -3.4560e+003 例9 求表面積為36而體積最大的長(zhǎng)方體的體積。解 設(shè)長(zhǎng)方體的三個(gè)棱長(zhǎng)分別為x1,x2,x3，則體積v=x1*x2*x3，且滿足約束條件2(x1*x2+x2*x3+x3*x1)=36。首先建立目標(biāo)函數(shù)的M文件myopt1.m，內(nèi)容如下：function f = myopt1(x)f =-x(1)*x(2)*x(3); % 負(fù)號(hào)表示最小值，即x(1)*

25、x(2)*x(3)的最大值再建立非線性約束函數(shù)的M文件mycon.m，內(nèi)容如下：function C,Ceq = mycon(x)C = ;Ceq = x(1)*x(2)+x(2)*x(3)+x(3)*x(1)-18;由于沒(méi)有非線性不等式，故C為。在命令窗口輸入：>> x,fval = fmincon('myopt1',1 1 1,0 0 0,'mycon')輸出結(jié)果為：x = 2.4495 2.4495 2.4495fval = -14.6969結(jié)果表明當(dāng)x1=x2=x3=2.4495時(shí)，體積最大值為14.6969。該問(wèn)題的精確解為x1=x2=x3

26、=時(shí)最大體積為。9.4 二次規(guī)劃二次規(guī)劃問(wèn)題（quadratic programming）的標(biāo)準(zhǔn)形式為：sub.to 其中，H、A、Aeq為矩陣，f、b、beq、lb、ub、x為向量。其它形式的二次規(guī)劃問(wèn)題都可轉(zhuǎn)化為標(biāo)準(zhǔn)形式。MATLAB中使用quadprog函數(shù)求解二次規(guī)劃，其調(diào)用格式有：x = quadprog(H,f,A,b) x = quadprog(H,f,A,b,Aeq,beq) x = quadprog(H,f,A,b,Aeq,beq,lb,ub) x = quadprog(H,f,A,b,Aeq,beq,lb,ub,x0)x = quadprog(H,f,A,b,Aeq,beq

27、,lb,ub,x0,options)x,fval = quadprog() x,fval,exitflag = quadprog() x,fval,exitflag,output = quadprog()x,fval,exitflag,output,lambda = quadprog()例10 求解下面二次規(guī)劃問(wèn)題解：則，在命令窗口中輸入：>> H = 1 -1; -1 2 ;f = -2; -6;A = 1 1; -1 2; 2 1;b = 2; 2; 3;lb = zeros(2,1);x,fval,exitflag,output,lambda = quadprog(H,f,A

28、,b, , ,lb)輸出結(jié)果為：x = %最優(yōu)解 0.6667 1.3333fval = %最優(yōu)值 -8.2222exitflag = %收斂 1output = iterations: 3 algorithm: 'medium-scale: active-set' firstorderopt: cgiterations: lambda = lower: 2x1 double upper: 2x1 double eqlin: 0x1 double ineqlin: 3x1 double>> lambda.ineqlinans = 3.1111 0.4444 0>

29、;> lambda.lowerans = 0 0此結(jié)果說(shuō)明第1、2個(gè)約束條件有效，其余無(wú)效。9.5 多目標(biāo)規(guī)劃多目標(biāo)規(guī)劃是指在一組約束下，對(duì)多個(gè)不同目標(biāo)函數(shù)進(jìn)行優(yōu)化。它的一般形式為其中：x、b、beq、lb、ub是向量；A、Aeq為矩陣；C(x)、Ceq(x)和F(x)是返回向量的函數(shù)；F(x)、C(x)、Ceq(x)可以是非線性函數(shù)；weight為權(quán)值系數(shù)向量，用于控制對(duì)應(yīng)的目標(biāo)函數(shù)與用戶定義的目標(biāo)函數(shù)值的接近程度；goal為用戶設(shè)計(jì)的與目標(biāo)函數(shù)相應(yīng)的目標(biāo)函數(shù)值向量；為一個(gè)松弛因子標(biāo)量；F(x)為多目標(biāo)規(guī)劃中的目標(biāo)函數(shù)向量。MATLAB中求解多目標(biāo)規(guī)劃的函數(shù)為fgoalattain，其

30、調(diào)用格式有：x = fgoalattain(fun,x0,goal,weight)x = fgoalattain(fun,x0,goal,weight,A,b)x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq)x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub)x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon)x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon,options)x,fval = fgoalattain()x,fval,attainfactor = fgoalattain()x,fval,attainfactor,exitflag = fgoalattain()x,fval,attainfactor,exitflag,output = fgoalattain()x,fval,attainfactor,exitflag,output,lambda = fgoalattain()輸入?yún)?shù)x0為初始解向量；fun為多目標(biāo)函數(shù)的文件名字符