大连理工优化方法大作业MATLAB编程.docx
《大连理工优化方法大作业MATLAB编程.docx》由会员分享，可在线阅读，更多相关《大连理工优化方法大作业MATLAB编程.docx（19页珍藏版）》请在冰豆网上搜索。
大连理工优化方法大作业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
升级会员 