梯度投影法MATLAB程序可执行

1、function x,minf=minRosen(f,A,b,x0,var,eps)%目标函数：f;%约束矩阵：A;%约束右端力量：b;%初始可行点：x0;%自变量向量：var;%精度：eps;%目标函数取最小值时的自变量值：x;%目标函数的最小值：minf;format long;if nargin = 5 eps=1.0e-6;endsyms l;x00=transpose(x0);n=length(var);sz=size(A);m=sz(1);gf=jacobian(f,var);bConti=1;while bConti k=0; s=0; A1=A; A2=A; b1=b; b2=

2、b;for i=1:m dfun=A(i,:)*x00-b(i); if abs(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); endendif k>0 A1=A1(1:k,:); b1=b1(1:k,:);endif s>0 A2=A2(1:s,:); b2=b2(1:s,:);endwhile 1 P=eye(n,n); if k>0 tM=transpose(A1); P=P-tM*inv(A1*tM)*A1; end

3、 gv=Funval(gf,var,x0); gv=transpose(gv); d=-P*gv ; if d=0 if k=0 x=x0; bConti=0; break; else w=inv(A1*tM)*A1*gv; if w>=0 x=x0; bConti=0; break; else u,index=min(w); sA1=size(A1); if sA1(1)=1 k=0; else k=sA1(2); A1=A1(1:(index-1),:);A1(index+1):sA1(2),:); %去掉A1对应的行 end end end else break; endendd1

4、=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,0.1); else lm=inf; for i=1:length(dd) if dd(i)<0 if bb(i)/dd(i)<lm lm=bb(i)/dd(i); end end end end xm,minf=minHJ(tmpf,0,lm,1.0e-14); tol=norm(xm*d); if tol<eps x=x0; break; end x0=x0+xm

5、*d1;% disp('x0'); x0 endminf=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 = s2 for i=0: (s1-1)/3) k = findstr(varc,var(3*i+1); index = (k-1)/3; varv(i+1) = var

6、val(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)format long;if nargin = 3 eps = 1.0e-6;end l = a + 0.382*(b-a);u = a + 0.618*(b-a);k=1;tol = b-a; while tol>eps && k<100000 fl = su

7、bs(f , findsym(f), l); fu = subs(f , findsym(f), u); if fl > 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);endif k = 100000 disp('找不到最小值！'); x = NaN; minf = NaN; return;endx = (a+b)/2;minf = subs(f, findsym(f),x);format short

8、;function minx,maxx = minJT(f,x0,h0,eps)format long;if nargin = 3 eps = 1.0e-6;end x1 = x0;k = 0;h = h0;while 1 x4 = x1 + h; k = k+1; f4 = subs(f, findsym(f),x4); f1 = subs(f, findsym(f),x1); if f4 < f1 x2 = x1; x1 = x4; f2 = f1; f1 = f4; h = 2*h; else if k=1 h = -h; x2 = x4; f2 = f4; else x3 = x

9、2; x2 = x1; x1 = x4; break; end endendminx = min(x1,x3);maxx = x1+x3 - minx;format short;% syms x1 x2 x3 ;% f=x12+x1*x2+2*x22+2*x32+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=x13+x22-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)

10、% 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*x12+2*x22-2*x1*x23-4*x17-6*x2; % x,mf=minRosen(f,1 1;1 5;-1 0;0 -1,2;5;0;0,-1 -1,x1 x2)-Main.msyms x1 x2 x3; % f=2*x12+2*x22-2*x1*x23-4*x17-6*x2; % var=x1,x2; % valst=-1,-1; % A=1 1;1 5;-1 0;0 -1

11、; % b=2 5 0 0' % f=x13+x22-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) MinRosenGrad

12、ientProjectionMethod.mfunction 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

13、 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 = 5 eps=1.0e-6; end syms l; x00=transpose(x0); n=length(var); sz=size(A); m=sz(1);% m is the nu

14、mber of rows of A 行数gf=jacobian(f,var);%calculate the jacobian matrix of the objective function 计算目标函数的雅可比矩阵bConti=1; while bConti k=0; s=0; A1=A; A2=A; b1=b; b2=b; for i=1:m dfun=A(i,:)*x00-b(i); %separate matrix A and b 分离矩阵A和b if abs(dfun)<0.0000001 %find matrixs that satisfy A1 x_k=b1 找到满足的矩阵

15、 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); end end if k>0 A1=A1(1:k,:); b1=b1(1:k,:); end if s>0 A2=A2(1:s,:); b2=b2(1:s,:); end while 1 P=eye(n,n); if k>0 tM=transpose(A1); P=P-tM*inv(A1*tM)*A1; %calculate

16、P; end gv=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=0 if k=0 x=x0; bConti=0; break; else w=-inv(A1*tM)*A1*gv; if w>=0 x

17、=x0; bConti=0; break; else u,index=min(w);%find the negative component in w sA1=size(A1); if sA1(1)=1 k=0; else k=sA1(2); A1=A1(1:(index-1),:);A1(index+1):sA1(1),:); %delete corresponding row in A1 删除对应的行A1 end end end else break; end end d1=transpose(d); y1=x0+l*d1; %new iteration variable 新的迭代变量 t

18、mpf=Funval(f,var,y1); bb=b2-A2*x00; dd=A2*d; if dd>=0 tmpI,lm= ForwardBackMethod(tmpf,0,0.001); %find the searching interval 找到搜索区间 else lm=inf; %find lambda_max for i=1:length(dd) % if(dd(i)>0) % % if dd(i)<0% % if bb(i)/dd(i)<lm lm=bb(i)/dd(i); end end end end if lm=inf lm=1e9; end xm,

19、minf=GoldenSectionSearch(tmpf,0,lm,1.0e-14); %guarantee lambda>0 保证> 0 %find the minimizer by one dimension searching method 找到由一维搜索方法得到目标 tol=norm(xm*d); if tol<eps x=x0; break; end x0=x0+xm*d1; disp('x0'); x0 end minf=Funval(f,var,x); GoldenSectionSearch.m 0.618搜索法确定局部最优值function

20、x,minf=GoldenSectionSearch(f,a,b,eps) %0.618 search method to find minimum value of function f 0.618搜索方法找到函数f的最小值format long; if nargin=3 eps=1.0e-6; end l=a+0.382*(b-a); u=a+0.618*(b-a); k=1; tol=b-a; while tol>eps&&k<100000 fl=subs(f ,findsym(f),l); fu=subs(f ,findsym(f),u); if fl>

21、;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 if k=100000 disp('找不到最小值！'); x=NaN; minf=NaN; return; end x=(a+b)/2; %return the minimizer 返回最小值minf=subs(f, findsym(f),x); format short; ForwardBackMethod.m 进退法确定搜索区间function left,right=ForwardBackMe

22、thod(f,x0,step) if nargin=2 step = 0.01 end if nargin=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; f