zoutendijk-可行方向法的matlab实现

2012-2013（1）专业课程实践论文可行方向法夏文春，0818180226，R数学08-2班姜皓缤，0618180213，R数学08-2班肖云龙，0818180214，R数学08-2班算法理论(一)、基本思想是：给定一个可行点)(kx之后，用某种方法确定一个改进的可行方向kd，然后沿方向kd，求解一个有约束的线搜索问题，得极小点k，令kkkkdxx)()1(，如果)1(kx不是最优解，则重复上述步骤。可行方向法就是利用线性规划方法来确定kd的。1)、线性约束问题：设x是问题nRfxeExb,Axxs.t.)(min的一个可行解，假定11bxA，22bxA，其中21AAA，21bbb则一个非零向量d是在点x点的一个可行方向，当且仅当0dA1，0Ed；如果0(Tx)df，则d是一改进方向。(二)、算法1)、线性不等式约束的Zoutendijk方法的计算步骤：1.求一初始可行解0x。，令k＝1，转2。2.对于可行点kx，设11kAxb，22kAxb，12(,)TTTAAA，12(,)TTTbbb求解问题11min.0011,1,,TjzfxdstAdPEddjn，得最优解kd，如果Tkfxd=0，计算结束，kx是K—T点；否则转3。3.求解线搜索问题maxmin().0kkfxdst(a)其中____max_22_2min|0,0,,0kkbdddbbAxdAdd设k为(a)式最优解，令1,1kkkkxxdkk，返回2。二、算法框图五流程图YYNYNbegin误差21,初始点0,0kX)(kXJ？12)(kXfoverkXmin()(),()11,1,2,...kTkTkjifXgXjJXdin2kover0:min(max(0,1,2,...kkkkkjfXDgXDjn））1:1kkkXXDkk算法代码1)、主程序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);%把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.1A1(i,:)=[];b1(i,:)=[];endendA1;A2;b1;b2;i2=rank(A2);AE=[A1;Aeq];[i1,j1]=size(AE);r=rank(AE);ifri1'行不满秩'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);ifdf0.1'最优解为：''@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'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,:);ifvvmaxvmax=v;endendendendendvmax;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,-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)ifz00.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+2*x2^2-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.利用可行方向法求解22121212121212min22246.25500fxxxxxxxstxxxxxx解：输入A=[11;15;-10;0-1]b=[2500]'x0=[00]Aeq=[]beq=[]label=1最后运行程序：zoutendijk(A,b,x0,Aeq,beq,label)得到结果：例2.利用可行方向法求解2212121212min4.115101200fxxxstxxxxxx解：输入A=[-1-1;-15-10;-10;0-1]b=[-1-200]'x0=[02]Aeq=[]beq=[]label=1;最后运行程序：zoutendijk(A,b,x0,Aeq,beq,label)得到结果：