梯度投影法 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)
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
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),:)];
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)
if bb(i)/dd(i)
lm=bb(i)/dd(i);
end %去掉A1对应的行
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
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;
k = k+1;
f4 = subs(f, findsym(f),x4);
f1 = subs(f, findsym(f),x1);
if f4
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;
end
end
end
minx = min(x1,x3);
maxx = x1+x3 - minx;
format short;
% syms x1 x2 x3 ;
% f=x1^2+x1*x2+2*x2^2+2*x3^2+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=x1^3+x2^2-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;72],[10 10 10],[x1 x2 x3])
% syms x1 x2;
%f=2*x1^2+2*x2^2-2*x1*x2^3-4*x1^7-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*x1^2+2*x2^2-2*x1*x2^3-4*x1^7-6*x2;
% var=[x1,x2];
% valst=[-1,-1];
% A=[1 1;1 5;-1 0;0 -1];
% b=[2 5 0 0]';
% f=x1^3+x2^2-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,fval]=fmincon('confun',valst,A,b)
MinRosenGradientProjectionMethod.m
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 vector 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);
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)
A1(k,:)=A(i,:);
b1(k,1)=b(i);
else%findmatrixs that satisfy A2 x_k
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; %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=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=transpose(d);
y1=x0+l*d1; %new iteration variable 新的迭代变量
tmpf=Funval(f,var,y1);
bb=b2-A2*x00;
dd=A2*d;
ifdd>=0
[tmpI,lm]= ForwardBackMethod(tmpf,0,0.001); %find the searching interval 找到搜索区间
else
lm=inf; %find lambda_max
fori=1:length(dd)
% if(dd(i)>0)
% %
ifdd(i)
% %
if bb(i)/dd(i)
lm=bb(i)/dd(i);
end
end
end
end
if lm==inf
lm=1e9;
end
[xm,minf]=GoldenSectionSearch(tmpf,0,lm,1.0e-14); %guarantee lambda>0 保证λ> 0 %find the minimizer by one dimension searching method 找到由一维搜索方法得到目标 tol=norm(xm*d);
iftol
x=x0;
break;
end
x0=x0+xm*d1;
disp('x0');
x0
end
minf=Funval(f,var,x);
GoldenSectionSearch.m 0.618搜索法确定局部最优值
function [x,minf]=GoldenSectionSearch(f,a,b,eps)
%0.618 search method to find minimum value of function f 0.618搜索方法找到函数f 的最小值
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
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; %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,findsym(f),{x1});
if(f1
step2=2*step;
x2=x1+step;
f2=subs(f,findsym(f),{x2});
while(f1>f2)
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});
while(f0>f2)
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=t^2+s+1
varc=findsym(varvec); %找出传递参数的变量[t s]中的 t s
s1=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);