andrewandrepowell/matlab_mcfp_cf_tdm_lp.m

## 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
	%% 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