matlab在科学计算中的应用7

null第七章 代数方程与最优化问题的求解第七章 代数方程与最优化问题的求解代数方程的求解 无约束最优化问题的计算机求解 有约束最优化问题的计算机求解 整数规划问题的计算机求解 7.1代数方程的求解 7.1.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]) % 同时绘制横轴null验证： >> 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-4null二元方程的图解法 例： >> 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)')null 方程的图解法 仅适用于一元、 二元方程的求根 问题。 7.1.2 多项式型方程的准解析解法7.1.2 多项式型方程的准解析解法例： >> ezplot('x^2+y^2-1'); hold on % 绘制第一方程的曲线 >> ezplot(‘0.75*x^3-y+0.9’) % 绘制第二方程 为关于x的6次多项式方程 应有6对根。图解法只能 显示求解方程的实根。null一般多项式方程的根可为实数，也可为复数。 可用MATLAB符号工具箱中的solve( )函数。 最简调用格式： S=solve(eqn1,eqn2,…,eqnn) （返回一个结构题型变量S，如S.x表示方程的根。） 直接得出根： （变量返回到MATLAB工作空间） [x,…]=solve(eqn1,eqn2,…,eqnn) 同上，并指定变量 [x,…]=solve(eqn1,eqn2,…,eqnn，’x,…’)null例： >> 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]null验证 >> [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的符号工具箱得出原始问题的高精度数值解，故称之为准解析解。null 例： >> [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 null验证： >> 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-028null例：求解 （含变量倒数） >> 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-030null例：求解 （带参数方程） >> 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]7.1.3 一般非线性方程数值解7.1.3 一般非线性方程数值解非线性方程的 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 形式为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的根。null例: 求 的根 解： >> fun='x^3-2*x-5'; >> z=fzero(fun,2) %初始估计值为2 z = 2.0946 >> format long >> opt=optimset('Tolx',1.0e-8); >> y=fzero('x^3-2*x-5',2,opt) y = 2.094551482117097.1.4 一般非线性方程组数值解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)null例： 自编函数： 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.0964nullc = 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]Tnull>> [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 初值改变有可能得出另外一组解。故初值的选择对解的影响很大，在某些初值下甚至无法搜索到方程的解。null例： 用solve( )函数求近似解析解 >> syms t; solve(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)* cos(2*t) - 0.5) ans = .67374570500134756702960220427474 ％不允许手工选择初值，只能获得这样的一个解。 可先用用图解法选取初值，再调用fsolve( )函数数值计算null>> 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-010null重新设置相关的控制变量，提高精度。 >> 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 = 0null例: 求 的根 解：化为标准形式 设初值点为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) null 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 = 0.56714303139736 0.56714303139736 fval = 1.0e-006 * -0.40590960570519 -0.40590960570519null例: 求矩阵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 = 17.2无约束最优化问题求解 7.2.1 解析解法和图解法7.2无约束最优化问题求解 7.2.1 解析解法和图解法null例： >> 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]) % 绘制出选定区间内 一阶导函数曲线null>> 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)单变量函数求最小值的(数值解法)单变量函数求最小值的(数值解法)标准形式为 : s.t. 函数 fminbnd 格式 x = fminbnd(fun,x1,x2) x = fminbnd(fun,x1,x2,options) [x,fval] = fminbnd(…) % fval为目标函数的最小值 [x,fval,exitflag] = fminbnd(…) exitflag为终止迭代的条件,若exitflag>0，收敛于x，exitflag=0，表示超过函数估计值或迭代的最大数字，exitflag<0表示函数不收敛于x； [x,fval,exitflag,output] = fminbnd(…) output为优化信息,若output=iterations表示迭代次数，output=funccount表示函数赋值次数，output=algorithm表示所使用的算法 null例: 计算下面函数在区间(0，1)内的最小值。 解： >> [x,fval,exitflag,output]=fminbnd('(x^3+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' null例:在[0,5]上求下面函数的最小值 解： function f = myfun(x) f = (x-3).^2 - 1; >> x=fminbnd(@myfun,0,5) x = 37.2.2 基于MATLAB的数值解法7.2.2 基于MATLAB的数值解法null例： >> 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.3056null >> 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 全局最优解与局部最优解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.5436null>> 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 7.2.4 利用梯度求解最优化问题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]) %测试算法的函数null>> 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.00000347608069null％求梯度向量(比较引入梯度的作用，但梯度的计算也费时间） >> 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.999999999999987.2.5非线性最小二乘 7.2.5非线性最小二乘 函数 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=[ ]。 null例: 求下面非线性最小二乘问题 初始解向量为x0=[0.3, 0.4]。 解：先建立函数文件，并保存为myfun.m. 由于lsqnonlin中的fun为向量形式而不是平方和形式，因此，myfun函数应由 建立： k=1,2,…,10 nullfunction 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.36227.3有约束最优化问题的计算机求解 7.3.1 约束条件与可行解区域7.3有约束最优化问题的计算机求解 7.3.1 约束条件与可行解区域有约束最优化问题的一般描述： 对于一般的一元问题和二元问题，可用图解法直接得出问题的最优解。null例：用图解方法求解： >> [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 7.3.2 线性规划问题的计算机求解7.3.2 线性规划问题的计算机求解null例：求解 >> 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)nullOptimization 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: [] null例：求解 >> 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; 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)nullResiduals: 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 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 二次型规划的求解7.3.3 二次型规划的求解null例：求解 >> 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)null例:求解下面二次规划问题 解： 则 null>>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: [ ]null例: 求二次规划的最优解 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 一般非线性规划问题的求解6.3.4 一般非线性规划问题的求解null例：求解 目标函数： 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 = [];null>> 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 = %目标函数调用的次数 108null简化非线约束函数 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 % 从略（计算结果同上一样）null例：利用梯度求解 梯度函数： >> 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)];null>> 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 6.3.5约束线性最小二乘 6.3.5约束线性最小二乘 有约束线性最小二乘的标准形式为: 其中：C、A、Aeq为矩阵；d、b、beq、lb、ub、x是向量。null函数 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乘子null例: 求解下面系统的最小二乘解 系统：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];nullb =[ 0.5251; 0.2026; 0.6721]; lb