2022年梯度投影法MATLAB程序可执行整理 .pdf

《2022年梯度投影法MATLAB程序可执行整理 .pdf》由会员分享，可在线阅读，更多相关《2022年梯度投影法MATLAB程序可执行整理 .pdf（10页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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

2、 A1=A; A2=A; b1=b; b2=b; fori=1:m dfun=A(i,:)*x00-b(i); if abs(dfun)0 A1=A1(1:k,:); b1=b1(1:k,:); end if s0 A2=A2(1:s,:); b2=b2(1:s,:); end 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 10 页 - - - - - - - - - while 1 P=eye(n,n); if k0 tM=transpose(A1); P=P-tM*

3、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 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

4、break; end end d1=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); else lm=inf; fori=1:length(dd) ifdd(i)0 if bb(i)/dd(i)lm lm=bb(i)/dd(i); end 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 10 页 - - - -

5、- - - - - end end end xm,minf=minHJ(tmpf,0,lm,1.0e-14); tol=norm(xm*d); iftoleps& kfu 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 10 页 - - - - - - - - - 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

6、); end if k = 100000 disp(找不到最小值！); x = NaN; minf = NaN; return; end x = (a+b)/2; minf = subs(f, findsym(f),x); format short; function minx,maxx = minJT(f,x0,h0,eps) format long; ifnargin = 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 = su

7、bs(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 = x2; x2 = x1; x1 = x4; break; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 10 页 - - - - - - - - - end end end minx = min(x1,x3); maxx =

8、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) % syms x1 x2 x3; % f=-x1*x2*x3; % x,mf=minRosen(f,-1,-2,

9、-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.m syms 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; % b=2 5 0 0; % f=x13+x22-2*x1-4*x2+6; % var=x1,x

10、2; % 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) MinRosenGradientProjectionMethod.m 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - -

11、- - - - 名师精心整理 - - - - - - - 第 5 页，共 10 页 - - - - - - - - - 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 vecto

12、r 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; ifnargin = 5 eps=1.0e-6; end syms l; x00=transpose(x0); n=length(var);

13、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; whilebConti k=0; s=0; A1=A; A2=A; b1=b; b2=b; fori=1:m dfun=A(i,:)*x00-b(i); %separate matrix A and b 分离矩阵A 和 b if abs(dfun)0.0000001 %find matrixs t

14、hat satisfy A1 x_k=b1 找到满足的矩阵k=k+1; A1(k,:)=A(i,:); b1(k,1)=b(i); else%findmatrixs that satisfy A2 x_k0 A1=A1(1:k,:); b1=b1(1:k,:); end if s0 A2=A2(1:s,:); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页，共 10 页 - - - - - - - - - b2=b2(1:s,:); end while 1 P=eye(n,n

15、); if k0 tM=transpose(A1); P=P-tM*inv(A1*tM)*A1; %calculate 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=

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

17、nspose(d); y1=x0+l*d1; %new iteration variable 新的迭代变量tmpf=Funval(f,var,y1); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页，共 10 页 - - - - - - - - - bb=b2-A2*x00; dd=A2*d; ifdd=0 tmpI,lm= ForwardBackMethod(tmpf,0,0.001); %find the searching interval 找到搜索区间else lm=

18、inf; %find lambda_max fori=1:length(dd) % if(dd(i)0) % % ifdd(i)0 % % if bb(i)/dd(i)0 保证 0 %find the minimizer by one dimension searching method 找到由一维搜索方法得到目标tol=norm(xm*d); iftoleps&kfu 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(找不到最小值！)

19、; 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=ForwardBackMethod(f,x0,step) ifnargin=2 step = 0.01 end ifnargin=1 x0=0; step = 0.01 end f0 =subs(f,findsym(f),x0); x1=x0+step; f1=subs(f,fi

20、ndsym(f),x1); if(f1f2) x0=x1; x1=x2; f0=f1; f1=f2; step=2*step; x2=x1+step; f2=subs(f,findsym(f),x2); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 9 页，共 10 页 - - - - - - - - - end left=min(x0, x2); right=max(x0, x2); else step=2*step; x2=x1-step; f2=subs(f,findsym(

21、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 left=min(x1,x2);%left end point right=max(x1,x2);%right end point end Funval.m functionfv=Funval(f,varvec,varval) var=findsym(f); %找出表达式包含的变量t,s f=t2+s+1 varc=findsym(varvec); %找出传递参数的变量t s中的t s s1

22、=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 如果变量的数量是不同的,专门处理它fori=0: (s1-1)/3) k=findstr(varc,var(3*i+1); index=(k-1)/3; varv(i+1) = varval(index+1); end fv=subs(f,var,varv); else fv=subs(f,varvec,varval); %return the value of the function 如果原来函数变量个数与传递参数中变量个数一致，调用subs函数，计算给定点处函数值end % disp(here) % f % varvec % varval %fv = subs(f,varvec,varval); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 10 页，共 10 页 - - - - - - - - -