【教学课件】第七章代数方程与最优化问题的求解

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,第七章,代数方程与最优化问题的求解,代数方程的求解,无约束最优化问题的计算机求解,有,约束最优化问题的计算机求解,整数规划问题的计算机求解,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 =,二元方程的图解法,例：, 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对根。图解法只能,显示求解方程的实根。,一般多项式方程的根可为实数，也可为复数。,可用MATLAB符号工具箱中的solve( )函数。,最简调用格式：,S=solve(eqn,1,eqn,2,eqn,n,),（,返回一个结构题型变量S，如S.x表示方程的根。）,直接得出根：,（变量返回到MATLAB工作空间）,x,=solve(eqn,1,eqn,2,eqn,n,),同上，并指定变量,x,=solve(eqn,1,eqn,2,eqn,n,，x,),例：, syms x y;, x,y=solve(x2+y2-1=0,75*x3/100-y+9/10=0),x =,y =,验证, 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),y(22),z(22),ans =,ans =,.37264668450644375527750811296627e-1,ans =,验证：, 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) % 将解代入求误差,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; 验证, norm(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+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 一般非线性方程数值解,非线性方程的标准形式为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的根。,例: 求 的根,解：, fun=x3-2*x-5;, z=fzero(fun,2) %初始估计值为2,z =,2.0946, format long, opt=optimset(Tolx,1.0e-8);, y=fzero(x3-2*x-5,2,opt),y =,2.09455148211709,7.1.4 一般非线性方程组数值解,格式：,最简求解语句,x=fsolve(Fun, x,0,),一般求解语句,x, f, flag, out=fsolve(Fun, x,0,opt, p,1, p,2,),若返回的flag 大于0，则表示求解成功，否则求解出现问题， opt 求解控制参数，结构体数据。,获得默认的常用变量,opt=optimset;,用这两种方法修改参数（解误差控制量）,opt.Tolx=1e-10; 或 set(opt.Tolx, 1e-10),例：,自编函数：,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,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, 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,初值改变有可能得出另外一组解。故初值的选择对解的影响很大，在某些初值下甚至无法搜索到方程的解。,例：,用solve( )函数求近似解析解, syms t; solve(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)* cos(2*t) - 0.5),ans =,不允许手工选择初值，只能获得这样的一个解。,可先用用图解法选取初值，再调用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-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,重新设置相关的控制变量，提高精度。, 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,例: 求 的根,解：化为标准形式,设初值点为x0=-5 -5。,function F = 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),Norm of First-order Trust-region,Iteration Func-count f(x) step optimality radius,0 3 47071.2 2.29e+004 1,1 6 12003.4 1 5.75e+003 1,2 9 3147.02 1 1.47e+003 1,3 12 854.452 1 388 1,4 15 239.527 1 107 1,5 18 67.0412 1 30.8 1,6 21 16.7042 1 9.05 1,7 24 2.42788 1 2.26 1,8 27 0.032658 0.759511 0.206 2.5,9 30 7.03149e-006 0.111927 0.00294 2.5,10 33 3.29525e-013 0.00169132 6.36e-007 2.5,Optimization terminated: first-order optimality is less than options.TolFun.,x =,fval =,1.0e-006 *,-0.40590960570519,-0.40590960570519,例: 求矩阵x使其满足方程，,并设初始解向量为x=1, 1; 1, 1。,解：function F = 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.1291 0.8602,1.2903 1.1612,Fval =,1.0e-003 *,0.1541 -0.1163,0.0109 -0.0243,exitflag =,1,7.2无约束最优化问题求解 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 =, y2=diff(y1); b=subs(y2,t,t0) % 并验证二阶导数为正,b =, 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),单变量函数求最小值的,标准形式为,: s.t.,函数,fminbnd,格式,x = fminbnd(fun,x1,x2),x = fminbnd(fun,x1,x2,options),x,fval = fminbnd() % fval为目标函数的最小值,x,fval,exitflag = fminbnd(),exitflag为终止迭代的条件,若,exitflag0,，收敛于,x,，,exitflag=0,，表示超过函数估计值或迭代的最大数字，,exitflag x,fval,exitflag,output=fminbnd(x3+cos(x)+x*log(x)/exp(x),0,1),x,=,0.5223,fval =,0.3974,exitflag =,1,output =,iterations: 9,funcCount: 11,algorithm: golden section search, parabolic interpolation,例:,在,0,5,上求下面函数的最小值,解：,function f = myfun(x),f = (x-3).2 - 1;, x=fminbnd(myfun,0,5),x =,3,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) 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, 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()函数，但,当所选函数高度不连续时，使用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);, 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(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,求梯度向量(比较引入梯度的作用，但梯度的计算也费时间）, 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*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,非线性最小二乘,函数,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= 。,例: 求下面非线性最小二乘问题,初始解向量为x0=0.3, 0.4。,解：先建立函数文件，并保存为myfun.m.,由于lsqnonlin中的fun为向量形式而不是平方和形式，因此，myfun函数应由 建立：,k=1,2,10,function F = myfun10(x),k = 1:10;,F = 2 + 2*k-exp(k*x(1)-exp(k*x(2);,然后调用优化程序：, x0 = 0.3 0.4;, x,resnorm = lsqnonlin(myfun10,x0),Optimization terminated: norm of the current step is less,than OPTIONS.TolX.,x =,0.2578 0.2578,resnorm =,124.3622,7.3有约束最优化问题的计算机求解 6.3.1 约束条件与可行解区域,有约束最优化问题的一般描述：,对于一般的一元问题和二元问题，可用图解法直接得出问题的最优解。,例：用图解方法求解：, 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); 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,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.7850,0.0000,3.3200,11.3850,2.5700,f_opt =,-89.5750,key =,1 求解成功,c =,iterations: 5,algorithm: medium-scale: activeset,cgiterations: ,例：求解, 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: 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,7.3.3 二次型规划的求解,例：求解, 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),Optimization terminated successfully.,x =,0.0000 0.6667 1.6667 2.6667,f_opt =,-23.6667 （解的目标值应为30（ -23.6667 ）6.3333),例:求解下面二次规划问题,解：,则,H = 1 -1; -1 2 ; f = -2; -6;,A = 1 1; -1 2; 2 1; b = 2; 2; 3;,lb = zeros(2,1);, x,fval,exitflag,output,lambda = quadprog(H,f,A,b, , ,lb),x = %最优解,0.6667,1.3333,fval = %最优值,-8.2222,exitflag = %收敛,1,output =,iterations: 3,algorithm: medium-scale: active-set,firstorderopt: ,cgiterations: ,例: 求二次规划的最优解,max f (x1, x2)=x1x2+3,sub.to x1+x2-2=0,解：化成标准形式：,H=0,-1;-1,0; f=0;0; Aeq=1 1;,x,fval,exitflag,output = quadprog(H,f, , ,Aeq),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)-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.5121,0.2170,3.5522,f_opt =,961.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 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 % 从略（计算结果同上一样）,例：利用梯度求解,梯度函数：, 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=; 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,约束线性最小二乘,有约束线性最小二乘的标准形式为:,其中：C、A、Aeq为矩阵；d、b、beq、lb、ub、x是向量。,函数,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乘子,例: 求解下面系统的最小二乘解,系统：Cx=d,约束:,先输入系统系数和x的上下界：,C = 0.9501 0.7620 0.6153 0.4057;,0.2311 0.4564 0.7919 0.9354;,0.6068 0.0185 0.9218 0.9169;,0.4859 0.8214 0.7382 0.4102;,0.8912 0.4447 0.1762 0.8936;,d = 0.0578; 0.3528; 0.8131; 0.0098; 0.1388;,A = 0.2027 0.2721 0.7467 0.4659;,0.1987 0.1988 0.4450 0.4186;,0.6037 0.0152 0.9318 0.8462;,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.1000,0.2152,0.3502,resnorm =,0.1672,residual =,0.0455,0.0764,-0.3562,0.1620,0.0784,exitflag =,1 %说明解x是收敛的,output =,iterations: 4,algorithm: medium-scale: active-set,firstorderopt: ,cgiterations: ,lambda =,lower: 4x1 double,upper: 4x1 double,eqlin: 0x1 double,ineqlin: 3x1 double,%通过lambda.ineqlin可查看非线性不等式约束是否有效。, lambda.ineqlin,ans =,0,0.2392,0,7.4整数规划问题的计算机求解,7.4.1 整数线性规划问题的求解,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,x4,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 -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,一般非线性整数规划问题与求解,err字符串为空，则返回最优解。,该函数尚有不完全之处，解往往不是精确整数，,可用下面语句将其化成整数。,例：,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,1,3;,errmsg,f,X=bnb20(c6miopt,x0,intlist,xm,xM,A,B,Aeq,Beq); X,X =,19.0000,0,3.8000,11.0000,3.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.00000000000000 25.00000001055334, if length(