zoutendijk 可行方向法的matlab实现

(完整)zoutendijk可行方向法的matlab实现(一)、基本思想是：给定一个可行点x(Q之后，用某种方法确定一个改进的可行方向％，然后沿方向1，求解一个有约束的线搜索问题，得极小点人，令x(k+1)=x(k)+人d，如果X(k+1)不是最优解，则重复上述步骤。可行方向法就是利用线性规划方法来确定dk的。1)、线性约束问题：设X是问题minf(x)<s.t.Ax<b,Ex=exWRn的一个可行解，假定Aix=«，A2x<b2，其中b1b-2」则一个非零向量d是在点x点的一个可行方向，当且仅当Ad<0，Ed=0；如果Vtf(x)d<0,则d是一改进方向。2)、非线性约束问题设x是问题minf(x)<s.t.g(x)<0,i=1,2,,mxWRn的一个可行解，令S=11xWRn,g(x)<。｝,/=4lg.(x)=0｝，即I是xWS点紧约束的指标集，设f和g(iwI)在x点可微，g(iWI)在x点连续，如果Vtf(x)d<0，Vtg(x)d<0(iwI)，则d是一改进的iii可行方向。(二)、算法(完整)zoutendijk可行方向法的matlab实现1)、线性不等式约束的Zoutendijk方法的计算步骤：1。求一初始可行解x。，令k=1,转2。2.对于可行,设Ax=b,Ax<b,At=(At,At),br=(br,br)求解问题1k12k21212minz=Vf(x>ds.tAd<01Ed=0-1<d<1,j=],•••,nj，得最优解d，如果Vf(x》d=0,计算结束，x是K—T点；否则

kkk转3。3。求解线搜索问题minf(x+Xd)s.t0<X<Xmax(a)其中XmaxminhjdId>°},d<w0bbadAd+3,d<0一一一--22k，--2k设x为(a)式最优解，令x=x+Xd,k=k+1,返回2.kk+1kkk2)、非线性不等式约束的Zoutendijk方法的计算步骤:1)选取允许误差81>0,82>0，求一初始可行点X(1)，令k=1,转2)。2)确定指标集I(x(k))=I|g(x(k))=°}。若I(x(k))=0，且Vf(x(k))<81，计算结束，取x*=X(k)；若I(x(k))=0，且Vf(x(k))>8］，令dk=-Vf(x(k))，转6);若I(x(k))。0，转4)。令x=x(k),求解线性规划问题(4-2)的最优解(Zk,dk);若zk<8，计算结束，取x=x(k)；否则令d=dk，转6)。6)求出线搜索问题Iminf（x（k）+Xd）[s.l.0<X<Xmax的最优解Xk，其中X=maxi|x（k）+XdeS};令x（k+1）=x（k）+Xd,k=k+1，返回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出的函数值formatlong;eps=1。0e_6;x0=transpose（x0）;%刚开始给的x0为行向量funcsz=length(x0);iflabel==1[m,n]=size(A);%把入分解为A1,A2，其中A1为起作用约束fork=1:1:100A1=A;A2=A；b1=b;b2=b；fori=m:一1:1ifabs(A2(i,：)*x0—b2(i,：))<0.1A2(i,:)=[];b2(i,:)=[];endendfori=m:一1:1ifabs(A1(i,：)*x0-b1(i,:))>=0.1endendA1；A2;b1；b2;i2=rank(A2)；AE=[A1；Aeq]；[i1,j1]=size(AE)；r=rank(AE)；ifr〈i1'行不满秩’returnendifi2==0，无效，returnend%求解线性规划问题得到可行下降方向d0s=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);ifdf〈0.1'最优解为：’(完整)zoutendijk可行方向法的matlab实现’@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'x0f_val=fval(x0)k'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'returnelse%进行一维搜索，求f(x(k+1))的最小值b_=b2-A2*x0；d_=A2大d0；[dh,dl]=size(d_);ul=1;fori=1:1:dhifd_(i,：)>=0u=1;elseu=0;endul=ul大u;endul;b_;d_;vmax=inf;iful==0vmax=inf;elsefori=1:1:dhifd_(i,:)〉0v=b_(i,:)/d_(i,:);ifv〈vmaxvmax=v;endendendendendendvmax;h=fmin(x0,d0，vmax);a=x0+h*d0；f_val=fval(a)；x0=x0+h大d0;’*大大大*大*大大大大**大大大'X0=x0f_val=fval(x0)k'大***大大*大大大大大***大'endendiflabel==0fork=1：1：100%确定指标集'f(x)梯度：’sf=diff_val(x0)sf=eval(sf)'f(x)梯度长度：'norm_s=norm(sf)GL=length(G);G_copy=G;fori=1：1：GLg=subs(G(i,:),x,x0)；G(i,:)=g;endG_zero=eval(G);fori=GL：-1：1ifabs(G_zero(i,:))>0.1G_zero(i,:)=[];G_copy(i,：)=口；I=length(G_zero)；endend’x0时为零的g(x):’G_copy'指标集I(x)：’Iadd=—ones(I，1)；%根据指标集确定不同情况下的搜索方向ifI==0ifnorm_s〈=10'最优解为：''@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'x0f_val=fval(x0)k’@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'returnelsed0=一sf;%线搜索问题vmax=100;d0=transpose(d0);h=fmin(x0,d0,vmax);x0=x0+h大d0;endelse%线性规划问题grad=jacobian(G_copy,x)；G_zero=subs(grad,x,x0)；G_zero=[G_zero,add];sf=[sf,—l];'线性规划问题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)ifz0<0。01'最优解为：''@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'x0f_val=fval(x0)k’@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'returnelsed0=d1;%线搜索问题vmax=10000;h=fmin(x0,d0,vmax);a=x0+h大d0;x0=x0+h大d0;endendkendend2)、回调函数func记录函数及其自变量信息，如：symsx1x2;f=2大x1"+2大x2”—2*x1大x2—4大x1—6*x2;x=[x1,x2];•fval计算函数在x0处得函数值functionf_val=fval(x0)x0二transpose(x0)；func；f_val=subs(f，x，x0)；diff_val(x0)计算函数在x0处得导数值functions=diff_val(x0)funcgrad二jacobian(f，x)；s=subs(grad,x,x0)；fmin(x0，d0,vmax)求函数在［0,vmax］上的最小值functionh=fmin(x0,d0,vmax)funcsymsh;a=x0+h大d0；f_val=inline(subs(f,x,a))；ifvmax==infmin_h=fminbnd(f_val,0,10000);elsemin_h=fminbnd(f_val,0,vmax)；endh=min_h;（四）、例题线性问题（以例1为例说明数据输入及其最后结果）例1minf（x）=2x2+2x2-2xx-4x-6xs.x1+x2<2<x+5x<5-1气<0-x2<0最优解x*=—，癸，f（x*）=-7.16"3131首先把func函数中相应例题的注释去掉；然后在matlab中输入如下数据：clearclcA=[11;15；-10;0-1]b=[2500]’x0=[00]Aeq=[]beq=口label=1最后运行程序：zoutendijk(A,b,x0,Aeq,beq,label)得到结果：kO二1.129032258064410.77419354838709f_v^l=-7.1612903225805^k=3f(x)=X2+4f(x)=X2+4X2mins.t15x+10x>12x>01x>0最优解x*=—,L,f(x*)=—"55J5A=[—1—1;—15-10;-10;0-1]b=[—1-1200],x0=[02]Aeq=[]beq=[]label=1;非线性问题

mins.tf（尤）=X2+xx+2X2—6x-2x-12x2