運籌學(xué)優(yōu)化試驗

運籌學(xué)優(yōu)化試驗一、用Excel求非線性規(guī)劃模型的命令格式和實例第2頁,共26頁，2024年2月25日，星期天例1求例2求第3頁,共26頁，2024年2月25日，星期天二、用Mathematica求非線性規(guī)劃模型的命令格式和實例1、用駐點法求函數(shù)f(x)在[a,b]區(qū)間上的最大和最小值1）先求駐點表zd=Solve[f[x]==0,x]2）再求函數(shù)在表zd中點處的{變量值、函數(shù)值}表與{端點值、端點函數(shù)值}并集表zdgh=Union[{x,f[x]}/.zdg，{{a，f(a)}，{b，f(b)}}]3）求表zdgh中的函數(shù)值的最大、最小值

fmax=Max[Transpose[zdgh][[2]]]fmin=Min[Transpose[zdgh][[2]]]4）求出函數(shù)取最大、最小值的點

x1=Position[zdgh,fmin]x2=Position[zdgh,fmax]xmin=zdgh[[x1[[1,1]]]]xmax=zdgh[[x2[[1,1]]]第4頁,共26頁，2024年2月25日，星期天注：學(xué)習(xí)編程后我們可把這種方法編成一個通用程序。

Transpose[zdgh]表示把表zdgh轉(zhuǎn)置

zdgh[[2]]表示獲取表zdgh的第2行

Transpose[zdgh][[2]]表示獲取表zdgh的第2列

Position[zdgh,fmin]表示表zdgh中與元素fmin相匹配的位置2、用改進法求函數(shù)f(x)在[a,b]區(qū)間上的最大和最小值1）繪制函數(shù)f(x)在[a，b]上的圖形，觀察函數(shù)f(x)在[a,b]區(qū)間上可能取最大和最小值的點為x2和x1附近。2）xmin=FindRoot[f’[x]==0,{x,x1}]xmax=FindRoot[f’[x]==0,{x,x2}]3)fmin=f[x]/.xminfmax=f[x]/.xmax注：該方法簡單直觀，但不適宜編通用程序。第5頁,共26頁，2024年2月25日，星期天3、用求極值命令求函數(shù)f(x)在[a,b]區(qū)間上的最大和最小值（P186）1）繪制函數(shù)f(x)在[a，b]上的圖形，觀察函數(shù)f(x)在[a,b]區(qū)間上可能取最大和最小值的點為x2和x1。2）xmin=FindMinmum[f[x],{x,x1}]xmax=FindMinmun[-f[x],{x,x2}]4、多元函數(shù)求極值命令xmin=FindMinmum[f[x,y,z],{x,x0}{y,y0},{z,z0}]表示求多元函數(shù)f(x,y,z)在點{xo,y0,z0}附近的極小值

xmax=FindMinmum[-f[x,y,z],{x,x0}{y,y0},{z,z0}]表示求多元函數(shù)f(x,y,z)在點{xo,y0,z0}附近的極大值5、有約束的非線性規(guī)劃求解命令NMinimize[{f(X),約束1，約束2,

…,約束m}，{決策變量}]

如：NMinimize[{E^(Sin[50x])+Sin[60E^y]+Sin[70Sin[x]]+Sin[Sin[80y]]-Sin[10(x+y)]+1/4(x^2+y^2),x^2+y^2￡1},{x,y}]第6頁,共26頁，2024年2月25日，星期天三、用MATLAB解優(yōu)化問題其中等式（3）、（4）、（5）的右邊可選用（1）或（2）的等式右邊.

函數(shù)fminbnd的算法基于黃金分割法和二次插值法，它要求目標函數(shù)必須是連續(xù)函數(shù)，并可能只給出局部最優(yōu)解.常用格式如下：（1）x=fminbnd(fun,x1,x2)（2）x=fminbnd(fun,x1,x2，options)（3）[x，fval]=fminbnd（…）（4）[x，fval，exitflag]=fminbnd（…）（5）[x，fval，exitflag，output]=fminbnd（…）第7頁,共26頁，2024年2月25日，星期天MATLAB(wliti1)主程序為wliti1.m:

f='2*exp(-x).*sin(x)';fplot(f,[0,8]);%作圖語句

[xmin,ymin]=fminbnd(f,0,8)f1='-2*exp(-x).*sin

(x)';[xmax,ymax]=fminbnd(f1,0,8)用Excel規(guī)劃求解工具求解用LINGO軟件求解第8頁,共26頁，2024年2月25日，星期天命令格式為:（1）x=fminunc（fun,X0）；或x=fminsearch（fun,X0）（2）x=fminunc（fun,X0，options）；或x=fminsearch（fun,X0，options）（3）[x，fval]=fminunc（...）；或[x，fval]=fminsearch（...）（4）[x，fval，exitflag]=fminunc（...）；或[x，fval，exitflag]=fminsearch（5）[x，fval，exitflag，output]=fminunc（...）；或[x，fval，exitflag，output]=fminsearch（...）

2.多元函數(shù)無約束優(yōu)化問題標準型為：min第9頁,共26頁，2024年2月25日，星期天[3]fminunc為中型優(yōu)化算法的步長一維搜索提供了兩種算法，由選項中參數(shù)LineSearchType控制：

LineSearchType=‘quadcubic’(缺省值)，混合的二次和三次多項式插值；

LineSearchType='cubicpoly'，三次多項式插使用fminunc和fminsearch可能會得到局部最優(yōu)解.說明:fminsearch是用單純形法尋優(yōu).fminunc算法見以下幾點說明：[1]fminunc為無約束優(yōu)化提供了大型優(yōu)化和中型優(yōu)化算法.由選項中的參數(shù)LargeScale控制：LargeScale=‘on'(默認值),使用大型算法LargeScale=‘off'(默認值),使用中型算法[2]fminunc為中型優(yōu)化算法的搜索方向提供了4種算法，由

選項中的參數(shù)HessUpdate控制：

HessUpdate='bfgs'（默認值），擬牛頓法的BFGS公式；HessUpdate='dfp'，擬牛頓法的DFP公式；HessUpdate='steepdesc'，最速下降法第10頁,共26頁，2024年2月25日，星期天例3minMATLAB(wliti3)

1.編寫M文件

fun1.m:function

f=fun1(x)f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);

2.輸入M文件wliti3.m如下:

x0=[-1,1];x=fminunc('fun1',x0);y=fun1(x)

3.運行結(jié)果:

x=0.5000-1.0000y=1.3029e-10=0用LINGO軟件求解用Excel規(guī)劃求解工具求解第11頁,共26頁，2024年2月25日，星期天MATLAB(wliti31)MATLAB(wliti32)第12頁,共26頁，2024年2月25日，星期天3.用fminsearch函數(shù)求解MATLAB(wliti41)輸入命令:

f='100*(x(2)-x(1)^2)^2+(1-x(1))^2';[x,fval,exitflag,output]=fminsearch(f,[-1.22])運行結(jié)果:

x=1.00001.0000fval=1.9151e-010exitflag=1output=iterations:108funcCount:202algorthm:'Nelder-Meadsimplexdirectsearch

'第13頁,共26頁，2024年2月25日，星期天4.用fminunc

函數(shù)MATLAB(wliti44)(1)建立M文件fun2.m

function

f=fun2(x)

f=100*(x(2)-x(1)^2)^2+(1-x(1))^2(2)主程序wliti44.m第14頁,共26頁，2024年2月25日，星期天

Rosenbrock函數(shù)不同算法的計算結(jié)果可以看出，最速下降法的結(jié)果最差.因為最速下降法特別不適合于從一狹長通道到達最優(yōu)解的情況.第15頁,共26頁，2024年2月25日，星期天四）、有約束優(yōu)化問題的數(shù)學(xué)軟件求解1.首先建立M文件fun.m,用來定義目標函數(shù)F（X）:functionf=fun(X);f=F(X);

其中X為n維變元向量，G(X)與Ceq(X)均為非線性函數(shù)組成的向量，MATLAB求解上述問題，基本步驟分三步：第16頁,共26頁，2024年2月25日，星期天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(…)輸出極值點M文件迭代的初值參數(shù)說明變量上下限第17頁,共26頁，2024年2月25日，星期天注意：[1]fmincon函數(shù)提供了大型優(yōu)化算法和中型優(yōu)化算法．默認時:若在fun函數(shù)中提供了梯度（options參數(shù)的GradObj設(shè)置為’on’），并且只有上下界存在或只有等式約束，fmincon函數(shù)將選擇大型算法．當(dāng)既有等式約束又有梯度約束時，使用中型算法．[2]fmincon函數(shù)的中型算法使用的是序列二次規(guī)劃法．在每一步迭代中求解二次規(guī)劃子問題，并用BFGS法更新拉格朗日Hesse矩陣．[3]fmincon函數(shù)可能會給出局部最優(yōu)解，這與初值X0的選取有關(guān)．第18頁,共26頁，2024年2月25日，星期天例1

min

f(x1,x2)=-2x1-6x2+x12-2x1x2+2x22s.t.

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

x1≥0,x2≥0MATLAB（youh1）1．寫成標準形式：

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．運算結(jié)果為：

x=0．66671．3333z=-8．2222s.t.第19頁,共26頁，2024年2月25日，星期天1．寫成標準形式：

s.t.

2x1+3x26

s.t.

x1+4x25

x1,x20例2第20頁,共26頁，2024年2月25日，星期天2．先建立M-文件fun3．m:

functionf=fun3(x);f=-x(1)-2*x(2)+(1/2)*x(1)^2+(1/2)*x(2)^2MATLAB(youh2)3．再建立主程序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．運算結(jié)果為：

x=0．76471．0588fval=-2．0294第21頁,共26頁，2024年2月25日，星期天1．先建立M文件fun4．m定義目標函數(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例3

2．再建立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];第22頁,共26頁，2024年2月25日，星期天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')MATLAB(youh3)4．運算結(jié)果為：

x=-1．22501