大连理工优化方法大作业MATLAB编程.docx
《大连理工优化方法大作业MATLAB编程.docx》由会员分享,可在线阅读,更多相关《大连理工优化方法大作业MATLAB编程.docx(19页珍藏版)》请在冰豆网上搜索。
![大连理工优化方法大作业MATLAB编程.docx](https://file1.bdocx.com/fileroot1/2023-1/26/d6b230f7-e295-4239-b58b-b1b972cb4e0e/d6b230f7-e295-4239-b58b-b1b972cb4e0e1.gif)
大连理工优化方法大作业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;
end
if(p>=0)&&(q>=0)
dk=d;
x=x+d*s;
flag=1;
end
k=k+1;
if(p>=0)&&(q<0)
a=d;
d=min{2*d,(d+b)/2};
end
end
%定义求函数值的函数fun,当输入为x0=(x1,x2)时,输出为ffunctionf=fun(x)
f=(x
(2)-x
(1)A2)A2+(1-x
(1)F2;
functiongf=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.00001.0000
dk=1.1102e-016
k=54
习题2
取初始=(0.0.0,0)r^'l用兵柜梯皮法求解下面无约東优化问题:
minf(x)=x孑—2x^X2十2x孑+x孑H-爲—X2天3十2xj+3|X2—*3,
其中步长g的选取可利用习題1戎精确一维披索.
注:
通过比习题验证共範梯度法求辉门无二次西数极小点至多需要“次迭代.
functionf=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);
endfunctiong=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,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(kg=feval(gfun,x0);%计算梯度itern=k-(n+1)*floor(k/(n+1));
itern=itern+1;
%计算搜索方向
if(itern==1)
d=-g;
else
beta=(g*g')/(g0*g0');
d=-g+beta*d0;
gd=g'*d;
if(gd>=0.0)
d=-g;
end
end
if(norm(g)break;
end
m=0;mk=0;
while(m<20)
if(feval(fun,xO+rhoAm*d)mk=m;break;
end
m=m+1;
end
x0=x0+rhoAmk*d;
val=feval(fun,x0);
g0=g;d0=d;
k=k+1;
end
x=x0;
val=feval(fun,x0);
end
结果:
>>x0=[0,0,0,0];
>>[x,val,k]=frcg('fun','gfun',x0)
x=
-4.0000-3.0000-1.00000
val=
-8.0000
k=
或者
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);
elseif(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);
end
end
f=fun(x);
end
functionf=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);
functiongf=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)';
functiondk=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;
end
if(p>=0)&&(q>=0)
dk=d;
flag=1;
end
if(p>=0)&&(q<0)
a=d;
d=min{2*d,(d+b)/2};
end
End
结果:
x=[0,0,0,0];
>>[x,f,k]=second(x)
-1.00900
x=-4.0147-3.0132
f=-7.9999
k=1
取初始点3=(0」)二考虑下面无约東优化问题:
minf(x)二冷+2x2+exp(xf+天孑),
其中歩长Qk的选取可別用习题1或精确一维搜索•搜索方向为一HNW
♦取垃=b
•取皿=R2f防)]"
9耳丈啟为BFG5公式亠
通过此习题体会上述三种算法的收敛速度.
function[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);
end
functionf=fun(x)
f=x
(1)+2*x
(2)A2+exp(x
(1)A2+x
(2)A2);
functiongf=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[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)';functiondk=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;
end
if(p>=0)&&(q>=0)
dk=d;
flag=1;
end
if(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.7729
x=-0.41960.0001
k=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);
end
functionf=fun(x)
f=x
(1)+2*x
(2)A2+exp(x(1F2+x
(2)A2);
functiongf=gfun(x)gf=[1+2*x
(1)*exp(x(1F2+x
(2)A2),4*x
(2)+2*x
(2)*(x(1F2+x
(2)A2)];
functionggf=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)';
functiondk=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;
end
if(p>=0)&&(q>=0)
dk=d;
flag=1;
end
if(p>=0)&&(q<0)
a=d;
if2*d>=(d+b)/2
d=(d+b)/2;
elsed=2*d;
end
end
End
结果:
x=[0,1];
[f,x,k]=third_2(x)
f=0.7729
x=-0.41930.0001
k=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};
end
functionf=fun(x)
f=x
(1)+2*x
(2)A2+exp(x(1F2+x
(2)A2);
functiongf=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)];
functionggf=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)';
functiondk=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;
end
if(p>=0)&&(q<0)
a=d;
if2*d>=(d+b)/2
d=(d+b)/2;
elsed=2*d;
end
end
end
结果:
x=[0,1];
[f,x,k]=third_3(x)
f=0.7729
0.0000
x=-0.4195
k=6
习题4
*U用有效集法求解下面勺勺二次规划问题:
mm
Si.
(XI一I)2+(x2一2.5)2
X1-2X2+2>0
-Xi—2>(2+6>0
-Xi+2X2+2>0
xi,x2>0
functioncallqpact
H=[20;02];
c=[-2-5]';
Ae=[];be=[];
Ai=[1-2;-1-2;-12;10;01];
bi=[-2-6-200]';
x0=[00]';
[x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,xO)
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;endwhile(k<=kmax)
Aee=[];
if(ne>0),Aee=Ae;end
for(j=1:
ni)
if(index(j)>0),Aee=[Aee;Ai(j,:
)];end
end
gk=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)));
end
if(y>=0)
exitflag=0;
else
exitflag=1;
for(i=1:
ni)
if(index(i)&(ne+sum(index(1:
i)))==jk)
index(i)=0;break;
end
end
k=k+1;
else
exitflag=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(tm1tm=tm1;ti=i;
end
end
end
alpha=min(alpha,tm);
x=x+alpha*dk;
if(tm<1),index(ti)=1;end
end
if(exitflag==0),break;end
k=k+1;
end
output.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);
else
x=-ginvH*c;
lambda=0;
end
结果
>>callqpact
x=
1.4000
1.7000
lambda=
0.8000
exitflag=
0
output=
fval:
-6.4500
iter:
7
function[x,mu,lambda,output]=multphr(fun,hf,gf,dfun,dhf,dgf,xO)
%功能:
用乘子法解一般约束问题:
minf(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=xO;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%调用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;
fori=1:
l
btak=btak+he(y2;
end
%更新乘子向量
fori=1:
m
temp=min(gi(i),lambda(i)/c);
btak=btak+tempA2;
end
btak=sqrt(btak);
ifbtak>epsilon
ifk>=2&&btak>theta*btaold
c=eta*c;
end
fori=1:
l
mu(i)=mu(i)-c*he(i);
end
lambda(i)=max(0,lambda(i)-c*gi(i));
end
k=k+1;
btaold=btak;
x0=x;
end
end
f=feval(fun,x);
output.fval=f;
output.iter=k;
%增广拉格朗日函数
functionpsi=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;
fori=1:
l
psi=psi-he(i)*mu(i);
s仁s1+he(y2;
psi=psi+0.5*c*s1;
s2=0;
fori=1:
m
s3=max(0,lambda(i)-c*gi(i));
s2=s2+s3A2-lambda(i)A2;
end
psi=psi+s2/(2*c);
%不等式约束函数文件g1.m
functiongi=g1(x)
gi=10*x
(1)-x
(1)A2+10*x
(2)-x
(2)A2-34;
%目标函数的梯度文件df1.m
functiong=df1(x)
g=[4,-2*x
(2)]';
%等式约束(向量)函数的Jacobi矩阵(转置)文件dh1.m
functiondhe=dh1(x)
dhe=[-2*x
(1),-2*x
(2)]'
%不等式约束(向量)函数的Jacobi矩阵(转置)文件dg1.m
functiondgi=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(kgk=feval(gfun,x0,varargin{:
});
if(norm(gk)break;
end
dk=-Bk\gk;
m=0;mk=0;
while(m<20)
newf=feval(fun,x0+rhoAm*dk,varargin{:
});
oldf=feval(fun,x0,varargin{:
});
if(newfbreak;
end
m=m+1;
end
x=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);
end
k=k+1;x0=x;
end
val=feval(fun,x0,varargin{:
});
结果
x=[22]';
[x,mu,lambda,output]=multphr('fun','hf','gf1','df','dh','dg',x0)
x=
1.0013
4.8987
mu=
0.7701
lambda=
0
0
0.9434
output=
fval:
-31.9923
iter:
4
习题6
利用序列二次规划方法求解习题5中的约束优化问题:
min4xi一好一12
s.t.25-x?
—x孑=Q
10xt一召+10旳-xj-34>0X1,X2>0
f=[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)
结果:
Optimizationterminated.
0.0000
0.5000
0.5000