梯度投影法MATLAB程序可执行

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;

A1=A;

A2=A;

b1=b;

b2=b;

fori=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);

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;

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 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=bb(i)/dd(i);

end

end

end

end [xm,minf]=minHJ(tmpf,0,lm,1.0e-14);

tol=norm(xm*d);

iftol

x=x0;

break;

end x0=x0+xm*d1;

%

disp('x0');

x0

end

minf=Funval(f,var,x)

functionfv = 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 fori=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)); end

fv = subs(f,var,varv);

else

fv = subs(f,varvec,varval);

end function [x,minf] = minHJ(f,a,b,eps) format long;

ifnargin == 3

eps = 1.0e-6;

end

l = a + 0.382*(b-a);

u = a + 0.618*(b-a);

k=1;

tol = b-a;

whiletol>eps&& k<100000

fl = subs(f , findsym(f), l);

fu = subs(f , findsym(f), u); iffl>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; 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;