bl = [40;60;70;100];
A1 =[1,2,6,2;
    2,6,9,4;
    5,7,1,5;
    4,8,9,8];
[A1_n, A1_m]=size(A1);
price=[4,9,10.5,6.5];
[A1_n, A1_m]=size(A1);
d = [0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 1, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 11, 11, 3, 7, 7, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 11, 11, 0, 0, 0, 11, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 11, 0, 0, 0, 3, 7, 7, 11, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 11, 11, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 3, 1, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 3, 3, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 3, 5, 5, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 11, 11, 0, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 7;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 7;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 3;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];

c=[0, 0, 0, 0, 0, 0, 10, 0, 0, 0, 102, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 10, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 10, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 13, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0, 9, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 10, 0, 21, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0, 11, 11, 102, 10, 60, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0, 11, 11, 0, 0, 0, 11, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 11, 0, 0, 0, 20, 70, 2, 11, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 11, 11, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 30, 100, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10, 17, 39, 31, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10, 36, 35, 25, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 11, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 31, 71;
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 78, 57;
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 47, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 72, 31;
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
a=[10,5];
[a_n, a_m]=size(a);

a=a(:);
[n m]=size(c);
nm=n*m;
Aeq = zeros(n, nm);
for i = 1:n
    for j = 1:m
        if d(i, j) > 0
            Aeq(i, (i - 1)*n + j) = 1;
        end
    end
    for k = 1:n
         if d(k, i) > 0
            if (i <= n - size(a))
                Aeq(i, (k - 1)*n + i) = -1;
            else
                Aeq(i, (k - 1)*n + i) = 1;
            end
         end
    end
end

Aeq1_1=[zeros(n,A1_m),Aeq];
Aeq1_1(1:A1_m,:)=0;
Aeq1_1(n-a_m+1:n,:)=0;
Aeq1_2=[-eye(A1_m),Aeq(1:A1_m,1:nm)];
Aeq1_3=[zeros(a_m,A1_m),Aeq(n-a_m+1:n,1:nm)];
Aeq1=[Aeq1_1;Aeq1_2;Aeq1_3];
%end =

% <=
A_1=[A1,zeros(A1_n,nm)];
A_3=[eye(A1_m),zeros(A1_m,nm)];
A=[A_1];
%end <=

b=bl;
beq=[zeros(n+A1_m,1);a];
ub=[ones(1,A1_m)*inf,reshape(d',1,[])];
lb=zeros(1,A1_m+nm);
intcon=1:(A1_m+nm);
ff=[-price';reshape(c',[],1)];

[x, fval, exitflag,lambda]=intlinprog(ff,intcon,A,b,Aeq1,beq,lb,ub);
res = x(A1_m+1:A1_m+nm,1);
u=x(1:A1_m,:);
res = reshape(res, n, []);

g = digraph;

for i = 1:n
    for j = 1:m
        if not(res(i,j) == 0 )
            g = addedge(g, j, i, res(i,j));
        end
    end
end
figure(1)
h = plot(g,'Layout','layered','EdgeLabel',g.Edges.Weight,'LineWidth', 2,'ArrowSize', 8, 'MarkerSize', 14, 'LineStyle', '-', 'Marker','o')
hs = struct(h);   %and ignore the warning message
for nhl = hs.NodeLabelHandles_
  nhl.Font.Size = 20;
end
for nhl = hs.EdgeLabelHandles_
  nhl.Font.Size = 25;
end
grid on
figure(2)
gg=digraph;
for i = 1:n
    for j = 1:m
        if not(d(i,j) == 0 )
            gg = addedge(gg, i, j, d(i,j));
        end
    end
end
hh = plot(gg,'Layout','layered','EdgeLabel',gg.Edges.Weight,'LineWidth', 2,'ArrowSize', 8, 'MarkerSize', 14, 'LineStyle', '-', 'Marker','o')
hhs = struct(hh);   %and ignore the warning message
for nhl = hhs.NodeLabelHandles_
  nhl.Font.Size = 20;
end
for nhl = hhs.EdgeLabelHandles_
  nhl.Font.Size = 25;
end

zapasi=bl-(sum((A1.*u')'))'
u
TC=sum(sum(res.*c'))
TR=price*u
Profit=TR-TC