clear all; clc;
% n power plants, m cities
n = 5;
m = 25;
% production capacity of the power plants:
w = 10*rand(n, 1);
% demands of the cities
d = rand(m, 1);
% price of sending energy from each plant to each city
C = rand(n, m);

% Make sure the total production capacity is at least the total demand
assert(sum(w) >= sum(d), 'Production capacity is too low to meet demand.');

%%
cvx_begin
    variable X(n, m)
    minimize sum(C(:).*X(:))
    subject to
        X >= 0; % can't send negative amounts of energy
        for i = 1 : n % for each power plant
            sum(X(i, :)) <= w(i); % can't exceed production capacity
        end
        for j = 1 : m % for each city
            sum(X(:, j)) >= d(j); % must meet the city's demand
        end
cvx_end

%%
clf;
imagesc(X);
colorbar; title('How much is sent from each plant to each city?');
xlabel('City #'); set(gca, 'XTick', 1:m);
ylabel('Plant #'); set(gca, 'YTick', 1:n);

% Compare demand to what we send
[sum(X, 1)', d]
% Compare what we send to production capacity
[sum(X, 2), w]