梯度投影法-MATLAB程序可执行

function[x,minf]=minRosen(f,A,b,x0,var,eps)%目标函数：f;%约束矩阵：A;%约束右端力量：b;%初始可行点：x0;%自变量向量：var;%精度：eps;%目标函数取最小值时的自变量值：x;%目标函数的最小值：minf;formatlong;ifnargin==5eps=1.0e-6;endsymsl;x00=transpose(x0);n=length(var);sz=size(A);m=sz(1);gf=jacobian(f,var);bConti=1;whilebContik=0;s=0;A1=A;A2=A;b1=b;b2=b;fori=1:mdfun=A(i,:)*x00-b(i);ifabs(dfun)0.000000001k=k+1;A1(k,:)=A(i,:);b1(k,1)=b(i);elses=s+1;A2(s,:)=A(i,:);b2(s,1)=b(i);endendifk0A1=A1(1:k,:);b1=b1(1:k,:);endifs0A2=A2(1:s,:);b2=b2(1:s,:);endwhile1P=eye(n,n);ifk0tM=transpose(A1);P=P-tM*inv(A1*tM)*A1;endgv=Funval(gf,var,x0);gv=transpose(gv);d=-P*gv;ifd==0ifk==0x=x0;bConti=0;break;elsew=inv(A1*tM)*A1*gv;ifw=0x=x0;bConti=0;break;else[u,index]=min(w);sA1=size(A1);ifsA1(1)==1k=0;elsek=sA1(2);A1=[A1(1:(index-1),:);A1((index+1):sA1(2),:)];%去掉A1对应的行endendendelsebreak;endendd1=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);elselm=inf;fori=1:length(dd)ifdd(i)0ifbb(i)/dd(i)lmlm=bb(i)/dd(i);endendendend[xm,minf]=minHJ(tmpf,0,lm,1.0e-14);tol=norm(xm*d);iftolepsx=x0;break;endx0=x0+xm*d1;%disp('x0');x0endminf=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~=s2fori=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));endfv=subs(f,var,varv);elsefv=subs(f,varvec,varval);endfunction[x,minf]=minHJ(f,a,b,eps)formatlong;ifnargin==3eps=1.0e-6;endl=a+0.382*(b-a);u=a+0.618*(b-a);k=1;tol=b-a;whiletoleps&&k100000fl=subs(f,findsym(f),l);fu=subs(f,findsym(f),u);ifflfua=l;l=u;u=a+0.618*(b-a);elseb=u;u=l;l=a+0.382*(b-a);endk=k+1;tol=abs(b-a);endifk==100000disp('找不到最小值！');x=NaN;minf=NaN;return;endx=(a+b)/2;minf=subs(f,findsym(f),x);formatshort;function[minx,maxx]=minJT(f,x0,h0,eps)formatlong;ifnargin==3eps=1.0e-6;endx1=x0;k=0;h=h0;while1x4=x1+h;k=k+1;f4=subs(f,findsym(f),x4);f1=subs(f,findsym(f),x1);iff4f1x2=x1;x1=x4;f2=f1;f1=f4;h=2*h;elseifk==1h=-h;x2=x4;f2=f4;elsex3=x2;x2=x1;x1=x4;break;endendendminx=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.msymsx1x2x3;%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.mfunction[x,minf]=MinRosenGradientProjectionMethod(f,A,b,x0,var,eps)%fistheobjectionfunction;%Aistheconstraintmatrix;约束矩阵%bistheright-hand-sidevectoroftheconstraints;%x0istheinitialfeasiblepoint;初始可行解%varisthevectorofindependentvariable;自变量向量%epsistheprecision;精度%xisthevalueoftheindependentvariablewhentheobjectivefunctionisminimum;自变量的值是当目标函数最小%minfistheminimumvalueoftheobjectivefunction;目标函数的最小值formatlong;ifnargin==5eps=1.0e-6;endsymsl;x00=transpose(x0);n=length(var);sz=size(A);m=sz(1);%misthenumberofrowsofA行数gf=jacobian(f,var);%calculatethejacobianmatrixoftheobjectivefunction计算目标函数的雅可比矩阵bConti=1;whilebContik=0;s=0;A1=A;A2=A;b1=b;b2=b;fori=1:mdfun=A(i,:)*x00-b(i);%separatematrixAandb分离矩阵A和bifabs(dfun)0.0000001%findmatrixsthatsatisfyA1x_k=b1找到满足的矩阵k=k+1;A1(k,:)=A(i,:);b1(k,1)=b(i);else%findmatrixsthatsatisfyA2x_kb2找到满足的矩阵s=s+1;A2(s,:)=A(i,:);b2(s,1)=b(i);endendifk0A1=A1(1:k,:);b1=b1(1:k,:);endifs0A2=A2(1:s,:);b2=b2(1:s,:);endwhile1P=eye(n,n);ifk0tM=transpose(A1);P=P-tM*inv(A1*tM)*A1;%calculateP;endgv=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==0ifk==0x=x0;bConti=0;break;elsew=-inv(A1*tM)*A1*gv;ifw=0x=x0;bConti=0;break;else[u,index]=min(w);%findthenegativecomponentinwsA1=size(A1);ifsA1(1)==1k=0;elsek=sA1(2);A1=[A1(1:(index-1),:);A1((index+1):sA1(1),:)];%deletecorrespondingrowinA1删除对应的行A1endendendelsebreak;endendd1=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找到搜索区