大连理工优化方法大作业MATLAB编程Word文档格式.docx
《大连理工优化方法大作业MATLAB编程Word文档格式.docx》由会员分享，可在线阅读，更多相关《大连理工优化方法大作业MATLAB编程Word文档格式.docx（23页珍藏版）》请在冰点文库上搜索。
functionf=fun（X）
%所求问题目标函数
f=X
（1）^2-2*X
（1）*X
（2）+2*X
（2）^2+X（3）^2+X（4）^2-X
（2）*X（3）+2*X
（1）+3*X
（2）-X（3）;
end
functiong=gfun（X）
%所求问题目标函数梯度
g=[2*X
（1）-2*X
（2）+2,-2*X
（1）+4*X
（2）-X（3）+3,2*X（3）-X
（2）-1,2*X（4）];
function[x,val,k]=frcg（fun,gfun,x0）
%功能：
用FR共轭梯度法求无约束问题最小值
%输入：
x0是初始点，fun和gfun分别是目标函数和梯度
%输出：
x、val分别是最优点和最优值，k是迭代次数
maxk=5000;
%最大迭代次数
rho=0.5;
sigma=0.4;
eps=10e-6;
n=length（x0）;
while（k<
maxk）
g=feval（gfun,x0）;
%计算梯度
itern=k-（n+1）*floor（k/（n+1））;
itern=itern+1;
%计算搜索方向
if（itern==1）
d=-g;
else
beta=（g*g'
）/（g0*g0'
）;
d=-g+beta*d0;
gd=g'
*d;
if（gd>
=0.0）
if（norm（g）<
eps）
break;
m=0;
mk=0;
while（m<
20）
if（feval（fun,x0+rho^m*d）<
feval（fun,x0）+sigma*rho^m*g'
*d）
mk=m;
break;
m=m+1;
x0=x0+rho^mk*d;
val=feval（fun,x0）;
g0=g;
d0=d;
k=k+1;
x=x0;
val=feval（fun,x0）;
>
x0=[0,0,0,0];
[x,val,k]=frcg（'
fun'
'
gfun'
x0）
x=
-4.0000-3.0000-1.00000
val=
-8.0000
k=
21
或者
function[x,f,k]=second（x）
dk=dfun（x）;
g0=gfun（x）;
s=-g0;
x=x+dk*s;
g1=gfun（x）;
while（norm（g1）>
=0.02）
if（k==3）
elseif（k<
3）
u=（（norm（g1））^2）/（norm（g0）^2）;
s=-g1+u*s;
g0=g1;
f=fun（x）;
f=x
（1）^2-2*x
（1）*x
（2）+2*x
（2）^2+x（3）^2+x（4）^2-x
（2）*x（3）+2*x
（1）+3*x
（2）-x（3）;
functiongf=gfun（x）
gf=[2*x
（1）-2*x
（2）+2,-2*x
（1）+4*x
（2）-x（3）+3,2*x（3）-x
（2）-1,2*x（4）];
function[p,q]=con（x,d）
ss=-gfun（x）;
p=fun（x）-fun（x+d*ss）+0.2*d*gfun（x）*（ss）'
q=gfun（x+d*ss）*（ss）'
-0.6*gfun（x）*（ss）'
functiondk=dfun（x）
[p,q]=con（x,d）;
=0）
End
x=[0,0,0,0];
[x,f,k]=second（x）
x=-4.0147-3.0132-1.00900
f=-7.9999
k=1
function[f,x,k]=third_1（x）
g=gfun（x）;
while（norm（g）>
=0.001）
s=-g;
dk=dfun（x,s）;
f=x
（1）+2*x
（2）^2+exp（x
（1）^2+x
（2）^2）;
gf=[1+2*x
（1）*exp（x
（1）^2+x
（2）^2）,4*x
（2）+2*x
（2）*（x
（1）^2+x
（2）^2）];
function[j_1,j_2]=con（x,d,s）
j_1=fun（x）-fun（x+d*s）+0.1*d*gfun（x）*（s）'
j_2=gfun（x+d*s）*（s）'
-0.5*gfun（x）*（s）'
functiondk=dfun（x,s）%获取步长
[p,q]=con（x,d,s）;
[f,x,k]=third_1（x）
f=0.7729
x=-0.41960.0001
k=8
（1）程序：
function[f,x,k]=third_2（x）
H=inv（ggfun（x））;
=0.001）
s=（-H*g'
）'
functionggf=ggfun（x）
ggf=[（4*x
（1）^2+2）*exp（x
（1）^2+x
（2）^2）,4*x
（1）*x
（2）*exp（x
（1）^2+x
（2）^2）;
4*x
（1）*x
（2）*exp（x
（1）^2+x
（2）^2）,4+（4*x
（2）^2+2）*exp（x
（1）^2+x
（2）^2）];
functiondk=dfun（x,s）%步长获取
b=10000;
if2*d>
=（d+b）/2
d=（d+b）/2;
elsed=2*d;
结果：
[f,x,k]=third_2（x）
x=-0.41930.0001
（2）程序：
function[f,x,k]=third_3（x）
X=cell
（2）;
g=cell
（2）;
X{1}=x;
H=eye
（2）;
g{1}=gfun（X{1}）;
s=（-H*g{1}'
dk=dfun（X{1},s）;
X{2}=X{1}+dk*s;
g{2}=gfun（X{2}）;
while（norm（g{2}）>
dx=X{2}-X{1};
dg=g{2}-g{1};
v=dx/（dx*dg'
）-（H*dg'
/（dg*H*dg'
h1=H*dg'
*dg*H/（dg*H*dg'
h2=dx'
*dx/（dx*dx'
h3=dg*H*dg'
*v'
*v;
H=H-h1+h2+h3;
X{1}=X{2};
norm（g{2}）;
x=X{2};
4*x
（1）*x
（2）*exp（x
（1）^2+x
（2）^2）,4+（4*x
（2）^2+2）*exp（x
（1）^2+x
（2）^2）;
function[p,q]=con（x,d,s）
p=fun（x）-fun（x+d*s）+0.1*d*gfun（x）*（s）'
q=gfun（x+d*s）*（s）'
functiondk=dfun（x,s）
b=d;
[f,x,k]=third_3（x）
x=-0.41950.0000
k=6
functioncallqpact
H=[20;
02];
c=[-2-5]'
Ae=[];
be=[];
Ai=[1-2;
-1-2;
-12;
10;
01];
bi=[-2-6-200]'
x0=[00]'
[x,lambda,exitflag,output]=qpact（H,c,Ae,be,Ai,bi,x0）
function[x,lamk,exitflag,output]=qpact（H,c,Ae,be,Ai,bi,x0）
epsilon=1.0e-9;
err=1.0e-6;
x=x0;
n=length（x）;
kmax=1.0e3;
ne=length（be）;
ni=length（bi）;
lamk=zeros（ne+ni,1）;
index=ones（ni,1）;
for（i=1:
ni）
if（Ai（i,:
）*x>
bi（i）+epsilon）,index（i）=0;
=kmax）
Aee=[];
if（ne>
0）,Aee=Ae;
for（j=1:
if（index（j）>
0）,Aee=[Aee;
Ai（j,:
）];
gk=H*x+c;
[m1,n1]=size（Aee）;
[dk,lamk]=qsubp（H,gk,Aee,zeros（m1,1））;
if（norm（dk）<
=err）
y=0.0;
if（length（lamk）>
ne）
[y,jk]=min（lamk（ne+1:
length（lamk）））;
if（y>
exitflag=0;
else
exitflag=1;
for（i=1:
if（index（i）&
（ne+sum（index（1:
i）））==jk）
index（i）=0;
exitflag=1;
alpha=1.0;
tm=1.0;
if（（index（i）==0）&
（Ai（i,:
）*dk<
0））
tm1=（bi（i）-Ai（i,:
）*x）/（Ai（i,:
）*dk）;
if（tm1<
tm）
tm=tm1;
ti=i;
alpha=min（alpha,tm）;
x=x+alpha*dk;
if（tm<
1）,index（ti）=1;
if（exitflag==0）,break;
output.fval=0.5*x'
*H*x+c'
*x;
output.iter=k;
function[x,lambda]=qsubp（H,c,Ae,be）
ginvH=pinv（H）;
[m,n]=size（Ae）;
if（m>
rb=Ae*ginvH*c+be;
lambda=pinv（Ae*ginvH*Ae'
）*rb;
x=ginvH*（Ae'
*lambda-c）;
x=-ginvH*c;
lambda=0;
结果
callqpact
1.4000
1.7000
lambda=
0.8000
exitflag=
0
output=
fval:
-6.4500
iter:
7
function[x,mu,lambda,output]=multphr（fun,hf,gf,dfun,dhf,dgf,x0）
%功能:
用乘子法解一般约束问题:
minf（x）,s.t.h（x）=0,g（x）.=0
%输入:
x0是初始点,fun,dfun分别是目标函数及其梯度；
%hf,dhf分别是等式约束（向量）函数及其Jacobi矩阵的转置；
%gf,dgf分别是不等式约束（向量）函数及其Jacobi矩阵的转置；
%输出:
x是近似最优点，mu,lambda分别是相应于等式约束和不等式约束的乘子向量;
%output是结构变量,输出近似极小值f,迭代次数,内迭代次数等
maxk=500;
c=2.0;
eta=2.0;
theta=0.8;
ink=0;
epsilon=0.00001;
he=feval（hf,x）;
gi=feval（gf,x）;
n=length（x）;
l=length（he）;
m=length（gi）;
mu=zeros（l,1）;
lambda=zeros（m,1）;
btak=10;
btaold=10;
while（btak>
epsilon&
k<
%调用BFGS算法程序求解无约束子问题[x,ival,ik]=bfgs（'
mpsi'
dmpsi'
x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c）;
ink=ink+ik;
he=feval（hf,x）;
btak=0;
fori=1:
l
btak=btak+he（i）^2;
%更新乘子向量
m
temp=min（gi（i）,lambda（i）/c）;
btak=btak+temp^2;
btak=sqrt（btak）;
ifbtak>
epsilon
ifk>
=2&
btak>
theta*btaold
c=eta*c;
mu（i）=mu（i）-c*he（i）;
lambda（i）=max（0,lambda（i）-c*gi（i））;
btaold=btak;
x0=x;
f=feval（fun,x）;
output.fval=f;
%增广拉格朗日函数
functionpsi=mpsi（x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c）
psi=f;
s1=0;
fori=1:
psi=psi-he（i）*mu（i）;
s1=s1+he（i）^2;
psi=psi+0.5*c*s1;
s2=0;
s3=max（0,lambda（i）-c*gi（i））;
s2=s2+s3^2-lambda（i）^2;
psi=psi+s2/（2*c）;
%不等式约束函数文件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*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）]'
function[x,val,k]=bfgs（fun,gfun,x0,varargin）
rho=0.55;
Bk=eye（n）;
gk=feval（gfun,x0,varargin{:
}）;
if（norm（gk）<
epsilon）
dk=-Bk\gk;
newf=feval（fun,x0+rho^m*dk,varargin{:
oldf=feval（fun,x0,varargin{:
if（newf<
oldf+sigma*rho^m*gk'
*dk）
end
x=x0+rho^mk*dk;
sk=x-x0;
yk=feval（gfun,x,varargin{:
}）-gk;
if（yk'
*sk>
Bk=Bk-（Bk*sk*sk'
*Bk）/（sk'
*Bk*sk）+（yk*yk'
）/（yk'
*sk）;
x0=x;
val=feval（fun,x0,varargin{:
x=[22]'
[x,mu,lambda,output]=multphr（'
hf'
gf1'
df'
dh'
dg'
x0）
1.0013
4.8987
mu=
0.7701
0.9434
-31.9923
4
f=[3,1,1];
A=[2,1,1;
1,-1,-1];
b=[2;
-1];
lb=[0,0,0];
x=linprog（f,A,b,zeros（3）,[0,0,0]'
lb）
Optimizationterminated.
0.0000
0.5000
升级会员 