梯度投影法matlab程序可执行.docx
《梯度投影法matlab程序可执行.docx》由会员分享，可在线阅读，更多相关《梯度投影法matlab程序可执行.docx（16页珍藏版）》请在冰豆网上搜索。
梯度投影法matlab程序可执行
function[x,minf]=minRosen（f,A,b,x0,var,eps）
%目标函数：
f;
%约束矩阵：
A;
%约束右端力量：
b;
%初始可行点：
x0;
%自变量向量：
var;
%精度：
eps;
%目标函数取最小值时的自变量值：
x;
%目标函数的最小值：
minf;
formatlong;
ifnargin==5
eps=;
end
symsl;
x00=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,:
）*x00-b（i）;
ifabs（dfun）<
k=k+1;
A1（k,:
）=A（i,:
）;
b1（k,1）=b（i）;
else
s=s+1;
A2（s,:
）=A（i,:
）;
b2（s,1）=b（i）;
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=Funval（gf,var,x0）;
gv=transpose（gv）;
d=-P*gv;
ifd==0
ifk==0
x=x0;
bConti=0;
break;
else
w=inv（A1*tM）*A1*gv;
ifw>=0
x=x0;
bConti=0;
break;
else
[u,index]=min（w）;
sA1=size（A1）;
ifsA1
（1）==1
k=0;
else
k=sA1
（2）;
A1=[A1（1:
（index-1）,:
）;A1（（index+1）:
sA1
（2）,:
）];%去掉A1对应的行
end
end
end
else
break;
end
end
d1=transpose（d）;
y1=x0+l*d1;
tmpf=Funval（f,var,y1）;
bb=b2-A2*x00;
dd=A2*d;
ifdd>=0
[tmpI,lm]=minJT（tmpf,0,;
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,;
tol=norm（xm*d）;
iftolx=x0;
break;
end
x0=x0+xm*d1;
%
disp（'x0'）;
x0
end
minf=Funval（f,var,x）
functionfv=Funval（f,varvec,varval）
var=findsym（f）;
varc=findsym（varvec）;
s1=length（var）;
s2=length（varc）;
m=floor（（s1-1）/3+1）;
varv=zeros（1,m）;
ifs1~=s2
fori=0:
（（s1-1）/3）
k=findstr（varc,var（3*i+1））;
index=（k-1）/3;
varv（i+1）=varval（index+1）;
%index（i+1）;
%varv（i+1）=varval（index（i+1））;
end
fv=subs（f,var,varv）;
else
fv=subs（f,varvec,varval）;
end
function[x,minf]=minHJ（f,a,b,eps）
formatlong;
ifnargin==3
eps=;
end
l=a+*（b-a）;
u=a+*（b-a）;
k=1;
tol=b-a;
whiletol>eps&&k<100000
fl=subs（f,findsym（f）,l）;
fu=subs（f,findsym（f）,u）;
iffl>fu
a=l;
l=u;
u=a+*（b-a）;
else
b=u;
u=l;
l=a+*（b-a）;
end
k=k+1;
tol=abs（b-a）;
end
ifk==100000
disp（'找不到最小值！
'）;
x=NaN;
minf=NaN;
return;
end
x=（a+b）/2;
minf=subs（f,findsym（f）,x）;
formatshort;
function[minx,maxx]=minJT（f,x0,h0,eps）
formatlong;
ifnargin==3
eps=;
end
x1=x0;
k=0;
h=h0;
while1
x4=x1+h;
k=k+1;
f4=subs（f,findsym（f）,x4）;
f1=subs（f,findsym（f）,x1）;
iff4x2=x1;
x1=x4;
f2=f1;
f1=f4;
h=2*h;
else
ifk==1
h=-h;
x2=x4;
f2=f4;
else
x3=x2;
x2=x1;
x1=x4;
break;
end
end
end
minx=min（x1,x3）;
maxx=x1+x3-minx;
formatshort;
%symsx1x2x3;
%f=x1^2+x1*x2+2*x2^2+2*x3^2+2*x2*x3+4*x1+6*x2+12*x3;
%[x,mf]=minRosen（f,[111;11-2],[6;-2],[113],[x1x2x3]）
%symsx1x2;
%f=x1^3+x2^2-2*x1-4*x2+6;
%[x,mf]=minRosen（f,[2,-1;1,1;-1,0;0,-1],[1;2;0;0],[12],[x1x2]）
%symsx1x2x3;
%f=-x1*x2*x3;
%[x,mf]=minRosen（f,[-1,-2,-2;1,2,2],[0;72],[101010],[x1x2x3]）
%symsx1x2;
%f=2*x1^2+2*x2^2-2*x1*x2^3-4*x1^7-6*x2;
%[x,mf]=minRosen（f,[11;15;-10;0-1],[2;5;0;0],[-1-1],[x1x2]）
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
symsx1x2x3;
%f=2*x1^2+2*x2^2-2*x1*x2^3-4*x1^7-6*x2;
%var=[x1,x2];
%valst=[-1,-1];
%A=[11;15;-10;0-1];
%b=[2500]';
%f=x1^3+x2^2-2*x1-4*x2+6;
%var=[x1,x2];
%valst=[00];
%A=[2,-1;1,1;-1,0;0,-1];
%b=[1200]';
var=[x1,x2,x3];
valst=[10,10,10];
f=-x1*x2*x3;
A=[-1,-2,-2;1,2,2];
b=[072]';
[x,mimfval]=MinRosenGradientProjectionMethod（f,A,b,valst,var）
[x2,fval]=fmincon（'confun',valst,A,b）
function[x,minf]=MinRosenGradientProjectionMethod（f,A,b,x0,var,eps）
%fistheobjectionfunction;
%Aistheconstraintmatrix;约束矩阵
%bistheright-hand-sidevectoroftheconstraints;
%x0istheinitialfeasiblepoint;初始可行解
%varisthevectorofindependentvariable;自变量向量
%epsistheprecision;精度
%xisthevalueoftheindependentvariablewhentheobjectivefunctionisminimum;自变量的值是当目标函数最小
%minfistheminimumvalueoftheobjectivefunction;目标函数的最小值
formatlong;
ifnargin==5
eps=;
end
symsl;
x00=transpose（x0）;
n=length（var）;
sz=size（A）;
m=sz
（1）;%misthenumberofrowsofA行数
gf=jacobian（f,var）;%calculatethejacobianmatrixoftheobjectivefunction计算目标函数的雅可比矩阵
bConti=1;
whilebConti
k=0;
s=0;
A1=A;
A2=A;
b1=b;
b2=b;
fori=1:
m
dfun=A（i,:
）*x00-b（i）;%separatemat
升级会员 