梯度投影法matlab程序可执行.docx
《梯度投影法matlab程序可执行.docx》由会员分享,可在线阅读,更多相关《梯度投影法matlab程序可执行.docx(18页珍藏版)》请在冰豆网上搜索。
![梯度投影法matlab程序可执行.docx](https://file1.bdocx.com/fileroot1/2023-7/22/e939e2d9-4471-4a19-ab9c-4aa5a6bbb832/e939e2d9-4471-4a19-ab9c-4aa5a6bbb8321.gif)
梯度投影法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);%separatematrixAandb分离矩阵A和b
ifabs(dfun)<%findmatrixsthatsatisfyA1x_k=b1找到满足的矩阵
k=k+1;
A1(k,:
)=A(i,:
);
b1(k,1)=b(i);
else%findmatrixsthatsatisfyA2x_ks=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;%calculateP;
end
gv=Funval(gf,var,x0);
gv=transpose(gv);
d=-P*gv;%calculatethesearchingdirection计算搜索方向
%flg=1;
%if(P==zeros(n))
%flg=0;
%end
%ifflg==1
%d=d/norm(d);%normorlizethesearchingdirection搜索方向
%end
%加入这部分会无止境的循环
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);%findthenegativecomponentinw
sA1=size(A1);
ifsA1
(1)==1
k=0;
else
k=sA1
(2);
A1=[A1(1:
(index-1),:
);A1((index+1):
sA1
(1),:
)];%deletecorrespondingrowinA1删除对应的行A1
end
end
end
else
break;
end
end
d1=transpose(d);
y1=x0+l*d1;%newiterationvariable新的迭代变量
tmpf=Funval(f,var,y1);
bb=b2-A2*x00;
dd=A2*d;
ifdd>=0
[tmpI,lm]=ForwardBackMethod(tmpf,0,;%findthesearchinginterval找到搜索区间
else
lm=inf;%findlambda_max
fori=1:
length(dd)
%if(dd(i)>0)
%%
ifdd(i)<0
%%
ifbb(i)/dd(i)lm=bb(i)/dd(i);
end
end
end
end
iflm==inf
lm=1e9;
end
[xm,minf]=GoldenSectionSearch(tmpf,0,lm,;%guaranteelambda>0保证λ>0
%findtheminimizerbyonedimensionsearchingmethod找到由一维搜索方法得到目标
tol=norm(xm*d);
iftolx=x0;
break;
end
x0=x0+xm*d1;
disp('x0');
x0
end
minf=Funval(f,var,x);
搜索法确定局部最优值
function[x,minf]=GoldenSectionSearch(f,a,b,eps)
%searchmethodtofindminimumvalueoffunctionf搜索方法找到函数f的最小值
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;%returntheminimizer返回最小值
minf=subs(f,findsym(f),x);
formatshort;
进退法确定搜索区间
function[left,right]=ForwardBackMethod(f,x0,step)
ifnargin==2
step=
end
ifnargin==1
x0=0;
step=
end
f0=subs(f,findsym(f),{x0});
x1=x0+step;
f1=subs(f,findsym(f),{x1});
if(f1<=f0)
step2=2*step;
x2=x1+step;
f2=subs(f,findsym(f),{x2});
while(f1>f2)
x0=x1;
x1=x2;
f0=f1;
f1=f2;
step=2*step;
x2=x1+step;
f2=subs(f,findsym(f),{x2});
end
left=min(x0,x2);
right=max(x0,x2);
else
step=2*step;
x2=x1-step;
f2=subs(f,findsym(f),{x2});
while(f0>f2)
x1=x0;
x0=x2;
f1=f0;
f0=f2;
step=2*step;
x2=x1-step;
f2=subs(f,findsym(f),{x2});
end
left=min(x1,x2);%leftendpoint
right=max(x1,x2);%rightendpoint
end
functionfv=Funval(f,varvec,varval)
var=findsym(f);%找出表达式包含的变量t,sf=t^2+s+1
varc=findsym(varvec);%找出传递参数的变量[ts]中的ts
s1=length(var);%函数的个数
s2=length(varc);%变量的个数
m=floor((s1-1)/3+1);%靠近左边的整数
varv=zeros(1,m);
ifs1~=s2%ifthenumberofvariableisdifferent,dealwithitspecially如果变量的数量是不同的,专门处理它
fori=0:
((s1-1)/3)
k=findstr(varc,var(3*i+1));
index=(k-1)/3;
varv(i+1)=varval(index+1);
end
fv=subs(f,var,varv);
else
fv=subs(f,varvec,varval);%returnthevalueofthefunction如果原来函数变量个数与传递参数中变量个数一致,调用subs函数,计算给定点处函数值
end
%disp('here')
%f
%varvec
%varval
%fv=subs(f,varvec,varval);