代数方程与最优化问题_第1页
代数方程与最优化问题_第2页
代数方程与最优化问题_第3页
代数方程与最优化问题_第4页
代数方程与最优化问题_第5页
已阅读5页,还剩81页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

代数方程与最优化问题第1页,课件共86页,创作于2023年2月7.1代数方程的求解

7.1.1代数方程的图解法一元方程的图解法例:>>ezplot('exp(-3*t)…*sin(4*t+2)+4*exp…(-0.5*t)*cos(2*t)-…0.5',[05])>>holdon,>>line([0,5],[0,0])%同时绘制横轴第2页,课件共86页,创作于2023年2月验证:>>symst;t=3.52028;vpa(exp(-3*t)*sin(4*t+2)+…4*exp(-0.5*t)*cos(2*t)-0.5)ans=-.19256654148425145223200161126442e-4第3页,课件共86页,创作于2023年2月二元方程的图解法例:>>ezplot('x^2*exp(-x*y^2/2)+exp(-x/2)*sin(x*y)')%第一个方程曲线>>holdon%保持当前坐标系>>ezplot(‘x^2*…cos(x+y^2)+…y^2*exp(x+y)')第4页,课件共86页,创作于2023年2月方程的图解法仅适用于一元、二元方程的求根问题。第5页,课件共86页,创作于2023年2月7.1.2多项式型方程的准解析解法例:>>ezplot('x^2+y^2-1');holdon%绘制第一方程的曲线>>ezplot(‘0.75*x^3-y+0.9’)%绘制第二方程为关于x的6次多项式方程应有6对根。图解法只能显示求解方程的实根。第6页,课件共86页,创作于2023年2月一般多项式方程的根可为实数,也可为复数。可用MATLAB符号工具箱中的solve()函数。最简调用格式:

S=solve(eqn1,eqn2,…,eqnn)

(返回一个结构题型变量S,如S.x表示方程的根。)直接得出根:(变量返回到MATLAB工作空间)

[x,…]=solve(eqn1,eqn2,…,eqnn)同上,并指定变量

[x,…]=solve(eqn1,eqn2,…,eqnn,’x,…’)第7页,课件共86页,创作于2023年2月例:>>symsxy;>>[x,y]=solve('x^2+y^2-1=0','75*x^3/100-y+9/10=0')x=[-.98170264842676789676449828873194][-.55395176056834560077984413882735-.35471976465080793456863789934944*i][-.55395176056834560077984413882735+.35471976465080793456863789934944*i][.35696997189122287798839037801365][.86631809883611811016789809418650-1.2153712664671427801318378544391*i][.86631809883611811016789809418650+1.2153712664671427801318378544391*i]y=[.19042035099187730240977756415289][.92933830226674362852985276677202-.21143822185895923615623381762210*i][.92933830226674362852985276677202+.21143822185895923615623381762210*i][.93411585960628007548796029415446][-1.4916064075658223174787216959259-.70588200721402267753918827138837*i][-1.4916064075658223174787216959259+.70588200721402267753918827138837*i]第8页,课件共86页,创作于2023年2月验证>>[eval('x.^2+y.^2-1')eval('75*x.^3/100-y+9/10')]ans=[0,0][0,0][0,0][-.1e-31,0][.5e-30+.1e-30*i,0][.5e-30-.1e-30*i,0]

由于方程阶次可能太高,不存在解析解。然而,可利用MATLAB的符号工具箱得出原始问题的高精度数值解,故称之为准解析解。第9页,课件共86页,创作于2023年2月例:>>[x,y,z]=solve('x+3*y^3+2*z^2=1/2','x^2+3*y+z^3=2','x^3+2*z+2*y^2=2/4');>>size(x)ans=271>>x(22),y(22),z(22)ans=-1.0869654762986136074917644096117ans=.37264668450644375527750811296627e-1ans=.89073290972562790151300874796949第10页,课件共86页,创作于2023年2月验证:>>err=[x+3*y.^3+2*z.^2-1/2,x.^2+3*y+z.^3-2,x.^3+2*z+2*y.^2-2/4];>>norm(double(eval(err)))ans=1.4998e-027多项式乘积形式也可,如把第三个方程替换一下。>>[x,y,z]=solve('x+3*y^3+2*z^2=1/2','x^2+3*y+z^3=2','x^3+2*z*y^2=2/4');>>err=[x+3*y.^3+2*z.^2-1/2,x.^2+3*y+z.^3-2,x.^3+2*z.*y.^2-2/4];>>norm(double(eval(err)))%将解代入求误差ans=5.4882e-028第11页,课件共86页,创作于2023年2月例:求解(含变量倒数)>>symsxy;>>[x,y]=solve('x^2/2+x+3/2+2/y+5/(2*y^2)+3/x^3=0',...'y/2+3/(2*x)+1/x^4+5*y^4','x,y');>>size(x)ans=261>>err=[x.^2/2+x+3/2+2./y+5./(2*y.^2)+3./x.^3,y/2+3./(2*x)+1./x.^4+5*y.^4];%验证>>norm(double(eval(err)))ans=8.9625e-030第12页,课件共86页,创作于2023年2月例:求解(带参数方程)>>symsabxy;>>[x,y]=solve('x^2+a*x^2+6*b+3*y^2=0','y=a+(x+3)','x,y')x=[1/2/(4+a)*(-6*a-18+2*(-21*a^2-45*a-27-24*b-6*a*b-3*a^3)^(1/2))][1/2/(4+a)*(-6*a-18-2*(-21*a^2-45*a-27-24*b-6*a*b-3*a^3)^(1/2))]y=[a+1/2/(4+a)*(-6*a-18+2*(-21*a^2-45*a-27-24*b-6*a*b-3*a^3)^(1/2))+3][a+1/2/(4+a)*(-6*a-18-2*(-21*a^2-45*a-27-24*b-6*a*b-3*a^3)^(1/2))+3]第13页,课件共86页,创作于2023年2月7.1.3一般非线性方程数值解非线性方程的标准形式为f(x)=0函数

fzero格式x=fzero(fun,x0)%用fun定义表达式f(x),x0为初始解。

x=fzero(fun,x0,options)

[x,fval]=fzero(…)%fval=f(x)

[x,fval,exitflag]=fzero(…)

[x,fval,exitflag,output]=fzero(…)

说明该函数采用数值解求方程f(x)=0的根。第14页,课件共86页,创作于2023年2月例:求的根解:>>fun='x^3-2*x-5';>>z=fzero(fun,2)%初始估计值为2z=2.0946>>formatlong>>opt=optimset('Tolx',1.0e-8);>>y=fzero('x^3-2*x-5',2,opt)y=2.09455148211709第15页,课件共86页,创作于2023年2月7.1.4一般非线性方程组数值解格式:最简求解语句

x=fsolve(Fun,x0)一般求解语句

[x,f,flag,out]=fsolve(Fun,x0,opt,p1,p2,…)若返回的flag大于0,则表示求解成功,否则求解出现问题,opt求解控制参数,结构体数据。获得默认的常用变量opt=optimset;用这两种方法修改参数(解误差控制量)opt.Tolx=1e-10;或set(opt.‘Tolx’,1e-10)第16页,课件共86页,创作于2023年2月例:自编函数:functionq=my2deq(p)q=[p(1)*p(1)+p(2)*p(2)-1;0.75*p(1)^3-p(2)+0.9];>>OPT=optimset;OPT.LargeScale='off';>>[x,Y,c,d]=fsolve('my2deq',[1;2],OPT)Optimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.x=0.35700.9341Y=1.0e-009*0.12150.0964第17页,课件共86页,创作于2023年2月c=1d=iterations:7funcCount:21algorithm:'trust-regiondogleg'firstorderopt:1.3061e-010%解回代的精度调用inline()函数:>>f=inline('[p(1)*p(1)+p(2)*p(2)-1;0.75*p(1)^3-p(2)+0.9]','p');>>[x,Y]=fsolve(f,[1;2],OPT);%结果和上述完全一致,从略。Optimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.若改变初始值x0=[-1,0]T第18页,课件共86页,创作于2023年2月>>[x,Y,c,d]=fsolve(f,[-1,0]',OPT);x,Y,kk=d.funcCountOptimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.x=-0.98170.1904Y=1.0e-010*0.5619-0.4500kk=15初值改变有可能得出另外一组解。故初值的选择对解的影响很大,在某些初值下甚至无法搜索到方程的解。第19页,课件共86页,创作于2023年2月例:用solve()函数求近似解析解>>symst;solve(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5)ans=.67374570500134756702960220427474%不允许手工选择初值,只能获得这样的一个解。可先用用图解法选取初值,再调用fsolve()函数数值计算第20页,课件共86页,创作于2023年2月>>formatlong>>y=inline('exp(-3*t).*sin(4*t+2)+4*exp(-0.5*t).*cos(2*t)-0.5','t');>>ff=optimset;ff.Display='iter';[t,f]=fsolve(y,3.5203,ff)NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius121.8634e-0095.16e-0051243.67694e-0193.61071e-0057.25e-0101Optimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.t=3.52026389294877f=-6.063776702980306e-010第21页,课件共86页,创作于2023年2月重新设置相关的控制变量,提高精度。>>ff=optimset;ff.TolX=1e-16;ff.TolFun=1e-30;>>ff.Display='iter';[t,f]=fsolve(y,3.5203,ff)NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius121.8634e-0095.16e-0051243.67694e-0193.61071e-0057.25e-01013605.07218e-01001Optimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.t=3.52026389244155f=0第22页,课件共86页,创作于2023年2月例:求的根解:化为标准形式设初值点为x0=[-5-5]。functionF=myfun(x)F=[2*x(1)-x(2)-exp(-x(1));-x(1)+2*x(2)-exp(-x(2))];>>x0=[-5;-5];%初始点>>options=optimset('Display','iter');%显示输出信息>>[x,fval]=fsolve(@myfun,x0,options)

第23页,课件共86页,创作于2023年2月NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius0347071.22.29e+00411612003.415.75e+0031293147.0211.47e+0031312854.45213881415239.5271107151867.0412130.8162116.704219.0517242.4278812.2618270.0326580.7595110.2062.59307.03149e-0060.1119270.002942.510333.29525e-0130.001691326.36e-0072.5Optimizationterminated:first-orderoptimalityislessthanoptions.TolFun.x=0.567143031397360.56714303139736fval=1.0e-006*-0.40590960570519-0.40590960570519第24页,课件共86页,创作于2023年2月例:求矩阵x使其满足方程,并设初始解向量为x=[1,1;1,1]。解:functionF=myfun(x)F=x*x*x-[1,2;3,4];>>x0=ones(2,2);%初始解向量>>options=optimset('Display','off');%不显示优化信息>>[x,Fval,exitflag]=fsolve(@myfun,x0,options)x=-0.12910.86021.29031.1612Fval=1.0e-003*0.1541-0.11630.0109-0.0243exitflag=1第25页,课件共86页,创作于2023年2月7.2无约束最优化问题求解

7.2.1解析解法和图解法第26页,课件共86页,创作于2023年2月例:>>symst;y=exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5;>>y1=diff(y,t);%求取一阶导函数>>ezplot(y1,[0,4])%绘制出选定区间内一阶导函数曲线第27页,课件共86页,创作于2023年2月>>t0=solve(y1)%求出一阶导数等于零的点t0=1.4528424981725411893375778048840>>y2=diff(y1);b=subs(y2,t,t0)%并验证二阶导数为正b=7.8553420253333601379464405534591>>t=0:0.4:4;y=exp(-3*t).*sin(4*t+2)+4*exp(-0.5*t).*cos(2*t)-0.5;>>plot(t,y)第28页,课件共86页,创作于2023年2月单变量函数求最小值的(数值解法)%求解区间优化问题标准形式为:s.t.函数

fminbnd格式x=fminbnd(fun,x1,x2)x=fminbnd(fun,x1,x2,options)[x,fval]=fminbnd(…)%fval为目标函数的最小值[x,fval,exitflag]=fminbnd(…)

exitflag为终止迭代的条件,若exitflag>0,收敛于x,exitflag=0,表示超过函数估计值或迭代的最大数字,exitflag<0表示函数不收敛于x;

[x,fval,exitflag,output]=fminbnd(…)

output为优化信息,若output=iterations表示迭代次数,output=funccount表示函数赋值次数,output=algorithm表示所使用的算法

第29页,课件共86页,创作于2023年2月例:

计算下面函数在区间(0,1)内的最小值。解:>>[x,fval,exitflag,output]=fminbnd('(x^3+cos(x)+x*log(x))/exp(x)',0,1)x=0.5223fval=0.3974exitflag=1output=iterations:9funcCount:11algorithm:'goldensectionsearch,parabolicinterpolation'第30页,课件共86页,创作于2023年2月例:在[0,5]上求下面函数的最小值

解:

functionf=myfun(x)f=(x-3).^2-1;>>x=fminbnd(@myfun,0,5)x=3第31页,课件共86页,创作于2023年2月7.2.2基于MATLAB的数值解法第32页,课件共86页,创作于2023年2月例:>>f=inline('(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2))','x');>>x0=[0;0];ff=optimset;ff.Display='iter';>>x=fminsearch(f,x0,ff)IterationFunc-countminf(x)Procedure13-0.000499937initial24-0.000499937reflect

72137-0.641424contractoutsideOptimizationterminatedsuccessfully:x=0.6111-0.3056第33页,课件共86页,创作于2023年2月

>>x=fminunc(f,[0;.0],ff)DirectionalIterationFunc-countf(x)Step-sizederivative12-2e-0080.001-429-0.5846690.3043530.343316-0.6413970.9503220.00191422-0.6414240.984028-1.45e-008x=0.6110-0.3055比较可知fminunc()函数效率高于fminsearch()函数,但当所选函数高度不连续时,使用fminsearch效果较好。故若安装了最优化工具箱则应调用fminunc()函数。第34页,课件共86页,创作于2023年2月7.2.3全局最优解与局部最优解例:>>f=inline('exp(-2*t).*cos(10*t)+exp(-3*(t+2)).*sin(2*t)','t');%目标函数>>t0=1;[t1,f1]=fminsearch(f,t0);[t1f1]ans=0.9228-0.1547>>t0=0.1;[t2,f2]=fminsearch(f,t0);[t2f2]ans=0.2945-0.5436第35页,课件共86页,创作于2023年2月>>symst;y=exp(-2*t)*cos(10*t)+exp(-3*(t+2))*sin(2*t);>>ezplot(y,[0,2.5]);set(gca,‘Ylim’,[-0.6,1])%在t[0,2.5]内的曲线>>ezplot(y,[-0.5,2.5]);set(gca,‘Ylim’,[-2,1.2])%在[-0.5,2.5]曲线>>t0=-0.2;[t3,f3]=fminsearch(f,t0);[t3f3]ans=-0.3340-1.9163第36页,课件共86页,创作于2023年2月7.2.4利用梯度求解最优化问题例:>>[x,y]=meshgrid…(0.5:0.01:1.5);…z=100*(y.^2-x).^2…+(1-x).^2;>>contour3(x,y,z,100),set(gca,'zlim',[0,310])%测试算法的函数第37页,课件共86页,创作于2023年2月>>f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2','x');>>ff=optimset;ff.TolX=1e-10;ff.TolFun=1e-20;ff.Display='iter';>>x=fminunc(f,[0;0],ff)Warning:Gradientmustbeprovidedfortrust-regionmethod;usingline-searchmethodinstead.DirectionalIterationFunc-countf(x)Step-sizederivative1210.5-4442023.01749e-0123.40397e-009-1.77e-013x=1.000001736959721.00000347608069第38页,课件共86页,创作于2023年2月%求梯度向量(比较引入梯度的作用,但梯度的计算也费时间)>>symsx1x2;f=100*(x2-x1^2)^2+(1-x1)^2;>>J=jacobian(f,[x1,x2])J=[-400*(x2-x1^2)*x1-2+2*x1,200*x2-200*x1^2]function[y,Gy]=c6fun3(x)%目标函数编写:y=100*(x(2)-x(1)^2)^2+(1-x(1))^2;%目标函数Gy=[-400*(x(2)-x(1)^2)*x(1)-2+2*x(1);200*x(2)-200*x(1)^2];%梯度>>ff.GradObj='on';x=fminunc('c6fun3',[0;0],ff)NormofFirst-orderIterationf(x)stepoptimalityCG-iterations196.38977e-0292.12977e-0121.6e-0141x=0.999999999999990.99999999999998第39页,课件共86页,创作于2023年2月7.2.5非线性最小二乘

函数

lsqnonlin格式x=lsqnonlin(fun,x0)x=lsqnonlin(fun,x0,lb,ub)x=lsqnonlin(fun,x0,lb,ub,options)%x0为初始解向量;fun为,i=1,2,…,m,

lb、ub定义x的下界和上界,options为指定优化参数,若x没有界,则lb=[],ub=[]。第40页,课件共86页,创作于2023年2月例:求下面非线性最小二乘问题初始解向量为x0=[0.3,0.4]。解:先建立函数文件,并保存为myfun.m.

由于lsqnonlin中的fun为向量形式而不是平方和形式,因此,myfun函数应由建立:

k=1,2,…,10

第41页,课件共86页,创作于2023年2月functionF=myfun10(x)k=1:10;F=2+2*k-exp(k*x(1))-exp(k*x(2));然后调用优化程序:

>>x0=[0.30.4];>>[x,resnorm]=lsqnonlin(@myfun10,x0)Optimizationterminated:normofthecurrentstepislessthanOPTIONS.TolX.x=0.25780.2578resnorm=124.3622第42页,课件共86页,创作于2023年2月7.3有约束最优化问题的计算机求解

7.3.1约束条件与可行解区域有约束最优化问题的一般描述:对于一般的一元问题和二元问题,可用图解法直接得出问题的最优解。第43页,课件共86页,创作于2023年2月例:用图解方法求解:>>[x1,x2]=meshgrid(-3:.1:3);%生成网格型矩阵>>z=-x1.^2-x2;%计算出矩阵上各点的高度>>i=find(x1.^2+x2.^2>9);z(i)=NaN;%找出x1^2+x2^2>9的坐标,并置函数值为NaN>>i=find(x1+x2>1);z(i)=NaN;%找出x1+x2>1的坐标,置为NaN>>surf(x1,x2,z);shadinginterp;>>max(z(:)),view(0,90)ans=3第44页,课件共86页,创作于2023年2月7.3.2线性规划问题的计算机求解第45页,课件共86页,创作于2023年2月例:求解>>f=-[21431]';A=[02142;345-1-1];>>B=[54;62];Ae=[];Be=[];xm=[0,0,3.32,0.678,2.57];>>ff=optimset;ff.LargeScale='off';%不使用大规模问题求解>>ff.TolX=1e-15;ff.TolFun=1e-20;ff.TolCon=1e-20;ff.Display='iter';>>[x,f_opt,key,c]=linprog(f,A,B,Ae,Be,xm,[],[],ff)第46页,课件共86页,创作于2023年2月Optimizationterminatedsuccessfully.x=19.78500.00003.320011.38502.5700f_opt=-89.5750key=1%求解成功c=iterations:5algorithm:'medium-scale:activeset'cgiterations:[]第47页,课件共86页,创作于2023年2月例:求解>>f=[-3/4,150,-1/50,6]';Aeq=[];Beq=[];>>A=[1/4,-60,-1/50,9;1/2,-90,-1/50,3];B=[0;0];>>xm=[-5;-5;-5;-5];xM=[Inf;Inf;1;Inf];ff=optimset;>>ff.TolX=1e-15;ff.TolFun=1e-20;ff.TolCon=1e-20;ff.Display='iter';>>[x,f_opt,key,c]=linprog(f,A,B,Aeq,Beq,xm,xM,[0;0;0;0],ff)第48页,课件共86页,创作于2023年2月Residuals:PrimalDualUpperDualityTotalInfeasInfeasBoundsGapRelA*x-bA'*y+z-w-f{x}+s-ubx'*z+s'*wError-------------------------------------------------------------Iter0:9.39e+0031.43e+0021.94e+0026.03e+0042.77e+001Iter10:0.00e+0006.15e-0260.00e+0001.70e-0244.10e-028Optimizationterminatedsuccessfully.x=-5.0000-0.19471.0000-5.0000f_opt=-55.4700key=1c=iterations:10cgiterations:0algorithm:'lipsol'第49页,课件共86页,创作于2023年2月7.3.3二次型规划的求解(线性)第50页,课件共86页,创作于2023年2月例:求解>>H=diag([2,2,2,2]);f=[-2,-4,-6,-8];>>OPT=optimset;OPT.LargeScale='off';%不使用大规模问题算法>>A=[1,1,1,1;3,3,2,1];B=[5;10];Aeq=[];Beq=[];LB=zeros(4,1);>>[x,f_opt]=quadprog(H,f,A,B,Aeq,Beq,LB,[],[],OPT)Optimizationterminatedsuccessfully.x=0.00000.66671.66672.6667f_opt=-23.6667%(解的目标值应为30+(-23.6667)=6.3333)第51页,课件共86页,创作于2023年2月例:求解下面二次规划问题解:则

第52页,课件共86页,创作于2023年2月>>H=[1-1;-12];f=[-2;-6];>>A=[11;-12;21];b=[2;2;3];>>lb=zeros(2,1);>>[x,fval,exitflag,output,lambda]=quadprog(H,f,A,b,[],[],lb)x=%最优解0.66671.3333fval=%最优值-8.2222exitflag=%收敛1output=iterations:3algorithm:'medium-scale:active-set'firstorderopt:[]cgiterations:[]第53页,课件共86页,创作于2023年2月例:求二次规划的最优解maxf(x1,x2)=x1x2+3sub.tox1+x2-2=0解:化成标准形式:>>H=[0,-1;-1,0];f=[0;0];Aeq=[11];>>[x,fval,exitflag,output]=quadprog(H,f,[],[],Aeq)第54页,课件共86页,创作于2023年2月6.3.4一般非线性规划问题的求解第55页,课件共86页,创作于2023年2月例:求解目标函数:functiony=opt_fun1(x)y=1000-x(1)*x(1)-2*x(2)*x(2)-x(3)*x(3)-x(1)*x(2)-x(1)*x(3);非线性约束函数:function[c,ceq]=opt_con1(x)ceq=[x(1)*x(1)+x(2)*x(2)+x(3)*x(3)-25;8*x(1)+14*x(2)+7*x(3)-56];c=[];第56页,课件共86页,创作于2023年2月>>ff=optimset;ff.LargeScale='off';ff.Display='iter';>>ff.TolFun=1e-30;ff.TolX=1e-15;ff.TolCon=1e-20;>>x0=[1;1;1];xm=[0;0;0];xM=[];A=[];B=[];Aeq=[];Beq=[];>>[x,f_opt,c,d]=fmincon('opt_fun1',x0,A,B,Aeq,Beq,xm,xM,'opt_con1',ff);>>x,f_opt,kk=d.funcCountx=3.51210.21703.5522f_opt=961.7152kk=%目标函数调用的次数108第57页,课件共86页,创作于2023年2月简化非线约束函数function[c,ceq]=opt_con2(x)ceq=x(1)*x(1)+x(2)*x(2)+x(3)*x(3)-25;c=[];求解:>>x0=[1;1;1];Aeq=[8,14,7];Beq=56;>>[x,f_opt,c,d]=fmincon('opt_fun1',x0,A,B,Aeq,Beq,xm,xM,'opt_con2',ff);maxDirectionalFirst-orderIterF-countf(x)constraintStep-sizederivativeoptimalityProcedure111955.33622.90.25-29518.3infeasible21116961.71501-6.3e-0156.97e-005HessianmodifiedtwiceOptimizationterminatedsuccessfully:Searchdirectionlessthan2*options.TolXandmaximumconstraintviolationislessthanoptions.TolConActiveConstraints:12>>x,f_opt,kk=d.funcCount%从略(计算结果同上一样)第58页,课件共86页,创作于2023年2月例:利用梯度求解梯度函数:>>symsx1x2x3;f=1000-x1*x1-2*x2*x2-x3*x3-x1*x2-x1*x3;>>J=jacobian(f,[x1,x2,x3])J=[-2*x1-x2-x3,-4*x2-x1,-2*x3-x1]改写目标函数:function[y,Gy]=opt_fun2(x)y=1000-x(1)*x(1)-2*x(2)*x(2)-x(3)*x(3)-x(1)*x(2)-x(1)*x(3);Gy=-[2*x(1)+x(2)+x(3);4*x(2)+x(1);2*x(3)+x(1)];第59页,课件共86页,创作于2023年2月>>x0=[1;1;1];xm=[0;0;0];xM=[];A=[];B=[];Aeq=[];Beq=[];>>ff=optimset;ff.GradObj=‘on’;%若知道梯度函数ff.Display='iter';ff.LargeScale='off';>>ff.TolFun=1e-30;ff.TolX=1e-15;ff.TolCon=1e-20;>>[x,f_opt,c,d]=fmincon('opt_fun2',x0,A,B,Aeq,Beq,xm,xM,'opt_con1',ff);>>x,f_opt,kk=d.funcCountx=3.51210.21703.5522f_opt=961.7152kk=95第60页,课件共86页,创作于2023年2月6.3.5约束线性最小二乘

有约束线性最小二乘的标准形式为:其中:C、A、Aeq为矩阵;d、b、beq、lb、ub、x是向量。第61页,课件共86页,创作于2023年2月函数lsqlin格式x=lsqlin(C,d,A,b)x=lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0,options)%求在约束条件下,方程Cx=d的最小二乘解x。Aeq、beq满足等式约束,若没有不等式约束,则设A=[],b=[]。lb、ub满足,若没有等式约束,则Aeq=[],beq=[]。x0为初始解向量,若x没有界,则lb=[],ub=[]。options为指定优化参数

[x,resnorm,residual,exitflag,output,lambda]=lsqlin(…)%resnorm=norm(C*x-d)^2,即2-范数。residual=C*x-d,即残差。exitflag为终止迭代的条件output表示输出优化信息lambda为解x的Lagrange乘子第62页,课件共86页,创作于2023年2月例:求解下面系统的最小二乘解系统:Cx=d约束:先输入系统系数和x的上下界:C=[0.95010.76200.61530.4057;…0.23110.45640.79190.9354;…0.60680.01850.92180.9169;…0.48590.82140.73820.4102;…0.89120.44470.17620.8936];d=[0.0578;0.3528;0.8131;0.0098;0.1388];A=[0.20270.27210.74670.4659;…0.19870.19880.44500.4186;…0.60370.01520.93180.8462];第63页,课件共86页,创作于2023年2月b=[0.5251;0.2026;0.6721];lb=-0.1*ones(4,1);ub=2*ones(4,1);然后调用最小二乘命令:

>>[x,resnorm,residual,exitflag,output,lambda]=lsqlin(C,d,A,b,[],[],lb,ub)x=-0.1000-0.10000.21520.3502resnorm=0.1672第64页,课件共86页,创作于2023年2月residual=0.04550.0764-0.35620.16200.0784exitflag=1%说明解x是收敛的output=iterations:4algorithm:'medium-scale:active-set'firstorderopt:[]cgiterations:[]第65页,课件共86页,创作于2023年2月lambda=lower:[4x1double]upper:[4x1double]eqlin:[0x1double]ineqlin:[3x1double]%通过lambda.ineqlin可查看不等式约束是否有效。>>lambda.ineqlinans=00.23920第66页,课件共86页,创作于2023年2月7.4整数规划问题的计算机求解7.4.1整数线性规划问题的求解(matlab2009a有问题)A、B线性等式和不等式约束,,约束式子由ctype变量控制,intlist为整数约束标示,how=0表示结果最优,2为无可行解,其余失败。第67页,课件共86页,创作于2023年2月例:>>f=-[21431]';A=[02142;345-1-1];intlist=ones(5,1);>>B=[54;62];ctype=[-1;-1];xm=[0,0,3.32,0.678,2.57];xM=inf*ones(5,1);>>[res,b]=ipslv_mex(f,A,B,intlist,xM,xm,ctype)res=1904105b=%因为返回的b=0,表示优化成功0第68页,课件共86页,创作于2023年2月>>[x1,x2,x3,x4,x5]=ndgrid(1:20,0:20,4:20,1:20,3:20);>>i=find((2*x2+x3+4*x4+2*x5<=54)&(3*x1+4*x2+5*x3-x4-x5<=62));>>f=-2*x1(i)-x2(i)-4*x3(i)-3*x4(i)-x5(i);[fmin,ii]=sort(f);%标号>>index=i(ii(1));x=[x1(index),x2(index),x3(index),x4(index),x5(index)]x=1904105%当把20换为30,一般计算机配置下实现不了,故求解整数规划时不适合采用穷举算法。第69页,课件共86页,创作于2023年2月次最优解>>fmin(1:10)'ans=-89-88-88-88-88-88-88-88-87-87>>in=i(ii(1:4));x=[x1(in),x2(in),x3(in),x4(in),x5(in),fmin(1:4)]x=1904105-891804113-881705104-881506104-88>>intlist=[1,0,0,1,1];%混合整数规划>>[res,b]=ipslv_mex(f,A,B,intlist,xM,xm,ctype)res=19.000003.800011.00003.0000b=0第70页,课件共86页,创作于2023年2月7.4.2一般非线性整数规划问题与求解(matlab2009a有问题)第71页,课件共86页,创作于2023年2月err字符串为空,则返回最优解。该函数尚有不完全之处,解往往不是精确整数,可用下面语句将其化成整数。第72页,课件共86页,创作于2023年2月例:functionf=c6miopt(x)f=-[21431]*x;>>A=[02142;345-1-1];intlist=ones(5,1);Aeq=[];Beq=[];>>B=[54;62];ctype=[-1;-1];>>xm=[0,0,4,1,3]';xM=20000*ones(5,1);x0=xm;>>[errmsg,f,X]=bnb20('c6miopt',x0,intlist,xm,xM,A,B,Aeq,Beq);X=X'X=19.000004.000010.00005.0000第73页,课件共86页,创作于2023年2月>>intlist=[1,0,0,1,1]';>>xm=[0,0,3.32,1,3]';>>[errmsg,f,X]=bnb20('c6miopt',x0,intlist,xm,xM,A,B,Aeq,Beq);XX=19.000003.800011.00003.0000第74页,课件共86页,创作于2023年2月例:编写函数:functiony=c7fun1(x)y=100*(x(2)-x(1)^2)^2+(4.5543-x(1))^2;>>x0=[1;1];xm=-100*[1;1];xM=100*[1;1];>>A=[];B=[];Aeq=[];Beq=[];intlist=[1,1]';>>[errmsg,f,x]=bnb20('c6fun1',x0,intlist,xm,xM,A,B,Aeq,Beq);x'ans=5.0000000000000025.00000001055334第75页,课件共86页,创作于2023年2月>>iflength(errmsg)==0,x=round(x),end;x=x';x=525%缩小范围。>>xm=-20*[1;1];

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论