梯度投影法MATLAB程序可执行.docx
- 文档编号:7219008
- 上传时间:2023-01-22
- 格式:DOCX
- 页数:16
- 大小:18.46KB
梯度投影法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=1.0e-6;
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)<0.000000001
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,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); iftol x=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=1.0e-6; end l=a+0.382*(b-a); u=a+0.618*(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+0.618*(b-a); else b=u; u=l; l=a+0.382*(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=1.0e-6; 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); iff4 x2=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]) ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ Main.m 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) MinRosenGradientProjectionMethod.m function[x,minf]=MinRosenGradientProjectionMethod(f,A,b,x0,var,eps) %fistheobjectionfunction; %Aistheconstraintmatrix;约束矩阵 %bistheright-hand-sidevectoroftheconstraints; %x0istheinitialfeasiblepoint;初始可行解 %varisthevectorofindependentvariable;自变量向量 %epsistheprecision;精度 %xisthevalueoftheindependentvariablewhentheobjectivefunctionisminimum;自变量的值是当目标函数最小 %minfistheminimumvalueoftheobjectivefunction;目标函数的最小值 formatlong; ifnargin==5 eps=1.0e-6; 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)<0.0000001%findmatrixsthatsatisfyA1x_k=b1找到满足的矩阵 k=k+1; A1(k,: )=A(i,: ); b1(k,1)=b(i); else%findmatrixsthatsatisfyA2x_k 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;%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,0.001);%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,1.0e-14);%guaranteelambda>0保证λ>0 %findtheminimizerbyonedimensionsearchingmethod找到由一维搜索方法得到目标 tol=norm(xm*d); iftol x=x0; break; end x0=x0+xm*d1; disp('x0'); x0 end minf=Funval(f,var,x); GoldenSectionSearch.m0.618搜索法确定局部最优值 function[x,minf]=GoldenSectionSearch(f,a,b,eps) %0.618searchmethodtofindminimumvalueoffunctionf0.618搜索方法找到函数f的最小值 formatlong; ifnargin==3 eps=1.0e-6; end l=a+0.382*(b-a); u=a+0.618*(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+0.618*(b-a); else b=u; u=l; l=a+0.382*(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; ForwardBackMethod.m进退法确定搜索区间 function[left,right]=ForwardBackMethod(f,x0,step) ifnargin==2 step=0.01 end ifnargin==1 x0=0; step=0.01 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 Funval.m 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); THANKS! ! ! 致力为企业和个人提供合同协议,策划案计划书,学习课件等等 打造全网一站式需求 欢迎您的下载,资料仅供参考
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 梯度 投影 MATLAB 程序 可执行
![提示](https://static.bdocx.com/images/bang_tan.gif)