Solves the Concurrent Flow variant of the Multi-Commodity Flow Problem in MATLAB with linprog.

mcfp_cf_script.m

%% Solve the Multi-Commodity Flow Problem in MATLAB.
%
% The goal of this script is to solve the Concurrent-Flow variant of the 
% Multi-Commodity Flow Problem ( MCFP ) with MATLAB's Linear Programming 
% ( LP ) tools. 
%
% The documentation generated from this script is a tex file. The PDF can 
% be generated from the tex, but \usagepackage{amsmath} needs to be added 
% to the tex file's preamble.

% It's important the workspace is cleared
% before executing this script!
clear all;
close all;

%% Define the Flow Network.
% 
% The flow network is defined as the following.
%
% <latex>
% \begin{align*}
% V & = \text{Node Set}\\
% E & = \text{Directed Edge ( Arc ) Set}\\
% e & = (u,v: u,v\in V) \\
% & \in E \\
% G & = \text{Directed Graph} \\
% & = (V,E) \\
% c & = \text{Capacity, } \forall e\in E \\
% \end{align*}
% </latex>

% Define the Nodes.
V = {'C','F','B'};

% Define the Arcs.
tuples = { ...
    { 'eAB' 'F' 'B' 1 }, ...
    { 'eBC' 'B' 'C' 1 }, ...
    { 'eCA' 'C' 'F' 1 } };
E = {};
for tuple=tuples
    e = tuple{1}{1};
    E{end+1} = e;
    E_s.(e) = tuple{1}{2};
    E_t.(e) = tuple{1}{3};
    E_c.(e) = tuple{1}{4};
end

% Plot the graph.
figure; 
plot(digraph(struct2cell(E_s),struct2cell(E_t)));
title( 'Digraph G' );

%% Define the Commodities.
%
% The commodities are each represented as a 3-tuple, containing a Source Node,
% a Target Node, and a Demand.
%
% <latex>
% \begin{align*}
% K & = \text{Commodity Set} \\
% & \subseteq V \times V \times \mathbf{R}_{\ge 0}\\
% s_k & = \text{Source Node of $k\in K$}\\
% t_k & = \text{Target Node of $k\in K$}\\
% d_k & = \text{Demand of $k\in K$}\\
% k & = (s_k,t_k,d_k) \\
% & \in K
% \end{align*}
% </latex>

% Define the commodities.
tuples = { 
    { 'k0' 'F' 'C' 0.5 },... 
    { 'k1' 'B' 'F' 0.5 } };
K = {};
for tuple=tuples
    k = tuple{1}{1};
    K{end+1} = k;
    K_s.(k) = tuple{1}{2};
    K_t.(k) = tuple{1}{3};
    K_d.(k) = tuple{1}{4};
end

%% Structure MCFP.
%
% In order to run MATLAB's linprog, the problem needs to be adapted to fit
% under the following form.
%
% <latex>
% \begin{align*}
% \text{Objective: } & \min_x o^T x \\
% \text{Constraints: } & A_{in}x \le b_{in},\\
% & A_{eq} x = b_{eq}, \\
% & b_l \le x \le b_h
% \end{align*}
% </latex>
%
% $o$, $b_{in}$, and $b_{eq}$ are all constant vectors. $A_{in}$ and $A_{eq}$ are
% constant matrices. The bounds $b_l$ and $b_u$ can include constants or
% infinity. Finally, $x$ is the vector of variables, most of which in MCFP
% are the flow of commodity over each edge.

% Declaration starting indices.
col = 1;
rowin = 1;
roweq = 1;

% The following anonymous function makes string comparison easier.
cs = @(s1,s2)strcmp(s1,s2);

%% MCFP Variables.
%
% The MCFP problem has the following MCFP variables.
%
% <latex>
% \begin{align*}
% f_{k,e} & = \text{Flow of Commodity $k$ at Edge $e$}, \forall k\in K, \forall e\in E \\
% & \ge 0 \\
% z & = \text{Objective Variable} \\
% & \ge 0
% \end{align*}
% </latex>

% Define the necessary objective variable with its constraint.
z = col;
bu(col) = Inf;
bl(col) = 0;
col = col+1;

% Define the flows with their constraints.
for k=K
    for e=E
        f.(k{1}).(e{1}) = col;
        bu(col) = Inf;
        bl(col) = 0;
        col = col+1;
    end
end

%% Define the Objective.
%
% Say Percentage of Injected Flow ( PIF ) describes the ratio of a 
% commodity's flow injected into the network over the commodity's demand.
% The objective of MCFP can then be described as maximizing the smallest
% PIF of all the commodities. 
%
% <latex>
% \begin{align*}
% z_k & = \text{PIF of Commodity $k\in K$} \\
% & = \frac{\sum_{v\in V} f_{k,(s_k,v)}}{d_k} \\
% \text{Objective } & : \text{ Maximize } \min\{ z_k : k \in K\}
% \end{align*}
% </latex>

% Define the objective constraints.
for k=K
    s = K_s.(k{1});
    d = K_d.(k{1});
    for v=V
        for e=E
            if ( cs(E_s.(e{1}),s) && cs(E_t.(e{1}),v{1}) )
                Ain(rowin,f.(k{1}).(e{1})) = -1/d;
            end
        end
    end
    Ain(rowin,z) = 1;
    rowin = rowin+1;
end

% Define the objective function.
o = zeros(1,col-1);
o(z) = -1;

%% Define the Capacity Constraint
%
% The Capacity Constraint can be summarized as the sum of flow at any edge
% must be less than or equal to the capacity of the edge.
%
% <latex>
% \begin{align*}
% c \ge \sum_{k \in K} f_{k,e}, \forall e \in E
% \end{align*}
% </latex>

% Define the capacity constraint.
for e=E
    for k=K
        Ain(rowin,f.(k{1}).(e{1})) = 1;
    end
    bin(rowin) = E_c.(e{1});
    rowin = rowin+1;
end

%% Define the Conservation Constraints.
%
% Say net flow describes for a commodity the sum of flow going into and out
% of a node. The conservation constraint can then be defined as net flow
% should equal zero if the node isn't any commodity's source or target.
% But, if the node is a source, the netflow should equal the demand. If a
% target, netflow should equal negative demand.
%
% <latex>
% \begin{align*}
% \Delta_{k,v} &= \text{Net Flow for Commodity $k$ of Edge $e$}\\
% & = \sum_{w\in V} f_{k,(v,w)} - \sum_{w\in V} f_{k,(w,v)}\\
% & = 0, \forall k \in K, \forall v \in \{w \in V: w\neq s_k, t_k\} \\
% & = d_k, \forall k \in K, \forall v \in \{w \in V: w = s_k\} \\
% & = -d_k, \forall k \in K, \forall v \in \{w \in V: w = t_k\} \\
% \end{align*}
% </latex>

% Define the conservation constraint.
for k=K
    s = K_s.(k{1});
    t = K_t.(k{1});
    d = K_d.(k{1});
    for v=V
        for w=V
            for e=E
                if ( cs(E_s.(e{1}),v{1}) && cs(E_t.(e{1}),w{1}) )
                    Aeq(roweq,f.(k{1}).(e{1})) = 1;
                end
                if ( cs(E_s.(e{1}),w{1}) && cs(E_t.(e{1}),v{1}) )
                    Aeq(roweq,f.(k{1}).(e{1})) = -1;
                end
            end
        end
        if ( ~cs(v{1},s) && ~cs(v{1},t) )
            beq(roweq) = 0;
        elseif ( v{1}==s )
            beq(roweq) = d;
        elseif ( v{1}==t )
            beq(roweq) = -d;
        end
        roweq = roweq+1;
    end
end

%% Optimized with MATLAB's linprog.
%
% Finally, the LP problem is solved with MATLAB's linprog.

% Apply MATLAB's linprog function to solve the LP problem.
[ x, fval, exitflag, output ] = ...
    linprog( o, Ain, bin, Aeq, beq, bl, bu );

%% Display the results.

% Build table with results.
Variables = {};
Variables{z} = 'z';
for k=K
    for e=E
        Variables{f.(k{1}).(e{1})} = ['f_',k{1},e{1}];
    end
end
T = table;
T.Variables = Variables';
T.Values = x;

% Display the results.
output
T