Solves a Time Division Multiplexing Concurrent Flow variant of the Multi Commodity Flow Problem in MATLAB.

matlab_mcfp_cf_tdm_lp.m

%% Solve the Time Division Multiplexing variant of the Multi Commodity Flow Problem in MATLAB.
%
% The goal of this script is to solve the Multi Commodity Flow Problem
% modified for Time Division Multiplexing ( TDM ) as seen in Contention
% Free Routing ( CFR ). The Multi Commodity Flow Problem is solved for the
% Concurrent Flow variant ( MCFP-CF ) with mixed-integer linear programming
% ( MILP ). 
%
% The documentation generated from this script is a tex file. The PDF can 
% be published 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*}
% G & = \text{Directed Graph} \\
% & = (V,E) \\
% V & = \text{Node Set} \\
% E & = \text{Arc Set} \\
% e & = (v_e\in V,w_e\in V) \\
% & \in E \\
% v_e & = \text{Source Node of Arc $e\in E$} \\
% w_e & = \text{Sink Node of Arc $e\in E$} \\
% \end{align*}
% </latex>


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

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

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

%% Define the Commodities.
%
% The commodities are defined as follows.
% 
% <latex>
% \begin{align*}
%     K & = \text{Commodity Set} \\
%     k & = (s_k\in V,t_k\in V,d_k\in \mathbf{Z}_{\ge 0} : s_k \neq t_k) \\
%     & \in K \\
%     s_k & = \text{Source Node of Commodity $k\in K$} \\
%     t_k & = \text{Sink Node of Commodity $k\in K$} \\
%     d_k & = \text{Demand of Commodity $k\in K$} \\
% \end{align*}
% </latex>

% Define the commodities.
tuples = { 
    { 'k1' 'A' 'C' 1 },... 
    { 'k2' 'B' 'A' 2 } };
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

%% Define the Slots.
%
% An important distinction between TDM MCFP-CF over regular MCFP-CF is that
% the arcs are divided into several slots.
%
% <latex>
% \begin{align*}
% L & = \text{Slot Sequence} \\
% & = (1,2,...,c) \\
% c & = \text{Capacity of each Arc $e\in E$ } \\
% & = \sum_{k\in K} d_k
% \end{align*}
% </latex>

% Define the Slots.
c = sum(struct2array(K_d));
L = 1:c;

%% Define variables needed for script.

% Variable related declaration.
col = 1;
intcont = [];

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

%% Define the Flows and Objective Function.
%
% The objective function has the following form. The goal is to maximize
% the smallest supply-demand ratio of all the commodities.
%
% <latex>
% \begin{align*}
% \text{Maximize} & : \min\left\{ \frac{\sum_{v\in V,l \in L}f_{k,(s_k,v),l}}{d_k}: k \in K \right\} \\ 
% f_{k,e,l} & = \text{Flow of Commodity $k\in K$ on Edge $e\in E$ for Slot $l\in L$} \\ 
% & \in \{0,1\} \\
% \end{align*}
% </latex>

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

% Define the variable related to the objective.
z = col;
bu(col) = Inf;
bl(col) = 0;
col = col+1;

% Constraints related declaration.
rowin = 1;
roweq = 1;
Ain = zeros(1,col-1);
bin = zeros(1,1);
Aeq = zeros(1,col-1);
beq = zeros(1,1);

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

%% Define the Slot Capacity Constraint.
%
% A slot can only be filled by a single commodity 
%
% <latex>
% \begin{align*}
% 1 & \ge \sum_{k\in K} f_{k,e,l} : \forall e\in E,\forall l \in L \\ 
% \end{align*}
% </latex>

% Define the Slot Capacity Constraint.
for e=E
    for l=L
        for k=K
            Ain(rowin,f.(k{1}).(e{1})(l)) = 1;
        end
        bin(rowin) = 1;
        rowin = rowin+1;
        
    end
end

%% Define the TDM Constraint.
%
% As defined in CFR, each node can buffer a single flow unit, thus the flow
% of a commodity should shift to an adjacent slot between edges who share the
% same node. However, this constraint should be broken for the commodity's
% source and sink nodes.
%
% <latex>
% \begin{align*}
% 0 & = \sum_{w \in V} f_{k,(w,v),l} - \sum_{w \in V} f_{k,(v,w),(l\text{ mod}|L|)+1} \\
%     & : \forall k\in K,\forall v\in \{w\in V: w\neq s_k,t_k \},\forall l \in L \\
% \end{align*}
% </latex>

% Define the TDM Constraint.
for k=K
    k_s = K_s.(k{1});
    k_t = K_t.(k{1});
    for v=V
        if ( cs(v{1},k_s) || cs(v{1},k_t) )
            continue;
        end
        for l=L
            for w=V
                for e=E
                    if ( cs(w{1},E_s.(e{1})) && cs(v{1},E_t.(e{1})) )
                        Aeq(roweq,f.(k{1}).(e{1})(l)) = 1;
                    end
                    if ( cs(v{1},E_s.(e{1})) && cs(w{1},E_t.(e{1})) )
                        Aeq(roweq,f.(k{1}).(e{1})(mod(l,c)+1)) = -1;
                    end
                end
            end
            beq(roweq) = 0;
            roweq = roweq+1;
        end
    end
end

%% Define the Supply Equals Demand Constraint.
%
% The supply of a commodity's source node should equal the demand of the
% commodity's sink node.
%
% <latex>
% \begin{align*}
% 0 & = \sum_{v\in V,l\in L} f_{k,(s_k,v),l} - \sum_{v\in V,l\in L} f_{k,(v,t_k),l} \\
%     & : \forall k\in K \\
% \end{align*}
% </latex>

% Define the Supply Equals Demand Constraint.
for k=K
    k_s = K_s.(k{1});
    k_t = K_t.(k{1});
    for e=E
        e_s = E_s.(e{1});
        e_t = E_t.(e{1});
        if ( cs(e_s,k_s) && cs(e_t,k_t) )
            Aeq(roweq,f.(k{1}).(e{1})(L(:))) = zeros(1,c);
        elseif ( cs(e_s,k_s) )
            Aeq(roweq,f.(k{1}).(e{1})(L(:))) = ones(1,c);
        elseif ( cs(e_t,k_t) )
            Aeq(roweq,f.(k{1}).(e{1})(L(:))) = -ones(1,c);
        end
    end
    beq(roweq) = 0;
    roweq = roweq+1;
end

%% Define the Demand Always Met Constraint.
%
% Clearly, the supply of each commodity must equal or exceed the
% commodity's demand.
%
% <latex>
% \begin{align*}
% d_k & \le \sum_{w\in V,l \in L} f_{k,(s_k,w),l} \\
%     & : \forall k \in K
% \end{align*}
% </latex>

% Define the Demand Always Met Constraint.
for k=K
    k_s = K_s.(k{1});
    k_d = K_d.(k{1});
    for e=E
        e_s = E_s.(e{1});
        if ( cs(e_s,k_s) )
            Ain(rowin,f.(k{1}).(e{1})(L(:))) = -ones(1,c);
        end
    end
    bin(rowin) = -k_d;
    rowin = rowin+1;
end

%% Solve with MATLAB's intlinprog.
%
% Optimize for the flows and objective variable with MATLAB's function for
% performing MILP.

% Solve for the variables with MILP.
[x,fval,exitflag,output] = intlinprog(o,intcont,Ain,bin,Aeq,beq,bl,bu);

%% Display the Table Results.

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

% Display the results.
output
T

%% Generate the Visualized Graph results.

% Declarations.
fn = 'temp.gif';

% Add an index value to the commodity tuple.
k_q = 1;
for k=K
    K_q.(k{1}) = k_q;
    k_q = k_q+1;
end

% Generate figure for each slot.
fsi = [];
for l=L
    
    % Create figure and grab necessary values.
    fh = figure;
    ph = plot(G);
    pxlim = fh.CurrentAxes.XAxis.Limits;
    pylim = fh.CurrentAxes.YAxis.Limits;
    title(['Digraph G for Time Slot ',num2str(l)]);
    
    % Highlight the flow of each commodity during the current time slot.
    for k=K
        k_q = K_q.(k{1});
        hls = {};
        for e=E
            if (x(f.(k{1}).(e{1})(l)))
                hls{end+1} = E_s.(e{1});
                hls{end+1} = E_t.(e{1});
            end
        end
        col = [0.5,k_q/numel(K),0];
        highlight(ph,hls,'EdgeColor',col,'LineWidth',2);
        
        % Add legend.
        text( ...
            pxlim(1), ...
            pylim(2)*.9-((pylim(2)-pylim(1))*(k_q-1)/numel(K)), ...
            ['---' k{1}],'Color',col);
    end
    axis square;
    
    % Update figure;
    drawnow
    if isempty(fsi)
        fsi = fh.Position;
    else
        fh.Position = fsi;
        refresh;
        drawnow;
    end
    
    % Build gid for animation.
    im = frame2im(getframe(fh));
    [imind,cm] = rgb2ind(im,256);
    if ( l==1 )
        imwrite(imind,cm,fn,'gif','Loopcount',inf);
    else
        imwrite(imind,cm,fn,'gif','WriteMode','append');
    end
end