大连理工大学优化作业程序.docx
《大连理工大学优化作业程序.docx》由会员分享，可在线阅读，更多相关《大连理工大学优化作业程序.docx（17页珍藏版）》请在冰点文库上搜索。
大连理工大学优化作业程序
1.1程序（Java）
publicclassWolfe_Powell{
publicstaticdoublegetFx（double[]x）{
doublex1=x[0];doublex2=x[1];
doubleFx=100*（x2-x1*x1）*（x2-x1*x1）+（1-x1）*（1-x1）;
returnFx;
}
publicstaticdouble[]getDeltFx（double[]x）{
doublex1=x[0];doublex2=x[1];
double[]deltFx=newdouble[2];
deltFx[0]=-400*（x2-x1*x1）*x1-2*（1-x1）;
deltFx[1]=200*（x2-x1*x1）;
returndeltFx;
}
publicstaticdoublegetDeltFx_Sk（double[]deltFx,double[]Sk）{
doublea=0;
for（inti=0;ia=a+deltFx[i]*Sk[i];
}
returna;
}
publicstaticdoublegetL（double[]x,double[]s）{
doublex1=x[0];doublex2=x[1];
doublec1=0.1,c2=0.5,a=0,b=1e8,L=1;
doubleFx0,Fx1，deltFx1_Sk,deltFx0_Sk,temp,temp2;
double[]deltFx0,deltFx1;
Fx0=getFx（x）;
deltFx0=getDeltFx（x）;
deltFx0_Sk=getDeltFx_Sk（deltFx0,s）;
temp=c2*getDeltFx_Sk（deltFx0,s）;
for（inti=0;i<1e8;i++）{
temp2=-c1*L*deltFx0_Sk;
x[0]=x1+L*s[0];
x[1]=x2+L*s[1];
Fx1=getFx（x）;
deltFx1=getDeltFx（x）;
deltFx1_Sk=getDeltFx_Sk（deltFx1,s）;
if（（Fx0-Fx1）>=temp2&&deltFx1_Sk>=temp）{
break;
}
elseif（（Fx0-Fx1）b=L;
L=（L+a）/2;
}
elseif（deltFx1_Ska=L;
L=（L+b）/2>=2*L?
（2*L）:
（L+b）/2;
}
}
System.out.println（"L="+L）;
System.out.println（"计算次数"+i）;
returnL;
}
publicstaticvoidmain（String[]args）{
Wolfe_Powelltemp=newWolfe_Powell（）;
double[]X={-1,1};double[]sk={1,1};temp.getL（X,sk）;
}
}
1.2实验结果
步长L=0.00390625x=[-0.9992,1.0324]计算次数8
2.1程序（Java）
publicclassGongE{
publicstaticdoublegetFx（double[]x）{
doublex1=x[0];
doublex2=x[1];
doubleFx=x1*x1-2*x1*x2+2*x2*x2+x3*x3-x2*x3+2*x1+3*x2-x3;
returnFx;
}
publicstaticdouble[]getDeltFx（double[]x）{
doublex1=x[0];
doublex2=x[1];
double[]deltFx=newdouble[x.length];
deltFx[0]=2*x1-2*x2+2;
deltFx[1]=-2*x1+4*x2-x3+3;
deltFx[2]=2*x3-x2-1;
returndeltFx;
}
publicstaticdouble[]getX（double[]x）{
double[]g0,g1;
double[]s0=newdouble[x.length];
double[]s1=newdouble[x.length];
doubleg0_L,g1_L,L,temp;
double[]x0=x;
intk=0;
g0=getDeltFx（x0）;
for（intj=0;js0[j]=-g0[j];
}
for（inti=0;i<2;i++,k++）{
g0=getDeltFx（x0）;
g0_L=getDeltFx_Sk（s0,s0）;
L=getL（x0,s0）;//例题一中的方法取得步长L
for（intj=0;jx0[j]=x0[j]+s0[j]*L;
}
g1=getDeltFx（x0）;
g1_L=getDeltFx_Sk（g1,g1）;
if（Math.sqrt（g1_L）<=1e-2）{
break;
}
else{
temp=g1_L/g0_L;
for（intj=0;js0[j]=-g1[j]+temp*s0[j];
}
}}
returnx0;
}
publicstaticvoidmain（String[]args）{
GongEtemp=newGongE（）;
double[]x={1,1};
double[]result=temp.getX（x）;
for（inti=0;iSystem.out.println（"result["+i+"]="+result[i]）;
}}}
2.2实验结果
最优点x*=[-4,-3,-1]最优解f*=-8
3.1公用程序（Java）
publicstaticdoublegetFx（double[]x）{//取得Fx值
doublex1=x[0];doublex2=x[1];
doubleFx=x1+2*x2*x2+Math.exp（x1*x1+x2*x2）;
returnFx;
}
publicstaticdouble[]getDeltFx（double[]x）{//取得Fx的梯度值
doublex1=x[0];doublex2=x[1];
double[]deltFx=newdouble[2];
deltFx[0]=1+2*x1*Math.exp（x1*x1+x2*x2）;
deltFx[1]=4*x2+2*x2*Math.exp（x1*x1+x2*x2）;
returndeltFx;
}
3.2.1最速下降法程序（Java）
publicclassFastWay{
publicstaticdouble[]getX（double[]x）{
double[]deltF0=newdouble[2];doubleL=0;
for（inti=0;i<1e1;i++）{
deltF0=getDeltFx（x）;
for（intj=0;jdeltF0[j]=-deltF0[j];
}
L=getL（x,deltF0）;//调用习题1的不精确搜索取得步长L
if（Math.sqrt（getDeltFx_Sk（deltF0,deltF0））<=1e-3）{
System.out.println（"最终计算次数"+i）;
System.out.println（"x1="+x[0]+"x2="+x[1]）;
break;
}
x[0]=x[0]+L*deltF0[0];x[1]=x[1]+L*deltF0[1];
}
returnx;
}
publicstaticvoidmain（String[]args）{
FastWaytemp=newFastWay（）;
double[]x0={2,2};temp.getX（x0）;
}
3.2.2最速下降法结果
最优点X*=[-0.41940]最优解f*=0.7729计算次数count=10
3.3.1牛顿法程序（Java）
publicstaticdouble[]getDeltFx（double[]x）{
doublex1=x[0];doublex2=x[1];
double[]one=newdouble[2];
doubleexp=Math.exp（Math.pow（x1,2）+Math.pow（x2,2））;
one[0]=1+2*x1*exp;one[1]=4*x2+2*x2*exp;
double[][]two=newdouble[2][2];
two[0][0]=2*exp*（1+2*Math.pow（x1,2））;
two[1][1]=2*exp*（1+2*Math.pow（x2,2））+4;
double[]deltFx=newdouble[2];
for（inti=0;i<2;i++）{
deltFx[0]=one[0]/two[0][0];
deltFx[1]=one[1]/two[1][1];
}
returndeltFx;
}
publicstaticvoidmain（String[]args）{
double[]x={1,0};
double[]DeltFx=newdouble[2];
for（inti=0;i<1e3;i++）{
DeltFx=getDeltFx（x）;
x[0]=x[0]-DeltFx[0];
x[1]=x[1]-DeltFx[1];
if（Math.sqrt（getDeltFx_Sk（DeltFx,DeltFx））<=1e-4）{
System.out.println（"计算次数为"+i）;
break;
}
}
System.out.println（"x1="+x[0]+"x2="+x[1]+"\n"）;
System.out.println（"Fx="+getFx（x））;
}
3.3.2牛顿法结果
最优点X*=[-0.4194,0]最优解f*=0.7729计算次数count=5
3.4.1BFGS法程序（matlab）
function[x,val,k]=bfgs（fun,gfun,x0）
maxk=1000;sigma=0.4;rho=0.55;epsion=1e-5;k=0;n=length（x0）;
Bk=eye（n）;%Bk=feval（'Hess',x0）;
while（kgk=feval（gfun,x0）;
if（norm（gk）dk=-Bk\gk;
m=0;mk=0;
while（m<20）
newf=feval（fun,x0+rho^m*dk）
oldf=feval（fun,x0）
if（newfmk=m;break;
end
m=m+1;
end
x=x0+rho^mk*dk;
sk=x-x0;
yk=feval（gfun,x）-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）;
3.4.2BFGS法结果
最优点X*=[-0.41940]最优解f*=0.7729计算次数count=4
4.1有效集法（matlab）
4.1.1主程序
function[x,Lagrange,exitflag,output]=TwoProg（H,c,Ae,be,Ai,bi,x0）
n=length（x0）;x=x0;ni=length（bi）;ne=length（be）;Lagrange=zeros（ne+ni,1）;index=ones（ni,1）;
for（i=1:
ni）
if（Ai（i,:
）*x>bi（i）+1e-9）,index（i）=0;end
end
%算法主程序
k=0;
while（k<=1e4）
%求解子问题
Temp=[];
if（ne>0）,Temp=Ae;end
for（j=1:
ni）
if（index（j）>0）,Temp=[Temp;Ai（j,:
）];end
end
gk=H*x+c;
[m1,n1]=size（Temp）;
[dk,Lagrange]=SubPro（H,gk,Temp,zeros（m1,1））;
if（norm（dk）<=1.0e-6）
y=0.0;
if（length（Lagrange）>ne）
[y,jk]=min（Lagrange（ne+1:
length（Lagrange）））;
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
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;
4.1.2目标函数
functionf=fun（x）
x1=x
（1）;x2=x
（2）;f=eval（'x1+2*x2^2+exp（x1^2+x2^2）'）;
4.1.3子问题函数
function[x,Lagrange]=SubPro（H,c,Ae,be）
[m,n]=size（Ae）;
ginvH=pinv（H）;
if（m>0）
rb=Ae*ginvH*c+be;
Lagrange=pinv（Ae*ginvH*Ae'）*rb;
x=ginvH*（Ae'*Lagrange-c）;
else
x=-ginvH*c;
Lagrange=0;
end
4.1.4运行函数
H=[2-2;-24];
c=[-2-6]';
Ae=[];be=[];
Ai=[1-2;-0.5-0.5;10;01];
bi=[-2-100]';
x0=[01]';
[x,lambda,exitflag,output]=qpact（H,c,Ae,be,Ai,bi,x0）
4.2有效集法结果
内部点初始点x0=[00]最优点X*=[0.81.2]最优解f*=-7.2迭代次数=10
边界点初始点x0=[11]最优点X*=[0.81.2]最优解f*=-7.2迭代次数=2
检验点初始点x0=[01]最优点X*=[0.81.2]最优解f*=-7.2迭代次数=7
5.1乘子法程序（matlab）
5.1.1chengZi程序---乘子法主程序
function[x,mu,Lagrange,output]=chengZi（fun,hf,gf,dfun,dhf,dgf,x0）
sigma=2.0;count=0;innerCount=0;
eta=2.0;θ=0.8;%PHR算法中的实参数θ
x=x0;he=feval（hf,x）;gi=feval（gf,x）;
n=length（x）;l=length（he）;m=length（gi）;
%选取乘子向量的初始值
mu=0.1*ones（l,1）;Lagrange=0.1*ones（m,1）;
btak=10;btaold=10;%用来检验终止条件的两个值
while（btak>1e-6&count<1e3）
%调用BFGS算法程序求解无约束子问题
[x,ival,ik]=bfgs（'Lagr','LagrTiDu',x0,fun,hf,gf,dfun,dhf,dgf,mu,Lagrange,sigma）;
innerCount=innerCount+ik;
he=feval（hf,x）;gi=feval（gf,x）;
btak=0.0;
for（i=1:
l）,btak=btak+he（i）^2;end
for（i=1:
m）
temp=min（gi（i）,Lagrange（i）/sigma）;
btak=btak+temp^2;
end
btak=sqrt（btak）;
ifbtak>1e-6
if（count>=2&btak>θ*btaold）
sigma=eta*sigma;
end
%更新乘子向量
for（i=1:
l）,mu（i）=mu（i）-sigma*he（i）;end
for（i=1:
m）
Lagrange（i）=max（0.0,Lagrange（i）-sigma*gi（i））;
end
end
count=count+1;
btaold=btak;
x0=x;
end
f=feval（fun,x）
output.inner_iter=innerCount;
output.iter=count;
output.bta=btak;
output.fval=f;
5.1.2f1程序---目标函数
functionf=f1（x）
f=4*x
（1）-x
（2）^2-12;
5.1.3h1程序---等式约束
functionhe=h1（x）
he=25-x
（1）^2-x
（2）^2;
5.1.4g1程序---不等式约束
functiongi=g1（x）
gi=10*x
（1）-x
（1）^2+10*x
（2）-x
（2）^2-34;
5.1.5df1程序---目标函数的梯度文件
functiong=df1（x）
g=[4,-2.0*x
（2）]';
5.1.6dhe程序---等式约束（向量）函数的Jacobi矩阵（转置）
functiondhe=dh1（x）
dhe=[-2*x
（1）,-2.0*x
（2）]';
5.1.7dgi程序---不等式约束（向量）函数的Jacobi矩阵（转置）
functiondgi=dg1（x）
dgi=[10-2*x
（1）,10-2*x
（2）;0,1;1,0]';
5.1.8LagrTiDu程序---增广拉格朗日函数的梯度程序
functionresult=LagrTiDu（x,fun,hf,gf,dfun,dhf,dgf,mu,Lagrange,sigma）
result=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）
result=result+（sigma*he（i）-mu（i））*dhe（:
i）;
end
for（i=1:
m）
result=result+（sigma*gi（i）-Lagrange（i））*dgi（:
i）;
end
5.1.9Lagr程序---增广拉格朗日函数程序
functionresult=Lagr（x,fun,hf,gf,dfun,dhf,dgf,mu,Lagrange,sigma）
f=feval（fun,x）;he=feval（hf,x）;gi=feval（gf,x）;
l=length（he）;m=length（gi）;
result=f;s1=0.0;
for（i=1:
l）
result=result-he（i）*mu（i）;
s1=s1+he（i）^2;
end
result=result+0.5
升级会员 