matlab与数学建模-第03章非线性

/第三章§1例1（投资决策问题）某企业有n个项目可供选择投资，并且至少要对其中一个 A元，投资于第i(i1,L,n)个项目需花 并预计可收益bi元。试选择最佳投资方案。解设投资决策变量为

决定投资第i

，i1,L,n0，决定不投资第i个项 则投资总额为aixi，投资总收益为bixi。因为该公司至少要对一个项目投资，i in0aixiixi(1xi)0,i1,L,结为总以及决策变量（取0或1）的限制条件下，极大化总收益和总投资之比。因此，nbinnmaxQnais.t.0aixiixi(1xi)0,i1,L,minf( hj(x)gi(x)

j1,L,qi1,L,p

xx1Lxn]T称为模型（NP）f称为目标函数，gi(i1,Lp)hjj1,Lq)称为约束函数。另外，gix)0(i1,L,p)称为等式约束，hjx)0(j1,Lq)称为不等式的约束。非线性规划的中非线性规划的数学模型写成以下形式minf(AxAeqxC(x) A,B,Aeq,Beq定义了线性约束A*XB,Aeq*XBeq，如果没有线性约束，则A=[],B=[],Aeq=[],Beq=[]；LB和UB是变量x的下界和上界，如果上界和下界没有约束，则LB=[]，UB=[]，如x无下界，则LB各分量都为-inf，如果x无上UB定义了优化参数，可以使用缺省的参数设置。 minf(x)x2x2x2 13 x2x2x213xx2x3 2x1x223x22x23x1,x2,x3解（i）Mfun1.m编写M文件fun2.mx(1)+x(2)^2+x(3)^3-20];%非线性不等式约束;%编写主程序文件example2.m如下：'fun2',options)记（NP）的可行域为K。x*Kf(x*)f( x则称x*是（NP）的整体最优解，f(x*)是(NP)的整体最优值。如果有f(x*)f( xK,xx*Kx*N(x*，使f(x*)f( xN(x*)IK则称x*是（NP）的局部最优解，f(x*)是(NP)的局部最优值。如果有f(x*)f( xN(x*)I则称x*是（NP）的严格局部最优解，f(x*)是(NP)的严格局部最想是：从一个选定的初始点x0Rn{xk}，使得当{xk是有穷点列时，其最后一个点是(NP)的最优解；当{xk是无穷点xkRn是某迭代方法的第kxk1Rn是第k1xk1xktk 这里tkR1pkRn,pk1pkxkxk1的方向。式（1）通常，把基本迭代格式（1）中的pk称为第k轮搜索方向，tk为沿pk方向的xRnp0，若存在0fxtpfxt(0,，称向量pfx处的下降方向。xRnp0，若存在t0xtpK向，称之为函数fx处关于K的可行下降方向。现在，给出用基本迭代格式（1）求解(NP)的一般步骤如下0°x0，令k:01°fxkK的可行下降方向作为搜索方向pk。2xkpk寻求适当的步长tk，使目标函数值3°求出下一个迭代点。按迭代格式（1）xk1xktkpk若xk1已满足某种终止条件，停fx为定义在n维欧氏空间E(n)中某个凸集R上的函数，若对任何实数(01)R中的任意x(1)x(2)f(x(1)(1)x(2))f(x(1))(1)f(x(2)若对每一个(01)x(1x(2Rf(x(1)(1)x(2))f(x(1))(1)f(x(2)

f(R{x|gj(x)0,j1,2,L,fx为凸函数，gjxj1,2,L,l为凸函数，这样的非线性规划称为 （2（at

f 若f(t)是[a,b]区间上的下单峰函数，介绍通过不断地缩短[a,b]的长度，来为了缩短区间[a,b]，逐步搜索得（2）的最优解t*的近似值，可以采用以下途径：在[ab中任取两个关于[ab是对称的点t1和t2（不妨设t2t1则必有t*[at]，因而[at]是缩短了的单峰区间；若f(t)f(t)，则有t*t*[tb

f(t)f(t，则[at

和[tb FibonacciF0F1FnFn2Fn1 n数列{FnFibonacciFn称为第nFibonacciFibonacci数之Fn1Fibonacci分数。FnFn1，其后各次分别为Fn2Fn3,LF1。由此，若t和t(tt是单峰区间[a

2 t1aFn1，t2aFn2ba ba taFn1(ba)

aFn2(b

F F 它们关于[a,b]确是对称的点过精度0，这就要求最后区间的长度不超过，即ba

应按照预先给定的精度，确定使（4）成立的最小整数n作为搜索次数，直到进行到第n个探索点时停止。f(t的单峰区间[ab的办法，来求得问题（2）的近似解，是Kiefer(1953年)，叫做Finbonacci法，具体步骤如下：1o选取初始数据，确定单峰区间[a0b0，给出搜索精度0，由（4）确定搜索次数n。2ok1,aabb，计算最初两个搜索点，按（3）计算t和t 3owhileknf1f(t1),f2f(t2iff1at;tt;taF(n1k)(b

2 1

F(nkF(n1k

bt;tt;tb (a1 2 F(nkkk4o当进行至kn1t1

1(ab) t22(ata(

)(b t1和t2这两点中，以函数值较小者为近似极小点，相应的函数值为近似极小值。并得最终区间[a,t1或[t2b]。3f(tt2t2的近似极小点，要求缩短后的区间不大于区间[1,3]的0.08倍。 0.618若01 称之为黄金分割数，其值为

510.6180339887L2黄金分割数FibonaccilimFn1n0.6182个探索点开始每增加一个探索点作一轮迭代以后，原单峰区间0.618倍。计算n个探索点的函数值可以把原区间[a0,b0]连续n1次，因为每次的缩短率均为，故最后的区间长度为

)n这就是说，当已知缩短的相对精度为时，可用下式计算探索点个数nn10.618minf( xE(n xk1xktk 告诉，点xk的负梯度方向pkf(xk)xkf下降最快的方向。为此，称负梯度方向fxkfxk处的时，用fxk0fxk)1°选取初始数据。选取初始点x0，给定终止误差，令k:02°求梯度向量。计算fxkf(xk3°pkf(xk)4°进行一维搜索。求tkf(xktpk)minf(xktpk

txk1xktkpkk:k1,2例 minf(x)x2 xx1x2)Tx02,2)T。解（i）fx)(2x150x2TM文件detaf.mf(xwhilenorm(g)>0.000001whilef>f0Newton考虑目标函数f在点xk处的二次近f(x)Q(x)

f(xk)f(xk)T(xxk)1(xxk)T2f(xk)(xxk)2f(xk 2f(xk)

x2f(xk) M

n2f(xk

)2f(xk)xx x

由于2fxk)正定，函数Q的驻点xk1是Qx)的极小点。为求Qxk1fxk2fxkxk1xk0，xk1xk[2f(xk)]1f(xk)对照基本迭代格式（1）可知从点xk出发沿搜索方向。pk[2f(xk)]1f(xktk1即可得Qx的最小点xk1。通常，把方向pk叫做从点xk出发的NewtonNewton方向并取步长为1的求解方法，称之为Newton法。其具体步骤如下：1x0，给定终止误差0，令k:02°求梯度向量。计算fxk

f(xk

3°构Newton方向。计算[2fxk)]1pk[2f(xk)]1f(xk)4°求下一迭代点。令xk1xkpk,k:k1，转2°。例5 用Newton法求解，minf(x)x425x4x2x1x02,2)T

11 100x32x2x112x2

124x2

12f 300x2122x2

1 1（ii）编写M文件nwfun.m如下：function[f,df,d2f]=nwfun(x);x,whilef>f0Newton法的优点是收敛速度快；缺点是有时不好用而需采取改进措施，此外，当维数较高时，计算[2f(xk)]1的工作量很大。变尺度法（VariableMetricAlgorithm）20多年来发展起来的，它不仅是求解尺度法—DFPDavidon1959年提出，后经Fletcher和Powell加以改进。已经知道，法的搜索方向是2fxk)]1fxk)，为了不计算二阶导数矩阵[2f(xk)]及其逆阵，设法构造另一个矩阵，用它来近二阶导数矩阵的逆阵[2f(xk)]1，这一类方法也称拟法（Quasi-NewtonMethod。下面研究如何构造这样的近似矩阵，并将它记H(k)。要求：每一步都能以矩阵最后应收敛于解点处的Hesse阵的逆阵。当fx)是二次函Hesse阵为常数阵Axkxk1处的梯度之f(xk1)f(xk)A(xk1xk或xk1xkA1[f(xk1)f(xkH(k1)满足xk1xkH(k1)[f(xk1)f(xk G(k)f(xk1)f(xkxkxk1

xkH(k1)G(k) 现假H(k)已知，用下H(k1)（设H(k)H(k1)均为对称正定阵H(k1)H(k)H(k 其中H(k)k次校正矩阵。显H(k1)Newton条件（9）即要xk(H(k)H(k))G(k或

H(k)G(k)xkH(k)G(kH(k)的一种比较简单的

H(k)xk(Q(k))TH(k)G(k)(W(k)其中Q(k)和W(k)为两个待将式（12）中的H(k)代入（11）xk(Q(k))TG(k)H(k)G(k)(W(k))TG(k)xkH(k)G(k(Q(k))TG(k)(W(k))TG(k)考虑到H(k)应为对称阵，最简单的办Q(k)kW(k)H(k)G(k

k(xk)TG(k)k(G(k))TH(k)G(k)若(xk)TGk)和(Gk)THkGk)

(xk)TG(k

(G(k))T

(G(k))TH(k)G(k

H(k)

xk(xk)T

H(k)G(k)(G(k))TH(k(G(k))TH(k)G(k

H(k

H(k

xk(xk)T(G(k))T

H(k)G(k)(G(k))TH(k(G(k))TH(k)G(k

上述矩阵称为尺度矩阵。通常，取第一个尺度矩阵H 为单位阵，以后的尺度)若H(k)为对称正定阵，则由式（18）产生的H(k1)现将DFP 1°给定初始点x0及梯度允许误差0。2fx0

H(0)I（单位矩阵p0H(0)f(x0p0方向进行一维搜索，确定最佳步长0minf(x0p0)x1x00

f(x00p04xk，算出fxkf(x0)则xk即为所求的近似解，停止迭代；否则，计算H(k)H(k

H(k

xk1(xk1 (G(k1))Txk

H(k1)G(k1)(G(k1))TH(k(G(k1))TH(k1)G(kpkH(k)fxkpk方向上进行一维搜索，得kkxk1xkk5°若xk1满足精度要求xk1即为所求的近似解，否则，转回4°，直到求理解，因而在实际应用中常为人们所采用。下面介绍Powell方法。°{p01pn1}。给定终止误差0，令k:02y0:xk，依次沿{p01pn1中的方向进行一维搜索。对应地得到辅助迭代点y1,y2,L,yn，即jyjyj1 pjf(yj1

pj1)

f(yj1tp j1,L, /3°pnyny0，若

xk1yn找出m

f(ym1)f(ym)max{f(yj1)f(yj)|1jf(y0)2f(yn)f(2yny0)2[f(ym1)f(ymxk1yn

pn:f(yntpn)minf(yntpn) t{p0,p1pn1}k1:{p0,L,pm1,pm1,L,pn1,pn}k:k12 minf(x)x其中x可以为标量或中fminunc的基本[X,FVAL]=FMINUNC(FUN,X0,OPTIONS,P1,P2,FUNfxFUN有三个返回例6 求函数f(x)100(xx2)2(1x)2的最小值。 解：编写Mfun2.m如下：function[f,g]=fun2(x);g=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-编写主函数文件example6.m如下：options=optimset('GradObj','on');function[f,df,d2f]=fun3(x); -编写主函数文件example62.moptions=optimset('GradObj','on','Hessian','on'); 例7 求函数f(x)sin(x)3取最小值时的x值。解编写f(x)的M文件fun4.m如下：functionf=fun4(x);编写主函数文件example7.m如下： 中二次规划的数学模型可表述如下：

1xTHxfTx,AxAeqx中求解二次规划令是[X,FVAL]= 看在指令中运行helpqurog后的帮助。8minf(x)2x24xx4x26x

x2

1 4x1x2 x1,x2解h=[4,-4;-f=[-6;-x x

Minfx)11.0250SUMT(SequentialUnconstrainedMinizationTechnique)。minf(gi(x)0,i1,L,s.t.hj(x)0,j1,L,k (x)0,m1,L,k取一个充分大的数M0 P(x,M)f(x)Mmax(gi(x),0)Mmin(hi(x),0)M|ki(x)

G(x)

H(x)

（或P(x,M)f(x)Msummax Msummin M||K(x) 这里G(x)g1 gr(x)，H(x)h1 hsK(x)k1(x) kt(t)， 增广目标函数PxM为目标函数的无约束极值问题例9求下列非线性规划minf(x)x2x2 x2x 2 x1

2 x,x 解(i)Mtest.mfunctiong=test(x); g=f-M*sum(min([x';zeros(1,2)]))-M*min(x(1)^2-M*abs(-x(1)- fminbndminf(x)x

x[x1,x2它的返回值是极小点x和函数的极小值。这里fun是用M文件定义的函数或 10求函数f(x)x3)21,x[0,5]的最小值。解编写Mfun5.mfunctionf=(x-3)^2- fseminf函数min{F(x)|C(x)0,Ceq(x)0,PHI(x,w)x

A*xAeq*x其中FUN用于定义目标函数F(x；X0x的初始值；NTHETA是半无穷约束PHI(xwSEMINFCON用于定义非线性不等式约束C(x)输入参量X和S，S是的取样步长，也许不被使用。例11 求函数f(x)(x10.5)2(x20.5)2(x30.5)2取最小值时的x值,K(x,w)sin(wx)cos(wx) (w50)2sin(wx)x 1 1 1 K(x,w)sin(w2x2)cos(w2x) (w50)2sin(wx)x

2 1w1100，1w2解（1）Mfun6.mfunctionf=fun6(x,s);编写Mfun7.mifs=[0.2,0;0.2%半无穷k1=sin(w1*x(1)).*cos(w1*x(2))-1/1000*(w1-50).^2-sin(w1*x(3))-x(3)-k2=sin(w2*x(2)).*cos(w2*x(1))-1/1000*(w2-50).^2-sin(w2*x(3))-x(3)- fminimax函数x A*xAeq*x

C(x)Ceq(x)F(xF1x),LFm(x)}。上述问题的命令为例 求函数族{f1(x),f2(x),f3(x),f4(x),f5(x)}取极大极小值时的x值其中f(x)2x2x248x40x f(x)x23x f3(x)x13x2f4(x)x1f5(x)x1x2解（1）Mfun8.mfunctionf=[2*x(1)^2+x(2)^2-48*x(1)--x(1)^2-x(1)+3*x(2)--x(1)-x(1)+x(2)-（2）调用函数fminimax例 已知函数f(x)ex1(4x22x24xx2x1)，且满足非线性约束1x1x2x1x2xx1求minfx

1 ex1(4x22x24xx8x6x ex1(4x1

1 解（1）Mfun9.mfunction[f,df]=fun9(x);编写Mfun10.mfunction[c,ceq,dc,dceq]=fun10(x);优化工具箱中的optimtool 1优化问题用户图形界面解法示意例 利用例1已经定义好的函数fun1和fun2。在命令窗口运行optimtool，就打后用鼠标点一下“start”按钮，就得到求解结果，再使用“file”菜单下的“ExporttoWorkspace…”选项，把计算结果输出到工作空间中去。§410，000m160km的正方形区域内，经常有若干架飞机作水平行计算（方向角误差不超过0.01度，要求飞机飞行方向角调整的幅度尽量小。1飞行记录数1345D为飞行管理区域的边长；为飞行管理区域，取直角坐标系使其为[0,D][0,D](x0y0为第i架飞机的初始 (xi(tyi(t为第i架飞机在tii为第i

0

为第

根据相对运动的观点在两架飞机i和j的飞行时，可以将飞机i视为不动而飞机j以相对速度vijvjvi(acosjacosi,asinjasiniv2asinji(sinji,cosji

2a

ji

),

ji

ii不妨设ji，此时相对飞行方向角为ij 2

22相对飞行方向(x0x0)2(x0x0)2(y0y0 0 则只要当相对飞行方向角ij02 2记0为调整前第j架飞机相对于第i架飞机的相对速度（矢量）与这两架线（j指向i的矢量）的夹角（以连线矢量为基准，逆时针方向为正，顺时针方向为01()

0mn＝相对速度vmn的幅角－从nm的连线矢量的幅0

ein（注意0表达式中的i表示虚数单位）这里 角度0（m,n1,2,L,6。6 01()0，i,j1,2,L,6，i i 130y0=[1408515550150q=[243236220.5159230xy0=[x0; form=1:6forif %往纯文本文件中写LINGO数据dlmwrite('txt1.txt',b0,'delimiter','\t','newline','PC','-append','roffset',求得0的值如2所示20的12345