梯度投影法matlab程序可执行

function [x,minf]=minRosen(f,A,b,x0,var,eps) %目标函数：f;%约束矩阵：A;%约束右端力量：b;%初始可行点：x0;%自变量向量：var;%精度：eps;%目标函数取最小值时的自变量值：x;%目标函数的最小值：minf;format long;if nargin == 5eps=;endsyms l;x00=transpose(x0);n=length(var);sz=size(A);m=sz(1);gf=jacobian(f,var);bConti=1;while bContik=0;s=0;A1=A;A2=A;b1=b;b2=b;for i=1:mdfun=A(i,:)*x00-b(i);if abs(dfun)<k=k+1;A1(k,:)=A(i,:);b1(k,1)=b(i);elses=s+1;A2(s,:)=A(i,:);b2(s,1)=b(i);endendif k>0A1=A1(1:k,:);b1=b1(1:k,:);endif s>0A2=A2(1:s,:);b2=b2(1:s,:);endwhile 1P=eye(n,n);if k>0tM=transpose(A1);P=P-tM*inv(A1*tM)*A1;endgv=Funval(gf,var,x0);gv=transpose(gv);d=-P*gv ;if d==0if k==0x=x0;bConti=0;break;elsew=inv(A1*tM)*A1*gv;if w>=0x=x0;bConti=0;break;else[u,index]=min(w);sA1=size(A1);if sA1(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;if dd>=0[tmpI,lm]=minJT(tmpf,0,;elselm=inf;for i=1:length(dd)if dd(i)<0if bb(i)/dd(i)<lmlm=bb(i)/dd(i);endendendend[xm,minf]=minHJ(tmpf,0,lm,;tol=norm(xm*d);if tol<epsx=x0;break;endx0=x0+xm*d1;%disp('x0');x0endminf=Funval(f,var,x)function fv = 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);if s1 ~= s2for i=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) format long;if nargin == 3eps = ;endl = a + *(b-a);u = a + *(b-a);k=1;tol = b-a;while tol>eps && k<100000fl = subs(f , findsym(f), l);fu = subs(f , findsym(f), u);if fl > fua = l;l = u;u = a + *(b - a);elseb = u;u = l;l = a + *(b-a);endk = k+1;tol = abs(b - a);endif k == 100000disp('找不到最小值！');x = NaN;minf = NaN;return;endx = (a+b)/2;minf = subs(f, findsym(f),x);format short;function [minx,maxx] = minJT(f,x0,h0,eps)format long;if nargin == 3eps = ;endx1 = x0;k = 0;h = h0;while 1x4 = x1 + h;k = k+1;f4 = subs(f, findsym(f),x4);f1 = subs(f, findsym(f),x1);if f4 < f1x2 = x1;x1 = x4;f2 = f1;f1 = f4;h = 2*h;elseif k==1h = -h;x2 = x4;f2 = f4;elsex3 = x2;x2 = x1;x1 = x4;break;endendendminx = min(x1,x3);maxx = x1+x3 - minx;format short;% syms x1 x2 x3 ;% f=x1^2+x1*x2+2*x2^2+2*x3^2+2*x2*x3+4*x1+6*x2+12*x3;% [x,mf]=minRosen(f,[1 1 1 ;1 1 -2],[6;-2],[1 1 3],[x1 x2 x3])% syms x1 x2;%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],[1 2],[x1 x2])% syms x1 x2 x3;% f=-x1*x2*x3;% [x,mf]=minRosen(f,[-1,-2,-2;1,2,2],[0;72],[10 10 10],[x1 x2 x3])% syms x1 x2;%f=2*x1^2+2*x2^2-2*x1*x2^3-4*x1^7-6*x2;% [x,mf]=minRosen(f,[1 1;1 5;-1 0;0 -1],[2;5;0;0],[-1 -1],[x1 x2])------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------syms x1 x2 x3;% f=2*x1^2+2*x2^2-2*x1*x2^3-4*x1^7-6*x2;% var=[x1,x2];% valst=[-1,-1];% A=[1 1;1 5;-1 0;0 -1];% b=[2 5 0 0]';% f=x1^3+x2^2-2*x1-4*x2+6;% var=[x1,x2];% valst=[0 0];% A=[2,-1;1,1;-1,0;0,-1];% b=[1 2 0 0]';var=[x1,x2,x3];valst=[10,10,10];f=-x1*x2*x3;A=[-1,-2,-2;1,2,2];b=[0 72]';[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)%f is the objection function;%A is the constraint matrix; 约束矩阵%b is the right-hand-side vector of the constraints;%x0 is the initial feasible point; 初始可行解%var is the vector of independent variable; 自变量向量%eps is the precision; 精度%x is the value of the independent variable when the objective function is minimum; 自变量的值是当目标函数最小%minf is the minimum value of the objective function; 目标函数的最小值format long;if nargin == 5eps=;endsyms l;x00=transpose(x0);n=length(var);sz=size(A);m=sz(1);% m is the number of rows of A 行数gf=jacobian(f,var);%calculate the jacobian matrix of the objective function 计算目标函数的雅可比矩阵bConti=1;while bContik=0;s=0;A1=A;A2=A;b1=b;b2=b;for i=1:mdfun=A(i,:)*x00-b(i); %separate matrix A and b 分离矩阵A和bif abs(dfun)< %find matrixs that satisfy A1 x_k=b1 找到满足的矩阵k=k+1;A1(k,:)=A(i,:);b1(k,1)=b(i);else%find matrixs that satisfy A2 x_k<b2 找到满足的矩阵s=s+1;A2(s,:)=A(i,:);b2(s,1)=b(i);endendif k>0A1=A1(1:k,:);b1=b1(1:k,:);endif s>0A2=A2(1:s,:);b2=b2(1:s,:);endwhile 1P=eye(n,n);if k>0tM=transpose(A1);P=P-tM*inv(A1*tM)*A1; %calculate P;endgv=Funval(gf,var,x0);gv=transpose(gv);d=-P*gv; %calculate the searching direction 计算搜索方向% flg=1;% if(P==zeros(n))% flg =0;% end% if flg==1% d=d/norm(d); %normorlize the searching direction 搜索方向% end% 加入这部分会无止境的循环if d==0if k==0x=x0;bConti=0;break;elsew=-inv(A1*tM)*A1*gv;if w>=0x=x0;bConti=0;break;else[u,index]=min(w);%find the negative component in wsA1=size(A1);if sA1(1)==1k=0;elsek=sA1(2);A1=[A1(1:(index-1),:);A1((index+1):sA1(1),:)]; %delete corresponding row in A1 删除对应的行A1endendendelsebreak;endendd1=transpose(d);y1=x0+l*d1; %new iteration variable 新的迭代变量tmpf=Funval(f,var,y1);bb=b2-A2*x00;dd=A2*d;if dd>=0[tmpI,lm]= ForwardBackMethod(tmpf,0,; %find the searching interval 找到搜索区间elselm=inf; %find lambda_maxfor i=1:length(dd)% if(dd(i)>0)% %if dd(i)<0% %if bb(i)/dd(i)<lmlm=bb(i)/dd(i);endendendendif lm==inflm=1e9;end[xm,minf]=GoldenSectionSearch(tmpf,0,lm,; %guarantee lambda>0 保证λ> 0%find the minimizer by one dimension searching method 找到由一维搜索方法得到目标tol=norm(xm*d);if tol<epsx=x0;break;endx0=x0+xm*d1;disp('x0');x0endminf=Funval(f,var,x);搜索法确定局部最优值function [x,minf]=GoldenSectionSearch(f,a,b,eps)% search method to find minimum value of function f 搜索方法找到函数f的最小值format long;if nargin==3eps=;endl=a+*(b-a);u=a+*(b-a);k=1;tol=b-a;while tol>eps&&k<100000fl=subs(f ,findsym(f),l);fu=subs(f ,findsym(f),u);if fl>fua=l;l=u;u=a+*(b-a);elseb=u;u=l;l=a+*(b-a);endk=k+1;tol=abs(b-a);endif k==100000disp('找不到最小值！');x=NaN;minf=NaN;return;endx=(a+b)/2; %return the minimizer 返回最小值minf=subs(f, findsym(f),x);format short;进退法确定搜索区间function [left,right]=ForwardBackMethod(f,x0,step)if nargin==2step =endif nargin==1x0=0;step =endf0 =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});endleft=min(x0, x2);right=max(x0, x2);elsestep=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});endleft=min(x1,x2);%left end pointright=max(x1,x2);%right end pointendfunction fv=Funval(f,varvec,varval)var=findsym(f); %找出表达式包含的变量t,s f=t^2+s+1varc=findsym(varvec); %找出传递参数的变量[t s]中的 t ss1=length(var); %函数的个数s2=length(varc); %变量的个数m=floor((s1-1)/3+1); %靠近左边的整数varv=zeros(1,m);if s1 ~= s2%if the number of variable is different, deal with it specially 如果变量的数量是不同的,专门处理它for i=0: ((s1-1)/3)k=findstr(varc,var(3*i+1));index=(k-1)/3;varv(i+1) = varval(index+1);endfv=subs(f,var,varv);elsefv=subs(f,varvec,varval); %return the value of the function 如果原来函数变量个数与传递参数中变量个数一致，调用subs函数，计算给定点处函数值end% disp('here')% f% varvec% varval%fv = subs(f,varvec,varval);。