约束优化方法与MATLAB实现.docx
《约束优化方法与MATLAB实现.docx》由会员分享，可在线阅读，更多相关《约束优化方法与MATLAB实现.docx（19页珍藏版）》请在冰点文库上搜索。
约束优化方法与MATLAB实现
例5-1MATLAB实现，用M函数文件形式求解：
symsst
f=（s-5）^2+4*（t-6）^2;
g=[s^2+t^2-64;s-t+10;10-s];
X=[6.27.88.2;11129.8];
[x,minf]=minconSimpSearch2（f,g,X,1.3,0.7,1,0.7,[st]）
复合形法minconSimpSearch2函数文件如下：
function[x,minf]=minconSimpSearch2（f,g,X,alpha,sita,gama,beta,var,eps）
%f:
目标函数
%g:
约束函数
%X:
初始复合形
%alpha:
反射系数
%sita:
紧缩系数
%gama:
扩展系数
%beta:
收缩系数
%var:
自变量向量
%eps:
精度
%x:
目标函数取最小值时的自变量
%minf:
目标函数的最小值
ifnargin==8%函数参量个数
eps=1.0e-6;
end
N=size（X）;
n=N
（2）;
FX=zeros（1,n）;
while1
fori=1:
n
FX（i）=subs（f,var,X（:
i））;
end
[XS,IX]=sort（FX）;
Xsorted=X（:
IX）;%按照IX的顺序重新排列X
px=sum（Xsorted（:
1:
（n-1））,2）/（n-1）;%sum（a,2）,a矩阵行相加Xsorted（:
1:
2）保留Xsort的1，2列。
中心点坐标。
Fpx=subs（f,var,px）;%中心点函数值
SumF=0;
fori=1:
n
SumF=SumF+（FX（IX（i））-Fpx）^2;%判断收敛end
SumF=sqrt（SumF/（n-1））;
ifSumF<=eps
x=Xsorted（:
1）;
break;
else
bcon_1=1;
cof_alpha=alpha;
whilebcon_1
x2=px+cof_alpha*（px-Xsorted（:
n））;%算反射点的坐标
gx2=subs（g,var,x2）;%看有没有出界
ifmin（gx2）>=0
bcon_1=0;
else
cof_alpha=0.7*（cof_alpha）;
end
end
fx2=subs（f,var,x2）;%反射点函数值
iffx2（1）
cof_gama=gama;
bcon_2=1;
whilebcon_2
x3=x2+cof_gama*（x2-px）;%扩张步骤，感觉应该用x2代贴第一部分px
gx3=subs（g,var,x3）;
fx3=subs（f,var,x3）;
ifmin（gx3）>=0
bcon_2=0;
iffx3（1）
count=1;
else
count=2;
end
else
bcon_2=0;
count=3;
end
end
ifcount==1
Xsorted（:
n）=x3;
X=Xsorted;
continue
else
Xsorted（:
n）=x2;
X=Xsorted;
continue
end
else
iffx2Xsorted（:
n）=x2;
X=Xsorted;
continue
else
iffx2Xsorted（:
n）=x2;
cof_beta=beta;
bcon_3=1;
whilebcon_3<4
x4=Xsorted（:
n）+cof_beta*（px-Xsorted（:
n））;
gx4=subs（g,var,x4）;
ifmin（gx4）>=0
bcon_3=5;
else
cof_beta=cof_beta/2;
bcon_3=bcon_3+1;
end
end
ifmin（gx4）>=0
fx4=subs（f,var,x4）;
FNnew=subs（f,var,Xsorted（:
n））;
iffx4Xsorted（:
n）=x4;
X=Xsorted;
continue
else
x0=Xsorted（:
1）;
fori=1:
n
Xsorted（:
i）=x0+sita*（Xsorted（:
i）-x0）;
end
end
else
x0=Xsorted（:
1）;
fori=1:
n
Xsorted（:
i）=x0+sita*（Xsorted（:
i）-x0）;
X=Xsorted;
continue
end
end
else
x0=Xsorted（:
1）;
fori=1:
n
Xsorted（:
i）=x0+sita*（Xsorted（:
i）-x0）;
X=Xsorted;
continue
end
end
end
end
end
X=Xsorted;
end
minf=subs（f,var,x）;
M函数文件的运行结果如下：
x=5.2186
6.0635
minf=0.0639
====================================================================
例5-2MATLAB实现，用M函数文件形式求解：
symsts;
f=（t-1）^2+（s-2）^2+1;
[x,mf]=minRosen（f,[-21;-1-1],[-1;-2],[0,0],[ts]）
Rosen梯度法投影法函数文件minRosen如下：
function[x,minf]=minRosen（f,A,b,x0,var,eps）
%目标函数：
f；
%约束矩阵：
A；
%约束右端常数向量：
b;
%初始可行点：
x0;
%自变量向量：
var;
%精度：
eps;
%目标函数最小值时的自变量值：
x;
%目标函数的最小值：
minf;
formatlong;
ifnargin==5
eps=1.0e-6;
end
symsl;
x0=transpose（x0）;
n=length（var）;
sz=size（A）;
m=sz
（1）;
gf=jacobian（f,var）;
bConti=1;
whilebConti
k=0;
s=0;
A1=A;
A2=A;
b1=b;
b2=b;
fori=1:
m
dfun=A（i,:
）*x0-b（i）;
ifabs（dfun）<0.000000001%对约束条件的系数矩阵和向量做分解
k=k+1;
A1（k,:
）=A（i,:
）;%A1代表等式约束系数矩阵
b1（k,1）=b（i）;%b1代表等式约束向量
else
s=s+1;
A2（s,:
）=A（i,:
）;%A2代表不等式约束系数矩阵
b2（s,1）=b（i）;%b2代表不等式约束系数矩阵
end
end
ifk>0
A1=A1（1:
k,:
）;
b1=b1（1:
k,:
）;
end
ifs>0
A2=A2（1:
s,:
）;
b2=b2（1:
s,:
）;
end
while1
P=eye（n,n）;
ifk>0
tM=transpose（A1）;
P=P-tM*inv（A1*tM）*A1;
end
gv=subs（gf,var,x0）;
gv=transpose（gv）;
d=-P*gv;%d搜索方向
ifd==0
ifk==0
x=x0;
bConti=0;
break;
else
w=inv（A1*tM）*A1*gv;
ifw>=0%W分量全为正
x=x0;
bConti=0;
break;
else
[u,index]=min（w）;
sA1=size（A1）;
ifsA1
（1）==1
k=0;
else
k=sA1
（2）;%选择W的一个负分量
A1=[A1（1:
（index-1）,:
）;A1（（index+1）:
sA1
（2）,:
）];%去掉A1对应的行
end
end
end
else
break;
end
end
yl=x0+l*d;
tmpf=subs（f,var,yl）;
bb=b2-A2*x0;
dd=A2*d;
ifdd>=0
[tmpI,lm]=minJT（tmpf,0,0.1）;%用进退法确定一维极值问题的极值区间
else
lm=inf;
fori=1:
length（dd）
ifdd（i）<0
ifbb（i）/dd（i）lm=bb（i）/dd（i）;
end
end
end
end
[xm,minf]=minHJ（tmpf,0,lm,1.0e-14）;%用黄金分割法求解一维极值问题
tol=norm（xm*d）;
iftolx=x0;
break;
end
x0=x0+xm*d;
end
minf=subs（f,var,x）;
M函数文件的运行结果如下：
x=0.5000
1.5000
mf=1.5000
===========================================================
例5-3MATLAB实现，用M函数文件形式求解：
symst;
a=4;b=3;
f=a*t;
g=[t-b];
[x,minf]=minNF（f,[5],g,10,0.5,[t]）
内点惩罚函数法文件minNF如下：
function[x,minf]=minNF（f,x0,g,u,v,var,eps）
%目标函数：
f；
%初始点：
x0；
%约束函数：
g；
%罚因子：
u；
%缩小系数：
v；
%自变量向量：
var；
%精度：
eps；
%目标函数取最小值时的自变量：
x；
%目标函数的最小值：
minf；
formatlong;
ifnargin==6
eps=1.0e-6;
end
k=0;
FE=0;
fori=1:
length（g）
FE=FE+1/g（i）;%构造罚函数
end
x1=transpose（x0）;
x2=inf;
while1
FF=u*FE;
SumF=f+FF;
[x2,minf]=minNT（SumF,transpose（x1）,var）;%用牛顿法求解无约束优化
Bx=subs（FE,var,x2）;
ifu*Bxifnorm（x2-x1）<=eps
x=x2;
break;
else
u=v*u;%修正参数
x1=x2;
end
else
ifnorm（x2-x1）<=eps
x=x2;
break;
else
u=v*u;%修正参数
x1=x2;
end
end
end
minf=subs（f,var,x）;
formatshort;
M函数文件的运行结果如下：
x=3.0000
minf=12.0000
=================================================================
例5-5用混合惩罚函数法求解下列优化问题
取初始点
。
解：
MATLAB实现，用M函数文件形式求解：
symsst;
f=s^2-s*t+1;
g=[s^2+t^2-15;s;t];h=[2*s+3*t-20];
[x,minf]=minMixFun（f,g,h,[33],2,0.5,[st]）
混合惩罚函数法文件minMixFun如下：
function[x,minf]=minMixFun（f,g,h,x0,r0,c,var,eps）
%目标函数：
f；
%不等式约束：
g；
%初始点：
x0；
%罚因子：
r0；
%缩小系数：
c；
%自变量向量：
var；
%精度：
eps；
%目标函数取最小值时的自变量:
x;
%目标函数的最小值：
minf;
gx0=subs（g,var,x0）;
ifgx0>=0;
else
disp（'初始点必须满足不等式约束！
'）;
x=NaN;
minf=NaN;
return;
end
ifr0<=0
disp（'初始障碍因子必须大于0！
'）;
x=NaN;
minf=NaN;
return;
end
ifc>=1||c<0
disp（'缩小系数必须大于0且小于1！
'）;
x=NaN;
minf=NaN;
return;
end
ifnargin==7
eps=1.0e-6;
end
FE=0;
fori=1:
length（g）
FE=FE+1/g（i）;
end
FH=transpose（h）*h;
x1=transpose（x0）;
x2=inf;
while1
FF=r0*FE+FH/sqrt（r0）;
SumF=f+FF;
[x2,minf]=minNT（SumF,transpose（x1）,var）;
ifnorm（x2-x1）<=eps
x=x2;
break;
else
r0=c*r0;
x1=x2;
end
end
minf=subs（f,var,x）;
M函数文件的运行结果如下：
x=2.0000
5.3333
minf=-5.6667
=====================================================================
例5-6用函数fmincon求解约束优化问题
s.t
解：
首先编制两个M文件，分别保存目标函数和约束函数
目标函数M文件
functionf=objfun（x）
f=5*x
（1）^2+2*x
（2）^2;
约束函数M文件
function[c,ceq]=confun（x）
c=1.5/x
（1）-x
（2）;
ceq=[];
调用优化工具箱fmincon解上述问题
[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN]=fmincon（@objfun,[12],[],[],[],[],[],[],@confun）
运行结果如下：
Warning:
Large-scale（trustregion）methoddoesnotcurrentlysolvethistypeofproblem,switchingtomedium-scale（linesearch）.
>Infminconat260
Optimizationterminated:
first-orderoptimalitymeasurelessthanoptions.TolFunandmaximumconstraintviolationislessthanoptions.TolCon.
Activeinequalities（towithinoptions.TolCon=1e-006）:
lowerupperineqlinineqnonlin
1
X=0.97401.5400
FVAL=9.4868
EXITFLAG=1
OUTPUT=
iterations:
7
funcCount:
33
stepsize:
1
algorithm:
'medium-scale:
SQP,Quasi-Newton,line-search'
firstorderopt:
1.8905e-007
cgiterations:
[]
message:
[1x144char]
LAMBDA=
lower:
[2x1double]
upper:
[2x1double]
eqlin:
[1x0double]
eqnonlin:
[1x0double]
ineqlin:
[1x0double]
ineqnonlin:
6.1601
GRAD=9.7400
6.1601
HESSIAN=24.6592-3.5437
-3.54371.6534
具体参数分析，fmincon函数介绍。
====================================================================
例5-7用函数fmincon求解约束优化问题
s.t
解：
首先建立目标函数M文件
functiony=optimFun（x）
y=x
（1）^4-4*x
（1）-8*x
（2）+15;
根据题目约束条件知，约束条件中既有线性约束，又有非线性约束，因此必须分别处理。
根据线性约束条件有
建立非线性约束M文件
function[c,ceq]=Fxxconfun（x）
c=9-x
（1）^2-x
（2）^2;
ceq=[];
调用优化工具箱fmincon函数求解
[x,fval,exitflag,output]=fmincon（@optimFun,[12],A,b,[],[],[],[],@fxxconfun）
运行结果如下：
Warning:
Large-scale（trustregion）methoddoesnotcurrentlysolvethistypeofproblem,switchingtomedium-scale（linesearch）.
>Infminconat260
Optimizationterminated:
Magnitudeofdirectionalderivativeinsearchdirectionlessthan2*options.TolFunandmaximumconstraintviolationislessthanoptions.TolCon.
Noactiveinequalities
x=1.2420-2.7308
fval=34.2582
exitflag=5
output=iterations:
17
funcCount:
74
stepsize:
1
algorithm:
'medium-scale:
SQP,Quasi-Newton,line-search'
firstorderopt:
0.0203
cgiterations:
[]
message:
[1x172char]
升级会员 