梯度投影法 MATLAB程序可執(zhí)行

1、.function x,minf=minRosen(f,A,b,x0,var,eps)%目標(biāo)函數(shù)：f;%約束矩陣：A;%約束右端力量：b;%初始可行點(diǎn)：x0;%自變量向量：var;%精度：eps;%目標(biāo)函數(shù)取最小值時(shí)的自變量值：x;%目標(biāo)函數(shù)的最小值：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 A1=A1(1:k,:); b1=b1(1:k,:);endif s0 A2=A2(1:s,:); b2=b2(1:s,:);endwhile 1 P=eye(n,n); if k0 tM=transpose(A1); P=P-tM*inv(A1*tM)*A1; end 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

3、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對應(yīng)的行 end end end else break; 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,0.1); else lm=inf;

4、 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 toleps & k 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 =

5、 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 = 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

6、= f4; h = 2*h; else if k=1 h = -h; x2 = x4; f2 = f4; else x3 = x2; 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

7、-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*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*x

8、23-4*x17-6*x2; % var=x1,x2; % valst=-1,-1; % A=1 1;1 5;-1 0;0 -1; % 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,f

9、val=fmincon(confun,valst,A,b) MinRosenGradientProjectionMethod.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

10、vector of independent variable; 自變量向量 %eps is the precision; 精度 %x is the value of the independent variable when the objective function is minimum; 自變量的值是當(dāng)目標(biāo)函數(shù)最小%minf is the minimum value of the objective function; 目標(biāo)函數(shù)的最小值 format long; if nargin = 5 eps=1.0e-6; end syms l; x00=transpose(x0); n=leng

11、th(var); sz=size(A); m=sz(1);% m is the number of rows of A 行數(shù)gf=jacobian(f,var);%calculate the jacobian matrix of the objective function 計(jì)算目標(biāo)函數(shù)的雅可比矩陣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

12、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_k0 A1=A1(1:k,:); b1=b1(1:k,:); end if s0 A2=A2(1:s,:); b2=b2(1:s,:); end while 1 P=eye(n,n); if k0 tM=transpose(A1); P=P-tM*inv(A1*tM)*A1; %calculate P; end gv=Funval(gf,var,x0); gv=transp

13、ose(gv); d=-P*gv; %calculate the searching direction 計(jì)算搜索方向% flg=1; % if(P=zeros(n) % flg =0; % end % if flg=1 % d=d/norm(d); %normorlize the searching direction 搜索方向% end % 加入這部分會無止境的循環(huán) 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);

14、%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 刪除對應(yīng)的行A1 end end end else break; end end d1=transpose(d); y1=x0+l*d1; %new iteration variable 新的迭代變量 tmpf=Funval(f,var,y1); bb=b2-A2*x00; dd=A2*

15、d; if dd=0 tmpI,lm= ForwardBackMethod(tmpf,0,0.001); %find the searching interval 找到搜索區(qū)間 else lm=inf; %find lambda_max for i=1:length(dd) % if(dd(i)0) % % if dd(i)0% % if bb(i)/dd(i)0 保證 0 %find the minimizer by one dimension searching method 找到由一維搜索方法得到目標(biāo) tol=norm(xm*d); if toleps&kfu a=l; l=u; u=a

16、+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 進(jìn)退法確定搜索區(qū)間function left,right=ForwardBackMethod(f,x0,step) if nargin=2

17、 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(f1f2) 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);

18、 while(f0f2) 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);%left end point right=max(x1,x2);%right end point end Funval.mfunction fv=Funval(f,varvec,varval) var=findsym(f); %找出表達(dá)式包含的變量t,s f=t2+s+1 varc=findsym(varvec); %找出傳遞參數(shù)的變量t s中的 t s s1=length(var); %函數(shù)的個(gè)數(shù) s2=length(varc); %變量的個(gè)數(shù) m=floor(s1-1)/3+1); %靠近左邊的整數(shù) varv=zeros(1,m); if s1 = s2%if the number of variable is different, deal with it specially 如果變量的數(shù)量是不同的,專門處