Matlab实现Zoutendijk编程例子.doc

用Zoutendijk方法求解下列问题：s.t. 取初始点，通过Matlab编程实现求解过程。公用函数如下：1、function f,x=func%设置目标函数syms x1 x2 x3;f=x12+2*x22+3*x32+x1*x2-2*x1*x3+x2*x3-4*x1-6*x2;x=x1,x2,x3; end2、function f_val=fval(x0)%求目标函数值x0=transpose(x0);f,x=func;f_val=subs(f,x,x0);end3、function s=diff_val(x0)%求目标函数梯度f,x=func;grad=jacobian(f,x);s=subs(grad,x,x0);end4、function h=fmin(x0,d0,vmax)%求函数最小值f,x=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);endh=min_h;endZoutendijk方法主函数function X0,f_val=zoutendijk(A,b,x0,Aeq,beq)%自定义函数diff_val(x0)作用是求所给函数在x0出的偏导数%自定义函数fval(x0)作用是求所给函数在x0出的函数值format long;eps=1.0e-6;x0=transpose(x0); %刚开始给的x0为行向量f,x=func;sz=length(x0);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 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 vvmax 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 *endend