數(shù)學(xué)建模講義-最優(yōu)化模型-非線性規(guī)劃

數(shù)學(xué)建模講義最優(yōu)化模型

---非線性規(guī)劃非線性規(guī)劃問(wèn)題1、非線性規(guī)劃模型2、非線性規(guī)劃的性質(zhì)3、非線性規(guī)劃的主要算法。4、用數(shù)學(xué)軟件包求解無(wú)約束最優(yōu)化問(wèn)題5、建模案例選講定義1:如果目標(biāo)函數(shù)或約束條件中至少有一個(gè)是非線性函數(shù)時(shí)的最優(yōu)化問(wèn)題就叫做非線性規(guī)劃問(wèn)題．1非線性規(guī)劃模型一般形式:其中是定義在En上的實(shí)值函數(shù)，簡(jiǎn)記:

定義1把滿足問(wèn)題(1)中條件的解稱為可行解(或可行點(diǎn))，所有可行點(diǎn)的集合稱為可行集（或可行域）．記為D．即問(wèn)題(1)可簡(jiǎn)記為．定義2

對(duì)于問(wèn)題(1),設(shè),若存在,使得對(duì)一切，且，都有，則稱x*是f(x)在D上的局部極小值點(diǎn)(局部最優(yōu)解)．特別地當(dāng)時(shí),若，則稱x*是f(x)在D上的嚴(yán)格局部極小值點(diǎn)（嚴(yán)格局部最優(yōu)解）．2非線性規(guī)劃的性質(zhì)定義3對(duì)于問(wèn)題(1),設(shè),對(duì)任意的,都有

.則稱x*是f(x)在D上的全局極小值點(diǎn)（全局最優(yōu)解）．特別地當(dāng)時(shí)，若，則稱x*是f(x)在D上的嚴(yán)格全局極小值點(diǎn)（嚴(yán)格全局最優(yōu)解）．2.2求解非線性規(guī)劃的基本算法罰函數(shù)法

內(nèi)點(diǎn)法，外點(diǎn)法，Lagrange乘子法序列線性規(guī)劃法序列二次規(guī)劃法信賴域法可行方向法

罰函數(shù)法

罰函數(shù)法基本思想是通過(guò)構(gòu)造罰函數(shù)把約束問(wèn)題轉(zhuǎn)化為一系列無(wú)約束最優(yōu)化問(wèn)題，進(jìn)而用無(wú)約束最優(yōu)化方法去求解．這類方法稱為序列無(wú)約束最小化方法．簡(jiǎn)稱為SUMT法．其一為SUMT外點(diǎn)法，其二為SUMT內(nèi)點(diǎn)法．SUTM外點(diǎn)法

其中T(x,M)稱為罰函數(shù)，M稱為罰因子，帶M的項(xiàng)稱為罰項(xiàng)，這里的罰函數(shù)只對(duì)不滿足約束條件的點(diǎn)實(shí)行懲罰:

當(dāng)時(shí)，滿足的約束條件，故罰項(xiàng)=0,不受懲罰．當(dāng)時(shí),必有，故罰項(xiàng)>0，要受懲罰．1、任意給定初始點(diǎn)x0，取M1>1，給定允許誤差，令k=1；2、求無(wú)約束極值問(wèn)題的最優(yōu)解，設(shè)為xk=x(Mk)，即；3、若存在,使，則取Mk>M()令k=k+1返回(2)，否則，停止迭代．得最優(yōu)解.計(jì)算時(shí)也可將收斂性判別準(zhǔn)則改為

.

SUTM外點(diǎn)法(罰函數(shù)法)的迭代步驟

罰函數(shù)法的缺點(diǎn)是：每個(gè)近似最優(yōu)解xk

往往不是可行解，而只能近似滿足約束，在實(shí)際問(wèn)題中這種結(jié)果可能不能使用；在解一系列無(wú)約束問(wèn)題中，計(jì)算量太大，特別是隨著Mk的增大，可能導(dǎo)致錯(cuò)誤．SUTM內(nèi)點(diǎn)法（障礙函數(shù)法）

內(nèi)點(diǎn)法的迭代步驟

近似規(guī)劃法的基本思想：將問(wèn)題(3)中的目標(biāo)函數(shù)和約束條件近似為線性函數(shù)，并對(duì)變量的取值范圍加以限制，從而得到一個(gè)近似線性規(guī)劃問(wèn)題，再用單純形法求解之，把其符合原始條件的最優(yōu)解作為(3)的解的近似．序列線性規(guī)劃法每得到一個(gè)近似解后，都從這點(diǎn)出發(fā)，重復(fù)以上步驟．

這樣，通過(guò)求解一系列線性規(guī)劃問(wèn)題，產(chǎn)生一個(gè)由線性規(guī)劃最優(yōu)解組成的序列，經(jīng)驗(yàn)表明，這樣的序列往往收斂于非線性規(guī)劃問(wèn)題的解。

序列線性規(guī)劃法的算法步驟如下

返回序列二次規(guī)劃法信賴域法1、二次規(guī)劃用MATLAB軟件求解,其輸入格式如下:1. x=quadprog(H,C,A,b);2. x=quadprog(H,C,A,b,Aeq,beq);3. x=quadprog(H,C,A,b,Aeq,beq,VLB,VUB);4. x=quadprog(H,C,A,b,Aeq,beq,VLB,VUB,x0);5. x=quadprog(H,C,A,b,Aeq,beq,VLB,VUB,x0,options);6. [x,fval]=quaprog(...);7. [x,fval,exitflag]=quaprog(...);8. [x,fval,exitflag,output]=quaprog(...);例1:

minf(x1,x2)=-2x1-6x2+x12-2x1x2+2x22

s.t.

x1+x2≤2-x1+2x2≤2

x1≥0,x2≥01、寫成標(biāo)準(zhǔn)形式：s.t.2、輸入命令：

H=[1-1;-12];

c=[-2;-6];

A=[11;-12];

b=[2;2];

Aeq=[];

beq=[];

VLB=[0;0];

VUB=[];[x,z]=quadprog(H,c,A,b,Aeq,beq,VLB,VUB)3、運(yùn)算結(jié)果為：

x=0.66671.3333

z=-8.22222、一般非線性規(guī)劃

其中x為n維變?cè)蛄?，G(x)與Ceq(x)均為非線性函數(shù)組成的向量，其它變量的含義與線性規(guī)劃、二次規(guī)劃中相同.用Matlab求解上述問(wèn)題，基本步驟分三步：

1.首先建立M文件fun.m,定義目標(biāo)函數(shù)F（X）:functionf=fun(x);f=f(x);3.建立主程序.非線性規(guī)劃求解的函數(shù)是fmincon,命令的基本格式如下：

(1)x=fmincon(‘fun’,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ō)明變量上下限[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法更新拉格朗日Hessian矩陣。大型算法采用的是較高級(jí)的信賴預(yù)算法.[3]fmincon函數(shù)可能會(huì)給出局部最優(yōu)解，這與初值x0的選取有關(guān)。注:1、寫成標(biāo)準(zhǔn)形式：

s.t.

2x1+3x26s.t

x1+4x25

x1,x20例22、先建立M-文件fun3.m:

functionf=fun3(x);f=-x(1)-2*x(2)+(1/2)*x(1)^2+(1/2)*x(2)^23、再建立主程序youh2.m：

x0=[1;1];A=[23;14];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=0.76471.0588fval=-2.02941．先建立M文件fun4.m,定義目標(biāo)函數(shù):

functionf=fun4(x);f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);x1+x2=0s.t.1.5+x1x2-x1-x20-x1x2–10

0例32．再建立M文件mycon.m定義非線性約束：

function[g,ceq]=mycon(x)g=[x(1)+x(2);1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10];3．主程序youh3.m為:x0=[-1;1];A=[];b=[];Aeq=[11];beq=[0];vlb=[];vub=[];[x,fval]=fmincon('fun4',x0,A,b,Aeq,beq,vlb,vub,'mycon')3.運(yùn)算結(jié)果為：

x=-1.22501.2250fval=1.8951

例4

1．先建立M-文件fun.m定義目標(biāo)函數(shù):functionf=fun(x);f=-2*x(1)-x(2);2．再建立M文件mycon2.m定義非線性約束：

function[g,ceq]=mycon2(x)g=[x(1)^2+x(2)^2-25;x(1)^2-x(2)^2-7];

3.主程序fxx.m為:x0=[3;2.5];VLB=[00];VUB=[510];[x,fval,exitflag,output]=fmincon('fun',x0,[],[],[],[],VLB,VUB,'mycon2')4.運(yùn)算結(jié)果為:x=4.00003.0000fval=-11.0000exitflag=1output=iterations:4funcCount:17stepsize:1algorithm:[1x44char]firstorderopt:[]cgiterations:[]注：NLP雖然可用現(xiàn)成的數(shù)學(xué)軟件求解(LINGO,MATLAB)，但是其結(jié)果常依賴于初值的選擇。而且所求得點(diǎn)也僅僅是一個(gè)局部最優(yōu)點(diǎn)（凸規(guī)劃除外），對(duì)于實(shí)際問(wèn)題往往不能得到滿意的結(jié)果.

這就要求建模的時(shí)候盡量不要建成一個(gè)很復(fù)雜的非線性規(guī)劃模型。適當(dāng)?shù)臅r(shí)候還要對(duì)模型進(jìn)行建化分析，或采用一些局部調(diào)整的方法已得到更好的近似解.應(yīng)用實(shí)例：供應(yīng)與選址

某公司有6個(gè)建筑工地要開工，每個(gè)工地的位置（用平面坐標(biāo)系a，b表示，距離單位：千米）及水泥日用量d(噸)由下表給出。目前有兩個(gè)臨時(shí)料場(chǎng)位于A(5,1)，B(2,7)，日儲(chǔ)量各有20噸。假設(shè)從料場(chǎng)到工地之間均有直線道路相連。（1）試制定每天的供應(yīng)計(jì)劃，即從A，B兩料場(chǎng)分別向各工地運(yùn)送多少噸水泥，使總的噸千米數(shù)最小。（2）為了進(jìn)一步減少噸千米數(shù)，打算舍棄兩個(gè)臨時(shí)料場(chǎng)，改建兩個(gè)新的，日儲(chǔ)量各為20噸，問(wèn)應(yīng)建在何處，節(jié)省的噸千米數(shù)有多大？（一）、建立模型

記工地的位置為(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，yj。（二）使用臨時(shí)料場(chǎng)的情形

使用兩個(gè)臨時(shí)料場(chǎng)A(5,1)，B(2,7).求從料場(chǎng)j向工地i的運(yùn)送量為Xij，在各工地用量必須滿足和各料場(chǎng)運(yùn)送量不超過(guò)日儲(chǔ)量的條件下，使總的噸千米數(shù)最小，這是線性規(guī)劃問(wèn)題.線性規(guī)劃模型為：設(shè)X11=X1,X21=X2,,X31=X3,X41=X4,X51=X5,,X61=X6X12=X7,X22=X8,,X32=X9,X42=X10,X52=X11,,X62=X12

編寫程序gying1.m計(jì)算結(jié)果為：x=[3.00005.00000.00007.00000.00001.00000.00000.00004.00000.00006.000010.0000]’fval=136.2275（三）改建兩個(gè)新料場(chǎng)的情形

改建兩個(gè)新料場(chǎng)，要同時(shí)確定料場(chǎng)的位置(xj,yj)和運(yùn)送量Xij，在同樣條件下使總噸千米數(shù)最小。這是非線性規(guī)劃問(wèn)題。非線性規(guī)劃模型為：設(shè)X11=X1,X21=X2,,X31=X3,X41=X4,X51=X5,,X61=X6X12=X7,X22=X8,,X32=X9,X42=X10,X52=X11,,X62=X12

x1=X13,y1=X14,x2=X15,y2=X16

（1）先編寫M文件liaoch.m定義目標(biāo)函數(shù)。MATLAB（liaoch）(2)取初值為線性規(guī)劃的計(jì)算結(jié)果及臨時(shí)料場(chǎng)的坐標(biāo):x0=[35070100406105127]';編寫主程序gying2.m.MATLAB（gying2）functionf=liaoch(x)a=[1.258.750.55.7537.25];b=[1.250.754.7556.57.75];d=[3547611];e=[2020];f1=0;fori=1:6s(i)=sqrt((x(13)-a(i))^2+(x(14)-b(i))^2);f1=s(i)*x(i)+f1;endf2=0;fori=7:12s(i)=sqrt((x(15)-a(i-6))^2+(x(16)-b(i-6))^2);f2=s(i)*x(i)+f2;endf=f1+f2;clear%x0=[35070100406105127]';%x0=[3.00005.00000.07077.000000.9293003.929306.000010.07076.38754.39435.75117.1867]';%x0=[3.00005.00000.30947.00000.01080.6798003.690605.989210.32025.53694.91945.82917.2852]';x0=[35471000005115.63484.86877.24797.7499]';A=[1111110000000000；0000001111110000];B=[20;20];Aeq=[100000100000000001000001000000000010000010000000000100000100000000001000001000000000010000010000];beq=[3547611]';vlb=[zeros(12,1);-inf;-inf;-inf;-inf];vub=[];[x,fval,exitflag]=fmincon('liaoch',x0,A,B,Aeq,beq,vlb,vub)(3)計(jì)算結(jié)果為：x=[3.00005.00000.07077.000000.9293003.929306.000010.07076.38754.39435.75117.1867]’fval=105.4626exitflag=1(4)若修改主程序gying2.m,取初值為上面的計(jì)算結(jié)果:x0=[3.00005.00000.07077.000000.9293003.929306.000010.07076.38754.39435