整理大连理工大学优化作业程序.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);
1)规划实施可能对相关区域、流域、海域生态系统产生的整体影响。
end
(4)预防或者减轻不良环境影响的对策和措施的合理性和有效性;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;
仍以