Matlab实现Zoutendijk编程例子.docxVIP

用 Zoutendijk 方法求解下列问题： min f (x) ? x2 ? 2x 2 ? 3x 2 ? x x 2x x x x 4x 6x s.t. 1 x ? 2x 1 2 2 3 1 2 ? x ? 4 3 1 3 2 3 1 2 x , x , x ? 0 1 2 3 取初始点x(1) ? (0,0,0) T ，通过 Matlab 编程实现求解过程。 公用函数如下： 1、function [f,x]=func %设置目标函数syms x1 x2 x3; f=x1^2+2*x2^2+3*x3^2+x1*x2-2*x1*x3+x2*x3-4*x1-6*x2; x=[x1,x2,x3]; end 2、function f_val=fval(x0) %求目标函数值x0=transpose(x0); [f,x]=func; f_val=subs(f,x,x0); end 3、function s=diff_val(x0) %求目标函数梯度[f,x]=func; grad=jacobian(f,x); s=subs(grad,x,x0); end 4、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); end h=min_h; end Zoutendijk 方法主函数 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 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 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 df0.1 最优解为 return else x0 f_val=fval(x0) k %进行一维搜索，求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 end  u=0; 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 **************** end end : · ]”J , 1心”吵 ， 刚日Rh th si \ d口I . 5 肛 Ill0 5 . 4. 1 ｀比id 6 令”“ 记： 5 1.心 让 d x ＾} ) . [