實驗數(shù)學(xué)六：MATLAB在優(yōu)化中的應(yīng)用.ppt

1,MATLAB在優(yōu)化中的應(yīng)用,2.1實驗?zāi)康?在研究與解決具體問題中，經(jīng)常遇到有關(guān)優(yōu)化問題，本實驗的目的是學(xué)會用MATLAB軟件求解一些優(yōu)化問題，包括求解線性規(guī)劃和非線性規(guī)劃等問題。,2.2實驗內(nèi)容,1、線性規(guī)劃求解,線性規(guī)劃是運籌學(xué)的一個重要分支，它起源于工業(yè)生產(chǎn)組織管理的決策問題。在數(shù)學(xué)上它用來確定多變量線性函數(shù)在變量滿足線性約束條件下的最優(yōu)值；隨著計算機的發(fā)展，出現(xiàn)了如單純形法等有效算法，它在工農(nóng)業(yè)、軍事、交通運輸、決策管理與規(guī)劃等領(lǐng)域中有廣泛的應(yīng)用。,2,問題一：,任務(wù)分配問題：某車間有甲、乙兩臺機床，可用于加工三種工件。假定這兩臺車床的可工作時間分別為800和900，三種工件的數(shù)量分別為400、600和500，且已知車床甲加工單位數(shù)量三種工件所需的時間和加工費分別為0.4、1.1、1和13、9、10，車床乙加工單位數(shù)量三種工件所需的時間和加工費分別為0.5、1.2、1.3和11、12、8。問怎樣分配車床的加工任務(wù)，才能既滿足加工工件的要求，又使加工費用最低？,3,解設(shè)在甲車床上加工工件1、2、3的數(shù)量分別為x1、x2、x3，在乙車床上加工工件1、2、3的數(shù)量分別為x4、x5、x6。,可建立以下線性規(guī)劃模型：,4,問題二：某廠每日8小時的產(chǎn)量不低于1800件。為了進行質(zhì)量控制，計劃聘請兩種不同水平的檢驗員。一級檢驗員的標(biāo)準(zhǔn)為：速度25件/小時，正確率98%，計時工資4元/小時；二級檢驗員的標(biāo)準(zhǔn)為：速度15件/小時，正確率95%，計時工資3元/小時。檢驗員每錯檢一次，工廠要損失2元。為使總檢驗費用最省，該工廠應(yīng)聘一級、二級檢驗員各幾名？,解設(shè)需要一級和二級檢驗員的人數(shù)分別為x1、x2人,則應(yīng)付檢驗員的工資為：,因檢驗員錯檢而造成的損失為：,故目標(biāo)函數(shù)為：,5,約束條件為：,線性規(guī)劃模型：,6,一般線性規(guī)劃問題的數(shù)學(xué)表達式：,max(min)f=,s.t,.,用矩陣向量符號表示：,max(min)f=cX,7,其中,max(min)f=cX,8,用MATLAB優(yōu)化工具箱解線性規(guī)劃的命令如下：,1、模型：minz=cXs.tAXb,命令：x=linprog（c，A，b）,2、模型：minz=cX,命令：x=linprog(c,A,b,Aeq,beq),如果沒有不等式：,存在，則令A(yù)=,b=.,9,3、模型：minz=cX,命令：x=linprog(c,A,b,Aeq,beq,vlb,vub),如果沒有等式約束：,則令A(yù)eq=,beq=.,存在，,4、命令：x=linprog(c,A,b,Aeq,beq,vlb,vub，x0),也用于求模型3，其中x0表示初始值。,5、命令：x,fval=linprog(),返回最優(yōu)解及處的目標(biāo)函數(shù)值fval.,10,解用命令3，編寫M文件xxgh1.m如下：,c=-0.4-0.28-0.32-0.72-0.64-0.6;,A=0.010.010.010.030.030.03;0.02000.0500;00.02000.050;000.03000.08;,b=850;700;100;900;,Aeq=;beq=;,vlb=0;0;0;0;0;0;vub=;,11,x,fval=linprog(c,A,b,Aeq,beq,vlb,vub),結(jié)果：,x=1.0e+004*3.50000.50003.00000.00000.00000.0000,fval=-2.5000e+004,12,例2,s.t,0,i=1,2j=1,2,3.,解用命令3，編寫M文件xxgh2.m如下：,c=10564812;,A=;,b=;,Aeq=111000;000111;100100;010010;001001;,13,beq=60;100;50;70;40;,vlb=0;0;0;0;0;0;vub=;,x,fval=linprog(c,A,b,Aeq,beq,vlb,vub),結(jié)果：,x=0.000020.000040.000050.000050.00000.0000,fval=940.0000,14,例3用MATLAB求問題一：,編寫M文件xxgh3.m如下:,c=1391011128;,A=0.41.11000；0000.51.21.3;,b=800;900;,Aeq=100100；010010；001001;,15,beq=400;600;500;,vlb=0;0;0;0;0;0;,vub=;,x,fval=linprog(c,A,b,Aeq,beq,vlb,vub),結(jié)果:,x=0.0000600.00000.0000400.00000.0000500.0000,fval=1.3800e+004,即在甲機床上加工600個工件2,在乙機床上加工400個工件1、500個工件3，可在滿足條件的情況下使總加工費最小為13800。,16,例4用MATLAB求解問題二：,編寫M文件xxgh4.m如下：,c=40;36;,A=-5-3;,b=-45;,Aeq=;,beq=;,vlb=0;0;,vub=9;15;,17,x,fval=linprog(c,A,b,Aeq,beq,vlb,vub),結(jié)果為：,x=9.00000.0000,fval=360,即只需聘用9個一級檢驗員。,注：本問題應(yīng)還有一個約束條件：x1、x2取整數(shù)。故它是一個整數(shù)線性規(guī)劃問題。這里把它當(dāng)成一個線性規(guī)劃來解，求得其最優(yōu)解剛好是整數(shù)：x1=9，x2=0，故它就是該整數(shù)規(guī)劃的最優(yōu)解。若用線性規(guī)劃解法求得的最優(yōu)解不是整數(shù)，將其取整后不一定是相應(yīng)整數(shù)規(guī)劃的最優(yōu)解，這樣的整數(shù)規(guī)劃應(yīng)用專門的方法求解。,18,例5廠址選擇問題,考慮A，B，C，三地，每地都出產(chǎn)一定數(shù)量的原材料，也消耗一定數(shù)量的產(chǎn)品（見下表6.1）。已知制成每噸產(chǎn)品需3噸原料，各地之間的距離為：A-B：150km，A-C:100km，B-C：200km.假定每萬噸原料運輸1km的運價是5000元，每萬噸產(chǎn)品運輸1km的運價是6000元。由于地區(qū)條件的差異，在不同地點設(shè)廠的生產(chǎn)費用也不同。問究竟在哪些地方設(shè)廠，規(guī)模多大，才使總費用最小？另外，由于其它條件的限制，在B處建廠的規(guī)模（生產(chǎn)的產(chǎn)品數(shù)量）不能超過5萬噸。,19,解：,根據(jù)題意，可以建立問題的數(shù)學(xué)模型（其中目標(biāo)函數(shù)包括原材料運輸費、產(chǎn)品運輸費和生產(chǎn)費）：,20,首先輸入下列系數(shù),c=75755050100100150240210120160220;,A=1-11-100330000-11001-100330000-11-110000330000000001100;,b=20;16;24;5;,Aeq=000000101010000000010101;,beq=7;13;,vlb=zeros(12,1);,x,fval=linprog(c,A,b,Aeq,beq,vlb,),21,x=0.00001.00000.00000.00000.00000.00007.00000.00000.00005.00000.00008.0000,fval=3.4850e+003,要使總費用最小，需要B地向A地運送1萬噸；A，B，C三地建廠規(guī)模分別為7萬噸、5萬噸和8萬噸；最小總費用為3485萬元。,22,2、無約束優(yōu)化問題求解,（1）一元函數(shù)無約束優(yōu)化問題求解,在實驗二中，我們已介紹過求一元函數(shù)的極值問題，實際上它是一元函數(shù)無約束優(yōu)化問題。,一元函數(shù)無約束優(yōu)化問題:minf（x）,常用格式如下：,1）x,fval=fminbnd(fun,x1,x2)：表示在x1,x2上函數(shù)fan的最小值點及最小值，fan為目標(biāo)函數(shù)的表達式字符串。,2）x,fval=fminbnd(fun,x1,x2，options)：其中options為指定優(yōu)化參數(shù)選項。,函數(shù)fminbnd的算法基于黃金分割法和二次插值法，要求目標(biāo)函數(shù)必須連續(xù)函數(shù)，并可能只給出局部最優(yōu)解。,23,例6求f=2,在0f=2*exp(-x)*sin(x);,xmin,ymin=fminbnd(f,0,8),xmin=3.9270ymin=-0.0279,畫出函數(shù)圖形:,fplot(f,0,8);,結(jié)果如圖,如何求最大值？,24,下面求函數(shù)最大值:,f1=-2*exp(-x)*sin(x);,xmin,ymin=fminbnd(f1,0,8),xmin=0.7854ymin=-0.6448,結(jié)果：xmin=3.9270ymin=-0.0279xmax=0.7854ymax=0.6448,25,(2)無約束多元函數(shù)優(yōu)化問題求解,標(biāo)準(zhǔn)型為：minf(X),其中X為向量。,命令格式為:,（1）x,fval=fminsearch（fun,x0）:其中x0為初始值，fan為目標(biāo)函數(shù)的表達式字符串。,（2）x,fval=fminsearch（fun,x0，options）：其中options為指定優(yōu)化參數(shù)選項。,例7,先作曲面的圖形及等高線圖,X,Y=meshgrid(-10:0.2:10,-10:0.2:10);,Z=4*X.2+25*Y.2-2*X+10*Y+2;,mesh(X,Y,Z),26,axis(-10,10,-10,10,-100,100);,contour(X,Y,Z,40);,grid,結(jié)果見圖6.1(a)(b),由圖可知，取x0=00,fun=inline(4*x(1)2+25*x(2)2-2*x(1)+10*x(2)+2,x);,27,X,fval=fminsearch(fun,0,0),X=0.2500-0.2000,fval=0.7500,28,3、二次規(guī)劃的求解,標(biāo)準(zhǔn)型為：,MinZ=,XTHX+cTXs.t.AX=b,VLBXVUB,用MATLAB軟件求解,其輸入格式如下:,1）x,fval=quadprog(H,C,A,b);,2）x,fval=quadprog(H,C,A,b,Aeq,beq);,3)x,fval=quadprog(H,C,A,b,Aeq,beq,VLB,VUB);,4)x,fval=quadprog(H,C,A,b,Aeq,beq,VLB,VUB,X0);,29,例8,寫成標(biāo)準(zhǔn)形式：,輸入命令：,H=2-2;-24;,30,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),運算結(jié)果為：,x=0.80001.2000,z=-7.2000,31,4、一般非線性規(guī)劃求解,非線性規(guī)劃標(biāo)準(zhǔn)型為：,minF(X)s.tAX=b,G(X),Ceq(X)=0VLB,X,VUB,其中X為n維變元向量，G(X)與Ceq(X)均為非線性函數(shù)組成的向量，其它變量的含義與線性規(guī)劃、二次規(guī)劃中相同.用MATLAB求解上述問題，基本步驟分三步：,1.首先建立M文件fun.m,定義目標(biāo)函數(shù)F（X）:,functionf=fun(X);,f=F(X);,32,2、若約束條件中有非線性約束:G(X),則建立M文件nonlcon.m定義函數(shù)G(X)與Ceq(X):,或Ceq(X)=0,functionG,Ceq=nonlcon(X),G=.,Ceq=.,3.建立主程序.非線性規(guī)劃求解的函數(shù)是fmincon,命令的基本格式如下：,(1)x,fval=fmincon(fun,X0,A,b),(2)x,fval=fmincon(fun,X0,A,b,Aeq,beq),(3)x,fval=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB),(4)x,fval=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon),33,(4)x,fval=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon),說明：fan為目標(biāo)函數(shù)的表達式字符串；如果沒有不等式：AXb存在，則令A(yù)=,b=；如果沒有等式約束：AeqX=beq存在，則令A(yù)eq=,beq=；如果沒有不等式：vlbxvub約束，則令vlb=,ulb=；nonlcon的作用是通過接收的向量x來計算非線性不等約束G(x)0和等式約束Ceq(x)=0分別在x處的估計G和Ceq,通過指定函數(shù)柄來使用；函數(shù)可能會出現(xiàn)局部最優(yōu)解，這與初值x0的選取有關(guān)。,34,例9,1）先建立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);,2）再建立M文件mycon.m定義非線性約束：,functiong,ceq=mycon(x),g=1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10;,ceq=;,35,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.運算結(jié)果為：,x=-1.22501.2250,fval=1.8951,36,例10供應(yīng)與選址,某公司有6個建筑工地要開工，每個工地的位置（用平面坐標(biāo)系a，b表示，距離單位：千米）及水泥日用量d(噸)由下表給出。目前有兩個臨時料場位于A(5,1)，B(2,7)，日儲量各有20噸。假設(shè)從料場到工地之間均有直線道路相連。（1）試制定每天的供應(yīng)計劃，即從A，B兩料場分別向各工地運送多少噸水泥，使總的噸千米數(shù)最小。（2）為了進一步減少噸千米數(shù)，打算舍棄兩個臨時料場，改建兩個新的，日儲量各為20噸，問應(yīng)建在何處，節(jié)省的噸千米數(shù)有多大？,37,（一）、建立模型,記工地的位置為(ai，bi)，水泥日用量為di，i=1,6;料場位置為(xj，yj)，日儲量為ej，j=1,2；從料場j向工地i的運送量為Xij。,目標(biāo)函數(shù)為：,約束條件為：,當(dāng)用臨時料場時決策變量為：Xij，,當(dāng)不用臨時料場時決策變量為：Xij，xj，yj。,38,（二）使用臨時料場的情形,使用兩個臨時料場A(5,1)，B(2,7).求從料場j向工地i的運送量為Xij，在各工地用量必須滿足和各料場運送量不超過日儲量的條件下，使總的噸千米數(shù)最小，這是線性規(guī)劃問題.線性規(guī)劃模型為：,其中,，i=1,2,6,j=1,2,為常數(shù)。,39,設(shè)X11=X1,X21=X2,X31=X3,X41=X4,X51=X5,X61=X6,X12=X7,X22=X8,X32=X9,X42=X10,X52=X11,X62=X12,編寫程序gying1.m,clear,a=1.258.750.55.7537.25;b=1.250.754.7556.57.75;d=3547611;x=52;y=17;e=2020;fori=1:6forj=1:2aa(i,j)=sqrt(x(j)-a(i)2+(y(j)-b(i)2);endend,40,CC=aa(:,1);aa(:,2);A=111111000000000000111111;B=20;20;Aeq=100000100000010000010000001000001000000100000100000010000010000001000001;beq=d(1);d(2);d(3);d(4);d(5);d(6);VLB=000000000000;VUB=;x0=123010010101;xx,fval=linprog(CC,A,B,Aeq,beq,VLB,VUB,x0),41,計算結(jié)果為：x=3.00005.00000.00007.00000.00001.00000.00000.00004.00000.00006.000010.0000fval=136.2275,即由料場A、B向6個工地運料方案為：,總的噸千米數(shù)為136.2275。,42,（三）改建兩個新料場的情形,改建兩個新料場，要同時確定料場的位置(xj,yj)和運送量Xij，在同樣條件下使總噸千米數(shù)最小。這是非線性規(guī)劃問題。非線性規(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,43,（1）先編寫M文件liaoch.m定義目標(biāo)函數(shù)。,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;,44,(2)取初值為線性規(guī)劃的計算結(jié)果及臨時料場的坐標(biāo):,x0=35070100406105127;,編寫主程序gying2.m.,clearx0=35070100406105127;,A=11111100000000000000001111110000;B=20;20;,Aeq=100000100000000001000001000000000010000010000000000100000100000000001000001000000000010000010000;,45,beq=3547611;v