zoutendijk 可行方向法的matlab实现

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

1) 、线性约束问题： 设x是问题







nRfxeExb,Axx s.t.)( min

的一个可行解，假定11bxA，22bxA，其中 21AAA，

21bbb

则一个非零向量d是在点x点的一个可行方向，当且仅当0dA1，0Ed；如果0(Tx)df，则d是一改进方向。 2) 、非线性约束问题 设x是问题







niRmifxxx ,,2,1 0,)(g s.t.)( min



的一个可行解，令0)(,|xxxingRS，0)(|xxigI，即I是Sx点紧约束的指标集，设f和)(Iigi在x点可微，)(Iigi在x点连续，如果0(Tx)df，)(0)(TIigidx，则d是一改进的可行方向。 (二) 、算法 1) 、线性不等式约束的Zoutendijk方法的计算步骤： 1.求一初始可行解0x。，令k＝1，转2。 2.对于可行点kx，设11kAxb，22kAxb，12(,)TTTAAA， 12(,)TTTbbb求解

问题11min . 0 0 11,1,,TjzfxdstAdPEddjn ，得最优解kd，如果Tkfxd=0，计算结束，kx是K—T点；否则转3。

3.求解线搜索问题

maxmin () . 0kkfxdst



(a)

其中 ____max_22_2min |0,0 ,,0kk

bdddbbAxdAdd







设k为(a)式最优解，令1,1kkkkxxdkk，返回2。 2) 、非线性不等式约束的Zoutendijk方法的计算步骤： 1) 选取允许误差01，02，求一初始可行点)1(x，令1k，转2)。 2）确定指标集0)(|)()()(kikgiIxx。 3) 若)()(kIx，且1)()(kfx，计算结束，取)(*kxx；若)()(kIx，且1)()(kfx，令)()(kkfxd，转6)；若)()(kIx，转4)。

4) 令)(kxx，求解线性规划问题(4-2)的最优解),(kkzd； 5) 若2kz，计算结束，取)(kxx；否则令kkdd＝，转6)。 6) 求出线搜索问题

max)(0 ..)( mintsfkkdx

的最优解k，其中Skkdx)(max|max；令kkkdxx)()1(，1kk，返回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'行不满秩' 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 vvmax=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)