模型见图片文件，以下代码不能运行，请高手赐教：

clear
X=[3696.50 290.90 610.17 63.86 823.02 851.40 354.75 404.42];  
Y=[25374.85 1025.83 1552.10 358.26 5359.02 5507.94 2566.74 1668.75];  
U=[12554.03 1092.84 2232.17 252.84 3307.42 2257.41 844.61 743.12];  
N=size(X,1);
K=size(X',1);
M=size(Y,1);
J=size(U,1);

lb=zeros(K+N+M+1+1,1);
lb(K+N+M+1:K+N+M+1+1)=-inf;
ub=[];
ub(K+N+M+1)=1;
A=zeros(1,K+N+M+1+1);
b=0;

epsilon=10^-10;
f=[zeros(1,K) -epsilon*ones(1,N+M) 0 1];
for i=1:K
Aeq=[X eye(N) zeros(N,M) -X(:,i) zeros(N,1)
     Y zeros(M,N) -eye(M) zeros(M,1) zeros(M,1)
     U zeros(J,N) zeros(J,M) zeros(J,1) -U(:,i)
     ones(1,K) zeros(1,N) zeros(1,M) -1 0];
beq=[zeros(N,1)
     Y(:,i)
     zeros(J,1)
     0];
w(:,i)=linprog(f,A,b,Aeq,beq,lb,ub);   
end

w
z=w(1:K,:)
s_minus=w(K+1:K+N,:)
s_plus=w(K+N+1:K+N+M,:)
beta=w(K+N+M+1,:)
lamda=w(K+N+M+1+1,:)