优化方法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. 加权法加权法是一种常用的解决多目标优化问题的方法。

Matlab在最优化问题中的应用举例

在企业生产和日常生活中，人们总是希望用最少的人力、物力、财力和时间去办更多的事，这就是所谓的最优化问题。

线性规划方法是解决最优化问题的有效方法之一，因此受到人们的普遍关注。

在企业生产过程中，生产计划安排直接影响到企业的经济效益，而生产计划本质就是在目标一定时，对于人力、时间和物质资源的优化配置问题。

1。

综述了最优化方法，归纳了最优化闯题中线性规划和非线性规划模型的解法，并给出了相应的matlab求解代码。

2。

提出了基于信息增益率的用电客户指标选择方法，根据信息增益率的大小选择对分类有贡献的指标。

关键词：Matlab，最优化方法，应用举例In enterprise production and daily life, people always hope with the least amount of human, material and financial resources and time to do more things, this is the so-called optimization problem. Linear programming method is to solve the optimal problem, so one of the effective method by people's attention. In enterprise production process, production plan directly affect the enterprise economic benefit, but in essence is the production plan for the target certain human, time and material resources optimization allocation problem.1·Studying the optimization,summing up the solutions ofoptimization problem for both linear and non-linear programming model and proposing the matlabcode．2·Proposing a new way based on information-gain-ratio to choose the powercustomer indices,selecting the indices which are more contributive to theclassification,in order to avoid over learning。

第8章 MATLAB优化设计

首先将原线性规划问题转换为线性规划的MATLAB标准型，如下所示：
MIN : Y f X 4 x1 5 x2 x3 MIN : Y C T X 3x1 2 x2 x3 17 AX b 2 x1 x2 9 s.t. x x x 10 s . t . Aeq X Beq 3 1 2 x1 , x2 , x3 0 lb X ub
options=optimset('TolX',1e-7, 'TolFun',1e-7, 'TolCon',1e-7);
%优化设置
[X,Y,exitflag,output,lambda] = linprog(C,A,b,Aeq,beq,lb,ub,X0,options)%解算
第8单元 MATLAB优化设计
第8单元 MATLAB优化设计
确定目标函数(总利润)：
f X 1.25 0.25 0.05 5 x1 1.25 0.25 0.03 7 x2 0.06 6 x3 0.11 4 x4 0.05 7 x5 2 0.35 0.05 10 x6 2 0.35 0.03 9 x7 0.06 8 x8 2.8 0.5 0.03 12 x9 0.11 11 x10 0.75 x1 0.79 x2 0.36 x3 0.44 x4 0.35 x5 1.15 x6 1.38 x7 0.48 x8 1.94 x9 1.21x10
第8单元 MATLAB优化设计
(2) 将原线性规划问题转换为线性规划的MATLAB标准型：

最优化方法matlab作业

实用最优化方法——matlab编程作业初值为[-1;1]其中g0、g1分别为不同x值下得导数，f0、f1为函数值MATLAB程序：x0=[-1;1];s0=[1;1];c1=0.1;c2=0.5;a=0;b=inf;d=1;n=0;x1=x0+d*s0;g0=[-400*(x0(2)-x0(1)^2)*x0(1)-2*(1-x0(1));200*(x0(2)-x0(1) ^2)];g1=[-400*(x1(2)-x1(1)^2)*x1(1)-2*(1-x1(1));200*(x1(2)-x1(1) ^2)];f1=100*(x1(2)-x1(1)^2)^2+(1-x1(1))^2;f0=100*(x0(2)-x0(1)^2)^2+(1-x0(1))^2;while((f0-f1<-c1*d*g0'*s0)||(g1'*s0<c2*g0'*s0))if ((f0-f1)<(-c1*d*g0'*s0))b=d;d=(d+a)/2;x1=x0+d*s0;g0=[-400*(x0(2)-x0(1)^2)*x0(1)-2*(1-x0(1));200*(x0 (2)-x0(1)^2)];g1=[-400*(x1(2)-x1(1)^2)*x1(1)-2*(1-x1(1));200*(x1 (2)-x1(1)^2)];f1=100*(x1(2)-x1(1)^2)^2+(1-x1(1))^2;f0=100*(x0(2)-x0(1)^2)^2+(1-x0(1))^2;elseif (((g1')*s0)<(c2*(g0')*s0))a=d;if(2*d<=(d+b)/2)d=2*d;elsed=(d+b)/2;endx1=x0+d*s0;g0=[-400*(x0(2)-x0(1)^2)*x0(1)-2*(1-x0(1));200*(x0(2) -x0(1)^2)];g1=[-400*(x1(2)-x1(1)^2)*x1(1)-2*(1-x1(1));200*(x1(2 )-x1(1)^2)];f1=100*(x1(2)-x1(1)^2)^2+(1-x1(1))^2;f0=100*(x0(2)-x0(1)^2)^2+(1-x0(1))^2;endx1df1=100*(x1(2)-x1(1)^2)^2+(1-x1(1))^2 计算结果：最优点：x1 = -0.99611.0039步长： d = 0.0039最优解：f1 = 3.9981目标函数：function f = fun2( x )f=x(1)^2-2*x(1)*x(2)+2*x(2)^2+x(3)^2-x(2)*x(3)+2*x(1)+3*x(2 )-x(3);end目标函数梯度：function g = gfun2( x )g=[2 -2 0;-2 4 -1;0 -1 2]*x+[2;3;-1];end源代码：（以初值为为（0；0；0））x0=[0;0;0]; %初始值eps=1.0e-5; %精度g0=gfun2(x0);s0=-g0;n=0;syms d1;while norm(g0)>epsif n<3g=gfun2(x0+d1*s0);d= double(solve(s0'*g));x1=x0+d*s0;g1=gfun2(x1);if norm(g1)<epsn=n+1;x0=x1;breakelses0=-g1+(norm(g1)^2/norm(g0)^2)*s0;x0=x1;g0=g1;endelseif n==3x0=x1;g0=gfun2(x0);s0=-g0;n=0;endn=n+1;x0nfun2(x0)计算结果：最优点：x0 = -4 -3 -1 迭代次数： n = 3 最优值： ans = -8（1）最速下降法：目标函数：function f= fun3_1(x )f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);end目标函数的梯度：function g= gfun3_1(x)g=[1+2*x(1)*exp(x(1)^2+x(2)^2);4*x(2)+2*x(2)*exp(x(1)^2+x(2 )^2)];end源代码（初值为（1；0））：x0=[1;0];%初始值eps=1.0e-5;%精度n=0;g0=gfun3(x0);syms d1;while norm(g0)>=epss0=-g0;g=gfun3(x0+d1*s0);d= double(solve(s0'*g));x1=x0+d*s0;g1=gfun3(x1);if( norm(g1)<eps)n=n+1;x0=x1;break;elsex0=x1;g0=gfun3(x0);endn=n+1;endf0=fun3(x0)x0n计算结果：最优值： f0 = 0.7729 最优点：x0 = -0.4194 0迭代次数：n = 1（2）牛顿法目标函数：function f= fun3_1(x )f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);end目标函数梯度：function g= gfun3_1(x)g=[1+2*x(1)*exp(x(1)^2+x(2)^2);4*x(2)+2*x(2)*exp(x(1)^2+x(2 )^2)];end目标函数的Hesse阵：function g2 = g2fun3(x)g2=[2*exp(x(1)^2+x(2)^2)+4*x(1)^2*exp(x(1)^2+x(2)^2),4*x(1) *x(2)*exp(x(1)^2+x(2)^2)4*x(1)*x(2)*exp(x(1)^2+x(2)^2),4+2*exp(x(1)^2+x(2)^2)+4*x(2 )^2*exp(x(1)^2+x(2)^2)];end源代码（初值为（1；0））：x0=[1;0];%初始值eps=1.0e-5;%精度n=0;g0=gfun3(x0);g20=g2fun3(x0);while norm(g0)>=epsd=-g20\g0;x1=x0+d;g1=gfun3(x1);if( norm(g1)<eps)n=n+1;x0=x1;break;elsex0=x1;g0=gfun3(x0);endn=n+1;endf0=fun3(x0)x0n计算结果：最优值：f0 = 0.7729最优点：x0 = -0.4194迭代次数：n = 63（3）利用BGFS法：目标函数：function f= fun3_1(x )f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);end目标函数梯度function g= gfun3_1(x)g=[1+2*x(1)*exp(x(1)^2+x(2)^2);4*x(2)+2*x(2)*exp(x(1)^2+x(2 )^2)];end源代码（初值为（1；0））：x0=[1;0];%初始值eps=1.0e-5;%精度n=0;g0=gfun3(x0);syms d1;h0=eye(2);while norm(g0)>=epss0=-h0*g0;g=gfun3(x0+d1*s0);d= double(solve(s0'*g));x1=x0+d*s0;g1=gfun3(x1);if( norm(g1)<eps)n=n+1;x0=x1;break;elseh0=h0-(h0*(g1-g0)*(g1-g0)'*h0)/((g1-g0)'*h0*(g1-g0))+...((x1-x0)*(x1-x0)')/((x1-x0)'*(g1-g0))+...((g1-g0)'*h0*(g1-g0))*((x1-x0)*(x1-x0)');x0=x1;g0=gfun3(x0);endn=n+1;endf0=fun3(x0)x0n计算结果：最优值：f0 = 0.7729最优点：x0 = -0.4194迭代次数：n = 1题四：求解子程序：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;endend源代码：H=[2 0;0 2];%目标函数的hesse阵c=[-2,-4];Ae=[0 0];be=[0;0];Ai=[1/2 0;-1 3];bi=[1;2];x0=[-1;0];%初始值内部点eps=1.0e-9; %µ±ax-b=epsÊ±µ±×öax-b=0´¦Àíerr=1.0e-6;k=0;x0(1)=x0(1)+x0(2);x=x0;n=length(x);max=1.0e3;ne=length(be);ni=length(bi);lamk=zeros(ne+ni,1);index=ones(ni,1);for i=(1:ni)if(bi(i)-Ai(i,:)*x>eps)index(i)=0;。

matlab中active-set最优化算法

matlab中active-set最优化算法下载温馨提示：该文档是我店铺精心编制而成，希望大家下载以后，能够帮助大家解决实际的问题。

文档下载后可定制随意修改，请根据实际需要进行相应的调整和使用，谢谢!并且，本店铺为大家提供各种各样类型的实用资料，如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等，如想了解不同资料格式和写法，敬请关注!Download tips: This document is carefully compiled by the editor. I hope that after you download them, they can help you solve practical problems. The document can be customized and modified after downloading, please adjust and use it according to actual needs, thank you!In addition, our shop provides you with various types of practical materials, such as educational essays, diary appreciation, sentence excerpts, ancient poems, classic articles, topic composition, work summary, word parsing, copy excerpts, other materials and so on, want to know different data formats and writing methods, please pay attention!Matlab中Active-set最优化算法是一种常用的数学优化方法，用于解决非线性最优化问题。

基于MATLAB与ANSYS的结构优化设计

当前比较成熟的各种遗传操作算子，借助它可以方便地完成各种问题的优化。为使遗传算法更高效的应用于结构优化设计，研究了在ＭＡＬＢ中调用ＡＳＳ的方法，ＴＡＮＹ实现了ＭＡＬＢ与ＡＳＳ的数据传递，ＴＡＮＹ并
用该方法对一钢框架结构进行了优化设计，验证了此方法的可行性。关键词：遗传算法；优化设计；ＮＹ；ＴＡＡＳＳＭＡＬＢ
Ａｂｓｒｃ：ｈｅｅｅｉａｇｒｔｍｓｏｅｋｎｆｉｔｌｇｎｐｉｚｄａｇｒｔｍｉｈｅｅｏｓｔａｔＴｇｎｔｃｌｏｈｉｎｉｄｏｎｅｌｅｔｏｔｍｉｅｌｏｈｗｈｃｄｖｌｐｉｉｉ
ｓｂｌｙｏｈｓｍｅｈｄ．ｉｉｉｆｔｉｔｏｔＫｅｒｓ：ｅｅｉｌｏｔｍ；ｐｉｙｗｏｄｇｎｔｃａｇｒｈｏｔｍｕｍｅｉｎ；ｉｄｓｇＡＮＳＹＳ；ＭＡＴＬＡＢ
伴随着数学、力学和计算机的发展，结构优化设计也逐渐发展、熟起来。ＡＳＳ是最早开发成ＮＹ
ｇａｕｌｎｒｃｎｅｒ．ＩｈｓｓｒｎｖｒｌｓｔａｉｎｓａｃｂｌｙＩａｓａｏｖｓｒｄａｌｉｅｅｔｙａｓｔａｔｇｏｅａｌｉｔｅｒｈａｉｔ．ｔｌｏｃｎｓｌｅｍｏｔｙｏｕｏｉｐｏｌｍｓｏｅｏｔｚｔｎｏｅｅｇｎｅｎ ¨－ｅｃｅｔｎｃｕａｅｙｔｒｕｈｔｅｃｍｂ－ｒｂｅｆｔｐｉａｉｆｈｎｉｅｒｇｈｍｉｏｔｉ３ｆｉｎｌａｄａｃｒｔｌｈｏｇｈｏｉｊｉｙ

大连理工大学秋季优化方法大作业

m=m+1; end x0=x0+rho^mk*d; val=feval(fun,x0); g0=g; d0=d; k=k+1; end x=x0; val=feval(fun,x);
//f(x)//
function f=fun(x) x1=[1 0]*x; x2=[0 1]*x; f=(1-x1)^2+100*(x2-x1^2)^2; //梯度函数// function g=gfun(x) x1=[1 0]*x; x2=[0 1]*x; g=[-2*(1-x1)-400*x1*(x2-x1^2); 200*(x2-x1^2)]; //运行过程// >> x0=[0 0]'
k=k+1; btaold=btak; x0=x; end f=feval(fun,x); output.fval=f; output.iter=k; output.inner_iter=ink; output.bta=btak;
//增广拉格朗日函数//
function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma) f=feval(fun,x); he=feval(hf,x); gi=feval(gf,x); l=length(he); m=length(gi); psi=f; s1=0.0; for(i=1:l) psi=psi-he(i)*mu(i); s1=s1+he(i)^2; end psi=psi+0.5*sigma*s1; s2=0.0; for(i=1:m) s3=max(0.0, lambda(i) - sigma*gi(i)); s2=s2+s3^2-lambda(i)^2; end psi=psi+s2/(2.0*sigma); //增广拉格朗日函数// function dpsi=dmpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma) dpsi=feval(dfun,x); he=feval(hf,x); gi=feval(gf,x); dhe=feval(dhf,x); dgi=feval(dgf,x); l=length(he); m=length(gi); for(i=1:l) dpsi=dpsi+(sigma*he(i)-mu(i))*dhe(:,i); end for(i=1:m) dpsi=dpsi+(sigma*gi(i)-lambda(i))*dgi(:,i); end //f(x)// function f=f1(x) f=4*x(1)-x(2)^2-12; //等式约束// function he=h1(x) he=25-x(1)^2-x(2)^2;

优化设计Matlab程序

进退法步骤：1. 给定初始点(0)x ，初始步长0h ，令(1)(0)0,,0h h x x k ===2.令(4)(1),1x x h k k =+=+ 3.若(4)(1)()()f x f x <，则转4，否则转5 4.(2)(1)(1)(4)(2)(1)(1)(4),,()(),()()x x x x f x f x f x f x ====，令h =2h ，转2 5.若k =1，则转6，否则转,7 6.令h =-h ，(2)(4)(2)(4),()()x x f x f x ==，转2 7. 令(3)(2)(2)(1)(1)(4),,x x x x x x ===，停止计算，极小点包含于区间(1)(3)[,]x x 或(3)(1)[,]x x%目标函数：f%初始点：x0%初始步长：h0%精度：eps%目标函数取包含极值的区间左端点：minx%目标函数取包含极值的区间右端点：maxxfunction [minx,maxx] = minJT(f,x0,h0,eps)format long ;if nargin == 3eps = 1.0e-6;endx1 = x0;k = 0;h = h0;while 1x4 = x1 + h;k = k+1;f4 = subs(f, findsym(f),x4); ! subs : Symbolic substitution in symbolic expression or matrixf1 = subs(f, findsym(f),x1); ! findsym : Determine variables in symbolic expression or matrixif f4 < f1x2 = x1;x1 = x4;f2 = f1;f1 = f4;h = 2*h;elseif k==1h = -h;x2 = x4;f2 = f4;elsex3 = x2;x2 = x1;x1 = x4;break;endendendminx = min(x1,x3); maxx = x1+x3 - minx; format short;syms t;f=t^4-t^2-2*t+5;[x1,x2]=minJT(f,0,0.1)黄金分割法：% 目标函数：f% 极值区间左端点：a% 极值区间右端点：b% 精度：eps% 目标函数取最小值时的自变量值：x % 目标函数的最小值：minffunction [x,minf] = minHJ(f,a,b,eps) format long;if nargin == 3eps = 1.0e-6;endl = a + 0.382*(b-a);u = a + 0.618*(b-a);k=1;tol = b-a;while tol>eps && k<100000fl = subs(f , findsym(f), l);fu = subs(f , findsym(f), u);if fl > fua = l;l = u;u = a + 0.618*(b - a);elseb = u;u = l;l = a + 0.382*(b-a);endk = k+1;tol = abs(b - a);endif k == 100000disp('找不到最小值');x = NaN;minf = NaN;return;endx = (a+b)/2;minf = subs(f, findsym(f),x);format short;抛物线法：% 目标函数：f% 极值区间左端点：a% 极值区间右端点：b% 精度：eps% 目标函数取最小值时的自变量值：x% 目标函数的最小值：minffunction [x,minf] = minPWX(f,a,b,eps)format long;if nargin == 3eps = 1.0e-6;endt0 = (a+b)/2;k = 0;tol = 1;while tol>epsfa = subs(f,findsym(f),a);fb = subs(f,findsym(f),b);ft0 = subs(f,findsym(f),t0);tu = fa*(b^2 - t0^2)+fb*(t0^2 - a^2)+ft0*(a^2 - b^2);td = fa*(b - t0)+fb*(t0 - a)+ft0*(a - b);t1 = tu/2/td;ft1 = subs(f,findsym(f),t1);tol = abs(t1 - t0);if ft1 <= ft0if t1<= t0b = t0;t0 = t1;elsea = t0;t0 = t1;endk = k+1;elseif t1<= t0a = t1;elseb = t1;endk = k+1;endendx = t1;minf = subs(f,findsym(f),x); format short;一维牛顿法：% 目标函数：f% 初始点：x0% 精度：eps% 目标函数取最小值时的自变量值：x % 目标函数的最小值：minffunction [x,minf] = minNewton(f,x0,eps) format long;if nargin == 2eps = 1.0e-6;enddf = diff(f);d2f = diff(df);k = 0;tol = 1;while tol>epsdfx = subs(df,findsym(df),x0);d2fx=subs(d2f,findsym(d2f),x0;x1=x0-dfx/d2fx;k = k + 1;tol = abs(dfx);x0 = x1;endx = x1;minf = subs(f,findsym(f),x);format short;最速下降法：%: 目标函数：f%: 初始点：x0%: 自变量向量：var%: 精度：eps%: 目标函数取最小值时的自变量值：x%: 目标函数的最小值：minffunction [x,minf] = minFD(f,x0,var,eps)format long;if nargin == 3eps = 1.0e-6;endsyms l;tol = 1;gradf = - jacobian(f,var);while tol>epsv = Funval(gradf,var,x0);!Objective function value of the current point tol = norm(v); ! Vector and matrix normsy = x0 + l*v;yf = Funval(f,var,y);[a,b] = minJT(yf,0,0.1); %初始点0，步长为0.1xm = minHJ(yf,a,b);x1 = x0 + xm*v;x0 = x1;endx = x1;minf = Funval(f,var,x);format short;用共轭梯度法求无约束问题minf(x),n x R ∈的算法步骤如下：1. 给定初始点(0)x ，及精度ε2. 若(0)()f x ε∇≤，停止，极小点为(0)x ，否则转3.3. 取(0)(0)()d f x =-∇，且置k=04. 用一维搜索法求k α，使得()()()()()min ()k k k k k f x d f x d αα+=+令(1)()()k k k x x d α+=+，转55. 若(1)()k f x ε+∇≤，停止，极小点为(1)k x +，否则转66. 若k+1=n ，令(0)()n x x =转3，否则转77. 令2(1)(1)(1)()2()()(),()k k k k k k k f x d f x d f x λλ+++∇=-∇+=∇，置1k k =+，转4程序举例：function [x,minf] = minGETD(f,x0,var,eps) format long ;if nargin == 3eps = 1.0e-6;endx0 = transpose(x0);n = length(var);syms l ;gradf = jacobian(f,var);v0 = Funval(gradf,var,x0);d = -transpose(v0);k = 0;while 1v = Funval(gradf,var,x0);tol = norm(v);if tol<=epsx = x0;break;endy = x0 + l*d;yf = Funval(f,var,y);[a,b] = minJT(yf,0,0.1);xm = minHJ(yf,a,b);x1 = x0 + xm*d;vk = Funval(gradf,var,x1);tol = norm(vk);if tol<=epsx = x1;break;endif k+1==nx0 = x1;continue;elselamda = dot(vk,vk)/dot(v,v);d = -transpose(vk) + lamda*d;k = k+1;x0 = x1;endendminf = Funval(f,var,x);format short;用牛顿法求无约束问题，步骤：1. 给定初始点(0)x ，及精度ε2. 若(0)(f x ε∇≤，停止，极小点为(0)x ，否则转3.3. 计算12()()k f x -⎡⎤∇⎣⎦，令2()1()[()]();k k k d f x f x -=-∇∇ 4. 令(1)()(),1k k k x x d k k +=+=+，转2程序举例：function [x,minf] = minNT(f,x0,var,eps) format long ;if nargin == 3eps = 1.0e-6;endtol = 1;x0 = transpose(x0);gradf = jacobian(f,var);jacf = jacobian(gradf,var);while tol>epsv = Funval(gradf,var,x0); tol = norm(v);pv = Funval(jacf,var,x0);p = -inv(pv)*transpose(v);p = double(p);x1 = x0 + p;x0 = x1;endx = x1;minf = Funval(f,var,x);format short;syms t s;f=(t-4)^2+(s+2)^2+1;[x, mf]=minNewton(f, [0 0],[t s])。

最优化方法的Matlab实现(公式(完整版))

第九章最优化方法的Matlab实现在生活和工作中，人们对于同一个问题往往会提出多个解决方案，并通过各方面的论证从中提取最佳方案。

最优化方法就是专门研究如何从多个方案中科学合理地提取出最佳方案的科学。

由于优化问题无所不在，目前最优化方法的应用和研究已经深入到了生产和科研的各个领域，如土木工程、机械工程、化学工程、运输调度、生产控制、经济规划、经济管理等，并取得了显著的经济效益和社会效益。

用最优化方法解决最优化问题的技术称为最优化技术，它包含两个方面的内容：1）建立数学模型即用数学语言来描述最优化问题。

模型中的数学关系式反映了最优化问题所要达到的目标和各种约束条件。

2）数学求解数学模型建好以后，选择合理的最优化方法进行求解。

最优化方法的发展很快，现在已经包含有多个分支，如线性规划、整数规划、非线性规划、动态规划、多目标规划等。

9.1 概述利用Matlab的优化工具箱，可以求解线性规划、非线性规划和多目标规划问题。

具体而言，包括线性、非线性最小化，最大最小化，二次规划，半无限问题，线性、非线性方程（组）的求解，线性、非线性的最小二乘问题。

另外，该工具箱还提供了线性、非线性最小化，方程求解，曲线拟合，二次规划等问题中大型课题的求解方法，为优化方法在工程中的实际应用提供了更方便快捷的途径。

9.1.1 优化工具箱中的函数优化工具箱中的函数包括下面几类：1．最小化函数表9-1 最小化函数表2．方程求解函数表9-2 方程求解函数表3．最小二乘（曲线拟合）函数表9-3 最小二乘函数表4．实用函数表9-4 实用函数表5．大型方法的演示函数表9-5 大型方法的演示函数表6．中型方法的演示函数表9-6 中型方法的演示函数表9.1.3 参数设置利用optimset函数，可以创建和编辑参数结构；利用optimget函数，可以获得o ptions优化参数。

● optimget函数功能：获得options优化参数。

语法：val = optimget(options,'param')val = optimget(options,'param',default)描述：val = optimget(options,'param') 返回优化参数options中指定的参数的值。

大连理工大学庞丽萍最优化方法MATLAB程序

班级：优化1班授课老师：庞丽萍姓名：学号：第二章12.（1）用修正单纯形法求解下列LP问题:>>clear>>A=[121100;123010;215001];[m,n]=size(A);b=[10;15;20];r=[-1-2-31];c=[-1-2-31];bs=[3:3];nbs=[1:4];a1=A(:,3);T=A(:,bs);a2=inv(T)*a1;b=inv(T)*b;A=[eye(m),a2];B=eye(m);xb=B\b;cb=c(bs);cn=c(nbs);con=1;M=zeros(1);while conM=M+1;t=cb/B;r=c-t*A;if all(r>=0)x(bs)=xb;x(nbs)=0;fx=cb*xb;disp(['当前解是最优解，minz=',num2str(fx)])disp('对应的最优解为,x=')disp(x)breakendrnbs=r(nbs);kk=find(rnbs==min(rnbs));k=kk(1);Anbs=A(:,nbs);yik=B\Anbs(:,k);xb=B\b;%yi0if all(yik<=0)disp('此LP问题无有限的最优解，计算结束',x)disp(xb)breakelsei=find(yik>0);w=abs(xb(i,1)./yik(i,1));l=find(w==min(w));rr=min(l);yrrk=yik(rr,1);Abs=A(:,bs);D=Anbs(:,k);Anbs(:,k)=Abs(:,rr);Abs(:,rr)=D;F=bs(rr);bs(rr)=nbs(k);nbs(k)=F;AA=[Anbs,Abs];EE=eye(m);EE(:,rr)=-yik./yrrk;Errk=EE;Errk(rr,rr)=1/yrrk;BB=Errk/B;B=inv(BB);cb=c(:,bs);xb=Errk*xb;x(bs)=xb;x(nbs)=0;fx=cb*xb;endif M>=1000disp('此问题无有限最优解')breakendend%结果当前解是最优解，minz=-15对应的最优解为,x=2.5000 2.5000 2.50000第三章30题DFP算法求函数极小点的计算程序function[x,val,k]=dfp(fun,gfun,x0)%功能:用DFP算法求解无约束问题:minf(x)%输入:x0是初始点,fun,gfun分别是目标函数及其梯度%输出:x,val分别是近似最优点和最优值,k是迭代次数.maxk=1e5;%给出最大迭代次数rho=0.55;sigma=0.4;epsilon=1e-5;k=0;n=length(x0);Hk=inv(feval('Hess',x0));%Hk=eye(n);while(k<maxk)gk=feval(gfun,x0);%计算梯度if(norm(gk)<epsilon),break;end%检验终止准则dk=-Hk*gk;%解方程组,计算搜索方向m=0;mk=0;while(m<20)%用Armijo搜索求步长if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk’*dk)mk=m;break;endm=m+1;end%DFP校正x=x0+rho^mk*dk;sk=x-x0;yk=feval(gfun,x)-gk;if(sk'*yk>0)Hk=Hk-(Hk*yk*yk'*Hk)/(yk'*Hk*yk)+(sk*sk')/(sk'*yk);endk=k+1;x0=x;endval=feval(fun,x0);%习题26的程序调用方式及结果：function y=fun(x)%UNTITLED Summary of this function goes here%Detailed explanation goes herey=(x(1)-1)^2+5*(x2-x(1)^2)^2endfunction y=gfun(x)%UNTITLED Summary of this function goes here%Detailed explanation goes herey=[diff(y,x1)diff(y,x2)]endx0=[20]’;[x,val,k]=dfp(fun,gfun,x0)%结果x=1.000001.00000val=k=6%习题27的程序调用方式及结果：function y=fun(x)%UNTITLED Summary of this function goes here %Detailed explanation goes herey=x1+2*x（2）^2+exp(x(1)^2+x(2)^2)endfunction y=gfun(x)%UNTITLED Summary of this function goes here %Detailed explanation goes herey=[diff(y,x1)diff(y,x2)]endx0=[10]’;[x,val,k]=dfp(fun,gfun,x0)%结果x=-0.419360val=0.77291k=536题编写Hooke-Jeeves方法求函数极小点的计算程序。

大连理工大学矩阵与数值分析MATLAB上机实验

二、解线性方程组 1．分别 Jacobi 迭代法和 Gauss-Seidel 迭代法求解线性方程组
3 1 0 0 x1 1 1 3 1 0 x2 0 , 0 1 2 1 x3 0 0 0 1 3 x4 0
Gauss 列主元消去法程序:
clc; clear; format long a=[2,4,3,1;8,2,0,0;5,0,4,0;9,0,0,5]; %系数矩阵 b=[12;6;23;16]; [n,m]=size(a); nb=length(b); det=1; for k=1:n-1 amax=0; for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k)); r=i; end end if amax<1e-10 return; end if r>k for j=k:n z=a(k,j); a(k,j)=a(r,j); a(r,j)=z; end z=b(k);
从小到大求和程序计算结果：
N 100 10000 1000000 从小到大求和程序得到的 SN 0.497512437810945 0.499975001249937 0.499999750000134 真实值�� = �� 2�� + 1 0.497512437810945 0.499975001249937 0.499999750000125 计算值有效位数 15 15 13
8
2
1 dx x
复化梯形公式程序
clc; clear; format long syms t m=int(1/t,2,8); %真实值 a=2; b=8; n=300; h=(b-a)/n; sum=0; f=inline('1/x'); for i=1:n-1 sum=sum+f(a+i.*h); end T=h/2*(f(a)+2*sum+f(b))

大连理工大学优化方法上机大作业

学院：专业：班级：学号：姓名：上机大作业1：1.最速下降法：function f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);endfunction x_star = steepest(x0,eps) gk = grad(x0);res = norm(gk);k = 0;while res > eps && k<=1000dk = -gk;ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + *ak*slopeak = ak/4;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;x0 = xk;gk = grad(xk);res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res); endx_star = xk;end>> clear>> x0=[0,0]';>> eps=1e-4;>> x=steepest(x0,eps)2.牛顿法：function f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; endfunction g = grad2(x)g = zeros(2,2);g(1,1)=2+400*(3*x(1)^2-x(2));g(1,2)=-400*x(1);g(2,1)=-400*x(1);g(2,2)=200;endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);endfunction x_star = newton(x0,eps)gk = grad(x0);bk = [grad2(x0)]^(-1);res = norm(gk);k = 0;while res > eps && k<=1000dk=-bk*gk;xk=x0+dk;k = k+1;x0 = xk;gk = grad(xk);bk = [grad2(xk)]^(-1);res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res); endx_star = xk;end>> clear>> x0=[0,0]';>> eps=1e-4;>> x1=newton(x0,eps)--The 1-th iter, the residual is--The 2-th iter, the residual isx1 =法：function f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);endfunction x_star = bfgs(x0,eps) g0 = grad(x0);gk=g0;res = norm(gk);Hk=eye(2);k = 0;while res > eps && k<=1000dk = -Hk*gk;ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + *ak*slopeak = ak/4;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;fa0=xk-x0;x0 = xk;go=gk;gk = grad(xk);y0=gk-g0;Hk=((eye(2)-fa0*(y0)')/((fa0)'*(y0)))*((eye(2)-(y0)*(fa0)')/((fa0)'*(y0)))+(fa0*(fa 0)')/((fa0)'*(y0));res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res);endx_star = xk;End>> clear>> x0=[0,0]';>> eps=1e-4;>> x=bfgs(x0,eps)4.共轭梯度法：function f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);endfunction x_star =CG(x0,eps) gk = grad(x0);res = norm(gk);k = 0;dk = -gk;while res > eps && k<=1000 ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + *ak*slope ak = ak/4;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;x0 = xk;g0=gk;gk = grad(xk);res = norm(gk);p=(gk/g0)^2;dk1=dk;dk=-gk+p*dk1;fprintf('--The %d-th iter, the residual is %f\n',k,res); endx_star = xk;end>> clear>> x0=[0,0]';>> eps=1e-4;>> x=CG(x0,eps)上机大作业2：function f= obj(x)f=4*x(1)-x(2)^2-12;endfunction [h,g] =constrains(x) h=x(1)^2+x(2)^2-25;g=zeros(3,1);g(1)=-10*x(1)+x(1)^2-10*x(2)+x(2)^2+34;g(2)=-x(1);g(3)=-x(2);endfunction f=alobj(x) %拉格朗日增广函数%N_equ等式约束个数?%N_inequ不等式约束个数N_equ=1;N_inequ=3;global r_al pena;%全局变量h_equ=0;h_inequ=0;[h,g]=constrains(x);%等式约束部分?for i=1:N_equh_equ=h_equ+h(i)*r_al(i)+(pena/2)*h(i).^2;end%不等式约束部分for i=1:N_inequh_inequ=h_inequ+pena)*(max(0,(r_al(i)+pena*g(i))).^2-r_al(i).^2); end%拉格朗日增广函数值f=obj(x)+h_equ+h_inequ;function f=compare(x)global r_al pena N_equ N_inequ;N_equ=1;N_inequ=3;h_inequ=zeros(3,1);[h,g]=constrains(x);%等式部分for i=1:1h_equ=abs(h(i));end%不等式部分for i=1:3h_inequ=abs(max(g(i),-r_al(i+1)/pena));endh1 = max(h_inequ);f= max(abs(h_equ),h1); %sqrt(h_equ+h_inequ);function [ x,fmin,k] =almain(x_al)%本程序为拉格朗日乘子算法示例算法%函数输入：% x_al:初始迭代点% r_al：初始拉格朗日乘子N-equ:等式约束个数N_inequ:不等式约束个数?%函数输出% X：最优函数点FVAL:最优函数值%============================程序开始================================ global r_al pena ; %参数（全局变量）pena=10; %惩罚系数r_al=[1,1,1,1];c_scale=2; %乘法系数乘数cta=; %下降标准系数e_al=1e-4; %误差控制范围max_itera=25;out_itera=1; %迭代次数%===========================算法迭代开始============================= while out_itera<max_iterax_al0=x_al;r_al0=r_al;%判断函数?compareFlag=compare(x_al0);%无约束的拟牛顿法BFGS[X,fmin]=fminunc(@alobj,x_al0);x_al=X; %得到新迭代点%判断停止条件?if compare(x_al)<e_aldisp('we get the opt point');breakend%c判断函数下降度?if compare(x_al)<cta*compareFlagpena=1*pena; %可以根据需要修改惩罚系数变量elsepena=min(1000,c_scale*pena); %%乘法系数最大1000disp('pena=2*pena');end%%?更新拉格朗日乘子[h,g]=constrains(x_al);for i=1:1%%等式约束部分r_al(i)= r_al0(i)+pena*h(i);endfor i=1:3%%不等式约束部分r_al(i+1)=max(0,(r_al0(i+1)+pena*g(i)));endout_itera=out_itera+1;end%+++++++++++++++++++++++++++迭代结束+++++++++++++++++++++++++++++++++ disp('the iteration number');k=out_itera;disp('the value of constrains'); compare(x_al)disp('the opt point');x=x_al;fmin=obj(X);>> clear>> x_al=[0,0];>> [x,fmin,k]=almain(x_al)上机大作业3： 1、>> clear alln=3; c=[-3,-1,-3]'; A=[2,1,1;1,2,3;2,2,1;-1,0,0;0,-1,0;0,0,-1];b=[2,5,6,0,0,0]';cvx_beginvariable x(n)minimize( c'*x)subject toA*x<=bcvx_endCalling SDPT3 : 6 variables, 3 equality constraints------------------------------------------------------------num. of constraints = 3dim. of linear var = 6*******************************************************************SDPT3: Infeasible path-following algorithms*******************************************************************version predcorr gam expon scale_dataNT 1 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime -------------------------------------------------------------------0|||+01|+00|+02|+01 +00| 0:0:00| chol 1 11|||||+01|+00 +01| 0:0:01| chol 1 12|||||+00|+00 +01| 0:0:01| chol 1 13|||||+00|+00 +00| 0:0:01| chol 1 14||||||+00 +00| 0:0:01| chol 1 15||||||+00 +00| 0:0:01| chol 1 16||||||+00 +00| 0:0:01| chol 1 17||||||+00 +00| 0:0:01| chol 1 18||||||+00 +00| 0:0:01|stop: max(relative gap, infeasibilities) <------------------------------------------------------------------- number of iterations = 8primal objective value = +00dual objective value = +00gap := trace(XZ) =relative gap =actual relative gap =rel. primal infeas (scaled problem) =rel. dual " " " =rel. primal infeas (unscaled problem) = +00rel. dual " " " = +00norm(X), norm(y), norm(Z) = +00, +00, +00norm(A), norm(b), norm(C) = +00, +00, +00Total CPU time (secs) =CPU time per iteration =termination code = 0DIMACS: +00 +00-------------------------------------------------------------------------------------------------------------------------------Status: SolvedOptimal value (cvx_optval):2、>> clear alln=2; c=[-2,-4]'; G=[,0;0,1]; A=[1,1;-1,0;0,-1]; b=[1,0,0]'; cvx_beginvariable x(n)minimize( x'*G*x+c'*x)subject toA*x<=bcvx_endCalling SDPT3 : 7 variables, 3 equality constraintsFor improved efficiency, SDPT3 is solving the dual problem.------------------------------------------------------------num. of constraints = 3dim. of socp var = 4, num. of socp blk = 1dim. of linear var = 3*******************************************************************SDPT3: Infeasible path-following algorithms*******************************************************************version predcorr gam expon scale_dataNT 1 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime -------------------------------------------------------------------0||||+00|+02| +01 +00| 0:0:00| chol 1 11|||||+01| +00 | 0:0:00| chol 1 12|||||+00| +00 | 0:0:00| chol 1 13|||||| | 0:0:00| chol 1 14|||||| | 0:0:00| chol 1 15|||||| | 0:0:00| chol 1 16|||||| | 0:0:00| chol 1 17|||||| | 0:0:00| chol 1 18|||||| | 0:0:00| chol 1 19|||||| | 0:0:00| chol 1 110|||||| | 0:0:00| chol 2 211|||||| | 0:0:00| chol 2 212|||||| | 0:0:00| chol 2 213|||||| | 0:0:00| chol 2 214|||||| | 0:0:00|stop: max(relative gap, infeasibilities) <------------------------------------------------------------------- number of iterations = 14primal objective value =dual objective value =gap := trace(XZ) =relative gap =actual relative gap =rel. primal infeas (scaled problem) =rel. dual " " " =rel. primal infeas (unscaled problem) = +00rel. dual " " " = +00norm(X), norm(y), norm(Z) = +00, +00, +00norm(A), norm(b), norm(C) = +00, +00, +00Total CPU time (secs) =CPU time per iteration =termination code = 0DIMACS: +00 +00-------------------------------------------------------------------------------------------------------------------------------Status: SolvedOptimal value (cvx_optval): -3。

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代入任意初始值，运算结果如下。