大连理工优化方法 增广拉格朗日方法MATLAB程序

合集下载

matlab解决凸优化和拉格朗日对偶方法

matlab解决凸优化和拉格朗日对偶方法

matlab解决凸优化和拉格朗日对偶方法
Matlab是一个强大的数值计算和科学编程工具,提供了丰富的函数和工
具箱来解决各种数学优化问题,包括凸优化和拉格朗日对偶方法。

在Matlab中,我们可以使用内置函数和工具箱来解决凸优化问题。

凸优
化是一类非常重要且广泛应用的数学优化问题,其目标是最小化或最大化凸
函数,在给定一些约束条件下,求解最优解。

Matlab中最常用的凸优化函数是"cvx"工具箱。

该工具箱提供了一套简洁
而强大的函数,可以轻松地定义凸优化问题,并使用内置的求解算法进行求解。

通过该工具箱,用户可以快速解决线性规划、二次规划、半定规划和凸
二次规划等问题。

除了凸优化,Matlab也提供了功能强大的函数来解决拉格朗日对偶方法。

拉格朗日对偶方法是一种用于解决约束优化问题的有效技术。

它通过将原问
题转化为拉格朗日函数,并通过求解对偶问题来近似求解原问题。

在Matlab中,我们可以使用"quadprog"函数来解决带约束的二次规划问题,其中可通过添加约束条件和求解问题的对偶问题来实现拉格朗日对偶方法。

此外,Matlab还提供了其他一些函数和工具箱,如"fmincon"和"linprog",这些函数可以用于解决不同类型的优化问题。

Matlab是一个功能强大的工具,可以通过其内置函数和工具箱来解决凸
优化和拉格朗日对偶方法。

无论是解决线性规划问题还是非线性优化问题,Matlab都提供了易于使用且高效的求解方法,可以帮助研究人员和工程师解
决复杂的数学优化问题。

如何用Matlab写拉格朗日函数?

如何用Matlab写拉格朗日函数?

如何用Matlab写拉格朗日函数?
谢邀。

首先拉格朗日函数具体公式如下:
编写一个名为lagrange.m的M文件,然后设n个节点数据以数组x0, y0输入(注意Matlab的数组下标从1开始),m个插值点以数组x输入,输出数组y 为m个插值。

图片内容如下:
纯文本内容如下(可直接复制使用):
function y=lagrange(x0,y0,x)
n=length(x0);m=length(x);
for i=1:m
z=x(i);
s=0.0;
for k=1:n
p=1.0;
for j=1:n
if j~=k
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=p*y0(k)+s;
end
y(i)=s;
end
保存后调用编写的程序,并运行。

在Matlab的命令窗口输入【lagrange (x,y,xh)】按【Enter】键即可得到拉格朗日插值函数计算的插值。

如果你对科学和科技内容感兴趣,欢迎订阅我的头条号。

我会在这里发布所有与科技、科学有关的有趣文章。

偶尔也回答有趣的问题,有问题可随时在评论区回复和讨论,看到即回。

(码字不易,若文章对你帮助可点赞支持~)。

大连理工优化方法大作业MATLAB编程

大连理工优化方法大作业MATLAB编程

function [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)^2)^2+(1-x(1))^2;function gf=gfun(x)gf=[-4*x(1)*(x(2)-x(1)^2)+2*(x(1)-1),2*(x(2)-x(1)^2)]; 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 =54function f= fun( X )%所求问题目标函数f=X(1)^2-2*X(1)*X(2)+2*X(2)^2+X(3)^2+ X(4)^2-X(2)*X(3)+2*X(1)+3*X(2)-X(3);endfunction 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)];endfunction [ x,val,k ] = frcg( fun,gfun,x0 )%功能:用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(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d) mk=m;break;endm=m+1;endx0=x0+rho^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 =21或者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))^2)/(norm(g0)^2); 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(1)^2-2*x(1)*x(2)+2*x(2)^2+x(3)^2+x(4)^2-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 = 1function [f,x,k]=third_1(x)k=0;g=gfun(x);while(norm(g)>=0.001)s=-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)^2+exp(x(1)^2+x(2)^2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)];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};endend结果: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)^2+exp(x(1)^2+x(2)^2); function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)]; function ggf=ggfun(x)ggf=[(4*x(1)^2+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+(4*x(2)^2+2)*exp(x(1)^2+x(2)^2)]; 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;else d=2*d;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)dx=X{2}-X{1};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)^2+exp(x(1)^2+x(2)^2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)];function ggf=ggfun(x)ggf=[(4*x(1)^2+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+(4*x(2)^2+2)*exp(x(1)^2+x(2)^2);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;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;else d=2*d;endendend结果:x=[0,1];[f,x,k]=third_3(x)f =0.7729x = -0.4195 0.0000 k=6function 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,x0) function [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) epsilon=1.0e-9; err=1.0e-6;k=0; x=x0; n=length(x); kmax=1.0e3;ne=length(be); ni=length(bi); lamk=zeros(ne+ni,1); index=ones(ni,1);for (i=1:ni)if(Ai(i,:)*x>bi(i)+epsilon), index(i)=0; endendwhile(k<=kmax)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;endendendk=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; end endif(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(fun,hf,gf,dfun,dhf,dgf,x0)%功能: 用乘子法解一般约束问题: 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;ink=0;epsilon=0.00001;x=x0;he=feval(hf,x);gi=feval(gf,x);n=length(x);l=length(he);m=length(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(i)^2;end%更新乘子向量for i=1:mtemp=min(gi(i),lambda(i)/c);btak=btak+temp^2;endbtak=sqrt(btak);if btak>epsilonif k>=2&&btak>theta*btaoldc=eta*c;endfor i=1:lmu(i)=mu(i)-c*he(i);endfor i=1:mlambda(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);s1=s1+he(i)^2;endpsi=psi+0.5*c*s1;s2=0;for i=1:ms3=max(0,lambda(i)-c*gi(i));s2=s2+s3^2-lambda(i)^2;endpsi=psi+s2/(2*c);%不等式约束函数文件g1.mfunction gi=g1(x)gi=10*x(1)-x(1)^2+10*x(2)-x(2)^2-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)newf=feval(fun,x0+rho^m*dk,varargin{:});oldf=feval(fun,x0,varargin{:});if(newf<oldf+sigma*rho^m*gk'*dk)mk=m;break;endm=m+1;endx=x0+rho^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{:});结果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: 4f=[3,1,1];A=[2,1,1;1,-1,-1];b=[2;-1];lb=[0,0,0];x=linprog(f,A,b,zeros(3),[0,0,0]',lb)结果:Optimization terminated.x =0.00000.50000.5000。

数据增广 matlab程序

数据增广 matlab程序

数据增广 matlab程序
在MATLAB中进行数据增广可以采用多种方法,包括图像处理、数据生成和增强等技术。

以下是一个简单的示例程序,用于对图像数据进行水平翻转和旋转增广:
matlab.
% 读取原始图像数据。

originalImage = imread('original.jpg');
% 水平翻转。

flippedImage = fliplr(originalImage);
% 旋转增广。

rotatedImage = imrotate(originalImage, 45);
% 保存增广后的图像数据。

imwrite(flippedImage, 'flipped.jpg');
imwrite(rotatedImage, 'rotated.jpg');
这个简单的示例程序展示了如何使用MATLAB对图像数据进行水平翻转和旋转增广。

当然,实际的数据增广可能涉及到更复杂的处理,比如添加噪声、改变亮度和对比度等。

针对不同的数据类型和增广需求,具体的程序会有所不同。

除了图像处理,对于其他类型的数据,比如文本数据或者时间序列数据,可以采用一些统计学方法或者模拟方法进行增广。

MATLAB提供了丰富的工具和函数,可以帮助你实现各种类型的数据增广操作。

例如,你可以使用数据生成函数来创建符合特定分布的数据,然后对其进行变换和扰动以进行增广。

总的来说,数据增广是一个非常灵活和多样化的过程,具体的实现方法会根据数据类型和增广目的的不同而有所差异。

希望这个简单的示例程序能够给你一些启发,让你在MATLAB中实现更复杂的数据增广操作。

优化方法作业第二版

优化方法作业第二版

优化方法上机大作业院系:化工与环境生命学部姓名:李翔宇学号:31607007 指导教师:肖现涛第一题:编写程序求解下述问题min y(z) = (1 —xi)2+ 100(x2 —xj)2.初始点取/ = 0,精度取s = le - 4,步长由Armijo线捜索生成,方向分别由下列方法生成:最速下降法BFGS方法共辄梯度法1.最速下降法源程序如下:function x_star = ZSXJ (x0,eps) gk = grad(xO);res = norm(g k);k = 0;while res > eps && k<=10000 dk = -gk;ak =1; f0 = fun(xO);f1 = fun(xO+ak*dk);slope = dot(gk,dk);while f1 > f0 + 0.0001*ak*slope ak = ak/2;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);end x_star = xk; end function f = fun(x)f = (1-x(1))A2 + 100*以(2)咲(1)八2)八2;end function g = grad(x) g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)A2-x(2)); g(2) = 200*(x(2)-x(1)A2);end 运行结果:>> x0=[0,0]';>> esp=1e-4;>> xk= ZSXJ(x0,eps)--The 1-th iter, the residual is 13.372079 --The 2-th iter, the residual is 12.079876 --The 3-th iter, the residual is 11.054105 --The 9144-th iter, the residual is 0.000105 --The 9145-th iter, the residual is 0.000102 --The 9146-th iter, the residual is0.000100 xk =0.99990.9998MATLAB 截屏:>> xO= [Qj 01’ :>> esp=le-4:>> xk=ZSXJ txOj eps)一The I-th it er?the residual is13.372075一The2-th it er j the residual is12,07&876一Iht3-th it亡「■the residual is11. 054105一The Wh it er j the residual is10. 421221一一Ih?5-th it er^the residual is10. 020369Tl_ -■:丄一一丄ll_ ______ I; Ji ___ '1J _一The9142-th iter,the residual is0.000111-- The9143-th iterj the residual is0.000108一一The9144-th iter>the匸esidual is D. 000105一-The9145-th iter f th?residual is0. 000102一Th?9143-th iter,the residual is0. 0001000.99990.99932. 牛顿法源程序如下:function x_star = NEWTON (x0,eps) gk = grad(x0);bk = [grad2(xO)]A(-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)]A(-1);res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res); end x_star = xk;endfunction f = fun(x)f = (1-x(1))A2 + 100*(x(2)-x(1)A2)A2;end function g = grad2(x)g = zeros(2,2); g(1,1)=2+400*(3*x(1)A2-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)A2-x(2));g(2) = 200*(x(2)-x(1)A2);end运行结果:>> xO=[O,O]';eps=1e-4;>> xk= NEWTON (x0,eps)--The 1-th iter, the residual is 447.213595--The 2-th iter, the residual is 0.000000xk =1.00001.0000MATALB 截屏;»x0= [0, DY :>> esp=1;»xk^NEirONtaO, eps)—The 1-th iter, the residual is 447,213595 一The 2-th iter, the residual is0.OOQOQOxk =L 0000|L 00003. BFGS 方法源程序如下:function 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 + 0.1*ak*slopeak = ak/2;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;fa0=xk-x0;x0 = xk;g0=gk;gk = grad(xk); y0=gk-g0;Hk=((eye(2)-fa0*(y0)')/((fa0)'*(y0)))*((eye(2)-(y0)*(fa0)')/((fa0)'*(y0)))+(fa0*(fa0)')/((fa0)'*(y0));res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res); end x_star = xk; end function f=fun(x)f=(1-x(1))A2 + 100*債(2)咲(1)八2)八2;endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)A2-x(2));g(2) = 200*(x(2)-x(1)A2);end运行结果:>> x0=[0,0]';>> esp=1e-4;>> xk= Bfgs(x0,eps)--The 1-th iter, the residual is 3.271712--The 2-th iter, the residual is 2.381565--The 3-th iter, the residual is 3.448742--The 1516-th iter, the residual is 0.000368 --The 1517-th iter, the residual is0.000099 xk =1.00011.0002MATLAB 截屏:x0= [0, 0]J:esp=le-4: xk=Bfgs (xOj eps)一The1th it erj the residual is3. 271712一The2th it erj the residual is2.381565一The3th it the residual is3.448742一The4th it er,the residual is3, 162431— The5th it erj the residual is2. 989084—The 1515-th iter, the residual is 0.000108—The 1516-th iter, the residual is 0.000368—The 1517-th iter, the residual is 0. 0000991.00011.00024. 共轭梯度法源程序如下:function x_star =Conj (x0,eps)gk = grad(x0);res = norm(gk);k = 0;dk = -gk;while res > eps && k<=1000ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + 0.1*ak*slopeak = ak/2;xk = x0 + ak*dk;f1 = fun(xk);endd0=dk;g0=gk;k=k+1;xO=xk;gk=grad(xk);f=(norm(gk)/norm(gO))A2;res=norm(gk);dk=-gk+f*d0;fprintf('--The %d-th iter, the residual is %f\n',k,res);endx_star = xk;endfunction f=fun(x)f=(1-x(1))A2+100*(x(2)-x(1)A2)A2;endfunction g=grad(x)g=zeros(2,1);g(1)=400*x(1)A3-400*x(1)*x(2)+2*x(1)-2;g(2)=-200*x(1)A2+200*x(2);end运行结果:>> xO=[O,O]';>> eps=1e-4;>> xk=Conj(xO,eps)--The 1-th iter, the residual is 3.271712--The 2-th iter, the residual is 1.380542--The 3-th iter, the residual is 4.527780--The 4-th iter, the residual is 0.850596--The 73-th iter, the residual is 0.001532 --The 74-th iter, the residual is 0.000402 --The 75-th iter, the residual is0.000134 --The 76-th iter, the residual is 0.000057 xk =0.99990.9999MATLAB 截屏:» ito= O]T:>> esp=l*-4 :>> xk=Conj eps)—The l-th it efj the residual is 3. 271712~Ihe2-th it erj the residual IS 1. 3805423-th it er^the residual IS 4. 527780—I he4-th it er^the residual IS0.8505&6—I he5-th it erj the residual is0.559005——Tin 尺一+ Vi■i十戶T十Ha■r 戶 uT H IT-al H QRR7J.il■■丄丄比U 3^111丄LBL,LILS丄M th UUUUZJ一The70-th iter,the residual is0. 001423一The71-th iter,the residual is0. 002841---The72-th iter^the residual is0. 002945一The73-th itetj the residual is0. 001532—Th?74-th iter,the residual is0. 000403—The75-th it the residual is0. 000134一Iht76-th iter,the r*sidual is0.0000570.99990-9999第二题:编写程序利用增广拉格朗日方法求解下述问题初始点取= 0,精度取s= lc-4.解:目标函数文件fl.mfunction f=f1(x)f=4*x(1)-x(2)A2-12;等式约束函数文件hl.mfunction he=h1(x)he=25-x(1)A2-x(2)A2;不等式约束函数文件gl.mfunction gi=g1(x)gi=10*x(1)-x(1)A2+10*x(2)-x(2)A2-34;目标函数的梯度文件dfl.mfunction g=df1(x)g = [4, 2.0*x(2)]';等式约束(向量)函数的Jacobi矩阵(转置)文件dhl.mfunction dhe=dh1(x)dhe = [-2*x(1), -2*x(2)]';不等式约束(向量)函数的Jacobi矩阵(转置)文件dg1.mfunction dgi=dg1(x)dgi = [10-2*x(1), 10-2*x(2)]';然后在Matlab命令窗口输入如下命令:x0=[0,0]';[X,mu,lambda,output]=multphr('f1','h1','g1','df1','dh1','dg1',x0);得到如下输出:x =4.898717426488211.00128197198571算法编程利用程序调用格式function [x, mu, lambda, out put ] =mul t phr (furij hf, gf, dfun, dhfj dgf, xO)paxk=500;si gma=2* 0; eta z2. 0; theta z O* 8;k=0; ink=0;epsilon=le-4;x=xO; he-feval (hf f x); gi=feval (gf> x);n=length(K); 1=1 eng th (he); nFlength(gi),rru=O> l*ones(lj 1); lairbda=O» l*ones(n A 1); btak=10; btaold=10;^hile(btak>epsilon & k<maxk)[爲ival, ik]=bfgsC 叩si"," dirpsi", xO* fun, hf, gf, dfun, dhf, dgf, m f 1 ambda, sigma); ink二ink+ik;he=feval (h£ x); gi=feval (gf, x);btak-0,0; 丹for (1=1:1)^ btak=btak+he(i) "2; endfor i=l:mteitp=min(gi (i), lairbda(i)/sigma); btak-btak+temp"2;end btak z sqrt(btak); if btak>epsilon if(k>=2 & btak > theta*btaold) sig^eta+sigma, end for iDu[i)z mu(i)-sigina*he(i); endfor (i=l:nOlambda(i)z max(0.0, lairbda(i)-signia*gi(i));endendk二k+1;btaold=btak;xO=x;endf=f eval (fun, x);output. fral=f; output. iter=k;output* inner_iter-ink;output,bta-btak;function [x, val, k]=bfgs (fun^ gfun, xO, varargin) maxk=500;rho=0. 55; sigmal=O. 4; epsi1onl = le—5;k=0; n=length(xO);Bk=eye (n) ; %Bk=feval (' Hess^ , xO);while(k<maxk)gk=feval (gfun, xO, varargin {: });if(norm(gk)<epsilonl), break; enddk=-Bk\gk;n»=0; ink=O;while(m<20) 亠newf=feval (fun, xO+rho^irF^dk, varacrgin {: }); oldf=feval (fun, xO, varargin {: });i f (newf <ol df+si gir)al*rho m*gk,*dk) rnk=m; break;endm=irrH;endx= xO+rho ^ink*dk;sk=x—xO; yk=feval (gfun^ x, varargin {: } ) -gk;i f (yk^ *sk>0)Bk=Bk—(Bk*sk*sk J *Bk)/(sk‘ *Bk*sk)+(yk*yk J )/(yk‘ *sk); end k=k+l;xO=x;endval = f eval (fun, xO, v ar argin {: });第三题:■下载安装 CVX h http: /cvx/■利用CVX 编写代码求解下述问题迟1 + 牝 一 1 < 0. xi > 0.工2 > 0 ■利用CVX 编写代码求解下述问题min —3J :I — X2 — 3帀S.t. 2淤]十念2 +②3 W 2x\ + 2x2 + 3x3 < 5 2叭 +2x 2 + 巧 W6 37 > 0.1.解:将目标函数改写为向量形式:x'*a*x-b*x程序代码:n=2;a=[0.5,0;0,1];b=[2 4];c=[1 1];cvx_beginvariable x(n)minimize( x'*a*x-b*x)subject toc * x <= 1x>=0cvx_end运算结果:Calling SDPT3 4.0: 7 variables, 3 equality constraintsFor improved efficiency, SDPT3 is solving the dual problem.nun 1-2—2x1 —■4J *2 subject tonum. 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 0.000 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime0|0.000|0.000|8.0e-001|6.5e+000|3.1e+002| 1.000000e+001 0.000000e+000| 0:0:00| chol 1 1 1|1.000|0.987|4.3e-007|1.5e-001|1.6e+001| 9.043148e+000 -2.714056e-001| 0:0:01| chol 1 1 2|1.000|1.000|2.6e-007|7.6e-003|1.4e+000| 1.234938e+000 -5.011630e-002| 0:0:01| chol 1 1 3|1.000|1.000|2.4e-007|7.6e-004|3.0e-001| 4.166959e-001 1.181563e-001| 0:0:01| chol4|0.892|0.877|6.4e-008|1.6e-004|5.2e-002| 2.773022e-001 2.265122e-001| 0:0:01| chol5|1.000|1.000|1.0e-008|7.6e-006|1.5e-002| 2.579468e-001 2.427203e-001| 0:0:01| chol6|0.905|0.904|3.1e-009|1.4e-006|2.3e-003| 2.511936e-001 2.488619e-001| 0:0:01| chol7|1.000|1.000|6.1e-009|7.7e-008|6.6e-004| 2.503336e-001 2.496718e-001| 0:0:01| chol8|0.903|0.903|1.8e-009|1.5e-008|1.0e-004| 2.500507e-001 2.499497e-001| 0:0:01| chol9|1.000|1.000|4.9e-010|3.5e-010|2.9e-005| 2.500143e-001 2.499857e-001| 0:0:01| chol10|0.904|0.904|5.7e-011|1.3e-010|4.4e-006| 2.500022e-001 2.499978e-001| 0:0:01| chol11|1.000|1.000|5.2e-013|1.1e-011|1.2e-006| 2.500006e-001 2.499994e-001| 0:0:01| chol12|1.000|1.000|5.9e-013|1.0e-012|1.8e-007| 2.500001e-001 2.499999e-001| 0:0:01| chol13|1.000|1.000|1.7e-012|1.0e-012|4.2e-008| 2.500000e-001 2.500000e-001| 0:0:01| chol14|1.000|1.000|2.3e-012|1.0e-012|7.3e-009| 2.500000e-001 2.500000e-001| 0:0:01|stop: max(relative gap, infeasibilities) < 1.49e-008number of iterations = 14primal objective value = 2.50000004e-001dual objective value = 2.49999996e-001gap := trace(XZ) = 7.29e-009relative gap = 4.86e-009actual relative gap = 4.86e-009rel. primal infeas (scaled problem) = 2.33e-012rel. dual " " " = 1.00e-012rel. primal infeas (unscaled problem) = 0.00e+000rel. dual " " " = 0.00e+000norm(X), norm(y), norm(Z) = 3.2e+000, 1.5e+000, 1.9e+000norm(A), norm(b), norm(C) = 3.9e+000, 4.2e+000, 2.6e+000 Total CPU time (secs) = 0.99CPU time per iteration = 0.07termination code = 0DIMACS: 3.3e-012 0.0e+000 1.3e-012 0.0e+000 4.9e-009 4.9e-009Status: SolvedOptimal value (cvx_optval): -32. 程序代码:n=3; a=[-3 -1 -3];b=[2;5;6];C=[2 1 1;1 2 3;2 2 1]; cvx_beginvariable x(n) minimize( a*x)subject toC * x <= bx>=0 cvx_end 运行结果:Calling SDPT3 4.0: 6 variables, 3 equality constraintsnum. of constraints = 3dim. of linear var = 6*******************************************************************SDPT3: Infeasible path-following algorithms*******************************************************************version predcorr gam expon scale_dataNT 1 0.000 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime0|0.000|0.000|1.1e+001|5.1e+000|6.0e+002|-7.000000e+001 0.000000e+000| 0:0:00| chol 1 1 1|0.912|1.000|9.4e-001|4.6e-002|6.5e+001|-5.606627e+000 -2.967567e+001| 0:0:00| chol 1 1 2|1.000|1.000|1.3e-007|4.6e-003|8.5e+000|-2.723981e+000 -1.113509e+001| 0:0:00| chol 1 1 3|1.000|0.961|2.3e-008|6.2e-004|1.8e+000|-4.348354e+000 -6.122853e+000| 0:0:00| chol 1 1 4|0.881|1.000|2.2e-008|4.6e-005|3.7e-001|-5.255152e+000 -5.622375e+000| 0:0:00| chol 1 1 5|0.995|0.962|1.6e-009|6.2e-006|1.5e-002|-5.394782e+000 -5.409213e+000| 0:0:00| chol 1 1 6|0.989|0.989|2.7e-010|5.2e-007|1.7e-004|-5.399940e+000 -5.400100e+000| 0:0:00| chol 1 1 7|0.989|0.989|5.3e-011|5.8e-009|1.8e-006|-5.399999e+000 -5.400001e+000| 0:0:00| chol 1 1 8|1.000|0.994|2.8e-013|4.3e-011|2.7e-008|-5.400000e+000 -5.400000e+000| 0:0:00|stop: max(relative gap, infeasibilities) < 1.49e-008 number of iterations = 8primal objective value = -5.39999999e+000 dual objective value = -5.40000002e+000 gap := trace(XZ) = 2.66e-008 relative gap = 2.26e-009 actual relative gap = 2.21e-009 rel. primal infeas (scaled problem) = 2.77e-013 rel. dual " " " = 4.31e-011 rel. primal infeas (unscaled problem) = 0.00e+000norm(X), norm(y), norm(Z) = 4.3e+000, 1.3e+000, 1.9e+000norm(A), norm(b), norm(C) = 6.7e+000, 9.1e+000, 5.4e+000Total CPU time (secs) = 0.11CPU time per iteration = 0.01termination code = 0DIMACS: 3.6e-013 0.0e+000 5.8e-011 0.0e+000 2.2e-009 2.3e-009Status: Solvedrel. dual = 0.00e+000Optimal value (cvx_optval): -5.4。

MATLAB中的优化算法及其使用方法

MATLAB中的优化算法及其使用方法

MATLAB中的优化算法及其使用方法1. 引言在科学与工程领域,优化问题是一类常见且重要的问题。

它涉及到在给定约束条件下,寻找最优解或使目标函数达到最小或最大值的问题。

在解决这类问题时,MATLAB是一个非常强大且常用的工具,它提供了多种优化算法和函数。

本文将介绍MATLAB中的部分常见优化算法及其使用方法。

2. 优化问题的形式化表示在应用优化算法之前,首先需要将优化问题进行形式化表示。

假设我们要解决一个优化问题,其中有一个目标函数f(x)和一组约束条件h(x) = 0和g(x) ≤ 0。

这里,x是一个n维向量,表示我们要优化的参数。

3. 无约束优化算法无约束优化算法用于解决没有约束条件的优化问题。

MATLAB中提供了多个无约束优化算法,常用的有fminunc和fminsearch。

3.1 fminunc函数fminunc函数是MATLAB中用于寻找无约束优化问题最小值的函数。

它基于梯度下降算法,通过迭代优化来逼近最优解。

使用fminunc函数,我们需要提供目标函数和初始解作为输入参数,并指定其他可选参数,如最大迭代次数和精度要求。

3.2 fminsearch函数fminsearch函数也是用于无约束优化问题的函数,但与fminunc不同的是,它使用了模拟退火算法来搜索最优解。

使用fminsearch函数,我们需要提供目标函数和初始解作为输入参数,并指定其他可选参数,如最大迭代次数和收敛容忍度。

4. 约束优化算法约束优化算法用于解决带有约束条件的优化问题。

MATLAB中提供了多个约束优化算法,常用的有fmincon和ga。

4.1 fmincon函数fmincon函数是MATLAB中用于求解约束优化问题的函数。

它基于拉格朗日乘子法,并使用内点法等技术来求解约束优化问题。

使用fmincon函数,我们需要提供目标函数、约束条件、初始解和约束类型等作为输入参数,并指定其他可选参数,如最大迭代次数和精度要求。

Matlab中的最优化算法详解

Matlab中的最优化算法详解

Matlab中的最优化算法详解最优化算法是数学和计算机科学中的一个重要研究领域,在实际问题求解中有着广泛的应用。

Matlab作为一个强大的数值计算和科学计算软件,提供了丰富的最优化算法工具箱,使得用户能够方便地进行各种最优化问题的求解。

本文将详细介绍Matlab中的最优化算法,包括基本概念、常用方法和应用示例。

1. 最优化问题的定义和基本概念在介绍最优化算法之前,我们首先需要了解最优化问题的定义和基本概念。

最优化问题可以定义为在给定条件下寻找使得目标函数达到最大或最小值的变量取值。

其中,目标函数是要最大化或最小化的函数,同时还需要考虑约束条件。

在Matlab中,最优化问题通常可以表示为以下形式:```minimize f(x)subject to c(x)≤0```其中,f(x)是目标函数,x是变量向量,c(x)是约束函数,≤0表示约束条件。

通过求解这个问题,我们可以得到最优解x*,使得目标函数f(x)取得最小值。

2. 常用的最优化算法Matlab提供了多种最优化算法,包括无约束优化算法和约束优化算法。

下面将介绍一些常用的最优化算法。

2.1 无约束优化算法无约束优化算法用于求解没有约束条件的最优化问题。

其中,最简单和常用的无约束优化算法是梯度下降法。

梯度下降法通过迭代的方式逐渐调整变量的取值,以使得目标函数逐渐减小。

Matlab中的fminunc函数是梯度下降法的实现,它可以非常方便地求解无约束优化问题。

用户只需要提供目标函数和初始变量值,fminunc就会自动进行迭代优化,最终得到最优解。

2.2 约束优化算法约束优化算法用于求解带有约束条件的最优化问题。

其中,最常用的约束优化算法是拉格朗日乘子法。

拉格朗日乘子法通过引入拉格朗日乘子来将约束条件转化为目标函数的一部分,进而将原始问题转化为无约束优化问题。

Matlab中的fmincon函数是拉格朗日乘子法的实现,它可以方便地求解约束优化问题。

用户需要提供目标函数、约束函数和初始变量值,fmincon会自动进行迭代优化,得到满足约束条件的最优解。

拉格朗日插值法matlab程序代码

拉格朗日插值法matlab程序代码

拉格朗日插值法matlab程序代码
使用拉格朗日插值法进行数据拟合是一种常见的数值计算方法。

在matlab中,我们可以使用polyfit函数来实现拉格朗日插值法。

下面是一个简单的matlab程序代码示例:
```matlab
% 定义原始数据
x = [1, 2, 3, 4, 5];
y = [2, 4, 6, 8, 10];
% 定义插值点
xi = 2.5;
% 使用拉格朗日插值法进行拟合
p = polyfit(x, y, length(x)-1);
yi = polyval(p, xi);
% 输出结果
fprintf('插值点 %f 的函数值为 %f\n', xi, yi);
```
在这个示例中,我们首先定义了原始数据x和y,然后定义了插值点xi。

接着,我们使用polyfit函数进行拉格朗日插值法拟合,其中length(x)-1表示使用n-1次多项式进行拟合,n为原始数据的长度。

最后,我们使用polyval函数计算插值点的函数值yi,并输出结果。

需要注意的是,拉格朗日插值法虽然可以很好地拟合数据,但在插值点附近的函数值可能会出现较大误差。

因此,在实际应用中,我们需要根据具体情况选择合适的插值方法。

优化方法MATLAB编程——大连理工大学

优化方法MATLAB编程——大连理工大学

优化方法上机大作业学院:姓名:学号:指导老师:肖现涛第一题源程序如下: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.最速下降法。

matlab程序优化的常用方法

matlab程序优化的常用方法

matlab程序优化的常用方法
1. 代码结构优化:合理的代码结构可以提高代码的可读性和可维护性。

2. 向量化:使用向量运算代替循环可以显著提高代码的效率。

3. 矩阵预分配:在使用循环更新矩阵时,预先为矩阵分配足够的内存空间可以提高效率。

4. 减少内存分配:避免在循环中频繁分配内存空间,可以显著提高效率。

5. 减少函数调用:函数调用的开销较大,应尽量避免不必要的函数调用。

6. 并行计算:利用Matlab 的并行计算能力,可以显著提高代码的效率。

7. JIT 编译:启用Matlab 的JIT 编译功能可以加速代码的执行。

8. 关闭debug 模式:在执行代码时关闭debug 模式可以加速代码的执行。

9. 使用矩阵运算代替逐元素运算:矩阵运算比逐元素运算效率更高。

10. 使用Mex 文件:使用Mex 文件可以显著提高代码的效率。

【设计】最优化方法及其Matlab程序设计

【设计】最优化方法及其Matlab程序设计

【关键字】设计最优化方法及其Matlab程序设计1.最优化方法概述在生活和工作中,人们对于同一个问题往往会提出多个解决方案,并通过各方面的论证,从中提取最佳方案。

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

最优化是每个人,每个单位所希望实现的事情。

对于产品设计者来说,是考虑如何用最少的材料,最大的性能价格比,设计出满足市场需要的产品。

对于企业的管理者来说,则是如何合理、充分使用现有的设备,减少库存,降低能耗,降低成本,以实现企业的最大利润。

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

用最优化方法解决最优化问题的技术称为最优化技术,它包含两个方面的内容:1)建立数学模型。

即用数学语言来描述最优化问题。

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

2)数学求解。

数学模型建好以后,选择合理的最优化算法进行求解。

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

2.最优化方法(算法)浅析最优化方法求解很大程度上依赖于最优化算法的选择。

这里,对最优化算法做一个简单的分类,并对一些比较常用的典型算法进行解析,旨在加深对一些最优化算法的理解。

最优化算法的分类方法很多,根据不同的分类依据可以得到不同的结果,这里根据优化算法对计算机技术的依赖程度,可以将最优化算法进行一个系统分类:线性规划与整数规划;非线性规划;智能优化方法;变分法与动态规划。

2.1 线性规划与整数规划线性规划在工业、农业、商业、交通运输、军事和科研的各个研究领域有广泛应用。

例如,在资源有限的情况下,如何合理使用人力、物力和资金等资源,以获取最大效益;如何组织生产、合理安排工艺流程或调制产品成分等,使所消耗的资源(人力、设备台时、资金、原始材料等)为最少等。

数值分析大作业(牛顿下山法,拉格朗日法,切比雪夫法)及Matlab程序

数值分析大作业(牛顿下山法,拉格朗日法,切比雪夫法)及Matlab程序

课程设计课程名称:数值分析设计题目:学号:姓名:完成时间:2014.11.18题目一: 解线性方程组的直接法 设方程组Ax b =,其中250002511125555111x x x x x x A x x x ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦, 矩阵中10.1(0,1,,5)k x k k =+=,b 由相应的矩阵元素计算,使解向量(1,1,,1)T x =。

(1) A 不变,对b 的元素6b 加一个扰动410-,求解方程组;(2) b 不变,对A 的元素22a 和66a 分别加一个扰动610-,求解方程组; (3) 对上述两种扰动方程组的解做误差分析。

一.数学原理:本计算采用直接法中的列主元高斯消元法,高斯列主元消元法原理如下: 1、设有n 元线性方程组如下:1111n n nn a a a a ⎛⎫ ⎪ ⎪ ⎪⎝⎭1nx x ⎛⎫ ⎪ ⎪ ⎪⎝⎭=1nb b ⎛⎫ ⎪ ⎪ ⎪⎝⎭2、第一步:如果a11!=0, 令l i1= ai1/a11, I= 2,3,……,n用(-li1)乘第一个方程加到第i 个方程上,得同解方程组:a (1)11 a (1)12 . . . a (1)1nx 1 b (1)1 a (1)21 a (1)22 . . . a (1)2n x 2 b (1)2 . . . . . . . = . a (1)n-11 a (1)n-12 . . a (1)n-1n x n-1 b (1)n-1 a (1)n1 a (1)n2 . . . a (1)nn x n b (1)n简记为:A (2) x = b (2) 其中a (2)ij = a (1)ij – l i1 * a (1)1j , I ,j = 2,3,..,nb (2)I = b (1)I – l i1 * b (1)1 , I = 2,3,...,n 第二步:如果a (2)22 != 0,令l i2= a (2)i2/a (2)22, I= 3,……,n依据同样的原理,对矩阵进行化间(省略),依次下去,直到完成!最后,得到上三角方程组:a(1)11 a(1)12. . . a(1)1nx1b(1)10 a(1)22 . . . a(1)2nx2b(1)2. . . . . . . = .0 0 . . a(n-1)n-1n xn-1b(n-1)n-10 0 . . . a(n)nn xnb(n)n简记为:A(n) x = b(n)最后从方程组的最后一个方程进行回代求解为:Xn = b(n) / a(n)nnXi = ( b(k)k- ∑ a(k)kj x j ) / a(k)kk二.解题过程:1.由题中所给条件可求出b。

大连理工优化方法大作业MATLAB编程

大连理工优化方法大作业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源程序

2000年11月9日%第二章第17题,修正单纯形法解习题11(1)clc;clear;A=[2 1 2;3 3 1];[m,n]=size(A);b=[4;3];r=[4 1 1];c=[4 1 1];bs=[2:3];nbs=[1:1]; a1=A(:,1);T=A(:,bs);a2=inv(T)*a1;b=inv(T)*b;A=[a2,eye(m)];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('当前解为')disp(x)breakendrnbs=r(nbs);kk=find(rnbs==min(rnbs));k=kk(1);Anbs=A(:,nbs);yik=B\Anbs(:,k);%yikxb=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%第二章第17题修正单纯型法解第11(2)A=[1 -1 2 -1;2 1 -3 1;1 1 1 1];[m,n]=size(A);b=[2;6;7];r=[2 1 -1 -1];c=[2 1 -1 -1];bs=[1:3];nbs=[4:4];a1=A(:,4);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%这是第三章十一题%f(x)=x^2+exp(-x);%x=[-2,3]%采用0.618法clc;clear;a=-2;b=3;%cout=1;MAX=100;k=0.618;e=0.00001;x1=a+(1-k)*(b-a);x2=a+k*(b-a);while(cout<MAX)if(b-a<e)xx=0.5*(b+a);disp('最优解x=')disp(xx)disp('最优值为f(x)=')disp(exp(-xx)+xx^2)break;elseif(b-a>e)if((x1^2+exp(-x1))>(x2^2+exp(-x2)))a=x1;b=b;x1=x2;x2=a+k*(b-a);cout=cout+1;elseif((x1^2+exp(-x1))<=(x2^2+exp(-x2))) a=a;b=x2;x2=x1;x1=a+(1-k)*(b-a);cout=cout+1;endendend%这是第三章第15题,要求采用不精确一维搜索算法二解问题clc;clear;syms x1 x2 ky=100*(x2-x1^2)^2+(1-x1)^2;v=[x1 x2];%X(i)g=jacobian(y,v);z=subs(y,{x1,x2},{x1+k,x2+k});%X(i)+k*S(i)z1=subs(z,{x1,x2},{-1,1});z2=diff(z1,k);k1=0;k2=100;k3=inf;cout=1;MAX=1000;c1=0.1;c2=0.5;T=0.1;s=[1;1];while(cout<MAX)y1=subs(y,{x1,x2},{-1,1});y2=subs(z1,k,k2);%f(x(i+1))g1=subs(g,{x1,x2},{-1,1});g3=subs(g,{x1,x2},{x1+k,x2+k});g2=subs(g3,{x1,x2,k},{-1,1,k2});if((y1-y2)>=(-c1*k2*(g1*s))&&(g2*s)>=c2*(g1*s))KK=k2;disp('最小步长=')disp(KK)break;elseif((y1-y2)<(-c1*k2*(g1*s))||(g2*s)<c2*(g1*s)) u1=subs(z1,k,k1);u2=subs(z1,k,k2);uu=subs(g,{x1,x2},{-1+k1,1+k1})*s;u=k1+0.5*(k2-k1)/(1+(u1-u2)/((k2-k1)*uu));%u= k1+(k2-k1)/((u2-u1)/((k2-k1)*uu)-1);if(u>k1&&u<k2)k3=k2;k2=u;elseif(u>=k2&&u<k3)k1=k2;k2=u;elseif(u<=k1)u=k1+T*(k3-k1);k2=u;elseif(u>=k3)u=k3-T*(k3-k1);k2=u;endm=[k1 k2 k3 u cout];cout=cout+1;endend%这是第三章第19题,要求采用阻尼牛顿法计算函数的极小值clc;clear;syms x1 x2 km=3;n=2;e=0.000000000000000000000001;MAX=10;cout=1;%y=x1^2+x2^2+x1*x2-3*x1;y=ff(x1,x2);y1=diff(y,x1);y2=diff(y,x2);y11=diff(y1,x1);y12=diff(y1,x2);y21=diff(y2,x1);y22=diff(y2,x2);while(cout<MAX)df1=subs(y1,{x1,x2},{m,n});df2=subs(y2,{x1,x2},{m,n});if ((double(df1))^2+(double(df2))^2>e)A=[y11 y12; y21 y22];B=[df1;df2];C=inv(A);S=-C*B;z=ff(x1+k*S(1,1),x2+k*S(2,1));z1=subs(diff(z,k),{x1,x2},{m,n});kk=solve(z1,k);S2=[m+kk*S(1,1);n+kk*S(2,1)];m=S2(1,1);n=S2(2,1);cout=cout+1;elseif((double(df1))^2+(double(df2))^2<=e) X=eye(2)*S2;disp('最优解,x=')disp(X);disp('最优值,f(x)=')disp(subs(y,{x1,x2},{X(1,1),X(2,1)})) cout=cout+1;break;endend%函数文件,单独成文件function y=ff(a,b)y=a^2+b^2+a*b-3*a;clc;clear;syms x1 x2 ke=0.000001;MAX=300;cout=1;m=2;n=0;X0=[m;n];temp=1;y=myBFGS(x1,x2);H0=eye(2);y1=diff(y,x1);y2=diff(y,x2);mmm=0;G0=subs([y1;y2],{x1,x2},{m,n});while(cout<MAX)df1=subs(y1,{x1,x2},{X0(1,1),X0(2,1)}); df2=subs(y2,{x1,x2},{X0(1,1),X0(2,1)}); m=X0(1,1);n=X0(2,1);if ((double(df1))^2+(double(df2))^2>e)S0=-H0*G0;z=myBFGS(x1+k*S0(1,1),x2+k*S0(2,1)); z1=diff(z,k);z2=subs(z1,{x1,x2},{m,n});k1=double(solve(z2,k));kk=1;for i=1:numel(k1)if(imag(k1(i))==0)if(k1(i)<1&&k1(i)>0)if(k1(i)<kk)kk=k1(i);endendendendX1=[m+kk*S0(1,1);n+kk*S0(2,1)];XX=subs(X1-X0,{x1,x2},{m,n});G1=subs([y1;y2],{x1,x2},{X1(1,1),X1(2,1)}); GG=G1-G0;u=1+GG'*H0*GG/(XX'*GG);H1=H0+(u*XX*XX'-H0*GG*XX'-XX*GG'*H0)/(XX'*GG);G0=G1;X0=X1;H0=H1;cout=cout+1;elseif((double(df1))^2+(double(df2))^2<=e)disp('最优解为,x=')disp(X0);break;endendfunction y=myBFGS(a,b)y=(a-1)^2+5*(b-a^2)^2;%这是26题的函数%第三章第31题,BFGS法解习题27clc;clear;syms x1 x2 ke=0.000001;MAX=300;cout=1;m=1;n=0;X0=[m;n];y=myBFGS(x1,x2);H0=eye(2);y1=diff(y,x1);y2=diff(y,x2);mmm=0;G0=subs([y1;y2],{x1,x2},{m,n});while(cout<MAX)df1=subs(y1,{x1,x2},{X0(1,1),X0(2,1)});df2=subs(y2,{x1,x2},{X0(1,1),X0(2,1)});m=X0(1,1);n=X0(2,1);if ((double(df1))^2+(double(df2))^2>e)S0=-H0*G0;z=myBFGS(x1+k*S0(1,1),x2+k*S0(2,1));z1=diff(z,k);z2=subs(z1,{x1,x2},{m,n});k1=double(solve(z2,k));kk=0;for i=1:numel(k1)if(imag(k1(i))==0)kk=k1(i);endendX1=[m+kk*S0(1,1);n+kk*S0(2,1)];XX=subs(X1-X0,{x1,x2},{m,n});G1=subs([y1;y2],{x1,x2},{X1(1,1),X1(2,1)});GG=G1-G0;u=1+GG'*H0*GG/(XX'*GG);H1=H0+(u*XX*XX'-H0*GG*XX'-XX*GG'*H0)/(XX'*GG);G0=G1;X0=X1;H0=H1;cout=cout+1;elseif((double(df1))^2+(double(df2))^2<=e)disp('最优解为,x=')disp(X0);break;endend%M函数文件,单独文件function y=myBFGS(a,b)y=a+2*b^2+exp(a^2+b^2);%这是27题的函数。

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

大连理工大学庞丽萍最优化方法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编程

大连理工优化方法大作业MATLAB编程

function [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)时,输出为ffunction f=fun(x)f=(x(2)-x(1)^2)^2+(1-x(1))^2;function gf=gfun(x)gf=[-4*x(1)*(x(2)-x(1)^2)+2*(x(1)-1),2*(x(2)-x(1)^2)];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';结果:s=[-1,1];[x,dk,k]=fjqx(x,s)x =-0.0000 1.0000dk =1.1102e-016k =54function f= fun( X )%所求问题目标函数f=X(1)^2-2*X(1)*X(2)+2*X(2)^2+X(3)^2+ X(4)^2- X(2)*X(3)+2*X(1)+3*X(2)-X(3);endfunction 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)];endfunction [ x,val,k ] = frcg( fun,gfun,x0 )%功能:用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;beta=(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(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d) mk=m;break;endm=m+1;endx0=x0+rho^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 =21或者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)k=0;g0=gfun(x);s=-g0;x=x+dk*s;g1=gfun(x);else if(k<3)u=((norm(g1))^2)/(norm(g0)^2);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(1)^2-2*x(1)*x(2)+2*x(2)^2+x(3)^2+x(4)^2-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};End结果:x=[0,0,0,0];>> [x,f,k]=second(x)x =-4.0147 -3.0132 -1.0090 0f = -7.9999k = 1function [f,x,k]=third_1(x)k=0;g=gfun(x);while(norm(g)>=0.001)s=-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)^2+exp(x(1)^2+x(2)^2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)]; 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;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};endend结果: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)^2+exp(x(1)^2+x(2)^2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)]; function ggf=ggfun(x)ggf=[(4*x(1)^2+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+(4*x(2)^2+2)*exp(x(1)^2+x(2)^2)];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;else d=2*d;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)dx=X{2}-X{1};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)^2+exp(x(1)^2+x(2)^2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)]; function ggf=ggfun(x)ggf=[(4*x(1)^2+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+(4*x(2)^2+2)*exp(x(1)^2+x(2)^ 2);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;endif(p>=0)&&(q>=0)dk=d;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.4195 0.0000k=6function 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,x0)function [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) epsilon=1.0e-9; err=1.0e-6;k=0; x=x0; n=length(x); kmax=1.0e3;ne=length(be); ni=length(bi); lamk=zeros(ne+ni,1);index=ones(ni,1);for (i=1:ni)if(Ai(i,:)*x>bi(i)+epsilon), index(i)=0; endendwhile(k<=kmax)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;endendendk=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(fun,hf,gf,dfun,dhf,dgf,x0)%功能: 用乘子法解一般约束问题: 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;ink=0;epsilon=0.00001;x=x0;he=feval(hf,x);gi=feval(gf,x);n=length(x);l=length(he);m=length(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(i)^2;end%更新乘子向量for i=1:mtemp=min(gi(i),lambda(i)/c);btak=btak+temp^2;endbtak=sqrt(btak);if btak>epsilonif k>=2&&btak>theta*btaoldc=eta*c;endfor i=1:lmu(i)=mu(i)-c*he(i);endfor i=1:mlambda(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);s1=s1+he(i)^2;endpsi=psi+0.5*c*s1;s2=0;for i=1:ms3=max(0,lambda(i)-c*gi(i));s2=s2+s3^2-lambda(i)^2;endpsi=psi+s2/(2*c);%不等式约束函数文件g1.mfunction gi=g1(x)gi=10*x(1)-x(1)^2+10*x(2)-x(2)^2-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.mfunction 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)newf=feval(fun,x0+rho^m*dk,varargin{:});oldf=feval(fun,x0,varargin{:});if(newf<oldf+sigma*rho^m*gk'*dk)mk=m;break;endm=m+1;endx=x0+rho^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{:});结果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: 4f=[3,1,1];A=[2,1,1;1,-1,-1];b=[2;-1];lb=[0,0,0];x=linprog(f,A,b,zeros(3),[0,0,0]',lb) 结果:Optimization terminated.x =0.00000.50000.5000。

优化方法MATLAB编程——大连理工大学

优化方法MATLAB编程——大连理工大学

The optimal solution is 0.000000.
The optimal "x" is "ans".
ans =
1.0000 1.0000 1.0000 1.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;
3
s0=-g0;
大连理工大学优化方法上机大作业
end end x=x0+lanmed*s0; x0=x; g0=g(x); s0=-g0; k=k+1; end end zy_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);
4
大连理工大学优化方法上机大作业
>> 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".
大连理工大学优化方法上机大作业

Matlab中的优化算法与应用

Matlab中的优化算法与应用

Matlab中的优化算法与应用引言随着科学技术的发展,优化问题在各个领域中的应用日益广泛。

优化算法作为解决这些问题的重要工具之一,也得到了越来越多的关注。

Matlab作为一种广泛使用的科学计算工具,不仅提供了丰富的优化算法库,也方便了用户在实际问题中的应用。

本文将介绍Matlab中一些常用的优化算法,并讨论它们在实际应用中的一些案例。

一、经典优化算法1.1 无约束优化算法无约束优化问题是指在没有约束条件的情况下,寻找一个最佳的解决方案。

常见的无约束优化算法有梯度下降法、共轭梯度法和拟牛顿法等。

梯度下降法是一种基于梯度信息的搜索算法,通过不断迭代寻找下降最快的方向来接近最优解。

它的优点是简单易实现,但在处理非凸函数或存在高度非线性的问题时可能陷入局部最优。

共轭梯度法是一种针对二次型无约束优化问题的优化算法,它利用了线性方程组的解的特性,迭代计算得到最优解。

相比于梯度下降法,共轭梯度法收敛速度更快。

拟牛顿法是一种通过构造目标函数的近似二次模型来逼近最优解的方法。

它综合了梯度下降法和牛顿法的优点,既考虑了梯度信息,又避免了计算Hessian矩阵的复杂性。

1.2 约束优化算法约束优化问题是在优化过程中需要满足一定的条件。

常见的约束优化算法有无约束优化算法的扩展版,如增广拉格朗日法、逐步二次规划法等。

增广拉格朗日法是一种通过引入拉格朗日乘子来处理约束条件的方法。

它将原始的优化问题转化为一个等价的无约束优化问题,从而简化了求解过程。

逐步二次规划法是一种通过不断迭代求解线性方程组和二次规划子问题的方法。

它在每一步中,通过求解一个二次规划问题来逼近最优解,直到满足约束条件为止。

二、Matlab中的优化工具箱Matlab提供了丰富的优化工具箱,包括优化算法、优化模型建模和求解、优化工具箱应用等。

2.1 优化算法Matlab中的优化工具箱提供了各种经典的优化算法,如梯度下降法、共轭梯度法、拟牛顿法等。

这些算法可以根据具体的问题选取合适的优化算法进行求解。

优化MATLAB程序的技巧

优化MATLAB程序的技巧

优化MATLAB程序的技巧已有 76 次阅读2012-10-19 22:57|个人分类:MATLAB|系统分类:科研笔记|关键词:的学习技巧其实是转载的,转载地址是:/s/blog_62983bd50100fyvn.html但是直接转载没法进行学习笔记的标注,因此现在拷贝到这里进行学习。

原文内容如下:这是薛定宇老师的建议,原文在大观园的首页:因为 MATLAB 语言是一种解释性语言,所以有时 MATLAB 程序的执行速度不是很理想。

这里将依照作者十多年的实际编程经验给出加快 MATLAB程序执行速度的一些建议和体会。

要点1:尽量避免使用循环:循环语句及循环体经常被认为是~MATLAB 编程的瓶颈问题。

改进这样的状况有两种方法:(1) 尽量用向量化的运算来代替循环操作。

我们将通过如下的例子来演示如何将一般的循环结构转换成向量化的语句。

〖例3.19〗考虑下面无穷级数求和问题:如果我们只求出其中前有限项,比如 100,000 项之和 (要精确地求出级数的和,无需求100000 项,几十项往往就能得出满意的精度。

这里主要是为了演示循环运算向量化的优越性。

),则可以采用下面的常规语句进行计算>>tic;s=0;for i=1:100000s=s+(1/2^i+1/3^i);endtocs =1.5000elapsed_time =1.9700如果采用向量化的方法,则可以得出下面结果。

可以看出,采取向量化的方法比常规循环运算效率要高得多。

>>tic;i=1:100000;s=sum(1./2.^i+1./3.^i);toc;s =1.5000elapsed_time =0.3800(2)在必须使用多重循环的情况下,如果两个循环执行的次数不同,则建议在循环的外循环执行循环次数少的,内环执行循环次数多的。

这样也可以显著提高速度。

〖例3.20〗考虑生成一个 5x10000 的 Hilbert 长方矩阵,该矩阵的定义是其第 i 行第 j 列元素为 h(i,j)=1/(i+j-1)。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

上机大作业II
定义目标函数fun
function f=fun(x)
x1=x(1);
x2=x(2);
f=4*x1-x2^2-12;
定义目标函数梯度函数dfun
function f=dfun(x)
x2=x(2);
f=[4;-2*x2];
定义等式约束函数hf
function qua=hf(x)
qua=25-x(1)^2-x(2)^2;
定义等式约束函数梯度函数dhf
function qua=dhf(x)
qua=[-2*x(1);-2*x(2)];
定义不等式约束函数gfun
function inq=gfun(x)
inq=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;
定义不等式约束梯度数dgf
function inq=dgf(x)
inq=[10-2*x(1);10-2*x(2)];
定义增广拉格朗日函数mpsi
function psi=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: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);
定义增广拉格朗日函数梯度函数dmpsi
function dpsi=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: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
定义BFGS法函数函数bfgs
function [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)
gk=feval(dmpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma);
if(norm(gk)<epsilon1)
break;
end
dk=-Bk\gk;
m=0;
mk=0;
while(m<20)
newf=feval(mpsi,x0+rho^m*dk,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma);
oldf=feval(mpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma);
if(newf<oldf+sigma1*rho^m*gk'*dk)
mk=m;
break;
end
m=m+1;
end
x=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);
end
k=k+1;
x0=x;
end
val=feval(mpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma);
定义增广拉格朗日乘子法函数multphr
function 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)
[x,v,ik]=bfgs('mpsi','dmpsi',x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma);
ink=ink+ik;
he=feval(hf,x);
gi=feval(gfun,x);
btak=0.0;
for i=1:l
btak=btak+he(i)^2;
end
for i=1:m
temp=min(gi(i),lambda(i)/sigma);
btak=btak+temp^2;
end
btak=sqrt(btak);
if btak>epsilon
if(k>=2&&btak > theta*btaold)
sigma=eta*sigma;
end
for i=1:l
mu(i)=mu(i)-sigma*he(i);
end
for i=1:m
lambda(i)=max(0.0,lambda(i)-sigma*gi(i));
end
end
k=k+1;
btaold=btak;
x0=x;
end
f=feval(fun,x);
x
f
mu
lambda
k
运行求解
>> x0=[0;0]
x0 =
>> multphr('fun','hf','gfun','dfun','dhf','dgf',x0)
x =
1.00128148956437
4.89871784708758
f =
-31.9923105871169
mu =
1.01559644571312
lambda =
0.754451167977228
k =
4。

相关文档
最新文档