1 / 61

第七章 代数方程与最优化问题的求解

第七章 代数方程与最优化问题的求解. 代数方程的求解 无约束最优化问题的计算机求解 有 约束最优化问题的计算机求解 整数规划问题的计算机求解. 7.1 代数方程的求解 7.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]) % 同时绘制横轴. 验证：

aysel
Download Presentation

第七章 代数方程与最优化问题的求解

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. 第七章代数方程与最优化问题的求解 • 代数方程的求解 • 无约束最优化问题的计算机求解 • 有约束最优化问题的计算机求解 • 整数规划问题的计算机求解

  2. 7.1代数方程的求解7.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]) % 同时绘制横轴

  3. 验证： >> 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

  4. 二元方程的图解法 例： >> ezplot('x^2*exp(-x*y^2/2)+exp(-x/2)*sin(x*y)') % 第一个方程曲线 >> hold on % 保护当前坐标系 >> ezplot(‘x^2 *… cos(x+y^2) +… y^2*exp(x+y)')

  5. 方程的图解法 仅适用于一元、 二元方程的求根 问题。

  6. 7.1.2 多项式型方程的准解析解法 • 例： >> ezplot('x^2+y^2-1'); hold on % 绘制第一方程的曲线 >> ezplot(‘0.75*x^3-y+0.9’) % 绘制第二方程 为关于x的6次多项式方程 应有6对根。图解法只能 显示求解方程的实根。

  7. 一般多项式方程的根可为实数，也可为复数。 可用MATLAB符号工具箱中的solve( )函数。 最简调用格式： S=solve(eqn1,eqn2,…,eqnn) （返回一个结构题型变量S，如S.x表示方程的根。） 直接得出根： （变量返回到MATLAB工作空间） [x,…]=solve(eqn1,eqn2,…,eqnn) 同上，并指定变量 [x,…]=solve(eqn1,eqn2,…,eqnn，’x,…’)

  8. 例： >> syms x y; >> [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]

  9. 验证 >> [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的符号工具箱得出原始问题的高精度数值解，故称之为准解析解。

  10. 例： >> [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 = 27 1 >> x(22),y(22),z(22) ans = -1.0869654762986136074917644096117 ans = .37264668450644375527750811296627e-1 ans = .89073290972562790151300874796949

  11. 验证： >> 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

  12. 例：求解 （含变量倒数） >> syms x y; >> [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 = 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

  13. 例：求解 （带参数方程） >> syms a b x y; >> [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]

  14. 7.1.3 一般非线性方程数值解 • 格式： 最简求解语句 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)

  15. 例： 自编函数： 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('my2deq',[1; 2],OPT) Optimization terminated successfully: First-order optimality is less than options.TolFun. x = 0.3570 0.9341 Y = 1.0e-009 * 0.1215 0.0964

  16. c = 1 d = 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)+0.9]','p'); >> [x,Y] = fsolve(f,[1; 2],OPT); % 结果和上述完全一致，从略。 Optimization terminated successfully: First-order optimality is less than options.TolFun. 若改变初始值 x0=[-1,0]T

  17. >> [x,Y,c,d]=fsolve(f,[-1,0]',OPT); x, Y, kk=d.funcCount Optimization terminated successfully: First-order optimality is less than options.TolFun. x = -0.9817 0.1904 Y = 1.0e-010 * 0.5619 -0.4500 kk = 15 初值改变有可能得出另外一组解。故初值的选择对解的影响很大，在某些初值下甚至无法搜索到方程的解。

  18. 例： 用solve( )函数求近似解析解 >> syms t; solve(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)* cos(2*t) - 0.5) ans = .67374570500134756702960220427474 ％不允许手工选择初值，只能获得这样的一个解。 可先用用图解法选取初值，再调用fsolve( )函数数值计算

  19. >> format long >> 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) 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 1 Optimization terminated successfully: First-order optimality is less than options.TolFun. t = 3.52026389294877 f = -6.063776702980306e-010

  20. 重新设置相关的控制变量，提高精度。 >> 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.67694e-019 3.61071e-005 7.25e-010 1 3 6 0 5.07218e-010 0 1 Optimization terminated successfully: First-order optimality is less than options.TolFun. t = 3.52026389244155 f = 0

  21. 7.2无约束最优化问题求解7.2.1 解析解法和图解法

  22. 例： >> 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]) % 绘制出选定区间内 一阶导函数曲线

  23. >> 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)

  24. 7.2.2 基于 MATLAB 的数值解法

  25. 例： >> 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) Iteration Func-count min f(x) Procedure 1 3 -0.000499937 initial 2 4 -0.000499937 reflect 72 137 -0.641424 contract outside Optimization terminated successfully: x = 0.6111 -0.3056

  26. >> x=fminunc(f,[0;.0],ff) Directional 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()函数。

  27. 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); [t1 f1] ans = 0.9228 -0.1547 >> t0=0.1; [t2,f2]=fminsearch(f,t0); [t2 f2] ans = 0.2945 -0.5436

  28. >> 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]) ％ 在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); [t3 f3] ans = -0.3340 -1.9163

  29. 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]) %测试算法的函数

  30. >> 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: 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

  31. ％求梯度向量(比较引入梯度的作用，但梯度的计算也费时间）％求梯度向量(比较引入梯度的作用，但梯度的计算也费时间） >> syms x1 x2; 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) Norm of First-order Iteration f(x) step optimality CG-iterations 19 6.38977e-029 2.12977e-012 1.6e-014 1 x = 0.99999999999999 0.99999999999998

  32. 7.3有约束最优化问题的计算机求解6.3.1 约束条件与可行解区域 • 有约束最优化问题的一般描述： 对于一般的一元问题和二元问题，可用图解法直接得出问题的最优解。

  33. 例：用图解方法求解： >> [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); shading interp; >> max(z(:)), view(0,90) ans = 3

  34. 7.3.2 线性规划问题的计算机求解

  35. 例：求解 >> 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='iter'; >> [x,f_opt,key,c]=linprog(f,A,B,Ae,Be,xm,[],[],ff)

  36. Optimization terminated successfully. x = 19.7850 0.0000 3.3200 11.3850 2.5700 f_opt = -89.5750 key = 1 ％求解成功 c = iterations: 5 algorithm: 'medium-scale: activeset' cgiterations: []

  37. 例：求解 >> 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; TolCon=1e-20; ff.Display='iter'; >> [x,f_opt,key,c]=linprog(f,A,B,Aeq,Beq,xm,xM, [0;0;0;0],ff)

  38. 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+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-028 Optimization terminated successfully. x = -5.0000 -0.1947 1.0000 -5.0000 f_opt = -55.4700 key = 1 c = iterations: 10 cgiterations: 0 algorithm: 'lipsol'

  39. 7.3.3 二次型规划的求解

  40. 例：求解 >> 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.6667 f_opt = -23.6667 ％（解的目标值应为30＋（ -23.6667 ）＝6.3333)

  41. 6.3.4 一般非线性规划问题的求解

  42. 例：求解 目标函数： 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 = [];

  43. >> 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.funcCount x = 3.5121 0.2170 3.5522 f_opt = 961.7152 kk = %目标函数调用的次数 108

  44. 简化非线约束函数 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); 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 infeasible 21 116 961.715 0 1 -6.3e-015 6.97e-005 Hessian modified twice Optimization terminated successfully: Search direction less than 2*options.TolX and maximum constraint violation is less than options.TolCon Active Constraints: 1 2 >> x, f_opt, kk=d.funcCount % 从略（计算结果同上一样）

  45. 例：利用梯度求解 梯度函数： >> 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 [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)];

  46. >> 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.funcCount x = 3.5121 0.2170 3.5522 f_opt = 961.7152 kk = 95

  47. 7.4整数规划问题的计算机求解7.4.1 整数线性规划问题的求解 A、B线性等式和不等式约束，约束式子由ctype变量控制，intlist 为整数约束标示，how＝0表示结果最优，2为无可行解，其余失败。

  48. 例： >> 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=[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 5 b = 0

  49. >> [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 = 19 0 4 10 5 ％当把20换为30，一般计算机配置下实现不了，故求解整数规划时不适合采用穷举算法。

  50. 次最优解 >> 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,0,1,1]; ％混合整数规划 >> [res,b]=ipslv_mex(f,A,B,intlist,xM,xm,ctype) res = 19.0000 0 3.8000 11.0000 3.0000 b = 0

More Related