优化方法作业第二版解读.docx
《优化方法作业第二版解读.docx》由会员分享，可在线阅读，更多相关《优化方法作业第二版解读.docx（19页珍藏版）》请在冰点文库上搜索。
优化方法作业第二版解读
优化方法上机大作业
院系：
化工与环境生命学部
姓名：
李翔宇
学号：
31607007
指导教师:
肖现涛
第一题：
1.最速下降法
源程序如下：
functionx_star=ZSXJ（x0,eps）
gk=grad（x0）;
res=norm（gk）;
k=0;
whileres>eps&&k<=10000
dk=-gk;
ak=1;f0=fun（x0）;
f1=fun（x0+ak*dk）;
slope=dot（gk,dk）;
whilef1>f0+0.0001*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]';
>>esp=1e-4;
>>xk=ZSXJ（x0,eps）
--The1-thiter,theresidualis13.372079
--The2-thiter,theresidualis12.079876
--The3-thiter,theresidualis11.054105
……………………………………………
--The9144-thiter,theresidualis0.000105
--The9145-thiter,theresidualis0.000102
--The9146-thiter,theresidualis0.000100
xk=
0.9999
0.9998
MATLAB截屏：
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
MATALB截屏;
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]';
>>esp=1e-4;
>>xk=Bfgs（x0,eps）
--The1-thiter,theresidualis3.271712
--The2-thiter,theresidualis2.381565
--The3-thiter,theresidualis3.448742
…………………………
--The1516-thiter,theresidualis0.000368
--The1517-thiter,theresidualis0.000099
xk=
1.0001
1.0002
MATLAB截屏：
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
MATLAB截屏：
第二题：
解：
目标函数文件f1.m
functionf=f1（x）
f=4*x
（1）-x
（2）^2-12;
等式约束函数文件h1.m
functionhe=h1（x）
he=25-x
（1）^2-x
（2）^2;
不等式约束函数文件g1.m
functiongi=g1（x）
gi=10*x
（1）-x
（1）^2+10*x
（2）-x
（2）^2-34;
目标函数的梯度文件df1.m
functiong=df1（x）
g=[4,2.0*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）]';
然后在Matlab命令窗口输入如下命令:
x0=[0,0]’;
[x,mu,lambda,output]=multphr（'f1','h1','g1','df1','dh1','dg1',x0）;
得到如下输出:
x=4.89871742648821
1.00128197198571
算法编程利用程序调用格式
第三题：
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
NT10.00010
itpstepdsteppinfeasdinfeasgapprim-objdual-objcputime
-------------------------------------------------------------------
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-011|1.2e-006|2.500006e-0012.499994e-001|0:
0:
01|chol22
12|1.000|1.000|5.9e-013|1.0e-012|1.8e-007|2.500001e-0012.499999e-001|0:
0:
01|chol22
13|1.000|1.000|1.7e-012|1.0e-012|4.2e-008|2.500000e-0012.500000e-001|0:
0:
01|chol22
14|1.000|1.000|2.3e-012|1.0e-012|7.3e-009|2.500000e-0012.500000e-001|0:
0:
01|
stop:
max（relativegap,infeasibilities）<1.49e-008
-------------------------------------------------------------------
numberofiterations=14
primalobjectivevalue=2.50000004e-001
dualobjectivevalue=2.49999996e-001
gap:
=trace（XZ）=7.29e-009
relativegap=4.86e-009
actualrelativegap=4.86e-009
rel.primalinfeas（scaledproblem）=2.33e-012
rel.dual"""=1.00e-012
rel.primalinfeas（unscaledproblem）=0.00e+000
rel.dual"""=0.00e+000
norm（X）,norm（y）,norm（Z）=3.2e+000,1.5e+000,1.9e+000
norm（A）,norm（b）,norm（C）=3.9e+000,4.2e+000,2.6e+000
TotalCPUtime（secs）=0.99
CPUtimeperiteration=0.07
terminationcode=0
DIMACS:
3.3e-0120.0e+0001.3e-0120.0e+0004.9e-0094.9e-009
-------------------------------------------------------------------
------------------------------------------------------------
Status:
Solved
Optimalvalue（cvx_optval）:
-3
2.程序代码：
n=3;
a=[-3-1-3];
b=[2;5;6];
C=[211;123;221];
cvx_begin
variablex（n）
minimize（a*x）
subjectto
C*x<=b
x>=0
cvx_end
运行结果：
CallingSDPT34.0:
6variables,3equalityconstraints
------------------------------------------------------------
num.ofconstraints=3
dim.oflinearvar=6
*******************************************************************
SDPT3:
Infeasiblepath-followingalgorithms
*******************************************************************
versionpredcorrgamexponscale_data
NT10.00010
itpstepdsteppinfeasdinfeasgapprim-objdual-objcputime
-------------------------------------------------------------------
0|0.000|0.000|1.1e+001|5.1e+000|6.0e+002|-7.000000e+0010.000000e+000|0:
0:
00|chol11
1|0.912|1.000|9.4e-001|4.6e-002|6.5e+001|-5.606627e+000-2.967567e+001|0:
0:
00|chol11
2|1.000|1.000|1.3e-007|4.6e-003|8.5e+000|-2.723981e+000-1.113509e+001|0:
0:
00|chol11
3|1.000|0.961|2.3e-008|6.2e-004|1.8e+000|-4.348354e+000-6.122853e+000|0:
0:
00|chol11
4|0.881|1.000|2.2e-008|4.6e-005|3.7e-001|-5.255152e+000-5.622375e+000|0:
0:
00|chol11
5|0.995|0.962|1.6e-009|6.2e-006|1.5e-002|-5.394782e+000-5.409213e+000|0:
0:
00|chol11
6|0.989|0.989|2.7e-010|5.2e-007|1.7e-004|-5.399940e+000-5.400100e+000|0:
0:
00|chol11
7|0.989|0.989|5.3e-011|5.8e-009|1.8e-006|-5.399999e+000-5.400001e+000|0:
0:
00|chol11
8|1.000|0.994|2.8e-013|4.3e-011|2.7e-008|-5.400000e+000-5.400000e+000|0:
0:
00|
stop:
max（relativegap,infeasibilities）<1.49e-008
-------------------------------------------------------------------
numberofiterations=8
primalobjectivevalue=-5.39999999e+000
dualobjectivevalue=-5.40000002e+000
gap:
=trace（XZ）=2.66e-008
relativegap=2.26e-009
actualrelativegap=2.21e-009
rel.primalinfeas（scaledproblem）=2.77e-013
rel.dual"""=4.31e-011
rel.primalinfeas（unscaledproblem）=0.00e+000
rel.dual"""=0.00e+000
norm（X）,norm（y）,norm（Z）=4.3e+000,1.3e+000,1.9e+000
norm（A）,norm（b）,n
升级会员 