Matlab实现Zoutendijk编程例子.doc

用Zoutendijk方法求解下列问题：

s.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 r<i1
'行不满秩'
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.1
'最优解为'
x0
f_val=fval(x0)
k
return
else
%进行一维搜索，求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
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 v<vmax
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