Matlab实现Zoutendijk编程例子

用Zoutendijk方法求解下列问题：22212312132312minf(x)x2x3xxx2xxxx4x6xs.t.123123x2xx4x,x,x0取初始点(1)Tx(0,0,0)，通过Matlab编程实现求解过程。公用函数如下：1、function[f,x]=func%设置目标函数symsx1x2x3;f=x1^2+2*x2^2+3*x3^2+x1*x2-2*x1*x3+x2*x3-4*x1-6*x2;x=[x1,x2,x3];end2、functionf_val=fval(x0)%求目标函数值x0=transpose(x0);[f,x]=func;f_val=subs(f,x,x0);end3、functions=diff_val(x0)%求目标函数梯度[f,x]=func;grad=jacobian(f,x);s=subs(grad,x,x0);end4、functionh=fmin(x0,d0,vmax)%求函数最小值[f,x]=func;symsh;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;endZoutendijk方法主函数function[X0,f_val]=zoutendijk(A,b,x0,Aeq,beq)%自定义函数diff_val(x0)作用是求所给函数在x0出的偏导数%自定义函数fval(x0)作用是求所给函数在x0出的函数值formatlong;eps=1.0e-6;x0=transpose(x0);%刚开始给的x0为行向量[f,x]=func;sz=length(x0);[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)kreturnelse%进行一维搜索，求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'****************'endend