梯度投影法-MATLAB程序可执行

。

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;

end

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

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;

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

end

end

end

end

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

tol=norm(xm*d);

。

if tol

x=x0;

break;

end

x0=x0+xm*d1;

%

disp('x0');

x0

end

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

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 > 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;

if nargin == 3

eps = 1.0e-6;