clc; clear;

function [X_opt, Z_min] = solve_transport(C, a, b)
    if sum(a) > sum(b)
        b = [b, sum(a) - sum(b)];
        C = [C, zeros(length(a), 1)];
    elseif sum(a) < sum(b)
        a = [a, sum(b) - sum(a)];
        C = [C; zeros(1, length(b))];
    end
    [m, n] = size(C);
    c_vec = reshape(C', m * n, 1);
    A_eq = zeros(m + n, m * n);
    for i = 1:m
        for j = 1:n
            A_eq(i, (i - 1) * n + j) = 1;
        end
    end
    for j = 1:n
        for i = 1:m
            A_eq(m + j, (i - 1) * n + j) = 1;
        end
    end
    b_eq = [a(:); b(:)];
    lb = zeros(m * n, 1);
    ctype = repmat("S", 1, m + n);
    vartype = repmat("C", 1, m * n);
    [x_vec, Z_min] = glpk(c_vec, A_eq, b_eq, lb, [], ctype, vartype, 1);
    X_opt = reshape(round(x_vec), n, m)';
end

C1 = [1, 3, 1, 4, 2;
      6, 2, 3, 1, 7;
      2, 4, 2, 2, 4];
a1 = [40, 150, 100];
b1 = [20, 80, 90, 60, 40];

[X1, Z1] = solve_transport(C1, a1, b1);
fprintf("Задача 1 (закрытая модель):\n");
fprintf("Минимальные затраты Zmin = %.0f руб.\n", Z1);
fprintf("Оптимальный план поставок X*:\n");
disp(X1);

C2 = [12, 16, 14, 16;
      18, 20, 18, 17;
      19, 14, 16, 12];
a2 = [30, 50, 45];
b2 = [40, 25, 35, 40];

[X2, Z2] = solve_transport(C2, a2, b2);
fprintf("\nЗадача 2 (открытая модель, с фиктивным поставщиком B4 = 15):\n");
fprintf("Минимальные расходы Zmin = %.0f тыс. руб.\n", Z2);
fprintf("Оптимальный план поставок X*:\n");
disp(X2);