优化方法.docx
《优化方法.docx》由会员分享，可在线阅读，更多相关《优化方法.docx（31页珍藏版）》请在冰豆网上搜索。
优化方法
优化方法上机大作业
上机大作业Ⅰ：
编写程序求解下述问题
2
+100（x
2
2
初
minxf（x）=（1-x1）
2–x1
）.
始点取x0=0,
精度取ε=1e-4,步长由Armijo
线搜索生成,
方向
分别由下列方法生成:
1
最速下降法
2
牛顿法
3BFGS方法
4
共轭梯度法
1.最速下降法
源程序如下：
functionx_star=steepest（x0,eps）
gk=grad（x0）;
res=norm（gk）;
k=0;
whileres>eps&&k<=1000
dk=-gk;
ak=1;f0=fun（x0）;
f1=fun（x0+ak*dk）;
slope=dot（gk,dk）;
whilef1>f0+0.1*ak*slope
ak=ak/2;
xk=x0+ak*dk;
f1=fun（xk）;
end
k=k+1;
x0=xk;
gk=grad（xk）;
res=norm（gk）;
fprintf（'--The%d-thiter,theresidualis%f\n',k,res）;
end
x_star=xk;
end
functionf=fun（x）
f=（1-x
（1））^2+100*（x
（2）-x
（1）^2）^2;
end
functiong=grad（x）
g=zeros（2,1）;
g
（1）=2*（x
（1）-1）+400*x
（1）*（x
（1）^2-x
（2））;
g
（2）=200*（x
（2）-x
（1）^2）;
end
运行结果：
>>x0=[0,0]';eps=1e-4eps=
1.0000e-004
>>xk=steepest（x0,eps）
--The1-thiter,theresidualis3.271712
--The2-thiter,theresidualis2.504194
--The3-thiter,theresidualis2.073282
⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯
--The998-thiter,theresidualis0.449447
--The999-thiter,theresidualis0.449447
--The1000-thiter,theresidualis0.449447
--The1001-thiter,theresidualis0.449447
xk=
0.6369
0.4038
2.牛顿法
源程序如下：
functionx_star=newton（x0,eps）
gk=grad（x0）;
bk=[grad2（x0）]^（-1）;
res=norm（gk）;
k=0;
whileres>eps&&k<=1000
dk=-bk*gk;
xk=x0+dk;
k=k+1;
x0=xk;
gk=grad（xk）;
bk=[grad2（xk）]^（-1）;
res=norm（gk）;
fprintf（'--The%d-thiter,theresidualis%f\n',k,res）;
end
x_star=xk;
end
functionf=fun（x）
f=（1-x
（1））^2+100*（x
（2）-x
（1）^2）^2;
end
functiong=grad2（x）
g=zeros（2,2）;
g（1,1）=2+400*（3*x
（1）^2-x
（2））;
g（1,2）=-400*x
（1）;
g（2,1）=-400*x
（1）;
g（2,2）=200;
end
functiong=grad（x）
g=zeros（2,1）;
g
（1）=2*（x
（1）-1）+400*x
（1）*（x
（1）^2-x
（2））;
g
（2）=200*（x
（2）-x
（1）^2）;
end
运行结果：
>>x0=[0,0]';eps=1e-4;
>>xk=newton（x0,eps）
--The1-thiter,theresidualis447.213595
--The2-thiter,theresidualis0.000000
xk=
1.0000
1.0000
3.BFGS方法
源程序如下：
functionx_star=bfgs（x0,eps）
g0=grad（x0）;
gk=g0;
res=norm（gk）;
Hk=eye
（2）;
k=0;
whileres>eps&&k<=1000
dk=-Hk*gk;
ak=1;f0=fun（x0）;
f1=fun（x0+ak*dk）;
slope=dot（gk,dk）;
whilef1>f0+0.1*ak*slope
ak=ak/2;
xk=x0+ak*dk;
f1=fun（xk）;
end
k=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-thiter,theresidualis%f\n',k,res）;
end
x_star=xk;
end
functionf=fun（x）
f=（1-x
（1））^2+100*（x
（2）-x
（1）^2）^2;
end
functiong=grad（x）
g=zeros（2,1）;
g
（1）=2*（x
（1）-1）+400*x
（1）*（x
（1）^2-x
（2））;
g
（2）=200*（x
（2）-x
（1）^2）;
end
运行结果：
>>x0=[0,0]';>>eps=1e-4;
>>xk=bfgs（x0,eps）
--The1-thiter,theresidualis3.271712
--The2-thiter,theresidualis2.381565
--The3-thiter,theresidualis3.448742
⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯
--The998-thiter,theresidualis0.004690
--The999-thiter,theresidualis0.008407
--The1000-thiter,theresidualis0.005314
--The1001-thiter,theresidualis0.010880
xk=
0.9955
0.9911
4.共轭梯度法
源程序如下：
functionx_star=conj（x0,eps）
gk=grad（x0）;
res=norm（gk）;
k=0;
dk=-gk;
whileres>eps&&k<=1000
ak=1;f0=fun（x0）;f1=fun（x0+ak*dk）;slope=dot（gk,dk）;
whilef1>f0+0.1*ak*slope
ak=ak/2;
xk=x0+ak*dk;
f1=fun（xk）;
end
d0=dk;g0=gk;
k=k+1;
x0=xk;gk=grad（xk）;f=（norm（gk）/norm（g0））^2;
res=norm（gk）;dk=-gk+f*d0;
fprintf（'--The%d-thiter,theresidualis%f\n',k,res）;
end
x_star=xk;
end
functionf=fun（x）
f=（1-x
（1））^2+100*（x
（2）-x
（1）^2）^2;
end
functiong=grad（x）
g=zeros（2,1）;
g
（1）=400*x
（1）^3-400*x
（1）*x
（2）+2*x
（1）-2;
g
（2）=-200*x
（1）^2+200*x
（2）;
end
运行结果：
>>x0=[0,0]';
>>eps=1e-4;
>>xk=Conj（x0,eps）
--The1-thiter,theresidualis3.271712
--The2-thiter,theresidualis1.380542
--The3-thiter,theresidualis4.527780
--The4-thiter,theresidualis0.850596
⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯
--The73-thiter,theresidualis0.001532
--The74-thiter,theresidualis0.000402
--The75-thiter,theresidualis0.000134
--The76-thiter,theresidualis0.000057
xk=
0.9999
0.9999
上机大作业Ⅱ：
编写程序利用增广拉格朗日方法求解下述问题
2
min4x1-x2-12
22
s.t.25-x1-x2=0
22
10x1–x1+10x2-x2-34≥0
x,x
≥0
1
2
初始点取x
0=0,精度取ε=1e-4.
主程序：
function[x,mu,lamda,output]=main（fun,hf,gf,dfun,dhf,dgf,x0）
maxk=2000;
theta=0.8;eta=2.0;
k=0;ink=0;
ep=1e-4;
sigma=0.4;
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）;lamda=0.1*ones（m,1）;
betak=10;betaold=10;
while（betak>ep&&k
[ik,x,val]=bfgs（'lagrang','dlagrang',x0,fun,hf,gf,dfun,dhf,dgf,mu,lamda,sigma）;
ink=ink+ik;
he=feval（hf,x）;gi=feval（gf,x）;
betak=sqrt（norm（he,2）^2+norm（min（gi,lamda/sigma）,2）^2）;
ifbetak>ep
mu=mu-sigma*he;
lamda=max（0.0,lamda-sigma*gi）;
if（k>=2&&betak>theta*betaold）
sigma=eta*sigma;
end
end
k=k+1;
betaold=betak;
x0=x;
end
f=feval（fun,x）;
output.fval=f;
output.iter=k;
output.inner_iter=ink;
output.beta=betak;
增广拉格朗日函数
functionlag=lagrang（x,fun,hf,gf,dfun,dhf,dgf,mu,lamda,sigma）
f=feval（fun,x）;he=feval（hf,x）;gi=feval（gf,x）;
l=length（he）;m=length（gi）;
lag=f;s1=0.0;
for（i=1:
l）
lag=lag-he（i）*mu（i）;
s1=s1+he（i）^2;
end
lag=lag+0.5*sigma*s1;
s2=0.0;
for（i=1:
m）
s3=max（0.0,lamda（i）-sigma*gi（i））;
s2=s2+s3^2-lamda（i）^2;
end
lag=lag+s2/（2.0*sigma）;
增广拉格朗日梯度函数
functiondlag=dlagrang（x,fun,hf,gf,dfun,dhf,dgf,mu,lamda,sigma）
dlag=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）
dlag=dlag+（sigma*he（i）-mu（i））*dhe（:
i）;
end
for（i=1:
m）
dlag=dlag+（sigma*gi（i）-lamda（i））*dgi（:
i）;
end
线搜索程序基于BFGS方法
function[k,x,val]=bfgs（fun,gfun,x0,varargin）
Max=1000;
ep=1.e-4;
beta=0.55;sigma1=0.4;
n=length（x0）;Bk=eye（n）;
k=0;
while（k
gk=feval（gfun,x0,varargin{:
}）;
if（norm（gk）
dk=-Bk\gk;
m=0;mk=0;
while（m<20）
newf=feval（fun,x0+beta^m*dk,varargin{:
}）;
oldf=feval（fun,x0,varargin{:
}）;
if（newf<=oldf+sigma1*beta^m*gk'*dk）
mk=m;break;
end
m=m+1;
end
x=x0+beta^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）;
end
k=k+1;
x0=x;
end
val=feval（fun,x0,varargin{:
}）;
目标函数文件
functionf=f1（x）
f=4*x
（1）-x
（2）^2-12;
等式约束文件
functionhe=h1（x）
he=25-x
（1）^2-x
（2）^2;
不等式约束
functiongi=g1（x）
gi=zeros（3,1）;
gi
（1）=10*x
（1）-x
（1）^2+10*x
（2）-x
（2）^2-34;
gi
（2）=x
（1）;
gi（3）=x
（2）;
目标函数梯度文件
functiong=df1（x）
g=[4;-2*x
（1）];
等式函数梯度文件
functiondhe=dh1（x）
dhe=[-2*x
（1）;-2*x
（2）];
不等式函数梯度文件
functiondgi=dg1（x）
dgi=[10-2*x
（1）,1,0
；10-2*x
（2）
，0,1];
输入数据
X0=[0;0]
[x,mu,lamda,output]=main（'f1','h1','g1','df1','dh1','dg1',x0）
输出数据
x=
1.0013
4.8987
mu=
0.0158
lamda=
0.5542
0
0
output=
fval:
-31.9924
iter:
5
inner_iter:
33
beta:
8.4937e-005
上机大作业Ⅲ：
1.解：
将目标函数改写为向量形式：
x'*a*x-b*x
程序代码：
n=2;
a=[0.5,0;0,1];
b=[24];
c=[11];
cvx_begin
variablex（n）
minimize（x'*a*x-b*x）
subjectto
c*x<=1
x>=0
cvx_end
运算结果:
CallingSDPT34.0:
7variables,3equalityconstraints
Forimprovedefficiency,SDPT3issolvingthedualproblem.
------------------------------------------------------------
num.ofconstraints=3
dim.ofsocpvar=4,num.ofsocpblk=1
dim.oflinearvar=3
*******************************************************************
SDPT3:
Infeasiblepath-followingalgorithms
*******************************************************************
versionpredcorrgamexponscale_data
NT
10.000
1
0
itpstepdsteppinfeasdinfeasgap
prim-obj
dual-obj
cputime
-------------------------------------------------------------------
0|0.000|0.000|8.0e-001|6.5e+000|3.1e+002|1.000000e+0010.000000e+000|0:
0:
00|chol11
1|1.000|0.987|4.3e-007|1.5e-001|1.6e+001|9.043148e+000-2.714056e-001|0:
0:
01|chol11
2|1.000|1.000|2.6e-007|7.6e-003|1.4e+000|1.234938e+000-5.011630e-002|0:
0:
01|chol11
3|1.000|1.000|2.4e-007|7.6e-004|3.0e-001|4.166959e-0011.181563e-001|0:
0:
01|chol11
4|0.892|0.877|6.4e-008|1.6e-004|5.2e-002|2.773022e-0012.265122e-001|0:
0:
01|chol11
5|1.000|1.000|1.0e-008|7.6e-006|1.5e-002|2.579468e-0012.427203e-001|0:
0:
01|chol11
6|0.905|0.904|3.1e-009|1.4e-006|2.3e-003|2.511936e-0012.488619e-001|0:
0:
01|chol11
7|1.000|1.000|6.1e-009|7.7e-008|6.6e-004|2.503336e-0012.496718e-001|0:
0:
01|chol11
8|0.903|0.903|1.8e-009|1.5e-008|1.0e-004|2.500507e-0012.499497e-001|0:
0:
01|chol11
9|1.000|1.000|4.9e-010|3.5e-010|2.9e-005|2.500143e-0012.499857e-001|0:
0:
01|chol11
10|0.904|0.904|5.7e-011|1.3e-010|4.4e-006|2.500022e-0012.499978e-001|0:
0:
01|chol22
11|1.000|1.000|5.2e-013|1.1e-0
升级会员 