梯度投影法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 = 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:mdf

2、un=A(i,:)*x00-b(i);if abs(dfun)0A1=A1(1:k,:);b1=b1(1:k,:);endif s0A2=A2(1:s,:);b2=b2(1:s,:);endwhile 1P=eye(n,n);if k0tM=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;elseu,index=m

3、in(w);sA1=size(A1);if sA1(1)=1k=0;elsek=sA1(2);去掉 A1 對應(yīng)的A1=A1(1:(index-1),:);A1(index+1):sA1(2),:);%行 endendendelsebreak;end endd1=transpose(d);y1=x0+l*d1;tmpf=Funval(f,var,y1);bb=b2-A2*x00;dd=A2*d;if dd=0tmpI,lm=minJT(tmpf,0,;elselm=inf;for i=1:length(dd)if dd(i)0if bb(i)/dd(i)lm lm=bb(i)/dd(i);end

4、endendendxm,minf=minHJ(tmpf,0,lm,;tol=norm(xm*d);if toleps & k 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,h

5、0,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

6、 - minx;format short;% syms x1 x2 x3 ;% f=xi2+x1*x2+2*x2A2+2*x3A2+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=x1A3+x2A2-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;

7、72,10 10 10,x1 x2 x3)% syms x1 x2;%f=2*x1A2+2*x2A2-2*x1*x2A3-4*x1A7-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*x1A2+2*x2A2-2*x1*x2A3-4*x1A7-6*x2;% var=x1,x2;% valst=-1,-1;% A=1 1;1 5;-1 0;0 -1;% b=2 5 0 0;% f=x1A3+x2A2-2*x1-4*x2+6;% var=x1,x2;% valst=0 0;% A=2,

8、-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 constrai

9、nt 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;自變量的值是當(dāng)目標(biāo)函數(shù)最小%minf is the minimum valu

10、e of the objective function;目標(biāo)函數(shù)的最小值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行數(shù)gf=jacobian(f,var);%calculate the jacobian matrix of the objective function計(jì)算目標(biāo)函數(shù)的雅可比矩陣bConti=1;while bContik=0;s=0;A1=A;A2=A;b1=b;b2=b;for i

11、=1:mdfun=A(i,:)*x00-b(i); %separate matrix A and bif abs(dfun) %find matrixs that satisfy A1 x_k=b1k=k+1;A1(k,:)=A(i,:);b1(k,1)=b(i);else%find matrixs that satisfy A2 x_k0A1=A1(1:k,:);b1=b1(1:k,:);endif s0A2=A2(1:s,:);b2=b2(1:s,:);endwhile 1找到滿足的矩陣分離矩陣A和b 找到滿足的矩陣P=eye(n,n);if k0tM=transpose(A1);P=P-

12、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% 加入這部分會無止境的循環(huán)計(jì)算搜索方向搜索方向if d=0if k=0x=x0;bConti=0;break;elsew=-inv(A1*tM)*A1*gv;if w=0x=x0

13、;bConti=0;break;elseu,index=min(w);%find the negative component in w sA1=size(A1);if sA1(1)=1 k=0;elsek=sA1(2);A1=A1(1:(index-1),:);A1(index+1):sA1(1),:);%deletecorresponding row in A1 刪除對應(yīng)的行A1endendendelse break;endendd1=transpose(d);y1=x0+l*d1; %new iteration variable新的迭代變量tmpf=Funval(f,var,y1);bb

14、=b2-A2*x00;dd=A2*d;if dd=0tmpI,lm= ForwardBackMethod(tmpf,0,; %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&kfua=l;l=u;u=a+*(b-

15、a);elseb=u;u=l;l=a+*(b-a);endk=k+1;tol=abs(b-a);endif k=100000disp( 找不到最小值！);x=NaN;minf=NaN;return;end返回最小值x=(a+b)/2; %return the minimizer minf=subs(f, findsym(f),x);format short;進(jìn)退法確定搜索區(qū)間function left,right=ForwardBackMethod(f,x0,step)if nargin=2step =endif nargin=1x0=0;step =endf0 =subs(f,findsym

16、(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);endleft=min(x0, x2);right=max(x0, x2);elsestep=2*step;x2=x1-step;f2=subs(f,findsym(f),x2);while(f0f2)x1=x0;x0=x2;f1=f0;f0=f2;step=2*step;x2=x1-step;f2=subs(f,findsym(f),x2);end

17、left=min(x1,x2);%left end pointright=max(x1,x2);%right end pointend function fv=Funval(f,varvec,varval)var=findsym；找出表達(dá)式包含的變量t,s f=tA2+s+1varc=findsym(varvec); % 找出傳遞參數(shù)的變量t s 中的 t ss1=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ù)量是不同的, 專門處理它for i=0