《matlab 代数方程与最优化问题的求解.pptx》由会员分享，可在线阅读，更多相关《matlab 代数方程与最优化问题的求解.pptx（60页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、代数方程的求解1 1 代数方程的图解法一元方程的图解法例：ezplot(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5,0 5)hold on,line(0,5,0,0)%同时绘制横轴第1页/共60页验证：syms t;t=3.52028;vpa(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5)ans=-.19256654148425145223200161126442e-4第2页/共60页二元方程的图解法例：ezplot(x2*exp(-x*y2/2)+exp(-x/2)*sin(x*y)%第一个
2、方程曲线 hold on%保护当前坐标系 ezplot(x2*cos(x+y2)+y2*exp(x+y)第3页/共60页方程的图解法仅适用于一元、二元方程的求根问题。第4页/共60页2 2 多项式型方程的准解析解法例：ezplot(x2+y2-1);hold on%绘制第一方程的曲线 ezplot(0.75*x3-y+0.9)%绘制第二方程为关于x的6次多项式方程应有6对根。图解法只能显示求解方程的实根。第5页/共60页一般多项式方程的根可为实数，也可为复数。可用MATLAB符号工具箱中的solve()函数。最简调用格式：S=solve(eqn1,eqn2,eqnn)（返回一个结构题型变量S，
3、如S.x表示方程的根。）直接得出根：（变量返回到MATLAB工作空间）x,=solve(eqn1,eqn2,eqnn)同上，并指定变量 x,=solve(eqn1,eqn2,eqnn，x,)第6页/共60页例：syms x y;x,y=solve(x2+y2-1=0,75*x3/100-y+9/10=0)x=-.98170264842676789676449828873194-.55395176056834560077984413882735-.35471976465080793456863789934944*i-.55395176056834560077984413882735+.35471
4、976465080793456863789934944*i .35696997189122287798839037801365 .86631809883611811016789809418650-1.2153712664671427801318378544391*i .86631809883611811016789809418650+1.2153712664671427801318378544391*iy=.19042035099187730240977756415289 .92933830226674362852985276677202-.21143822185895923615623381
5、762210*i .92933830226674362852985276677202+.21143822185895923615623381762210*i .93411585960628007548796029415446-1.4916064075658223174787216959259-.70588200721402267753918827138837*i-1.4916064075658223174787216959259+.70588200721402267753918827138837*i第7页/共60页验证 eval(x.2+y.2-1)eval(75*x.3/100-y+9/10
6、)ans=0,0 0,0 0,0 -.1e-31,0.5e-30+.1e-30*i,0.5e-30-.1e-30*i,0 由于方程阶次可能太高，不存在解析解。然而，可利用MATLAB的符号工具箱得出原始问题的高精度数值解，故称之为准解析解。第8页/共60页例：x,y,z=solve(x+3*y3+2*z2=1/2,x2+3*y+z3=2 ,x3+2*z+2*y2=2/4);size(x)ans=27 1 x(22),y(22),z(22)ans=-1.0869654762986136074917644096117 ans=.37264668450644375527750811296627e-1
7、ans=.89073290972562790151300874796949 第9页/共60页验证：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*y3+2*z2=1/2,x2+3*y+z3=2,x3+2*z*y2=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)%将解代入求误差an
8、s=5.4882e-028第10页/共60页例：求解（含变量倒数）syms x y;x,y=solve(x2/2+x+3/2+2/y+5/(2*y2)+3/x3=0,.y/2+3/(2*x)+1/x4+5*y4,x,y);size(x)ans=26 1 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第11页/共60页例：求解（带参数方程）syms a b x y;x,y=solve(x2+a*x2+6*b+3*y2=0,y=a+(x+
9、3),x,y)x=1/2/(4+a)*(-6*a-18+2*(-21*a2-45*a-27-24*b-6*a*b-3*a3)(1/2)1/2/(4+a)*(-6*a-18-2*(-21*a2-45*a-27-24*b-6*a*b-3*a3)(1/2)y=a+1/2/(4+a)*(-6*a-18+2*(-21*a2-45*a-27-24*b-6*a*b-3*a3)(1/2)+3 a+1/2/(4+a)*(-6*a-18-2*(-21*a2-45*a-27-24*b-6*a*b-3*a3)(1/2)+3第12页/共60页3 3 一般非线性方程数值解格式：最简求解语句 x=fsolve(Fun,x0
10、)一般求解语句 x,f,flag,out=fsolve(Fun,x0,opt,p1,p2,)若返回的flag 大于0，则表示求解成功，否则求解出现问题，opt 求解控制参数，结构体数据。获得默认的常用变量 opt=optimset;用这两种方法修改参数（解误差控制量）opt.Tolx=1e-10;或 set(opt.Tolx,1e-10)第13页/共60页例：自编函数：function q=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(my2
11、deq,1;2,OPT)Optimization terminated successfully:First-order optimality is less than options.TolFun.x=0.3570 0.9341Y=1.0e-009*0.1215 0.0964第14页/共60页c=1d=iterations:7 funcCount:21 algorithm:trust-region dogleg firstorderopt:1.3061e-010 解回代的精度调用inline()函数：f=inline(p(1)*p(1)+p(2)*p(2)-1;0.75*p(1)3-p(2)
12、+0.9,p);x,Y=fsolve(f,1;2,OPT);%结果和上述完全一致，从略。Optimization terminated successfully:First-order optimality is less than options.TolFun.若改变初始值 x0=-1,0T第15页/共60页 x,Y,c,d=fsolve(f,-1,0,OPT);x,Y,kk=d.funcCountOptimization terminated successfully:First-order optimality is less than options.TolFun.x=-0.9817
13、0.1904Y=1.0e-010*0.5619 -0.4500kk=15 初值改变有可能得出另外一组解。故初值的选择对解的影响很大，在某些初值下甚至无法搜索到方程的解。第16页/共60页例：用solve()函数求近似解析解 syms t;solve(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5)ans=.67374570500134756702960220427474不允许手工选择初值，只能获得这样的一个解。可先用用图解法选取初值，再调用fsolve()函数数值计算第17页/共60页 format long y=inline(exp(-3*t)
14、.*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)Norm of First-order Trust-region Iteration Func-count f(x)step optimality radius 1 2 1.8634e-009 5.16e-005 1 2 4 3.67694e-019 3.61071e-005 7.25e-010 1Optimization terminated successfully:First-order optima
15、lity is less than options.TolFun.t=3.52026389294877f=-6.063776702980306e-010第18页/共60页重新设置相关的控制变量，提高精度。ff=optimset;ff.TolX=1e-16;ff.TolFun=1e-30;ff.Display=iter;t,f=fsolve(y,3.5203,ff)Norm of First-order Trust-region Iteration Func-count f(x)step optimality radius 1 2 1.8634e-009 5.16e-005 1 2 4 3.67
16、694e-019 3.61071e-005 7.25e-010 1 3 6 0 5.07218e-010 0 1Optimization terminated successfully:First-order optimality is less than options.TolFun.t=3.52026389244155f=0第19页/共60页无约束最优化问题求解 1 1 解析解法和图解法第20页/共60页例：syms t;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
17、)%绘制出选定区间内一阶导函数曲线第21页/共60页 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)第22页/共60页2 2 基于 MATLAB MATLAB 的数值解法第23页/共60页例：f=inline(x(1)2-2*x(1)
18、*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)Iteration Func-count min f(x)Procedure 1 3 -0.000499937 initial 2 4 -0.000499937 reflect 72 137 -0.641424 contract outsideOptimization terminated successfully:x=0.6111 -0.3056第24页/共60页 x=fminunc(f,0;0,ff)Direction
19、al Iteration Func-count f(x)Step-size derivative 1 2 -2e-008 0.001 -4 2 9 -0.584669 0.304353 0.343 3 16 -0.641397 0.950322 0.00191 4 22 -0.641424 0.984028 -1.45e-008 x=0.6110 -0.3055 比较可知 fminunc()函数效率高于fminsearch()函数，故若安装了最优化工具箱则应调用fminunc()函数。第25页/共60页3 3 全局最优解与局部最优解例：f=inline(exp(-2*t).*cos(10*t)
20、+exp(-3*(t+2).*sin(2*t),t);%目标函数 t0=1;t1,f1=fminsearch(f,t0);t1 f1ans=0.9228 -0.1547 t0=0.1;t2,f2=fminsearch(f,t0);t2 f2ans=0.2945 -0.5436第26页/共60页 syms t;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)在t0,2.5内的曲线 ezplot(y,-0.5,2.5);set(gca,Ylim,-2,1.2)在-0.5,2.5曲线 t0
21、=-0.2;t3,f3=fminsearch(f,t0);t3 f3ans=-0.3340 -1.9163第27页/共60页4 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)%测试算法的函数第28页/共60页 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,
22、ff)Warning:Gradient must be provided for trust-region method;using line-search method instead.Directional Iteration Func-count f(x)Step-size derivative 1 2 1 0.5 -4 44 202 3.01749e-012 3.40397e-009 -1.77e-013 x=1.00000173695972 1.00000347608069第29页/共60页求梯度向量(比较引入梯度的作用，但梯度的计算也费时间）syms x1 x2;f=100*(x2
23、-x12)2+(1-x1)2;J=jacobian(f,x1,x2)J=-400*(x2-x12)*x1-2+2*x1,200*x2-200*x12function 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)Norm of First-order Iteration f(x)step optimality CG-iterations19
24、 6.38977e-029 2.12977e-012 1.6e-014 1x=0.99999999999999 0.99999999999998第30页/共60页有约束最优化问题的计算机求解 1 1 约束条件与可行解区域有约束最优化问题的一般描述：对于一般的一元问题和二元问题，可用图解法直接得出问题的最优解。第31页/共60页例：用图解方法求解：x1,x2=meshgrid(-3:.1:3);%生成网格型矩阵 z=-x1.2-x2;%计算出矩阵上各点的高度 i=find(x1.2+x2.29);z(i)=NaN;%找出 x12+x229 的坐标，并置函数值为 NaN i=find(x1+x21
25、);z(i)=NaN;%找出 x1+x21的坐标,置为 NaN surf(x1,x2,z);shading interp;max(z(:),view(0,90)ans=3第32页/共60页2 2 线性规划问题的计算机求解第33页/共60页例：求解 f=-2 1 4 3 1;A=0 2 1 4 2;3 4 5-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=ite
26、r;x,f_opt,key,c=linprog(f,A,B,Ae,Be,xm,ff)第34页/共60页Optimization terminated successfully.x=19.7850 0.0000 3.3200 11.3850 2.5700f_opt=-89.5750key=1 求解成功c=iterations:5 algorithm:medium-scale:activeset cgiterations:第35页/共60页例：求解 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;
27、-5;-5;-5;xM=Inf;Inf;1;Inf;ff=optimset;ff.TolX=1e-15;ff.TolFun=1e-20;TolCon=1e-20;ff.Display=iter;x,f_opt,key,c=linprog(f,A,B,Aeq,Beq,xm,xM,0;0;0;0,ff)第36页/共60页Residuals:Primal Dual Upper Duality Total Infeas Infeas Bounds Gap Rel A*x-b A*y+z-w-f x+s-ub x*z+s*w Error -Iter 0:9.39e+003 1.43e+002 1.94e
28、+002 6.03e+004 2.77e+001 Iter 1:6.38e-012 1.21e+001 0.00e+000 3.50e+003 1.78e+000 Iter 10:0.00e+000 6.15e-026 0.00e+000 1.70e-024 4.10e-028Optimization terminated successfully.x=-5.0000 -0.1947 1.0000 -5.0000f_opt=-55.4700key=1c=iterations:10 cgiterations:0 algorithm:lipsol第37页/共60页3 3 二次型规划的求解第38页/
29、共60页例：求解 f=-2,-4,-6,-8;H=diag(2,2,2,2);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)Optimization terminated successfully.x=0.0000 0.6667 1.6667 2.6667f_opt=-23.6667 （解的目标值应为30（-23.6667）6.3333)第39页/共60页4 4 一般非线性规划
30、问题的求解第40页/共60页例：求解目标函数：function y=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=;第41页/共60页 ff=optimset;ff.LargeScale=off;ff.Display=iter;ff.TolFun=1e-30;ff.TolX=1e-15;ff.TolCon=1
31、e-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.5121 0.2170 3.5522f_opt=961.7152kk=%目标函数调用的次数 108第42页/共60页简化非线约束函数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_op
32、t,c,d=fmincon(opt_fun1,x0,A,B,Aeq,Beq,xm,xM,opt_con2,ff);max Directional First-order Iter F-count f(x)constraint Step-size derivative optimality Procedure 1 11 955.336 22.9 0.25 -295 18.3 infeasible21 116 961.715 0 1 -6.3e-015 6.97e-005 Hessian modified twice Optimization terminated successfully:Sea
33、rch direction less than 2*options.TolX and maximum constraint violation is less than options.TolConActive Constraints:1 2 x,f_opt,kk=d.funcCount%从略（计算结果同上一样）第43页/共60页例：利用梯度求解梯度函数：syms x1 x2 x3;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
34、 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);第44页/共60页 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=fminc
35、on(opt_fun2,x0,A,B,Aeq,Beq,xm,xM,opt_con1,ff);x,f_opt,kk=d.funcCountx=3.5121 0.2170 3.5522f_opt=961.7152kk=95 第45页/共60页整数规划问题的计算机求解 1 1整数线性规划问题的求解A、B线性等式和不等式约束，约束式子由ctype变量控制，intlist为整数约束标示，how0表示结果最优，2为无可行解，其余失败。第46页/共60页例：f=-2 1 4 3 1;A=0 2 1 4 2;3 4 5-1-1;intlist=ones(5,1);B=54;62;ctype=-1;-1;xm=
36、0,0,3.32,0.678,2.57;xM=inf*ones(5,1);res,b=ipslv_mex(f,A,B,intlist,xM,xm,ctype)%因为返回的 b=0，表示优化成功res=19 0 4 10 5b=0第47页/共60页 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 f=-2*x1(i)-x2(i)-4*x3(i)-3*x4(i)-x5(i);fmin,ii=sort(f);index=i(ii(1);x=x1(ind
37、ex),x2(index),x3(index),x4(index),x5(index)x=19 0 4 10 5当把20换为30，一般计算机配置下实现不了，故求解整数规划时不适合采用穷举算法。第48页/共60页次最优解 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=19 0 4 10 5 -89 18 0 4 11 3 -88 17 0 5 10 4 -88 15 0 6 10 4 -88 intlist=1,0
38、,0,1,1;混合整数规划 res,b=ipslv_mex(f,A,B,intlist,xM,xm,ctype)res=19.0000 0 3.8000 11.0000 3.0000b=0第49页/共60页2 2一般非线性整数规划问题与求解第50页/共60页err字符串为空，则返回最优解。该函数尚有不完全之处，解往往不是精确整数，可用下面语句将其化成整数。第51页/共60页例：function f=c6miopt(x)f=-2 1 4 3 1*x;A=0 2 1 4 2;3 4 5-1-1;intlist=ones(5,1);Aeq=;Beq=;B=54;62;ctype=-1;-1;xm=0
39、,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=XX=19.0000 0 4.0000 10.0000 5.0000第52页/共60页 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.0000 0 3.8000 11.0000 3.0000第53页/共60页例：编写函数：function y=c7fun1(x)y=100*
40、(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);xans=5.00000000000000 25.00000001055334第54页/共60页 if length(errmsg)=0,x=round(x),end;x=x;x=5 25缩小范围。xm=-20*1;1;xM=20*1;1;errmsg,f,x=bnb20(c6fun1,x0,intlist,xm,xM
41、,A,B,Aeq,Beq);xans=4 16扩大范围用穷举法得出相同的解。x1,x2=meshgrid(-200:200);f=100*(x2-x1.2).2+(4.5543-x1).2;fmin,i=sort(f(:);x=x1(i(1),x2(i(1)x=5 25第55页/共60页3 0-13 0-1规划问题求解MATLAB 7.0 版本提供的 0-1 线性规划问题当然也可以用前面的函数求解第56页/共60页例：f=-3,2,-5;A=1 2-1;1 4 1;1 1 0;0 4 1;B=2;4;3;6;x=bintprog(f,A,B,)第57页/共60页 x1,x2,x3=meshgr
42、id(0,1);i=find(x1+2*x2-x3=2)&(x1+4*x2+x3=4)&(x1+x2=3)&(4*x1+x3 f=-3*x1(i)+2*x2(i)-5*x3(i);fmin,ii=sort(f);index=i(ii(1);x=x1(index),x2(index),x3(index)x=1 0 1%还可以列出所有的可行解 x=x1(i(ii),x2(i(ii),x3(i(ii);x fmin ans=1 0 1 -8 0 0 1 -5 1 0 0 -3 0 0 0 0 0 1 0 2第58页/共60页例：f=-3,2,-5;xm=0;0;0;xM=1;1;1;intlist=1;1;1;A=1 2-1;1 4 1;1 1 0;0 4 1;B=2;4;5;6;ctype=-1*ones(4,1);res,b=ipslv_mex(f,A,B,intlist,xM,xm,ctype)res=1 0 1b=0第59页/共60页谢谢您的观看！第60页/共60页