clear
X=[7	7	7
5	9	7
4	6	5
5	9	8
6	8	5];
Y=[4	4
7	7
5	7
6	2
3	6];
n=size(X,1);
m=size(X',1);
s=size(Y',1);
for i=1:n
    Aeq=[X(i,:) zeros(1,s)];
    beq=1;
    f=[zeros(1,m)  -Y(i,:)];
    A=[-X Y];
    b=zeros(n,1);
    LB=zeros(m+s,1);
    w(:,i)=LINPROG(f,A,b,Aeq,beq,LB);
    c(i)=(w(m+1:m+s,i))'*Y(i,:)';
end

for i=1:n
for j=1:n
    E(i,j)=(w(m+1:m+s,i)'*Y(j,:)')/((w(1:m,i))'*X(j,:)');
end
end
for j=1:n
    a(j)=sum(E(:,j))/n;
end

for num=1:1000
for j=1:n
    for d=1:n
        Aeq=[X(j,:) zeros(1,s)];
        beq=1;
        f=[zeros(1,m)  -Y(j,:)];
        A1=[-X Y];
        A2=[a(d)*X(d,:)  -Y(d,:)];
        A=[A1;A2];
        b=zeros(n+1,1);
        LB=zeros(m+s,1);
        w1(:,d,j)=LINPROG(f,A,b,Aeq,beq,LB);
    end
end

for j=1:n
   for d=1:n
   E1(d,j)=w1(m+1:m+s,d,j)'*Y(j,:)';
   end
   a1(j)=sum(E1(:,j))/n; 
end 

if max(abs(a1-a))<0.0001
  a1
  break
else  
  a=a1;    
end

end