代数方程与最优化问题

1、代数方程与最优化问题第1页，共86页，2022年，5月20日，13点31分，星期日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)% 同时绘制横轴第2页，共86页，2022年，5月20日，13点31分，星期日验证： syms t ; t=3.52028; vpa(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5) ans =第3页，共86页，2022年，5月20日，

2、13点31分，星期日二元方程的图解法例： ezplot(x2*exp(-x*y2/2)+exp(-x/2)*sin(x*y) % 第一个方程曲线 hold on % 保持当前坐标系 ezplot(x2 *cos(x+y2) +y2*exp(x+y)第4页，共86页，2022年，5月20日，13点31分，星期日方程的图解法仅适用于一元、二元方程的求根问题。第5页，共86页，2022年，5月20日，13点31分，星期日7.1.2 多项式型方程的准解析解法例： ezplot(x2+y2-1); hold on % 绘制第一方程的曲线 ezplot(0.75*x3-y+0.9) % 绘制第二方程为关于

3、x的6次多项式方程应有6对根。图解法只能显示求解方程的实根。第6页，共86页，2022年，5月20日，13点31分，星期日一般多项式方程的根可为实数，也可为复数。可用MATLAB符号工具箱中的solve( )函数。最简调用格式： S=solve(eqn1,eqn2,eqnn) （返回一个结构题型变量S，如S.x表示方程的根。）直接得出根： （变量返回到MATLAB工作空间） x,=solve(eqn1,eqn2,eqnn)同上，并指定变量 x,=solve(eqn1,eqn2,eqnn，x,)第7页，共86页，2022年，5月20日，13点31分，星期日例： syms x y; x,y=sol

4、ve(x2+y2-1=0,75*x3/100-y+9/10=0)x =y =第8页，共86页，2022年，5月20日，13点31分，星期日验证 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页，2022年，5月20日，13点31分，星期日例： x,y,z=solve(x+3*y3+2*z2=1/

5、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-1ans =第10页，共86页，2022年，5月20日，13点31分，星期日验证： 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*z

6、2=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第11页，共86页，2022年，5月20日，13点31分，星期日例：求解（含变量倒数） 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.

7、/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页，2022年，5月20日，13点31分，星期日例：求解（带参数方程） 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*

8、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第13页，共86页，2022年，5月20日，13点31分，星期日7.1.3 一般非线性方程数值解非线性方程的标准形式为f(x)=0 函数 fzero 格式 x = fzero (fun,x0) %用fun定义表达式f(x)，x0为初始解。 x = fzero (fun,x0,options) x,fval = fze

9、ro() %fval=f(x) x,fval,exitflag = fzero() x,fval,exitflag,output = fzero() 说明 该函数采用数值解求方程f(x)=0的根。第14页，共86页，2022年，5月20日，13点31分，星期日例: 求 的根解： fun=x3-2*x-5; z=fzero(fun,2) %初始估计值为2z = 2.0946 format long opt=optimset(Tolx,1.0e-8); y=fzero(x3-2*x-5,2,opt)y = 2.09455148211709第15页，共86页，2022年，5月20日，13点31分，星

10、期日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页，2022年，5月20日，13点31分，星期日例：自编函数：function q = my2deq(p) q=p(1)*p(

11、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.9341Y = 1.0e-009 * 0.1215 0.0964第17页，共86页，2022年，5月20日，13点31分，星期日c = 1d = iterations: 7 funcCou

12、nt: 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,0T第18页，共86页，202

13、2年，5月20日，13点31分，星期日 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 0.1904Y = 1.0e-010 * 0.5619 -0.4500kk = 15 初值改变有可能得出另外一组解。故初值的选择对解的影响很大，在某些初值下甚至无法搜索到方程的解。第19页，共86页，2022年，5月20日，13点31分，星期日例：用solve(

14、 )函数求近似解析解 syms t; solve(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)* cos(2*t) - 0.5)ans =不允许手工选择初值，只能获得这样的一个解。可先用用图解法选取初值，再调用fsolve( )函数数值计算第20页，共86页，2022年，5月20日，13点31分，星期日 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 o

15、f 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 optimality is less than options.TolFun.t = 3.52026389294877f = -6.063776702980306e-010第21页，共86页，2022年，5月2

16、0日，13点31分，星期日重新设置相关的控制变量，提高精度。 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 1Optimizat

17、ion terminated successfully: First-order optimality is less than options.TolFun.t = 3.52026389244155f = 0第22页，共86页，2022年，5月20日，13点31分，星期日例: 求 的根解：化为标准形式设初值点为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); % 显示

18、输出信息 x,fval = fsolve(myfun,x0,options) 第23页，共86页，2022年，5月20日，13点31分，星期日 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 1

19、6.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.5Optimization terminated: first-order optimality is less than options.TolFun.x =fval = 1.0e-006 * -0.40590960570519 -0.40590960570519第24页，共86页，2022年

20、，5月20日，13点31分，星期日例: 求矩阵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.1612Fval = 1.0e-003 * 0.1541 -0.1163 0.0109 -0.0243exitflag = 1第25页，共86页，

21、2022年，5月20日，13点31分，星期日7.2无约束最优化问题求解 7.2.1 解析解法和图解法第26页，共86页，2022年，5月20日，13点31分，星期日例： 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) % 绘制出选定区间内一阶导函数曲线第27页，共86页，2022年，5月20日，13点31分，星期日 t0=solve(y1) % 求出一阶导数等于零的点t0 = y2=diff(y1); b=subs(y2,t,t0) % 并验证二阶

22、导数为正 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)第28页，共86页，2022年，5月20日，13点31分，星期日单变量函数求最小值的(数值解法)%求解区间优化问题标准形式为 : s.t.函数 fminbnd格式 x = fminbnd(fun,x1,x2) x = fminbnd(fun,x1,x2,options) x,fval = fminbnd() % fval为目标函数的最小值 x,fval,exitflag = fminbnd() exitflag为终止迭代的条件,若ex

23、itflag0，收敛于x，exitflag=0，表示超过函数估计值或迭代的最大数字，exitflag x,fval,exitflag,output=fminbnd(x3+cos(x)+x*log(x)/exp(x),0,1)x = 0.5223fval = 0.3974exitflag = 1output = iterations: 9 funcCount: 11 algorithm: golden section search, parabolic interpolation 第30页，共86页，2022年，5月20日，13点31分，星期日例:在0,5上求下面函数的最小值 解： functi

24、on f = myfun(x) f = (x-3).2 - 1; x=fminbnd(myfun,0,5)x = 3第31页，共86页，2022年，5月20日，13点31分，星期日7.2.2 基于MATLAB的数值解法第32页，共86页，2022年，5月20日，13点31分，星期日例： 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

25、3 -0.000499937 initial 2 4 -0.000499937 reflect 72 137 -0.641424 contract outsideOptimization terminated successfully: x = 0.6111 -0.3056第33页，共86页，2022年，5月20日，13点31分，星期日 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

26、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()函数。第34页，共86页，2022年，5月20日，13点31分，星期日7.2.3 全局最优解与局部最优解例： f=inline(exp(-2*t).*cos(10*t)+exp(-3*(t+2).*sin(2*t),t); % 目标函数 t0=1;

27、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第35页，共86页，2022年，5月20日，13点31分，星期日 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=-

28、0.2; t3,f3=fminsearch(f,t0); t3 f3ans = -0.3340 -1.9163第36页，共86页，2022年，5月20日，13点31分，星期日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页，2022年，5月20日，13点31分，星期日 f=inline(100*(x(2)-x(1)2)2+(1-x(1)2,x); ff=optimset; ff.Tol

29、X=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.00000

30、347608069第38页，共86页，2022年，5月20日，13点31分，星期日求梯度向量(比较引入梯度的作用，但梯度的计算也费时间） 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*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.

31、GradObj=on; x=fminunc(c6fun3,0;0,ff) Norm of First-order Iteration f(x) step optimality CG-iterations19 6.38977e-029 2.12977e-012 1.6e-014 1x = 0.99999999999999 0.99999999999998第39页，共86页，2022年，5月20日，13点31分，星期日7.2.5非线性最小二乘 函数 lsqnonlin 格式 x = lsqnonlin(fun,x0) x = lsqnonlin(fun,x0,lb,ub) x = lsqnonli

32、n(fun,x0,lb,ub,options) %x0为初始解向量；fun为 ，i=1,2,m， lb、ub定义x的下界和上界, options为指定优化参数，若x没有界，则lb= ，ub= 。第40页，共86页，2022年，5月20日，13点31分，星期日例: 求下面非线性最小二乘问题初始解向量为x0=0.3, 0.4。解：先建立函数文件，并保存为myfun.m. 由于lsqnonlin中的fun为向量形式而不是平方和形式，因此，myfun函数应由 建立： k=1,2,10 第41页，共86页，2022年，5月20日，13点31分，星期日function F = myfun10(x) k =

33、 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.2578resnorm = 124.3622第42页，共86页，2022年，5月20日，13点31分，星期日7.3有约束最优化问题的计算机求解 7.3.1 约束条件与可行解区域有约束最优化问题的一般描述： 对于一般的一元问题和

34、二元问题，可用图解法直接得出问题的最优解。第43页，共86页，2022年，5月20日，13点31分，星期日例：用图解方法求解： 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第44页，共86页，

35、2022年，5月20日，13点31分，星期日7.3.2 线性规划问题的计算机求解第45页，共86页，2022年，5月20日，13点31分，星期日例：求解 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,

36、Ae,Be,xm,ff)第46页，共86页，2022年，5月20日，13点31分，星期日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: 第47页，共86页，2022年，5月20日，13点31分，星期日例：求解 f=-3/4,150,-1/50,6; Aeq=; Beq=; A=1/4,-60,-1/50

37、,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页，2022年，5月20日，13点31分，星期日Residuals: Primal Dual Upper Duality Total Infeas Infeas Bounds Gap Rel A

38、*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 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

39、第49页，共86页，2022年，5月20日，13点31分，星期日7.3.3 二次型规划的求解（线性）第50页，共86页，2022年，5月20日，13点31分，星期日例：求解 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

40、= 0.0000 0.6667 1.6667 2.6667f_opt = -23.6667 （解的目标值应为30（ -23.6667 ）6.3333)第51页，共86页，2022年，5月20日，13点31分，星期日例:求解下面二次规划问题解：则 第52页，共86页，2022年，5月20日，13点31分，星期日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.6

41、667 1.3333fval = %最优值 -8.2222exitflag = %收敛 1output = iterations: 3 algorithm: medium-scale: active-set firstorderopt: cgiterations: 第53页，共86页，2022年，5月20日，13点31分，星期日例: 求二次规划的最优解max f (x1, x2)=x1x2+3sub.to x1+x2-2=0解：化成标准形式：H=0,-1;-1,0; f=0;0; Aeq=1 1; x,fval,exitflag,output = quadprog(H,f, , ,Aeq) 第

42、54页，共86页，2022年，5月20日，13点31分，星期日6.3.4 一般非线性规划问题的求解第55页，共86页，2022年，5月20日，13点31分，星期日例：求解目标函数：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 = ;第56页，共86页，2022年，5月20日，1

43、3点31分，星期日 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.5121 0.2170 3.5522f_opt = 961.7152kk = %目标函数调用的次数 10

44、8第57页，共86页，2022年，5月20日，13点31分，星期日简化非线约束函数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 P

45、rocedure 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: Search direction less than 2*options.TolX and maximum constraint violation is less than options.TolConActive Constraints: 1 2 x, f_opt, kk=d.funcCou

46、nt % 从略（计算结果同上一样）第58页，共86页，2022年，5月20日，13点31分，星期日例：利用梯度求解梯度函数： 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

47、(1); 2*x(3)+x(1);第59页，共86页，2022年，5月20日，13点31分，星期日 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.funcCo

48、untx = 3.5121 0.2170 3.5522f_opt = 961.7152kk = 95 第60页，共86页，2022年，5月20日，13点31分，星期日约束线性最小二乘 有约束线性最小二乘的标准形式为:其中：C、A、Aeq为矩阵；d、b、beq、lb、ub、x是向量。第61页，共86页，2022年，5月20日，13点31分，星期日函数 lsqlin 格式 x = lsqlin(C,d,A,b) x = lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0,options) %求在约束条件下，方程Cx = d的最小二乘解x。 Aeq、beq满足等式约束，若没有不等式约束，

49、则设 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页，2022年，5月20日，13点31分，星期日例: 求解下面系统的最小二乘解系统：Cx=d约束:先输入系

50、统系数和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;第63页，共86页

51、，2022年，5月20日，13点31分，星期日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.3502resnorm = 0.1672第64页，共86页，2022年，5月20日，13点31分，星期日residual = 0.0455 0.0764 -0.3562 0.1620 0.0784exitf

52、lag = 1 %说明解x是收敛的output = iterations: 4 algorithm: medium-scale: active-set firstorderopt: cgiterations: 第65页，共86页，2022年，5月20日，13点31分，星期日lambda = lower: 4x1 double upper: 4x1 double eqlin: 0 x1 double ineqlin: 3x1 double%通过lambda.ineqlin可查看不等式约束是否有效。 lambda.ineqlinans = 0 0.2392 0第66页，共86页，2022年，5月2

53、0日，13点31分，星期日7.4整数规划问题的计算机求解7.4.1 整数线性规划问题的求解(matlab2009a 有问题）A、B线性等式和不等式约束,，约束式子由ctype变量控制，intlist为整数约束标示，how0表示结果最优，2为无可行解，其余失败。第67页，共86页，2022年，5月20日，13点31分，星期日例： 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=ipsl

54、v_mex(f,A,B,intlist,xM,xm,ctype) res = 19 0 4 10 5b = % 因为返回的 b=0，表示优化成功 0第68页，共86页，2022年，5月20日，13点31分，星期日 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(index),x2(i

55、ndex),x3(index),x4(index),x5(index)x = 19 0 4 10 5当把20换为30，一般计算机配置下实现不了，故求解整数规划时不适合采用穷举算法。第69页，共86页，2022年，5月20日，13点31分，星期日次最优解 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

56、 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.0000b = 0第70页，共86页，2022年，5月20日，13点31分，星期日7.4.2一般非线性整数规划问题与求解(matlab2009a 有问题）第71页，共86页，2022年，5月20日，13点31分，星期日err字符串为空，则返回最优解。该函数尚有不完全之处，解往往不是精确整数，可用下面语句将其化成整数。第72页，共86页，2022年，5月20日，13点31分，星

57、期日例：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=XX = 19.0000 0 4.0000 10.0000 5.0000第73页，共86页，2022年，5月20日，13点31分，星期日 intlist=1,

58、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第74页，共86页，2022年，5月20日，13点31分，星期日例：编写函数：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

59、,x0,intlist,xm,xM,A,B, Aeq,Beq); xans = 5.00000000000000 25.00000001055334第75页，共86页，2022年，5月20日，13点31分，星期日 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,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第76页，共86页，2022年