1、第七章代数方程与最优化问题的求解,代数方程的求解 无约束最优化问题的计算机求解 有约束最优化问题的计算机求解 整数规划问题的计算机求解,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) % 同时绘制横轴,验证： syms t ; t=3.52028; vpa(exp(-3*t)*sin(4*t+2)+ 4*exp(-0.5*t)*cos(2*t)-0.5)ans = -.19256654148425
2、145223200161126442e-4,二元方程的图解法 例： ezplot(x2*exp(-x*y2/2)+exp(-x/2)*sin(x*y) % 第一个方程曲线 hold on % 保护当前坐标系 ezplot(x2 * cos(x+y2) + y2*exp(x+y),方程的图解法 仅适用于一元、 二元方程的求根 问题。,7.1.2 多项式型方程的准解析解法,例： ezplot(x2+y2-1); hold on % 绘制第一方程的曲线 ezplot(0.75*x3-y+0.9) % 绘制第二方程为关于x的6次多项式方程 应有6对根。图解法只能 显示求解方程的实根。,一般多项式方程的
3、根可为实数，也可为复数。 可用MATLAB符号工具箱中的solve( )函数。最简调用格式：S=solve(eqn1,eqn2,eqnn) （返回一个结构题型变量S，如S.x表示方程的根。）直接得出根： （变量返回到MATLAB工作空间）x,=solve(eqn1,eqn2,eqnn)同上，并指定变量x,=solve(eqn1,eqn2,eqnn，x,),例： syms x y; x,y=solve(x2+y2-1=0,75*x3/100-y+9/10=0) x = -.98170264842676789676449828873194 -.553951760568345600779844138
4、82735-.35471976465080793456863789934944*i -.55395176056834560077984413882735+.35471976465080793456863789934944*i .35696997189122287798839037801365 .86631809883611811016789809418650-1.2153712664671427801318378544391*i .86631809883611811016789809418650+1.2153712664671427801318378544391*i y = .19042035
5、099187730240977756415289 .92933830226674362852985276677202-.21143822185895923615623381762210*i .92933830226674362852985276677202+.21143822185895923615623381762210*i .93411585960628007548796029415446 -1.4916064075658223174787216959259-.70588200721402267753918827138837*i -1.491606407565822317478721695
6、9259+.70588200721402267753918827138837*i,验证 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的符号工具箱得出原始问题的高精度数值解，故称之为准解析解。,例： 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
7、),y(22),z(22) ans =-1.0869654762986136074917644096117ans =.37264668450644375527750811296627e-1 ans =.89073290972562790151300874796949,验证： 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+
8、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) % 将解代入求误差 ans =5.4882e-028,例：求解 （含变量倒数） 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; 验证 n
9、orm(double(eval(err) ans =8.9625e-030,例：求解 （带参数方程） syms a b x y; x,y=solve(x2+a*x2+6*b+3*y2=0,y=a+(x+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+
10、1/2/(4+a)*(-6*a-18-2*(-21*a2-45*a-27-24*b-6*a*b-3*a3)(1/2)+3,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),例：自编函数： function
11、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.35700.9341 Y =1.0e-009 *0.12150.0964,c =1 d = iterations: 7funcCount: 21algorit
12、hm: trust-region doglegfirstorderopt: 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,0T, x,Y,c,d=fsolve(f,-1,0,OPT
13、); x, Y, kk=d.funcCount Optimization terminated successfully:First-order optimality is less than options.TolFun. x =-0.98170.1904 Y =1.0e-010 *0.5619-0.4500 kk =15初值改变有可能得出另外一组解。故初值的选择对解的影响很大，在某些初值下甚至无法搜索到方程的解。,例：用solve( )函数求近似解析解 syms t; solve(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)* cos(2*t) - 0.5) ans
14、 = .67374570500134756702960220427474 不允许手工选择初值，只能获得这样的一个解。可先用用图解法选取初值，再调用fsolve( )函数数值计算, 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-regionIteration Func-count f(x) step optimality radius1
15、2 1.8634e-009 5.16e-005 12 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,重新设置相关的控制变量，提高精度。 ff=optimset; ff.TolX=1e-16; ff.TolFun=1e-30; ff.Display=iter; t,f=fsolve(y,3.52
16、03,ff)Norm of First-order Trust-regionIteration Func-count f(x) step optimality radius1 2 1.8634e-009 5.16e-005 12 4 3.67694e-019 3.61071e-005 7.25e-010 13 6 0 5.07218e-010 0 1 Optimization terminated successfully:First-order optimality is less than options.TolFun. t =3.52026389244155 f =0,7.2无约束最优化
17、问题求解 7.2.1 解析解法和图解法,例： 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) % 绘制出选定区间内 一阶导函数曲线, t0=solve(y1) % 求出一阶导数等于零的点 t0 =1.4528424981725411893375778048840 y2=diff(y1); b=subs(y2,t,t0) % 并验证二阶导数为正b = 7.8553420253333601379464405534591 t=0:0.4:4;y=exp(
18、-3*t).*sin(4*t+2)+4*exp(-0.5*t).*cos(2*t)-0.5; plot(t,y),7.2.2 基于 MATLAB 的数值解法,例： 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) Procedure1 3 -0.000499937 initial2 4 -0.000499937 reflect72 137 -0.641424 c
19、ontract outside Optimization terminated successfully:x =0.6111-0.3056, 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()函数效
20、率高于fminsearch()函数，故若安装了最优化工具箱则应调用fminunc()函数。,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, syms t; y=exp(-2*t)*cos(10*t)+exp(-3*(t+2)*sin(2*t);
21、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=-0.2; t3,f3=fminsearch(f,t0); t3 f3 ans =-0.3340 -1.9163,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) %测试算法的函数, f=inline(10
22、0*(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.40
23、397e-009 -1.77e-013 x =1.000001736959721.00000347608069,求梯度向量(比较引入梯度的作用，但梯度的计算也费时间） syms x1 x2; f=100*(x2-x12)2+(1-x1)2; J=jacobian(f,x1,x2) J = -400*(x2-x12)*x1-2+2*x1, 200*x2-200*x12 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*
24、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.999999999999990.99999999999998,7.3有约束最优化问题的计算机求解 6.3.1 约束条件与可行解区域,有约束最优化问题的一般描述：对于一般的一元问题和二元问题，可用图解法直接得出问题的最优解。,例：用图解方法求解： x1,x2=meshgrid(-
25、3:.1:3); % 生成网格型矩阵 z=-x1.2-x2; % 计算出矩阵上各点的高度 i=find(x1.2+x2.29); z(i)=NaN; % 找出 x12+x229 的坐标，并置函数值为 NaN i=find(x1+x21); z(i)=NaN; % 找出 x1+x21的坐标,置为 NaN surf(x1,x2,z); shading interp; max(z(:), view(0,90) ans =3,7.3.2 线性规划问题的计算机求解,例：求解 f=-2 1 4 3 1; A=0 2 1 4 2; 3 4 5 -1 -1; B=54; 62; Ae=; Be=; xm=0,
26、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),Optimization terminated successfully. x =19.78500.00003.320011.38502.5700 f_opt =-89.5750 key =1 求解成功 c = iterations: 5algorithm:
27、 medium-scale: activesetcgiterations: ,例：求解 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),Residuals: Prima
28、l Dual Upper Duality TotalInfeas Infeas Bounds Gap RelA*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+001Iter 1: 6.38e-012 1.21e+001 0.00e+000 3.50e+003 1.78e+000Iter 10: 0.00e+000 6.15e-026 0.00e+000 1.70e-024 4.10e-028 Optimization terminated successfully
29、. x = -5.0000 -0.1947 1.0000 -5.0000 f_opt =-55.4700 key =1 c = iterations: 10cgiterations: 0algorithm: lipsol,7.3.3 二次型规划的求解,例：求解 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
30、,OPT) Optimization terminated successfully. x =0.0000 0.6667 1.6667 2.6667 f_opt =-23.6667 （解的目标值应为30（ -23.6667 ）6.3333),6.3.4 一般非线性规划问题的求解,例：求解目标函数： 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)-
31、25; 8*x(1)+14*x(2)+7*x(3)-56;c = ;, 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.51210.21703.5522 f_opt =96
32、1.7152 kk = %目标函数调用的次数108,简化非线约束函数 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 Procedur
33、e 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 andmaximum constraint violation is less than options.TolCon Active Constraints:12 x, f_opt, kk=d.funcCount % 从略（计
34、算结果同上一样）,例：利用梯度求解梯度函数： 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);, x0=1;1;1; xm=0;0;0; xM=
35、; 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.51210.21703.5522 f_opt =961.7152 kk =95,7.4整数规划问题的计算机求解7.4.1 整数线性
36、规划问题的求解,A、B线性等式和不等式约束，约束式子由ctype变量控制，intlist 为整数约束标示，how0表示结果最优，2为无可行解，其余失败。,例： 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, x1,x2,x3,x
37、4,x5=ndgrid(1:20,0:20,4:20,1:20,3:20); i=find(2*x2+x3+4*x4+2*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(index),x2(index),x3(index),x4(index),x5(index) x =19 0 4 10 5 当把20换为30，一般计算机配置下实现不了，故求解整数规划时不适合采用穷举算法。,次最优解 fmin(1:10) ans =-89 -88 -88 -88 -88 -88 -88 -88 -87
38、 -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 -8918 0 4 11 3 -8817 0 5 10 4 -8815 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,7.4.2一般非线性整数规划问题与求解,err字符串为空，则返回最优解。 该函数尚有不完全之处，解往往不是精确整数，可用
39、下面语句将其化成整数。,例：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,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.0000 0 4.0000 10.0000 5.0000, intlist=1,0,0,1,1; xm=0,0,3.32
40、,1,3; errmsg,f,X=bnb20(c6miopt,x0,intlist,xm,xM,A,B,Aeq,Beq); X X =19.000003.800011.00003.0000,例：编写函数： function y=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.0000000000
41、0000 25.00000001055334, if length(errmsg)=0, x=round(x), end; x=x; x =525 缩小范围。 xm=-20*1;1; xM=20*1;1; errmsg,f,x=bnb20(c6fun1,x0,intlist,xm,xM,A,B,Aeq,Beq); x ans =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,7.4.3 0-1规划
42、问题求解,MATLAB 7.0 版本提供的 0-1 线性规划问题,当然也可以用前面的函数求解,例： f=-3,2,-5; A=1 2 -1; 1 4 1; 1 1 0; 0 4 1; B=2;4;5;6; x=bintprog(f,A,B,), x1,x2,x3=meshgrid(0,1); i=find(x1+2*x2-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 -80 0 1 -51 0 0 -30 0 0 00 1 0 2,例： 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 =101 b =0,
