优化方法MATLAB编程——大连理工大学
完整版优化设计Matlab编程作业
化设计hl4HU©0⑥ 3 hlu 凹内r d X1州fci-rU-fFF卢F ♦ 忡下¥为+1 —*— S-ll-« F41：Si —MATLABoftiHMirjirCfiffliiiiJ PHI■1**■ 温不平?」11,・—喜M - 〜FT 文词一时y 片 34ml 3F*L9TR0i. Jill!-LkftLgWf 1S1CSI掰f 1 ■ >A A A »W I % ：k Dnfl w I ■ J k^lXMprfaMk tjn nn Alflhw初选 x0=[1,1] 程序：Step 1: Write an Mfle objfunl.m.function f1=objfun1(x)f1=x(1)人2+2*x(2)入2-2*x(1)*x(2)-4*x(1);Step 2: Invoke one of the unconstrained optimization routinesx0=[1,1]；>> options = 0Ptimset('LargeScale','off);>> [x,fval,exitflag,output] = fminunc(@objfun1,x0,options)运行结果： x =4.0000 2.0000 fval = -8.0000exitflag =1 output = iterations: 3 funcCount: 12 stepsize: 1 firstorderopt: 2.3842e-007algorithm: 'medium-scale: Quasi-Newton line search message: [1x85 char]非线性有约束优化1. Min f(x)=3 x ： + x 2+2 x 1-3 x 2+5 Subject to:g 2(x)=5 X 1-3 X 2 -25 < 0 g (x)=13 X -41 X 2 < 0 3 12g 4(x)=14 < X 1 < 130无约束优化 min f(x)=X 2 + x 2-2 x 1 x 2-4 x 1g5 (x)=2 < X 2 < 57初选x0=[10,10]Step 1: Write an M-file objfun2.mfunction f2=objfun2(x)f2=3*x(1)人2+x(2)人2+2*x(1)-3*x(2)+5;Step 2: Write an M-file confunl.m for the constraints. function [c,ceq]=confun1(x) % Nonlinear inequality constraints c=[x(1)+x(2)+18;5*x(1)-3*x(2)-25;13*x(1)-41*x(2)人2;14-x(1);x(1)-130;2-x(2);x(2)-57];% Nonlinear inequality constraints ceq=[];Step 3: Invoke constrained optimization routinex0=[10,10]; % Make a starting guess at the solution>> options = optimset('LargeScale','off);>> [x, fval]=...fmincon(@objfun2,x0,[],[],[],[],[],[],@confun1,options)运行结果：x =3.6755 -7.0744 fval =124.14952.min f (x) =4x2 + 5x2s.t. g 1(x) = 2X] + 3x2- 6 < 0g (x) = x x +1 > 0初选x0=[1,1]Step 1: Write an M-file objfun3.m function f=objfun3(x) f=4*x(1)人2 + 5*x(2)人2Step 2: Write an M-file confun3.m for the constraints. function [c,ceq]=confun3(x) %Nonlinear inequality constraints c=[2*x(1)+3*x(2)-6;-x(1)*x(2)-1];% Nonlinear equality constraints ceq口；Step 3: Invoke constrained optimization routinex0=[1,1];% Make a starting guess at the solution>> options = optimset('LargeScale','off);>> [x, fval]=...fmincon(@objfun,x0,[],[],[],[],[],[],@confun,options)运行结果：Optimization terminated: no feasible solution found. Magnitude of search direction less than2*options.TolX but constraints are not satisfied.x =11fval =-13实例：螺栓连接的优化设计图示为一压气机气缸与缸盖连接的示意图。
提高matlab运算速度的几种方法
由于matlab是一种解释性语言，所以在matlab程序中最忌讳直接使用循环语句，如果不得已要使用for循环，可以采用以下方法提高速度。
1、使用6.5以上版本，对循环已作优化；2、尽可能转化为矩阵运算；3、转化为二进制执行文件运算，如使用matlab内带的编译系统或matcom以及com组件技术。
其中com组件技术最方便的就是利用com builder来实现，这里重点介绍。
com builder是matlab6.5才有的，也是mathworks公司推荐使用于混合编程的，这些日子进行了全方位的摸索，感觉是爽呆了，下面我们一起来揭开它的神秘面纱。
此系列分为以下几块：1.matlab下做com组件2.vb,c#.net实现调用3.vc实现调用4.打包5.优缺点评注其中2,3部分可以选择一个看matlab下做com组件com是component object module的简称，它是一种通用的对象接口，任何语言只要按照这种接口标准，就可以实现调用它。
matlab6.5新推出来的combuilder就是把matlab下的程序做成com组件，供其他语言调用。
我们先准备两个测试文件，并copy一个图片到c盘下，起名叫1.jpg(这些你都可以改的，我这儿是为了程序方便)：第一个叫im_test.m如下:function im_test %这个文件不带输入与输出I=imread('c:\1.jpg'); %因为以前带有imshow的程序用mcc编成dll后用不%了，所以想试combuilder是否imshow(I); %能支持这些函数第二个叫split2rgb.m，就是前些日子Zosco发给我的那个程序，因为它用mcc编成dll后有问题，所以我在这儿继续将它进行测试，而且它也带有多个输入及输出参数，所以也正好拿来测试在matlab的workspace下打comtool，就打开了matlab com builder，点击file-new project，新建一个工程，在component name里填上comtest，Class name里填上一个sgltest(并将自动生成classes里的comtest remove掉),compliecode in选c或c++都无所谓，将Complier options里的Use Handle Graphics library的复选框画上，点击ok就行了。
Matlab中的非线性优化和非线性方程求解技巧
Matlab中的非线性优化和非线性方程求解技巧在科学和工程领域中，我们经常会遇到一些复杂的非线性问题，例如最优化问题和方程求解问题。
解决这些问题的方法主要分为线性和非线性等，其中非线性问题是相对复杂的。
作为一种强大的数值计算工具，Matlab提供了许多专门用于解决非线性优化和非线性方程求解的函数和方法。
本文将介绍一些常用的Matlab中的非线性优化和非线性方程求解技巧。
非线性优化是指在给定一些约束条件下，寻找目标函数的最优解的问题。
在实际应用中，往往需要根据实际情况给出一些约束条件，如等式约束和不等式约束。
Matlab中的fmincon函数可以用于求解具有约束条件的非线性优化问题。
其基本语法如下：[x,fval] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub)其中，fun是目标函数，x0是初始值，A、b是不等式约束矩阵和向量，Aeq、beq是等式约束矩阵和向量，lb、ub是变量的上下边界。
x表示最优解，而fval表示最优解对应的目标函数值。
另外，非线性方程求解是指寻找使得方程等式成立的变量值的问题。
Matlab中提供的fsolve函数可以用于求解非线性方程。
其基本语法如下：x = fsolve(fun,x0)其中，fun是方程函数，x0是初始值，x表示方程的解。
除了fmincon和fsolve函数之外，Matlab还提供了一些其他的非线性优化和非线性方程求解函数，例如lsqnonlin、fminunc等，这些函数分别适用于无约束非线性优化问题和带约束非线性方程求解问题。
除了直接调用这些函数外，Matlab还提供了一些可视化工具和辅助函数来帮助我们更好地理解和解决非线性问题。
例如，使用Matlab的优化工具箱可以实现对非线性优化问题的求解过程可视化，从而更直观地观察到优化算法的收敛过程。
此外，Matlab还提供了一些用于计算梯度、雅可比矩阵和海塞矩阵的函数，这些函数在求解非线性问题时非常有用。
大连理工优化方法-增广拉格朗日方法MATLAB程序
大连理工优化方法-增广拉格朗日方法MATLAB程序上机大作业II定义目标函数funfunction f=fun(x)x1=x(1);x2=x(2);f=4*x1-x2^2-12;定义目标函数梯度函数dfunfunction f=dfun(x)x2=x(2);f=[4;-2*x2];定义等式约束函数hffunction qua=hf(x)qua=25-x(1)^2-x(2)^2;定义等式约束函数梯度函数dhffunction qua=dhf(x)qua=[-2*x(1);-2*x(2)];定义不等式约束函数gfunfunction inq=gfun(x)inq=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;定义不等式约束梯度数dgffunction inq=dgf(x)inq=[10-2*x(1);10-2*x(2)];定义增广拉格朗日函数mpsifunctionpsi=mpsi(x,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma)f=feval(fun,x);he=feval(hf,x);gi=feval(gfun,x);l=length(he);m=length(gi);psi=f;s1=0;for i=1:lpsi=psi-he(i)*mu(i);s1=s1+he(i)^2;endpsi=psi+0.5*sigma*s1;s2=0.0;for i=1:ms3=max(0.0, lambda(i) - sigma*gi(i));s2=s2+s3^2-lambda(i)^2;endpsi=psi+s2/(2.0*sigma);定义增广拉格朗日函数梯度函数dmpsifunctiondpsi=dmpsi(x,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma) dpsi=feval(dfun,x);he=feval(hf,x);gi=feval(gfun,x);dhe=feval(dhf,x);dgi=feval(dgf,x);l=length(he);m=length(gi);for i=1:ldpsi=dpsi+(sigma*he(i)-mu(i))*dhe(:,i);endfor i=1:mdpsi=dpsi+(sigma*gi(i)-lambda(i))*dgi(:,i);end定义BFGS法函数函数bfgsfunction[x,val,k]=bfgs(mpsi,dmpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda ,sigma) maxk=1000;rho=0.5;sigma1=0.4;epsilon1=1e-4;k=0;n=length(x0);Bk=eye(n);while(k<maxk)< p="">gk=feval(dmpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigm a);if(norm(gk)<epsilon1)< p="">break;enddk=-Bk\gk;m=0;mk=0;while(m<20)newf=feval(mpsi,x0+rho^m*dk,fun,hf,gfun,dfun,dhf,dgf,mu,l ambda,sigma);oldf=feval(mpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma);if(newf<oldf+sigma1*rho^m*gk'*dk)< p="">mk=m;break;endm=m+1;endx=x0+rho^mk*dk;sk=x-x0;yk=feval(dmpsi,x,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma) -gk;if(yk'*sk>0)Bk=Bk-((Bk*sk)*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);endk=k+1;x0=x;endval=feval(mpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma );定义增广拉格朗日乘子法函数multphrfunction answer=multphr(fun,hf,gfun,dfun,dhf,dgf,x0)maxk=5000;sigma=2.0;eta=2.0;theta=0.8;k=0;ink=0;epsilon=1e-4;x=x0;he=feval(hf,x);gi=feval(gfun,x);l=length(he);m=length(gi);mu=0.1*ones(l,1);lambda=0.1*ones(m,1);btak=10;btaold=10;while(btak>epsilon&&k<maxk)< p="">[x,v,ik]=bfgs('mpsi','dmpsi',x0,fun,hf,gfun,dfun,dhf,dgf,mu,la mbda,sigma); ink=ink+ik;he=feval(hf,x);gi=feval(gfun,x);btak=0.0;for i=1:lbtak=btak+he(i)^2;endfor i=1:mtemp=min(gi(i),lambda(i)/sigma);btak=btak+temp^2;endbtak=sqrt(btak);if btak>epsilonif(k>=2&&btak > theta*btaold)sigma=eta*sigma;endfor i=1:lmu(i)=mu(i)-sigma*he(i);endfor i=1:mlambda(i)=max(0.0,lambda(i)-sigma*gi(i)); endendk=k+1;btaold=btak;x0=x;endf=feval(fun,x);xfmulambdak运行求解>> x0=[0;0]x0 =>> multphr('fun','hf','gfun','dfun','dhf','dgf',x0) x = 1.001281489564374.89871784708758f =-31.9923105871169mu =1.01559644571312lambda =0.754451167977228k =4</maxk)<></oldf+sigma1*rho^m*gk'*dk)<></epsilon1)<></maxk)<>。
大连理工大学MATLAB概述
命令窗口
历史 命令
命令提示符
工作路径
MATLAB变量命名规则
变量名、函数名对字母的大小写是敏感的。 变量名第一个字母必须是英文字母。 变量名可以包含英文字母、下划线和数字。 变量名不能包含空格、标点。
MATLAB预定义变量
预定义变量名 ans eps pi
Inf或inf i或j
数据类型 int8, int16, int32, int64 uint8, uint16, uint32, uint64 single double logical char cell struct funcion_handle
说明 有符号整数 无符号整数 单精度浮点型 双精度浮点型 逻辑型 字符型 单元数组型 结构体型 函数句柄型
MATLAB数值表示
缺省的数据类型为双精度浮点型
例如：3 -10 0.001 1.3e10 1.256e-6 基本操作
ceil( ), floor(), round() %取整
single( )
%单精度浮点型
double( )
%双精度浮点型
MATLAB四则运算符
运算 加 减 乘 除 幂
例：计算sin(45ْ ) >>sin(45*pi/180)
Matalb中正弦函数sin就是常见的正弦函数。 它的参数值是以“弧度”为单位的。 Matlab对字母大小写是敏感的。
例：计算 2ex0.5 1 的值，其中x=4.92。
>>sqrt(2*exp(4.92+0.5)+1)
Matalb中开平方—sqrt(x)，是英文square root的缩写 。 Matalb中指数函数exp(x)，常见的表达方式。
大连理工大学优化方法上机作业
大连理工大学优化方法上机作业本页仅作为文档页封面，使用时可以删除This document is for reference only-rar21year.March优化方法上机大作业学院：电子信息与电气工程学部姓名：学号：指导老师：上机大作业（一）%目标函数function f=fun(x)f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;end%目标函数梯度function gf=gfun(x)gf=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)]; End%目标函数Hess矩阵function He=Hess(x)He=[1200*x(1)^2-400*x(2)+2,-400*x(1);-400*x(1), 200;];end%线搜索步长function mk=armijo(xk,dk)beta=0.5; sigma=0.2;m=0; maxm=20;while (m<=maxm)if(fun(xk+beta^m*dk)<=fun(xk)+sigma*beta^m*gfun(xk)'*dk) mk=m; break;endm=m+1;endalpha=beta^mknewxk=xk+alpha*dkfk=fun(xk)newfk=fun(newxk)%最速下降法function [k,x,val]=grad(fun,gfun,x0,epsilon)%功能：梯度法求解无约束优化问题：minf(x)%输入：fun，gfun分别是目标函数及其梯度，x0是初始点，% epsilon为容许误差%输出：k是迭代次数，x，val分别是近似最优点和最优值maxk=5000; %最大迭代次数beta=0.5; sigma=0.4;k=0;while(k<maxk)gk=feval(gfun,x0); %计算梯度dk=-gk; %计算搜索方向if(norm(gk)<epsilon), break;end%检验终止准则m=0;mk=0;while(m<20) %用Armijo搜索步长if(feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk) mk=m;break;endm=m+1;endx0=x0+beta^mk*dk;k=k+1;endx=x0;val=feval(fun,x0);>> x0=[0;0];>> [k,x,val]=grad('fun','gfun',x0,1e-4)迭代次数：k =1033x =0.99990.9998val =1.2390e-008%牛顿法x0=[0;0];ep=1e-4;maxk=10;k=0;while(k<maxk)gk=gfun(x0);if(norm(gk)<ep)x=x0miny=fun(x)k0=kbreak;elseH=inv(Hess(x0));x0=x0-H*gk;k=k+1;endendx =1.00001.0000miny =4.9304e-030迭代次数k0 =2%BFGS方法function [k,x,val]=bfgs(fun,gfun,x0,varargin) %功能：梯度法求解无约束优化问题：minf(x)%输入：fun，gfun分别是目标函数及其梯度，x0是初始点，% epsilon为容许误差%输出：k是迭代次数，x，val分别是近似最优点和最优值N=1000;epsilon=1e-4;beta=0.55;sigma=0.4;n=length(x0);Bk=eye(n);k=0;while(k<N)gk=feval(gfun,x0,varargin{:});if(norm(gk)<epsilon), break;enddk=-Bk\gk;m=0;mk=0;while(m<20)newf=feval(fun,x0+beta^m*dk,varargin{:});oldf=feval(fun,x0,varargin{:});if(newf<=oldf+sigma*beta^m*gk'*dk)mk=m;break;endm=m+1;endx=x0+beta^mk*dk;sk=x-x0;yk=feval(gfun,x,varargin{:})-gk;if(yk'*sk>0)Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);endk=k+1;x0=x;endval=feval(fun,x0,varargin{:});>> x0=[0;0];>> [k,x,val]=bfgs('fun','gfun',x0)k =20x =1.00001.0000val =2.2005e-011%共轭梯度法function [k,x,val]=frcg(fun,gfun,x0,epsilon,N)if nargin<5,N=1000;endif nargin<4, epsilon=1e-4;endbeta=0.6;sigma=0.4;n=length(x0);k=0;while(k<N)gk=feval(gfun,x0);itern=k-(n+1)*floor(k/(n+1));itern=itern+1;if(itern==1)dk=-gk;elsebetak=(gk'*gk)/(g0'*g0);dk=-gk+betak*d0; gd=gk'*dk;if(gd>=0),dk=-gk;endendif(norm(gk)<epsilon),break;endm=0;mk=0;while(m<20)if(feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk) mk=m;break;endm=m+1;endx=x0+beta^m*dk;g0=gk; d0=dk;x0=x;k=k+1;endval=feval(fun,x);>> x0=[0;0];[k,x,val]=frcg('fun','gfun',x0,1e-4,1000)k =122x =1.00011.0002val =7.2372e-009上机大作业（二）%目标函数function f_x=fun(x)f_x=4*x(1)-x(2)^2-12;%等式约束条件function he=hf(x)he=25-x(1)^2-x(2)^2;end%不等式约束条件function gi_x=gi(x,i)switch icase 1gi_x=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;case 2gi_x=x(1);case 3gi_x=x(2);otherwiseend%求目标函数的梯度function L_grad=grad(x,lambda,cigma)d_f=[4;2*x(2)];d_g(:,1)=[-2*x(1);-2*x(2)];d_g(:,2)=[10-2*x(1);10-2*x(2)];d_g(:,3)=[1;0];d_g(:,4)=[0;1];L_grad=d_f+(lambda(1)+cigma*hf(x))*d_g(:,1);for i=1:3if lambda(i+1)+cigma*gi(x,i)<0L_grad=L_grad+(lambda(i+1)+cigma*gi(x,i))*d_g(:,i+1);continueendend%增广拉格朗日函数function LA=lag(x,lambda,cee)LA=fun(x)+lambda(1)*hf(x)+0.5*cee*hf(x)^2;for i=1:3LA=LA+1/(2*cee)*(min(0,lambda(i+1)+cee*gi(x,i))^2-lambda(i+1)^2); endfunction xk=BFGS(x0,eps,lambda,cigma)gk=grad(x0,lambda,cigma);res_B=norm(gk);k_B=0;a_=1e-4;rho=0.5;c=1e-4;length_x=length(x0);I=eye(length_x);Hk=I;while res_B>eps&&k_B<=10000dk=-Hk*gk;m=0;while m<=5000if lag(x0+a_*rho^m*dk,lambda,cigma)-lag(x0,lambda,cigma)<=c*a_*rho^m*gk'*dkmk=m;break;endm=m+1;endak=a_*rho^mk;xk=x0+ak*dk;delta=xk-x0;y=grad(xk,lambda,cigma)-gk;Hk=(I-(delta*y')/(delta'*y))*Hk*(I-(y*delta')/(delta'*y))+(delta*delta')/(delta'*y);k_B=k_B+1;x0=xk;gk=y+gk;res_B=norm(gk);end%增广拉格朗日法function val_min=ALM(x0,eps)lambda=zeros(4,1);cigma=5;alpha=10;k=1;res=[abs(hf(x0)),0,0,0];for i=1:3res(1,i+1)=norm(min(gi(x0,i),-lambda(i+1)/cigma)); endres=max(res);while res>eps&&k<1000xk=BFGS(x0,eps,lambda,cigma);lambda(1)=lambda(1)+cigma*hf(xk);for i=1:3lambda(i+1)=lambda(i+1)+min(0,lambda(i+1)+gi(x0,1)); endk=k+1;cigma=alpha*cigma;x0=xk;res=[norm(hf(x0)),0,0,0];for i=1:3res(1,i+1)=norm(min(gi(x0,i),-lambda(i+1)/cigma)); endres=max(res);endval_min=fun(xk);fprintf('k=%d\n',k);fprintf('fmin=%.4f\n',val_min);fprintf('x=[%.4f;%.4f]\n',xk(1),xk(2));>> x0=[0;0];>> val_min=ALM(x0,1e-4)k=10fmin=-31.4003x=[1.0984;4.8779]val_min =-31.4003上机大作业（三）A=[1 1;-1 0;0 -1];n=2;b=[1;0;0];G=[0.5 0;0 2];c=[2 4];cvx_solver sdpt3cvx_beginvariable x(n)minimize (x'*G*x-c*x)subject toA*x<=bcvx_enddisp(x)Status: SolvedOptimal value (cvx_optval): -2.40.40000.6000A=[2 1 1;1 2 3;2 2 1;-1 0 0;0 -1 0;0 0 -1]; n=3;b=[2;5;6;0;0;0];C=[-3 -1 -3];cvx_solver sdpt3cvx_beginvariable x(n)minimize (C*x)subject toA*x<=bcvx_enddisp(x)Status: SolvedOptimal value (cvx_optval): -5.40.20000.00001.600011。
大连理工优化方法大作业MATLAB编程
fun ctio n [x,dk,k]=fjqx(x,s) flag=0;a=0;b=0;k=0;d=1;while (flag==0)[p,q]=getpq(x,d,s);if (P<0)b=d;d=(d+a)/2;endif(p>=0) &&( q>=0)dk=d;x=x+d*s;flag=1;endk=k+1;if (p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2};endend%定义求函数值的函数 fun ，当输入为 x0= (x1 ， x2 )时，输出为 f function f=fun(x)f=(x(2)-x(1)A2)A2+(1-x(1)F2;function gf=gfun(x)gf=[-4*x(1)*(x (2) -x(1)A2)+2*(x(1)-1),2*(x(2)-x(1)A2)];function [p,q]=getpq(x,d,s)p=fun(x)-fun(x+d*s)+0.20*d*gfun(x)*s';q=gfun(x+d*s)*s'-0.60*gfun(x)*s';结果：x=[0,1];s=[-1,1];[x,dk,k]=fjqx(x,s)x =-0.0000 1.0000dk =1.1102e-016k =54取初始= (0.0. 0,0)r^'l用兵柜梯皮法求解下面无约東优化问题：min f (x) = x孑—2x^X2 十2x孑 + x孑H-爲—X2天3 十 2xj + 3|X2 —*3,其中步长g的选取可利用习題1戎精确一维披索.注:通过比习题验证共範梯度法求辉门无二次西数极小点至多需要“次迭代.fun ctio n f= fun( X )%所求问题目标函数f=X(1)A2-2*X(1)*X (2)+2*X(2)A2+X(3)A2+ X(4) A2-X( 2)*X(3)+2*X(1)+3*X(2)-X(3);end function g= gfun( X )%所求问题目标函数梯度g=[2*X(1)-2*X(2)+2,-2*X(1)+4*X(2)-X(3)+3,2*X (3) -X (2)-1,2*X(4)];end function [ x,val,k ] = frcg( fun,gfun,xO )%功能：用FR共轭梯度法求无约束问题最小值%输入：x0是初始点，fun和gfun分别是目标函数和梯度%输出：x、val分别是最优点和最优值，k是迭代次数maxk=5000; %最大迭代次数rho=0.5;sigma=0.4;k=0;eps=10e-6;n=length(x0);while (k<maxk)g=feval(gfun,x0); % 计算梯度 itern=k-(n+1)*floor(k/(n+1));itern=itern+1;%计算搜索方向if (itern==1)d=-g;elsebeta=(g*g')/(g0*g0');d=-g+beta*d0;gd=g'*d;if (gd>=0.0)d=-g;endendif (norm(g)<eps)break ;endm=0;mk=0;while (m<20)if(feval(fu n,xO+rhoAm*d)<feval(fu n,xO)+sigma*rhoAm*g'*d) mk=m; break ;endm=m+1;endx0=x0+rho A mk*d;val=feval(fun,x0);g0=g;d0=d;k=k+1;endx=x0;val=feval(fun,x0);end结果：>> x0=[0,0,0,0];>> [ x,val,k ] = frcg( 'fun','gfun',x0 ) x =-4.0000 -3.0000 -1.0000 0val =-8.0000k =或者function [x,f,k]=second(x)k=0;dk=dfun(x);g0=gfun(x);s=-g0;x=x+dk*s;g1=gfun(x);while (norm(g1)>=0.02)if (k==3)k=0;g0=gfun(x);s=-g0;x=x+dk*s;g1=gfun(x);else if (k<3)u=(( norm(g1))A2)/( norm(gO)A2); s=-g1+u*s;k=k+1;g0=g1;dk=dfun(x);x=x+dk*s;g1=gfun(x);endendf=fun(x);endfunction f=fun(x)f=x(1F2-2*x(1)*x (2)+2*x (2)A2+x(3)A2+x(4)A2-x (2) *x (3)+2*x(1)+3*x(2)-x(3); function gf=gfun(x)gf=[2*x(1)-2*x(2)+2,-2*x(1)+4*x(2)-x(3)+3,2*x(3)-x(2)-1,2*x(4)];function [p,q]=con(x,d)ss=-gfun(x);p=fun(x)-fun(x+d*ss)+0.2*d*gfun(x)*(ss)';q=gfun(x+d*ss)*(ss)'-0.6*gfun(x)*(ss)';function dk=dfun(x)flag=0;a=0;d=1;while (flag==0)[p,q]=con(x,d);if (p<0)b=d;d=(d+a)/2;endif (p>=0)&&(q>=0)dk=d;flag=1;endif (p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2};endEnd结果： x=[0,0,0,0];>> [x,f,k]=second(x)x =-4.0147 -3.0132-1.0090 0 f = -7.9999k = 1取初始点3 = （0」）二考虑下面无约東优化问题：min f（x）二冷 + 2x2 + exp（xf + 天孑），其中歩长Qk的选取可別用习题1或精确一维搜索•搜索方向为一HNW ♦取垃=b•取皿=R2f防）］"9耳丈啟为BFG5公式亠通过此习题体会上述三种算法的收敛速度.fun ctio n [f,x,k]=third_1(x) k=0;g=gfu n(x);while (norm(g)>=0.001) s=-g;dk=dfu n( x,s);x=x+dk*s;k=k+1;g=gfu n(x);f=fun( x);endfun ctio n f=fun(x)f=x(1)+2*x(2)A2+exp(x(1)A2+x(2)A2);fun ctio n gf=gfu n(x)gf=[1+2*x(1)*exp(x(1)A2+x(2)A2),4*x(2)+2*x(2)*(x(1)A2+x(2)A2)]; function[j_1,j_2]=con(x,d,s)j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)'; j_2=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)'; function dk=dfun(x,s) % 获取步长 flag=0;a=0;d=1;while (flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;endif (p>=0)&&(q>=0)dk=d;flag=1;endif (p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2}; end结果：x=[0,1];[f,x,k]=third_1(x)f =0.7729x = -0.4196 0.0001k =8(1 ) 程序：function [f,x,k]=third_2(x)k=0;H=inv(ggfun(x));g=gfun(x);while (norm(g)>=0.001)s=(-H*g')';dk=dfun(x,s);x=x+dk*s;k=k+1;g=gfun(x);f=fun(x);endfunction f=fun(x)f=x(1)+2*x(2)A2+exp(x(1F2+x(2)A2);function gf=gfun(x) gf=[1+2*x(1)*exp(x(1F2+x(2)A2),4*x(2)+2*x(2)*(x(1F2+x(2)A2)]; function ggf=ggfun(x)ggf=[(4*x(1)A2+2)*exp(x(1)A2+x (2) A2),4*x(1)*x (2) *exp(x(1)A2+x(2)A2);4*x(1)*x(2)*exp(x(1)A2+x(2)A2),4+(4*x(2)A2+2)*exp(x(1)A2+x(2)A2)];function [j_1,j_2]=con(x,d,s)j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)';j_2=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)'; function dk=dfun(x,s) % 步长获取flag=0;a=0;d=1;b=10000;while (flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;flag=1;endif (p>=0)&&(q<0)a=d;if 2*d>=(d+b)/2d=(d+b)/2;endendEnd结果：x=[0,1];[f,x,k]=third_2(x)f =0.7729x = -0.4193 0.0001k =8(2) 程序：function [f,x,k]=third_3(x) k=0;X=cell(2);g=cell(2);X{1}=x;H=eye(2);g{1}=gfun(X{1});s=(-H*g{1}')';dk=dfun(X{1},s);X{2}=X{1}+dk*s;g{2}=gfun(X{2});while (norm(g{2})>=0.001)dg=g{2}-g{1};v=dx/(dx*dg')-(H*dg')'/(dg*H*dg');h1=H*dg'*dg*H/(dg*H*dg');h2=dx'*dx/(dx*dx');h3=dg*H*dg'*v'*v;H=H-h1+h2+h3;k=k+1;X{1}=X{2};g{1}=gfun(X{1});s=(-H*g{1}')';dk=dfun(X{1},s);X{2}=X{1}+dk*s;g{2}=gfun(X{2});norm(g{2});f=fun(x);x=X{2};endfunction f=fun(x)f=x(1)+2*x(2)A2+exp(x(1F2+x(2)A2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)A2+x(2)A2),4*x(2)+2*x(2)*(x(1)A2+x(2)A2)];function ggf=ggfun(x)ggf=[(4*x(1)A2+2)*exp(x(1)A2+x(2)A2),4*x(1)*x(2)*exp(x(1)A2+x(2)A2);4*x(1)*x(2)* exp(x(1)A2+x(2)A2),4+(4*x(2)A2+2)*exp(x(1)A2+x(2)A2);function [p,q]=con(x,d,s)p=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)';q=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)';function dk=dfun(x,s)flag=0;a=0;d=1;b=10000;while (flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;if (p>=0)&&(q>=0)dk=d;flag=1;endif (p>=0)&&(q<0)a=d;if 2*d>=(d+b)/2d=(d+b)/2;else d=2*d;endendend结果：x=[0,1];[f,x,k]=third_3(x)f =0.7729x = -0.41950.0000 k=6*U 用有效集法求解下面勺勺二次规划问题:(XI 一 I)2 + (x 2 一 2.5)2 X1 - 2X2 + 2 > 0-Xi — 2>(2 + 6 > 0-Xi + 2X2 + 2 > 0xi,x 2 > 0function callqpactH=[2 0; 0 2];c=[-2 -5]';Ae=[ ]; be=[];Ai=[1 -2; -1 -2; -1 2;1 0;0 1];bi=[-2 -6 -2 0 0]';x0=[0 0]';[x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,xO)fun ctio n [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) epsilo n=1.0e-9; err=1.0e-6;k=0; x=x0; n=len gth(x); kmax=1.0e3;n e=le ngth(be); ni=le ngth(bi); lamk=zeros( ne+n i,1); in dex=ones(n i,1);for (i=1:ni)if(Ai(i,:)*x>bi(i)+epsil on), i ndex(i)=0; end while (k<=kmax)mmSi.Aee=[];if (ne>0), Aee=Ae; endfor (j=1:ni)if (index(j)>0), Aee=[Aee; Ai(j,:)]; end endgk=H*x+c;[m1,n1] = size(Aee);[dk,lamk]=qsubp(H,gk,Aee,zeros(m1,1)); if (norm(dk)<=err)y=0.0;if (length(lamk)>ne)[y,jk]=min(lamk(ne+1:length(lamk))); endif (y>=0)exitflag=0;elseexitflag=1;for (i=1:ni)if (index(i)&(ne+sum(index(1:i)))==jk) index(i)=0; break ;endendk=k+1;elseexitflag=1;alpha=1.0; tm=1.0;for (i=1:ni)if ((index(i)==0)&(Ai(i,:)*dk<0))tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk);if (tm1<tm)tm=tm1; ti=i;endendendalpha=min(alpha,tm);x=x+alpha*dk;if (tm<1), index(ti)=1; endendif (exitflag==0), break ; endk=k+1;endoutput.fval=0.5*x'*H*x+c'*x;output.iter=k;function [x,lambda]=qsubp(H,c,Ae,be) ginvH=pinv(H); [m,n]=size(Ae);if (m>0)rb=Ae*ginvH*c + be;lambda=pinv(Ae*ginvH*Ae')*rb; x=ginvH*(Ae'*lambda-c);elsex=-ginvH*c;lambda=0;end结果>>callqpactx =1.40001.7000lambda =0.8000exitflag =output =fval: -6.4500iter: 7function [x,mu,lambda,output]=multphr(fu n, hf,gf,dfu n, dhf,dgf,xO)%功能：用乘子法解一般约束问题：min f(x), s.t. h(x)=0, g(x).=0%输入：x0是初始点，fun, dfun分别是目标函数及其梯度；% hf, dhf分别是等式约束(向量)函数及其 Jacobi矩阵的转置；% gf, dgf分别是不等式约束(向量)函数及其 Jacobi矩阵的转置；%输出：x是近似最优点，mu, lambda分别是相应于等式约束和不等式约束的乘子向量% output是结构变量，输出近似极小值f,迭代次数，内迭代次数等maxk=500;c=2.0;eta=2.0;theta=0.8;k=0;i nk=0;epsilo n=0.00001;x=xO;he=feval(hf,x);gi=feval(gf,x);n=len gth(x);l=le ngth(he);m=le ngth(gi);mu=zeros(l,1);lambda=zeros(m,1);btak=10;btaold=10;while (btak>epsilon&&k<maxk)%调用BFGS算法程序求解无约束子问题[x,ival,ik]=bfgs( 'mpsi' ,'dmpsi' ,x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c);ink=ink+ik;he=feval(hf,x);gi=feval(gf,x);btak=0;for i=1:lbtak=btak+he(y2;end% 更新乘子向量for i=1:mtemp=min(gi(i),lambda(i)/c);btak=btak+temp A2;endbtak=sqrt(btak);if btak>epsilonif k>=2&&btak>theta*btaoldc=eta*c;endfor i=1:lmu(i)=mu(i)-c*he(i);endlambda(i)=max(0,lambda(i)-c*gi(i));endk=k+1;btaold=btak;x0=x;endendf=feval(fun,x);output.fval=f;output.iter=k;%增广拉格朗日函数function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c) f=feval(fun,x);he=feval(hf,x);gi=feval(gf,x);l=length(he);m=length(gi);psi=f;s1=0;for i=1:lpsi=psi-he(i)*mu(i);s仁 s1+he(y2;psi=psi+0.5*c*s1;s2=0;for i=1:ms3=max(0,lambda(i)-c*gi(i));s2=s2+s3A2-lambda(i)A2;endpsi=psi+s2/(2*c);% 不等式约束函数文件 g1.mfunction gi=g1(x)gi=10*x(1)-x(1)A2+10*x(2)-x(2)A2-34;% 目标函数的梯度文件df1.mfunction g=df1(x)g=[4, -2*x(2)]';% 等式约束(向量)函数的Jacobi 矩阵(转置)文件 dh1.m function dhe=dh1(x)dhe=[-2*x(1), -2*x(2)]'% 不等式约束(向量)函数的Jacobi 矩阵(转置)文件 dg1.m function dgi=dg1(x)dgi=[10-2*x(1), 10-2*x(2)]';function [x,val,k]=bfgs(fun,gfun,x0,varargin) maxk=500; rho=0.55;sigma=0.4;epsilon=0.00001;k=0;n=length(x0);Bk=eye(n);while (k<maxk)gk=feval(gfun,x0,varargin{:});if (norm(gk)<epsilon)break ;enddk=-Bk\gk;m=0;mk=0;while (m<20)n ewf=feval(fu n, x0+rho A m*dk,vararg in {:});oldf=feval(fun,x0,varargin{:});if(newf<oldf+sigma*rhoAm*gk'*dk) mk=m;break ;endm=m+1;endx=x0+rhoAmk*dk;sk=x-x0;yk=feval(gfun,x,varargin{:})-gk;if (yk'*sk>0)Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);endk=k+1;x0=x;endval=feval(fun,x0,varargin{:});结果x=[2 2]';[x,mu,lambda,output]=multphr( 'fun' ,'hf' ,'gf1' ,'df' ,'dh' ,'dg' ,x0) x =1.00134.8987mu =0.7701lambda =0.9434output =fval: -31.9923iter: 4利用序列二次规划方法求解习题5中的约束优化问题:min 4xi 一好一 12s.t. 25 - x? —x孑=Q10x一召 + 10旳-xj - 34 > 0 X1,X2 > 0tf=[3,1,1]；A=[2,1,1;1,-1,-1];b=[2;-1];lb=[0,0,0]； x=li nprog(f,A,b,zeros(3),[0,0,0]',lb)结果：Optimization terminated.0.00000.50000.5000。
使用Matlab进行神经网络优化问题求解的方法
使用Matlab进行神经网络优化问题求解的方法一、引言在当今信息时代，神经网络已经成为解决复杂问题的重要工具。
随着计算能力的提升，神经网络优化问题的求解变得越来越重要。
而Matlab作为一种强大的科学计算软件，能够提供丰富的工具和函数来解决神经网络优化问题。
本文将介绍如何使用Matlab来解决神经网络优化问题。
二、神经网络优化问题的建模在使用Matlab解决神经网络优化问题之前，首先需要对问题进行建模。
通常来说，神经网络优化问题可以分为两类：单目标优化问题和多目标优化问题。
单目标优化问题是指希望优化网络的某个特定输出，常见的问题有回归问题和分类问题。
而多目标优化问题则是希望在多个指标上获得最优解，常见的问题有多目标分类和多目标回归问题。
在建模过程中，需要确定网络的结构和参数。
神经网络的结构通常由输入层、隐藏层和输出层组成。
输入层接受原始数据，隐藏层进行特征提取，输出层给出最终的结果。
而参数则包括权重和偏置，这些参数需要进行调整以达到最优解。
三、使用Matlab解决单目标优化问题1. 数据准备在解决单目标优化问题之前，首先需要准备好数据集。
数据集应该包含输入值和对应的目标值。
2. 网络训练使用Matlab的神经网络工具箱，可以方便地进行网络训练。
首先，需要创建一个神经网络对象，并设置好网络的结构和参数。
然后，使用训练函数对网络进行训练，常见的训练函数有Levenberg-Marquardt算法和梯度下降算法。
通过训练函数，可以不断调整网络的权重和偏置，直到达到最优解。
3. 网络评估训练完网络后，需要对网络进行评估。
可以使用测试数据集来评估网络的性能，通常采用预测误差、准确率等指标来评估网络的表现。
四、使用Matlab解决多目标优化问题解决多目标优化问题与解决单目标优化问题的方法类似，只是目标变成了多个。
可以使用多种方法来解决多目标优化问题，如加权法、约束法和分级法等。
1. 加权法加权法是一种常用的解决多目标优化问题的方法。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
优化方法上机大作业学院：姓名：学号：指导老师：肖现涛第一题源程序如下：function zy_x = di1ti(x)%di1ti是用来求解优化作业第一题的函数。
x0=x; yimuxulong=0.000001;g0=g(x0);s0=-g0;A=2*ones(100,100);k=0;while k<100lanmed=-(g0)'*s0/(s0'*A*s0);x=x0+lanmed*s0;g=g(x);k=k+1;if norm(g)<yimuxulongzy_x=x;fprintf('After %d iterations,obtain the optimal solution.\n \n The optimal solution is \n %f.\n\nThe optimal "x" is "ans".',k,f(x) )break;endmiu=norm(g)^2/norm(g0)^2;s=-g+miu*s0;g0=g; s0=s;x0=x;endfunction f=f(x)f=(x'*ones(100,1))^2-x'*ones(100,1);function g=g(x)g=(2*x'*ones(100,1))*ones(100,1)-ones(100,1);代入x0,运行结果如下：>> x=zeros(100,1);>> di1ti(x)After 1 iterations,obtain the optimal solution.The optimal solution is-0.250000.The optimal "x" is "ans".ans =0.005*ones(100,1).第二题1.最速下降法。
源程序如下：function zy_x=di2titidu(x)%该函数用来解大作业第二题。
x0=x; yimuxulong=1e-5; k=0;g0=g(x0); s0=-g0;while k>=0if norm(g0)<yimuxulongbreak;elselanmed=10;c=0.1;i=0;while i>=0&i<100x=x0+lanmed*s0;if f(x)>(f(x0)+c*lanmed*g0'*s0) lanmed=lanmed/2;i=i+1;elsebreak;endendx=x0+lanmed*s0;x0=x;g0=g(x);s0=-g0;k=k+1;endendzy_x=x;zyj=f(x);fprintf('after %d iterations,obtain the optimal solution.\n\nThe optimal solution is %f.\n\n The optimal "x" is "ans".\n',k,zyj);function f=f(x)x1=[1 0 0 0]*x;x2=[0 1 0 0]*x;x3=[0 0 1 0]*x;x4=[0 0 0 1]*x;f=(x1-1)^2+(x3-1)^2+100*(x2-x1^2)^2+100*(x4-x3^2)^2;function g=g(x)x1=[1 0 0 0]*x;x2=[0 1 0 0]*x;x3=[0 0 1 0]*x;x4=[0 0 0 1]*x;g=[2*(x1-1)-400*x1*(x2-x1^2);200*(x2-x1^2);2*(x3-1)-400*x3*(x4-x3^2);200*(x 4-x3^2)];>> x=[-1.2 1 -1.2 1]';>> di2titidu(x)after 5945 iterations,obtain the optimal solution. The optimal solution is 0.000000.The optimal "x" is "ans".ans =1.00001.00001.00001.00002.牛顿法源程序如下：function zy_x=di2tinewton(x)%该函数用来解大作业第二题。
x0=x; yimuxulong=1e-5; k=0;g0=g(x0); h0=h(x0);s0=-inv(h0)*g0;while k>=0&k<1000if norm(g0)<yimuxulongbreak;elsex=x0+s0;x0=x;g0=g(x);h0=h(x);s0=-inv(h0)*g0;k=k+1;endendzy_x=x;zyj=f(x);fprintf('after %d iterations,obtain the optimal solution.\n\nThe optimal solution is %f.\n\n The optimal "x" is "ans".\n',k,zyj);function f=f(x)x1=[1 0 0 0]*x;x2=[0 1 0 0]*x;x3=[0 0 1 0]*x;x4=[0 0 0 1]*x;f=(x1-1)^2+(x3-1)^2+100*(x2-x1^2)^2+100*(x4-x3^2)^2;function g=g(x)x1=[1 0 0 0]*x;x2=[0 1 0 0]*x;x3=[0 0 1 0]*x;x4=[0 0 0 1]*x;g=[2*(x1-1)-400*x1*(x2-x1^2);200*(x2-x1^2);2*(x3-1)-400*x3*(x4-x3^2);200*(x 4-x3^2)];function h=h(x)x1=[1 0 0 0]*x;x2=[0 1 0 0]*x;x3=[0 0 1 0]*x;x4=[0 0 0 1]*x;h=[2+1200*x1^2-400*x2 -400*x1 0 0;-400*x1 200 0 0;0 0 2+1200*x3^2-400*x4 -400*x3;0 0 -400*x3 200];代入初始值，运行结果如下：>> x=[-1.2 1 -1.2 1]';>> di2tinewton(x)after 6 iterations,obtain the optimal solution.The optimal solution is 0.000000.The optimal "x" is "ans".ans =1.00001.00001.00001.0000可以看出，用Newton法经过6次迭代就能求出最优解。
3.BFGS法源程序如下：function zy_x=di2tiBFGS(x)%该函数用来解大作业第二题。
x0=x; yimuxulong=1e-5; k=0;g0=g(x0); H0=eye(4);s0=-H0*g0;while k>=0&k<100if norm(g0)<yimuxulongbreak;elselanmed=10;c=0.1;i=0;while i>=0&i<100x=x0+lanmed*s0;if f(x)>(f(x0)+c*lanmed*g0'*s0)lanmed=lanmed/2;i=i+1;elsebreak;endendx=x0+lanmed*s0;dete_x=x-x0;dete_g=g(x)-g0;miu=1+dete_g'*H0*dete_g/(dete_x'*dete_g);H=H0+(miu*dete_x*dete_x'-H0*dete_g*dete_x'-dete_x*dete_g'*H0)/(dete_x'*dete _g);s=-H*g(x);x0=x;s0=s;H0=H;g0=g(x);k=k+1;endendzy_x=x;zyj=f(x);fprintf('after %d iterations,obtain the optimal solution.\n\nThe optimal solution is %f.\n\n The optimal "x" is "ans".\n',k,zyj);function f=f(x)x1=[1 0 0 0]*x;x2=[0 1 0 0]*x;x3=[0 0 1 0]*x;x4=[0 0 0 1]*x;f=(x1-1)^2+(x3-1)^2+100*(x2-x1^2)^2+100*(x4-x3^2)^2;function g=g(x)x1=[1 0 0 0]*x;x2=[0 1 0 0]*x;x3=[0 0 1 0]*x;x4=[0 0 0 1]*x;g=[2*(x1-1)-400*x1*(x2-x1^2);200*(x2-x1^2);2*(x3-1)-400*x3*(x4-x3^2);200*(x 4-x3^2)];代入初始值，计算结果如下：>> x=[-1.2 1 -1.2 1]';>> di2tiBFGS(x)after 53 iterations,obtain the optimal solution.The optimal solution is 0.000000.The optimal "x" is "ans".ans =1.00001.00001.00001.0000第三题1.惩罚函数法源程序如下：function zy_x=di3ti(x)%该函数用来解大作业第三题。
x0=x; M=100; c=4; m=1;while m>0g0=g(x0,M); yimuxulong=1e-5;k=0;s0=-inv(H(x0,M))*g0; while k>=0if norm(g0)<yimuxulongbreak;elsex=x0+s0; %牛顿法；x0=x;g0=g(x,M);s0=-inv(H(x0,M))*g0;k=k+1;endendif max([abs(h(x)),g1(x),g2(x),g3(x)])<0.5break;elseM=M*c;m=m+1;endendzy_x=x;zyj=f(x);fprintf('after %d iterations,obtain the optimal solution.\n\nThe optimal solution is %f.\n\n The optimal "x" is "ans".\n',m,zyj);function F=F(x,M)x1=[1 0]*x;x2=[0 1]*x;F=4*x1-x2^2-12+M*(h^2+g1^2+g2^2+g3^2);function g=g(x,M)x1=[1 0]*x;x2=[0 1]*x;g=[4+M*(-4*(25-x1^2-x2^2)*x1+2*(10*x1-x1^2+10*x2-x2^2-34)*(10-2*x1)+2*x1);-2*x2+M*(-4*(25-x1^2-x2^2)*x2+2*(10*x1-x1^2+10*x2-x2^2-34)*(10-2*x2)+2*x2)];function H=H(x,M)x1=[1 0]*x;x2=[0 1]*x;H=[M*(24*x1^2-120*x1+8*x2^2-40*x2+238),M*(16*x1*x2-40*x1-40*x2+200);M*(16*x 1*x2-40*x1-40*x2+200),-2+M*(24*x2^2-120*x2+8*x1^2-40*x1+238)];function f=f(x)x1=[1 0]*x;x2=[0 1]*x;f=4*x1-x2^2-12;function h=h(x)x1=[1 0]*x;x2=[0 1]*x;h=25-x1^2-x2^2;function g1=g1(x)x1=[1 0]*x;x2=[0 1]*x;g=10*x1-x1^2+10*x2-x2^2-34; if g<0g1=g;elseg1=0;endfunction g2=g2(x)x1=[1 0]*x;x2=[0 1]*x;if x1>=0g2=0;elseg2=x1;endfunction g3=g3(x)x1=[1 0]*x;x2=[0 1]*x;if x2>=0g3=0;elseg3=x2;end代入任意初始值，运算结果如下。