大连理工大学优化方法上机大作业.docx

大连理工大学优化方法上机大作业.docx

《大连理工大学优化方法上机大作业.docx》由会员分享，可在线阅读，更多相关《大连理工大学优化方法上机大作业.docx（27页珍藏版）》请在冰点文库上搜索。

大连理工大学优化方法上机大作业

2016年大连理工大学优化方法上机大作业

学院：

专业：

班级：

学号：

姓名：

上机大作业1：

1.最速下降法：

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

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/4;

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

>>clear

>>x0=[0,0]';

>>eps=1e-4;

>>x=steepest（x0,eps）

2.牛顿法：

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

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

>>clear

>>x0=[0,0]';

>>eps=1e-4;

>>x1=newton（x0,eps）

--The1-thiter,theresidualis447.213595

--The2-thiter,theresidualis0.000000

x1=

1.0000

1.0000

3.BFGS法：

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

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/4;

xk=x0+ak*dk;

f1=fun（xk）;

end

k=k+1;

fa0=xk-x0;

x0=xk;

go=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

>>clear

>>x0=[0,0]';

>>eps=1e-4;

>>x=bfgs（x0,eps）

4.共轭梯度法：

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

functionx_star=CG（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/4;

xk=x0+ak*dk;

f1=fun（xk）;

end

k=k+1;

x0=xk;

g0=gk;

gk=grad（xk）;

res=norm（gk）;

p=（gk/g0）^2;

dk1=dk;

dk=-gk+p*dk1;

fprintf（'--The%d-thiter,theresidualis%f\n',k,res）;

end

x_star=xk;

end

>>clear

>>x0=[0,0]';

>>eps=1e-4;

>>x=CG（x0,eps）

上机大作业2：

functionf=obj（x）

f=4*x

（1）-x

（2）^2-12;

end

function[h,g]=constrains（x）

h=x

（1）^2+x

（2）^2-25;

g=zeros（3,1）;

g

（1）=-10*x

（1）+x

（1）^2-10*x

（2）+x

（2）^2+34;

g

（2）=-x

（1）;

g（3）=-x

（2）;

end

functionf=alobj（x）%拉格朗日增广函数

%N_equ等式约束个数?

%N_inequ不等式约束个数

N_equ=1;

N_inequ=3;

globalr_alpena;%全局变量

h_equ=0;

h_inequ=0;

[h,g]=constrains（x）;

%等式约束部分?

fori=1:

N_equ

h_equ=h_equ+h（i）*r_al（i）+（pena/2）*h（i）.^2;

end

%不等式约束部分

fori=1:

N_inequ

h_inequ=h_inequ+（0.5/pena）*（max（0,（r_al（i）+pena*g（i）））.^2-r_al（i）.^2）;

end

%拉格朗日增广函数值

f=obj（x）+h_equ+h_inequ;

functionf=compare（x）

globalr_alpenaN_equN_inequ;

N_equ=1;

N_inequ=3;

h_inequ=zeros（3,1）;

[h,g]=constrains（x）;

%等式部分

fori=1:

1

h_equ=abs（h（i））;

end

%不等式部分

fori=1:

3

h_inequ=abs（max（g（i）,-r_al（i+1）/pena））;

end

h1=max（h_inequ）;

f=max（abs（h_equ）,h1）;%sqrt（h_equ+h_inequ）;

function[x,fmin,k]=almain（x_al）

%本程序为拉格朗日乘子算法示例算法%函数输入：

%x_al:

初始迭代点

%r_al：

初始拉格朗日乘子N-equ:

等式约束个数N_inequ:

不等式约束个数?

%函数输出

%X：

最优函数点FVAL:

最优函数值

%============================程序开始================================

globalr_alpena;%参数（全局变量）

pena=10;%惩罚系数

r_al=[1,1,1,1];

c_scale=2;%乘法系数乘数

cta=0.5;%下降标准系数

e_al=1e-4;%误差控制范围

max_itera=25;

out_itera=1;%迭代次数

%===========================算法迭代开始=============================

whileout_itera

x_al0=x_al;

r_al0=r_al;

%判断函数?

compareFlag=compare（x_al0）;

%无约束的拟牛顿法BFGS

[X,fmin]=fminunc（@alobj,x_al0）;

x_al=X;%得到新迭代点

%判断停止条件?

ifcompare（x_al）

disp（'wegettheoptpoint'）;

break

end

%c判断函数下降度?

ifcompare（x_al）

pena=1*pena;%可以根据需要修改惩罚系数变量

else

pena=min（1000,c_scale*pena）;%%乘法系数最大1000

disp（'pena=2*pena'）;

end

%%?

更新拉格朗日乘子

[h,g]=constrains（x_al）;

fori=1:

1

%%等式约束部分

r_al（i）=r_al0（i）+pena*h（i）;

end

fori=1:

3

%%不等式约束部分

r_al（i+1）=max（0,（r_al0（i+1）+pena*g（i）））;

end

out_itera=out_itera+1;

end

%+++++++++++++++++++++++++++迭代结束+++++++++++++++++++++++++++++++++

disp（'theiterationnumber'）;

k=out_itera;

disp（'thevalueofconstrains'）;

compare（x_al）

disp（'theoptpoint'）;

x=x_al;

fmin=obj（X）;

>>clear

>>x_al=[0,0];

>>[x,fmin,k]=almain（x_al）

上机大作业3：

1、

>>clearall

n=3;c=[-3,-1,-3]';A=[2,1,1;1,2,3;2,2,1;-1,0,0;0,-1,0;0,0,-1];b=[2,5,6,0,0,0]';

cvx_begin

variablex（n）

minimize（c'*x）

subjectto

A*x<=b

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+01|5.1e+00|6.0e+02|-7.000000e+010.000000e+00|0:

0:

00|chol11

1|0.912|1.000|9.4e-01|4.6e-02|6.5e+01|-5.606627e+00-2.967567e+01|0:

0:

01|chol11

2|1.000|1.000|1.3e-07|4.6e-03|8.5e+00|-2.723981e+00-1.113509e+01|0:

0:

01|chol11

3|1.000|0.961|2.3e-08|6.2e-04|1.8e+00|-4.348354e+00-6.122853e+00|0:

0:

01|chol11

4|0.881|1.000|2.2e-08|4.6e-05|3.7e-01|-5.255152e+00-5.622375e+00|0:

0:

01|chol11

5|0.995|0.962|1.6e-09|6.2e-06|1.5e-02|-5.394782e+00-5.409213e+00|0:

0:

01|chol11

6|0.989|0.989|2.7e-10|5.2e-07|1.7e-04|-5.399940e+00-5.400100e+00|0:

0:

01|chol11

7|0.989|0.989|5.3e-11|5.8e-09|1.8e-06|-5.399999e+00-5.400001e+00|0:

0:

01|chol11

8|1.000|0.994|2.8e-13|4.3e-11|2.7e-08|-5.400000e+00-5.400000e+00|0:

0:

01|

stop:

max（relativegap,infeasibilities）<1.49e-08

-------------------------------------------------------------------

numberofiterations=8

primalobjectivevalue=-5.39999999e+00

dualobjectivevalue=-5.40000002e+00

gap:

=trace（XZ）=2.66e-08

relativegap=2.26e-09

actualrelativegap=2.21e-09

rel.primalinfeas（scaledproblem）=2.77e-13

rel.dual"""=4.31e-11

rel.primalinfeas（unscaledproblem）=0.00e+00

rel.dual"""=0.00e+00

norm（X）,norm（y）,norm（Z）=4.3e+00,1.3e+00,1.9e+00

norm（A）,norm（b）,norm（C）=6.7e+00,9.1e+00,5.4e+00

TotalCPUtime（secs）=0.71

CPUtimeperiteration=0.09

terminationcode=0

DIMACS:

3.6e-130.0e+005.8e-110.0e+002.2e-092.3e-09

-------------------------------------------------------------------

------------------------------------------------------------

Status:

Solved

Optimalvalue（cvx_optval）:

-5.4

2、

>>clearall

n=2;c=[-2,-4]';G=[0.5,0;0,1];

A=[1,1;-1,0;0,-1];b=[1,0,0]';

cvx_begin

variablex（n）

minimize（x'*G*x+c'*x）

subjectto

A*x<=b

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-01|6.5e+00|3.1e+02|1.000000e+010.000000e+00|0:

0:

00|chol11

1|1.000|0.987|4.3e-07|1.5e-01|1.6e+01|9.043148e+00-2.714056e-01|0:

0:

00|chol11

2|1.000|1.000|2.6e-07|7.6e-03|1.4e+00|1.234938e+00-5.011630e-02|0:

0:

00|chol11

3|1.000|1.000|2.4e-07|7.6e-04|3.0e-01|4.166959e-011.181563e-01|0:

0:

00|chol11

4|0.892|0.877|6.4e-08|1.6e-04|5.2e-02|2.773022e-012.265122e-01|0:

0:

00|chol11

5|1.000|1.000|1.0e-08|7.6e-06|1.5e-02|2.579468e-012.427203e-01|0:

0:

00|chol11

6|0.905|0.904|3.1e-09|1.4e-06|2.3e-03|2.511936e-012.488619e-01|0:

0:

00|chol11

7|1.000|1.000|6.1e-09|7.7e-08|6.6e-04|2.503336e-012.496718e-01|0:

0:

00|chol11

8|0.903|0.903|1.8e-09|1.5e-08|1.0e-04|2.500507e-012.499497e-01|0:

0:

00|chol11

9|1.000|1.000|4.9e-10|3.5e-10|2.9e-05|2.500143e-012.499857e-01|0:

0:

00|chol11

10|0.904|0.904|4.7e-11|1.3e-10|4.4e-06|2.500022e-012.499978e-01|0:

0:

00|chol22

11|1.000|1.000|2.3e-12|9.4e-12|1.2e-06|2.500006e-012.499994e-01|0:

0:

00|chol22

12|1.000|1.000|4.7e-13|1.0e-12|1.8e-07|2.500001e-012.499999e-01|0:

0:

00|chol22

13|1.000|1.000|2.0e-12|1.0e-12|4.2e-08|2.500000e-012.500000e-01|0:

0:

00|chol22

14|1.000|1.000|2.6e-12|1.0e-12|7.3e-09|2.500000e-012.500000e-01|0:

0:

00|

stop:

max（relativegap,infeasibilities）<1.49e-08

-------