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

(一) 、基本思想是：给定一个可行点之后，用某种方法确定一个改进的可行方向，然后沿方向，求解一个有约束的线搜索问题，得极小点，令，如果不是最优解，则重复上述步骤。可行方向法就是利用线性规划方法来确定的。1) 、线性约束问题：设是问题的一个可行解，假定，其中，则一个非零向量是在点点的一个可行方向，当且仅当，；如果，则是一改进方向。2) 、非线性约束问题设是问题的一个可行解，令，即是点紧约束的指标集，设和在点可微，在点连续，如果，则是一改进的可行方向。(二) 、算法1) 、线性不等式约束的Zoutendijk方法的计算步骤：1.求一初始可行解。，令k1，转2。 2.对于可行点，设， 求解问题 ，得最优解，如果=0，计算结束，是KT点；否则转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为行向量funcsz=length(x0);if label=1m,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 A1(i,:)=; b1(i,:)=; end end A1; A2; b1; b2; i2=rank(A2); AE=A1;Aeq; i1,j1=size(AE); r=rank(AE); if ri1 行不满秩 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 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 v0.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 z00.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 endend2) 、回调函数l func 记录函数及其自变量信息，如：syms x1 x2;f=2*x12+2*x22-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)funcgrad=jacobian(f,x);s=subs(grad,x,x0);l fmin(x0,d0,vmax)求函数在0,vmax上的最小值function h=fmin(x0,d0,vmax)funcsyms 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);endh=min_h;(四) 、例题l 线性问题（以例1为例说明数据输入及其最后结果）例1最优解首先把func函数中相应例题的注释去掉；然后在matlab中输入如下数据：clearclcA=1 1;1 5;-1 0;0 -1b=2 5 0 0x0=0 0Aeq=beq=label=1最后运行程序：zoutendijk(A,b,x0,Aeq,beq,label)得到结果：例2最优解A=-1 -1;-15 -10;-1 0;0 -1b=-1 -12 0 0x0