大连理工大学优化方法上机大作业.docx
《大连理工大学优化方法上机大作业.docx》由会员分享，可在线阅读，更多相关《大连理工大学优化方法上机大作业.docx（24页珍藏版）》请在冰豆网上搜索。
大连理工大学优化方法上机大作业
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+*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,theresidualis
--The2-thiter,theresidualis
x1=
法：
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+*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+*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+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=;%下降标准系数
e_al=1e-4;%误差控制范围
max_itera=25;
out_itera=1;%迭代次数
%===========================算法迭代开始=============================
whileout_iterax_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+
升级会员 