數(shù)學(xué)建模第6講 非線(xiàn)性規(guī)劃

1、數(shù)學(xué)建模與數(shù)學(xué)實(shí)驗(yàn)非線(xiàn)性規(guī)劃 8/8/20221數(shù)學(xué)建模實(shí)驗(yàn)?zāi)康膶?shí)驗(yàn)內(nèi)容2 掌握用數(shù)學(xué)軟件求解優(yōu)化問(wèn)題1 直觀(guān)了解非線(xiàn)性規(guī)劃的基本內(nèi)容1非線(xiàn)性規(guī)劃的基本理論4實(shí)驗(yàn)作業(yè)2 用數(shù)學(xué)軟件求解非線(xiàn)性規(guī)劃3 鋼管訂購(gòu)及運(yùn)輸優(yōu)化模型8/8/20222數(shù)學(xué)建模*非線(xiàn)性規(guī)劃的基本解法非線(xiàn)性規(guī)劃的基本概念非線(xiàn)性規(guī)劃 返回8/8/20223數(shù)學(xué)建模 定義 如果目標(biāo)函數(shù)或約束條件中至少有一個(gè)是非線(xiàn)性函數(shù),則最優(yōu)化問(wèn)題就叫做非線(xiàn)性規(guī)劃問(wèn)題非現(xiàn)性規(guī)劃的基本概念 一般形式: （1） 其中 ， 是定義在 Rn 上的實(shí)值函數(shù)，簡(jiǎn)記: 其它情況: 求目標(biāo)函數(shù)的最大值，或約束條件小于等于零兩種情況，都可通過(guò)取其相反數(shù)化為上述一般

2、形式1nj1ni1nR :h ,R :g ,R :RRRf()nTnRxxxX=,21L()()=.,.,2,1 0 m;1,2,., 0. ljXhiXgtsji8/8/20224數(shù)學(xué)建模 定義1 把滿(mǎn)足問(wèn)題（1）中條件的解 稱(chēng)為可行解（或可行點(diǎn)），所有可行點(diǎn)的集合稱(chēng)為可行集（或可行域）記為D即 問(wèn)題(1)可簡(jiǎn)記為 定義2 對(duì)于問(wèn)題(1),設(shè) ,若存在 ,使得對(duì)一切 ,且 ,都有 ,則稱(chēng)X*是f(X)在D上的局部極小值點(diǎn)（局部最優(yōu)解）特別地,當(dāng) 時(shí)，若 ,則稱(chēng)X*是f(X)在D上的嚴(yán)格局部極小值點(diǎn)（嚴(yán)格局部最優(yōu)解）定義3 對(duì)于問(wèn)題(1),設(shè) ,若對(duì)任意的 ,都有則稱(chēng)X*是f(X)在D上的全局

3、極小值點(diǎn)（全局最優(yōu)解）特別地,當(dāng) 時(shí)，若 ,則稱(chēng)X*是f(X)在D上的嚴(yán)格全局極小值點(diǎn)（嚴(yán)格全局最優(yōu)解） 返回)(nRX()()njiRXXhXg XD = =,0,0|()(),Xf Xf *8/8/20225數(shù)學(xué)建模非線(xiàn)性規(guī)劃的基本解法SUTM外點(diǎn)法SUTM內(nèi)點(diǎn)法（障礙罰函數(shù)法）1 罰函數(shù)法2 近似規(guī)劃法 返回8/8/20226數(shù)學(xué)建模 罰函數(shù)法 罰函數(shù)法基本思想是通過(guò)構(gòu)造罰函數(shù)把約束問(wèn)題轉(zhuǎn)化為一系列無(wú)約束最優(yōu)化問(wèn)題，進(jìn)而用無(wú)約束最優(yōu)化方法去求解這類(lèi)方法稱(chēng)為序列無(wú)約束最小化方法簡(jiǎn)稱(chēng)為SUMT法 其一為SUMT外點(diǎn)法，其二為SUMT內(nèi)點(diǎn)法8/8/20227數(shù)學(xué)建模 其中T(X,M)稱(chēng)為罰函數(shù)

4、，M稱(chēng)為罰因子，帶M的項(xiàng)稱(chēng)為罰項(xiàng)，這里的罰函數(shù)只對(duì)不滿(mǎn)足約束條件的點(diǎn)實(shí)行懲罰:當(dāng) 時(shí)，滿(mǎn)足各 ，故罰項(xiàng)為0，不受懲罰當(dāng) 時(shí),必有約束條件 ，故罰項(xiàng)大于0，要受懲罰SUTM外點(diǎn)法8/8/20228數(shù)學(xué)建模 罰函數(shù)法的缺點(diǎn)：每個(gè)近似最優(yōu)解Xk往往不是容許解，而只能近似滿(mǎn)足約束，在實(shí)際問(wèn)題中這種結(jié)果可能不能使用；在解一系列無(wú)約束問(wèn)題中，計(jì)算量太大，特別是隨著Mk的增大，可能導(dǎo)致錯(cuò)誤1任意給定初始點(diǎn) X0，取M11，給定允許誤差 ，令k=1；2求無(wú)約束極值問(wèn)題 的最優(yōu)解，設(shè)Xk=X(Mk)，即 ；3若存在 ，使 ,則取MkM( ),令k=k+1返回（2），否則，停止迭代得最優(yōu)解 計(jì)算時(shí)也可將收斂性判別

5、準(zhǔn)則 改為 SUTM外點(diǎn)法(罰函數(shù)法)的迭代步驟8/8/20229數(shù)學(xué)建模SUTM內(nèi)點(diǎn)法（障礙函數(shù)法）()()()()()()()為障礙因子.為障礙項(xiàng)，或其中稱(chēng)或 ：構(gòu)造障礙函數(shù)rXgrXgrXgrXfrXIXgrXfrXIrXImiimiimiimii=+=+=11111 ln1)(),( ln, 8/8/202210數(shù)學(xué)建模 內(nèi)點(diǎn)法的迭代步驟8/8/202211數(shù)學(xué)建模 近似規(guī)劃法的基本思想：將問(wèn)題(3)中的目標(biāo)函數(shù) 和約束條件 近似為線(xiàn)性函數(shù)，并對(duì)變量的取值范圍加以限制，從而得到一個(gè)近似線(xiàn)性規(guī)劃問(wèn)題，再用單純形法求解之，把其符合原始條件的最優(yōu)解作為(3)的解的近似近似規(guī)劃法每得到一個(gè)近似

6、解，都從這點(diǎn)出發(fā)，重復(fù)以上步驟 這樣，通過(guò)求解一系列線(xiàn)性規(guī)劃問(wèn)題，產(chǎn)生一個(gè)由線(xiàn)性規(guī)劃最優(yōu)解組成的序列，經(jīng)驗(yàn)表明，這樣的序列往往收斂于非線(xiàn)性規(guī)劃問(wèn)題的解8/8/202212數(shù)學(xué)建模 近似規(guī)劃法的算法步驟如下：8/8/202213數(shù)學(xué)建模 返回8/8/202214數(shù)學(xué)建模 用MATLAB軟件求解,其輸入格式如下: 1x=quadprog(H,C,A,b); 2x=quadprog(H,C,A,b,Aeq,beq); 3x=quadprog(H,C,A,b,Aeq,beq,VLB,VUB); 4x=quadprog(H,C,A,b, Aeq,beq ,VLB,VUB,X0); 5x=quadprog

7、(H,C,A,b,Aeq,beq,VLB,VUB,X0,options); 6x,fval=quaprog(); 7x,fval,exitflag=quaprog(); 8x,fval,exitflag,output=quaprog();1二次規(guī)劃8/8/202215數(shù)學(xué)建模例1 min f(x1,x2)=-2x1-6x2+x12-2x1x2+2x22 s.t. x1+x22 -x1+2x22 x10, x20 MATLAB（youh1）1寫(xiě)成標(biāo)準(zhǔn)形式： 2輸入命令： H=1 -1; -1 2; c=-2 ;-6;A=1 1; -1 2;b=2;2; Aeq=;beq=; VLB=0;0;VU

8、B=; x,z=quadprog(H,c,A,b,Aeq,beq,VLB,VUB)3運(yùn)算結(jié)果為： x =06667 13333 z = -82222s.t.8/8/202216數(shù)學(xué)建模 1 首先建立M文件fun.m,用來(lái)定義目標(biāo)函數(shù)F（X）:function f=fun(X);f=F(X);2一般非線(xiàn)性規(guī)劃 其中X為n維變?cè)蛄?，G(X)與Ceq(X)均為非線(xiàn)性函數(shù)組成的向量，其他變量的含義與線(xiàn)性規(guī)劃、二次規(guī)劃中相同用MATLAB求解上述問(wèn)題，基本步驟分三步：8/8/202217數(shù)學(xué)建模3 建立主程序.求解非線(xiàn)性規(guī)劃的函數(shù)是fmincon,命令的基本格式如下： (1) x=fmincon(fu

9、n,X0,A,b) (2) x=fmincon(fun,X0,A,b,Aeq,beq) (3) x=fmincon(fun,X0,A,b, Aeq,beq,VLB,VUB) (4) x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon)(5)x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon,options) (6) x,fval= fmincon() (7) x,fval,exitflag= fmincon() (8)x,fval,exitflag,output= fmincon()輸出極值點(diǎn)M文件迭代的初值參數(shù)說(shuō)明變

10、量上下限8/8/202218數(shù)學(xué)建模注意：1 fmincon函數(shù)提供了大型優(yōu)化算法和中型優(yōu)化算法默認(rèn)時(shí): 若在fun函數(shù)中提供了梯度（options參數(shù)的GradObj設(shè)置為on），并且只有上下界存在或只有等式約束，fmincon函數(shù)將選擇大型算法當(dāng)既有等式約束又有梯度約束時(shí)，使用中型算法2 fmincon函數(shù)的中型算法使用的是序列二次規(guī)劃法在每一步迭代中求解二次規(guī)劃子問(wèn)題，并用BFGS法更新拉格朗日Hesse矩陣3 fmincon函數(shù)可能會(huì)給出局部最優(yōu)解，這與初值X0的選取有關(guān)8/8/202219數(shù)學(xué)建模1寫(xiě)成標(biāo)準(zhǔn)形式： s.t. 2x1+3x2 6 s.t. x1+4x2 5 x1,x2

11、0例28/8/202220數(shù)學(xué)建模2先建立M-文件 fun3m: function f=fun3(x); f=-x(1)-2*x(2)+(1/2)*x(1)2+(1/2)*x(2)2MATLAB(youh2)3再建立主程序youh2m： x0=1;1; A=2 3 ;1 4; b=6;5; Aeq=;beq=; VLB=0;0; VUB=; x,fval=fmincon(fun3,x0,A,b,Aeq,beq,VLB,VUB)4運(yùn)算結(jié)果為： x = 07647 10588 fval = -202948/8/202221數(shù)學(xué)建模1先建立M文件fun4m定義目標(biāo)函數(shù): function f=fun

12、4(x); f=exp(x(1) *(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2)+1); x1+x2=0 s.t. 1.5+x1x2 - x1 - x2 0 -x1x2 10 0例3 2再建立M文件myconm定義非線(xiàn)性約束： function g,ceq=mycon(x) g=x(1)+x(2);15+x(1)*x(2)-x(1)-x(2); -x(1)*x(2)-10;8/8/202222數(shù)學(xué)建模3主程序youh3m為:x0=-1;1;A=;b=;Aeq=1 1;beq=0;vlb=;vub=;x,fval=fmincon(fun4,x0,A,b,Aeq,beq,

13、vlb, vub,mycon)MATLAB(youh3)4 運(yùn)算結(jié)果為： x = -12250 12250 fval = 189518/8/202223數(shù)學(xué)建模 例4 1先建立M文件funm定義目標(biāo)函數(shù): function f=fun(x); f=-2*x(1)-x(2);2再建立M文件mycon2m定義非線(xiàn)性約束：function g,ceq=mycon2(x)g=x(1)2+x(2)2-25;x(1)2-x(2)2-7;8/8/202224數(shù)學(xué)建模3 主程序fxxm為: x0=3;25; VLB=0 0;VUB=5 10; x,fval,exitflag,output =fmincon(f

14、un,x0, VLB,VUB,mycon2)MATLAB(fxx(fun)8/8/202225數(shù)學(xué)建模4 運(yùn)算結(jié)果為: x = 40000 30000fval =-110000exitflag = 1output = iterations: 4 funcCount: 17 stepsize: 1 algorithm: 1x44 char firstorderopt: cgiterations: 返回8/8/202226數(shù)學(xué)建模應(yīng)用實(shí)例： 供應(yīng)與選址 某公司有6個(gè)建筑工地要開(kāi)工，每個(gè)工地的位置（用平面坐標(biāo)系a，b表示，距離單位：km）及水泥日用量d(t)由下表給出目前有兩個(gè)臨時(shí)料場(chǎng)位于A(5,1

15、)，B(2,7)，日儲(chǔ)量各有20t假設(shè)從料場(chǎng)到工地之間均有直線(xiàn)道路相連 （1）試制定每天的供應(yīng)計(jì)劃，即從A，B兩料場(chǎng)分別向各工地運(yùn)送多少水泥，可使總的噸千米數(shù)最小 （2）為了進(jìn)一步減少?lài)嵡讛?shù)，打算舍棄兩個(gè)臨時(shí)料場(chǎng)，改建兩個(gè)新的，日儲(chǔ)量各為20t，問(wèn)應(yīng)建在何處，節(jié)省的噸千米數(shù)有多大？8/8/202227數(shù)學(xué)建模（一）建立模型 記工地的位置為(ai，bi)，水泥日用量為di，i=1,6;料場(chǎng)位置為(xj，yj)，日儲(chǔ)量為ej，j=1,2；料場(chǎng)j向工地i的運(yùn)送量為Xij當(dāng)用臨時(shí)料場(chǎng)時(shí)決策變量為：Xij，當(dāng)不用臨時(shí)料場(chǎng)時(shí)決策變量為：Xij，xj，yj8/8/202228數(shù)學(xué)建模（二）使用臨時(shí)料場(chǎng)的情

16、形 使用兩個(gè)臨時(shí)料場(chǎng)A(5,1)，B(2,7)求從料場(chǎng)j向工地i的運(yùn)送量Xij . 在各工地用量必須滿(mǎn)足和各料場(chǎng)運(yùn)送量不超過(guò)日儲(chǔ)量的條件下，使總的噸千米數(shù)最小，這是線(xiàn)性規(guī)劃問(wèn)題 線(xiàn)性規(guī)劃模型為：設(shè) X11=X1, X21= X 2, X31= X 3, X41= X 4, X51= X 5, X61= X 6X12= X 7, X22= X 8, X32= X 9, X42= X 10, X52= X 11, X62= X 12 編寫(xiě)程序gying1mMATLAB（gying1）8/8/202229數(shù)學(xué)建模計(jì)算結(jié)果為：x = 30000 50000 00000 70000 00000 1000

17、0 00000 00000 40000 00000 60000 100000fval = 13622758/8/202230數(shù)學(xué)建模（三）改建兩個(gè)新料場(chǎng)的情形 改建兩個(gè)新料場(chǎng)，要同時(shí)確定料場(chǎng)的位置(xj,yj)和運(yùn)送量Xij，在同樣條件下使總噸千米數(shù)最小這是非線(xiàn)性規(guī)劃問(wèn)題非線(xiàn)性規(guī)劃模型為：8/8/202231數(shù)學(xué)建模設(shè) X11=X1, X21= X 2, X31= X 3, X41= X 4, X51= X 5, X61= X 6 X12= X 7, X22= X 8, X32= X 9, X42= X 10, X52= X 11, X62= X 12 x1=X13, y1=X14, x2=X

18、15, y2=X16 （1）先編寫(xiě)M文件liaochm定義目標(biāo)函數(shù)MATLAB（liaoch）(2) 取初值為線(xiàn)性規(guī)劃的計(jì)算結(jié)果及臨時(shí)料場(chǎng)的坐標(biāo): x0=3 5 0 7 0 1 0 0 4 0 6 10 5 1 2 7;編寫(xiě)主程序gying2mMATLAB（gying2）8/8/202232數(shù)學(xué)建模(3) 計(jì)算結(jié)果為：x= 30000 50000 00707 70000 0 09293 0 0 39293 0 60000 100707 63875 43943 57511 71867fval = 1054626exitflag = 18/8/202233數(shù)學(xué)建模(4) 若修改主程序gying2m

19、, 取初值為上面的計(jì)算結(jié)果:x0= 30000 50000 00707 70000 0 09293 0 0 39293 0 60000 100707 63875 43943 57511 71867 則得結(jié)果為:x=30000 50000 03094 70000 00108 06798 0 0 36906 0 59892 103202 55369 49194 58291 72852fval =1034760exitflag = 1總的噸千米數(shù)比上面結(jié)果略?xún)?yōu) (5) 若再取剛得出的結(jié)果為初值, 卻計(jì)算不出最優(yōu)解MATLAB（gying2）MATLAB（gying2）8/8/202234數(shù)學(xué)建模(6) 若取初值為: x0=3 5 4 7 1 0 0 0 0 0 5 11 56348 48687 72479 77499, 則計(jì)算結(jié)果為:x=30000 50000 40000 70000