zoutendijk可行方向法的matlab实现.doc

(一) 、基本思想是：
给定一个可行点之后，用某种方法确定一个改进的可行方向，然后沿方向，求解一个有约束的线搜索问题，得极小点，令，如果不是最优解，则重复上述步骤。可行方向法就是利用线性规划方法来确定的。
1) 、线性约束问题：
设是问题

的一个可行解，假定，，其中
，
则一个非零向量是在点点的一个可行方向，当且仅当，；如果，则是一改进方向。
2) 、非线性约束问题
设是问题

的一个可行解，令，，即是点紧约束的指标集，设和在点可微，在点连续，如果，，则是一改进的可行方向。
(二) 、算法
1) 、线性不等式约束的Zoutendijk方法的计算步骤：
1.求一初始可行解。，令k＝1，转2。
2.对于可行点，设，，， 求解问题 ，得最优解，如果=0，计算结束，是K—T点；否则转3。
3.求解线搜索问题
(a)
其中

设为(a)式最优解，令，返回2。
2) 、非线性不等式约束的Zoutendijk方法的计算步骤：
1) 选取允许误差，，求一初始可行点，令，转2)。
2）确定指标集。
3) 若，且，计算结束，取；若，且，令，转6)；若，转4)。
4) 令，求解线性规划问题(4-2)的最优解；
5) 若，计算结束，取；否则令，转6)。
6) 求出线搜索问题

的最优解，其中；令，，返回2)。
(三) 、程序源码
1) 、主程序
简单说明：此程序可以处理线性和非线性问题，程序主要由label得值来判断，当label=1时运行线性约束部分，label=0时运行非线性约束部分
function [X0,f_val]=zoutendijk(A,b,x0,Aeq,beq,label)
%自定义函数diff_val(x0)作用是求所给函数在x0出的偏导数
%自定义函数fval(x0)作用是求所给函数在x0出的函数值
format long;
eps=1.0e-6;
x0=transpose(x0);%刚开始给的x0为行向量

func
sz=length(x0);
if label==1
[m,n]=size(A);
%把A分解为A1,A2，其中A1为起作用约束
for k=1:1:100
A1=A;
A2=A;
b1=b;
b2=b;
for i=m:-1:1
if abs(A2(i,:)*x0-b2(i,:)) < 0.1
A2(i,:)=[];
b2(i,:)=[];
end
end
for i=m:-1:1
if abs(A1(i,:)*x0-b1(i,:))>=0.1
A1(i,:)=[];
b1(i,:)=[];
end
end
A1;
A2;
b1;
b2;

i2=rank(A2);
AE=[A1;Aeq];
[i1,j1]=size(AE);
r=rank(AE);
if r<i1
'行不满秩'
return
end

if i2==0
'无效'
return
end
%求解线性规划问题得到可行下降方向d0
s=diff_val(x0);
c=double(s);
lb=-1*ones(sz,1);
ub=ones(sz,1);
k1=length(b1);
k2=length(beq);
p=zeros(k1,1);
q=zeros(k2,1);
[d0,mn,m1,m2,m3]=linprog(c,A1,p,Aeq,q,lb,ub);
d0;mn;

df=abs(s*d0);
if df<0.1
'最优解为：'
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
x0
f_val=fval(x0)
k
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
return
else
%进行一维搜索，求f(x(k+1))的最小值
b_=b2-A2*x0;
d_=A2*d0;
[dh,dl]=size(d_);
ul=1;
for i=1:1:dh
if d_(i,:)>=0
u=1;
else
u=0;
end
ul=ul*u;
end
ul;b_;d_;
vmax=inf;
if ul==0
vmax=inf;
else
for i=1:1:dh
if d_(i,:)>0
v=b_(i,:)/d_(i,:);
if v<vmax
vmax=v;
end
end
end
end
end
vmax;
h=fmin(x0,d0,vmax);
a=x0+h*d0;
f_val=fval(a);
x0=x0+h*d0;
'****************'
X0=x0
f_val=fval(x0)
k
'****************'
end
end
if label==0
for k=1:1:100
%确定指标集
'f(x)梯度：'
sf=diff_val(x0)
sf=eval(sf)

'f(x)梯度长度：'
norm_s=norm(sf)
GL=length(G);
G_copy=G;
for i=1:1:GL
g=subs(G(i,:),x,x0);
G(i,:)=g;
end
G_zero=eval(G);
for i=GL:-1:1
if abs(G_zero(i,:))>0.1
G_zero(i,:)=[];
G_copy(i,:)=[];
I=length(G_zero);
end
end
'x0时为零的g(x)：'
G_copy
'指标集I(x)：'
I
add=-ones(I,1);
%根据指标集确定不同情况下的搜索方向
if I==0
if norm_s<=10
'最优解为：'
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
x0
f_val=fval(x0)
k
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
return
else
d0=-sf;
%线搜索问题
vmax=100;
d0=transpose(d0);
h=fmin(x0,d0,vmax);
x0=x0+h*d0;
end
else
%线性规划问题
grad=jacobian(G_copy,x);
G_zero=subs(grad,x,x0);
G_zero=[G_zero,add];
sf=[sf,-1];
'线性规划问题A矩阵：'
Ac=[sf;G_zero]
lb=-1*ones(sz,1);
ub=ones(sz,1);
p=zeros(I+1,1);
c=zeros(1,sz);
c=[c,1];
[dz,mn,m1,m2,m3]=linprog(c,Ac,p,[],[],lb,ub);
dz;
mn;
sd=length(dz);
d1=dz(1:sd-1,1:1);
z0=dz(sd,1);
z0=abs(z0)
if z0<0.01
'最优解为：'
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
x0
f_val=fval(x0)
k
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
return
else
d0=d1;
%线搜索问题
vmax=10000;
h=fmin(x0,d0,vmax);
a=x0+h*d0;
x0=x0+h*d0;
end
end
k
end
end

2) 、回调函数

l func 记录函数及其自变量信息，如：
syms x1 x2;
f=2*x1^2+2*x2^2-2*x1*x2-4*x1-6*x2;
x=[x1,x2];

l fval 计算函数在x0处得函数值
function f_val=fval(x0)
x0=transpose(x0);
func;
f_val=subs(f,x,x0);

l diff_val(x0)计算函数在x0处得导数值
function s=diff_val(x0)
func
grad=jacobian(f,x);
s=subs(grad,x,x0);

l fmin(x0,d0,vmax)求函数在[0,vmax]上的最小值
function h=fmin(x0,d0,vmax)
func
syms h;
a=x0+h*d0;
f_val=inline(subs(f,x,a));
if vmax==inf
min_h=fminbnd(f_val,0,10000);
else
min_h=fminbnd(f_val,0,vmax);
end
h=min_h;
(四) 、例题
l 线性问题（以例1为例说明数据输入及其最后结果）
例1

最优解
首先把func函数中相应例题的注释去掉；
然后在matlab中输入如下数据：
clear
clc
A=[1 1;1 5;-1 0;0 -1]
b=[2 5 0 0]'
x0=[0 0]
Aeq=[]
beq=[]
label=1
最后运行程序：zoutendijk(A,b,x0,Aeq,beq,label)
得到结果：

例2

最优解
A=[-1 -1;-15 -10;-1 0;0 -1]
b=[-1 -12 0 0]'
x0=[0 2]
Aeq=[]
beq=[]
label=1;
l 非线性问题
例

首先把func函数中相应例题的注释去掉；
然后在matlab中输入如下参数：
clear
clc
A=[]
b=[]
x0=[0 0 0]
Aeq=[]
beq=[]
label=0
最后运行程序：zoutendijk(A,b,x0,Aeq,beq,label)
得到结果：