%%% libraryFR_optimalFlowAssignment
%
% Usage: [exitFlag flowTable flowRoutingMatrix cost] =
% libraryFR_optimalFlowAssignment(traff_trafficMatrix, linkTable , linkCapacities , linkCostPerGbps)
%
% Abstract: This function solves the Traffic Flow Routing over a network
% topology. It implements an optimal Flow Assignment problem.
% The problem is established as linear programming problem,
% and solved with the linprog solver
% The function provides the solution which minimizes a given cost function.
% The cost function employed is the sum of the costs in each link.
% The cost in each link is give by the product of the Gbps
% carried and cost per Gbps in the link.
%
% Arguments:
% o In:
% traff_trafficMatrix(NxN): Average traffic flow offered between node
% pairs. The Traffic Matrix is a two-dimensional matrix with N (N:
% number of nodes) rows and N columns. An entry(s,d) means the average
% traffic flow from node 's' to node 'd', expressed in Gbps. The main
% diagonal is full of 0s.
%
% . linkTable(M,2): M-by-2 integer matrix. Each row is a link in the graph.
% first column is the ingress node (1...N), second the egress node
% (1...N)
%
% . linkCapacities(M,1): Capacity in Gbps of each link. This is the
% maximum traffic that the link can carry.
%
% . linkCostPerGbps(M,1): Cost of carrying 1 Gbps in each link. If the
% vector is constant, the average number of hops is minimized. If the
% vector contains a measure of link distance, the average distance traversed
% by the carried traffic is minimized
%
% o Out:
% . exitFlag:
% 0, if LINPROG converged to a solution X.(It means that all the flows
% can be routed)
% 1, if maximum number of iterations reached.
% 2, if the traffic flow routing is not FEASIBLE.
%
% . flowTable(F,3): F-by-3 matrix. First column de ingress node, second
% the egress node, third the traffic offered in Gbps.
%
% . flowRoutingMatrix (F,L): F-by-L integer matrix where F is the number of
% flows and L is the number of lightpaths. Each row is a flow 'f' and each
% column is a lightpath 'l'. If a flow 'f' uses a lightpath 'l', the
% entry (f,l) is equal to the value of the flow 'f' carried by the
% lightpath 'f'. If no lightpath is used by the flow 'f', the entry is
% equal to '0'. If a row is a row of zeros, the the traffic flows was
% not routed
%
% cost: value of the cost function at optimum
function [exitFlag flowTable flowRoutingMatrix cost] = libraryFR_optimalFlowAssignment(traff_trafficMatrix, linkTable , linkCapacities , linkCostPerGbps, varargin)
NUMBERNODES = size (traff_trafficMatrix,1);
NUMBERLINKS = size (linkTable,1);
% Create the flowTable from the trafficMatrix
numberOfFlows = nnz(traff_trafficMatrix);
flowTable = zeros (numberOfFlows,3);
contFlow = 0;
for ingressNode=1:NUMBERNODES
for egressNode=1:NUMBERNODES
if (ingressNode ~= egressNode) && (traff_trafficMatrix(ingressNode,egressNode) ~= 0),
contFlow = contFlow + 1;
flowTable (contFlow,:) = [ingressNode egressNode traff_trafficMatrix(ingressNode,egressNode)];
end
end
end
NUMBERFLOWS = contFlow;
NUMBERDECISIONVARIABLES = NUMBERFLOWS * NUMBERLINKS;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Non-negativity constraints
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = -eye ( NUMBERDECISIONVARIABLES);
b = zeros ( NUMBERDECISIONVARIABLES , 1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Flow conservation constraints (one per flow and node)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for nodeID=1:NUMBERNODES
% inputLinks are the incoming edges
inputLinks = find (linkTable (:,2) == nodeID);
% outputLinks are the outgoing edges
outputLinks = find (linkTable (:,1) == nodeID);
for flowId=1:NUMBERFLOWS
flowOriginID = flowTable (flowId , 1);
flowDestID = flowTable (flowId , 2);
flowGbps = flowTable (flowId , 3);
constraint = zeros (1 , NUMBERDECISIONVARIABLES);
for contOutputLink=1:length (outputLinks)
linkId = outputLinks (contOutputLink);
columnA = columnMatrixA (linkId , flowId , NUMBERFLOWS);
constraint (columnA) = 1;
end
for contInputLink=1:length (inputLinks)
linkId = inputLinks (contInputLink);
columnA = columnMatrixA (linkId , flowId , NUMBERFLOWS);
constraint (columnA) = -1;
end
A = [A ; constraint ; -constraint];
if (nodeID == flowOriginID)
b = [b ; flowGbps ; -flowGbps]; % node is ingress of the flow
elseif (nodeID == flowDestID)
b = [b ; -flowGbps ; flowGbps]; % node is ingress of the flow
else
b = [b ; 0 ; 0]; % node is intermediate to the flow
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constraint carried traffic in each link limited by link capacity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for linkId = 1:NUMBERLINKS
constraint = zeros (1, NUMBERDECISIONVARIABLES);
for flowId=1:NUMBERFLOWS
columnA = columnMatrixA (linkId , flowId , NUMBERFLOWS);
constraint (columnA) = 1;
end
A = [A ; constraint];
b = [b ; linkCapacities(linkId)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cost function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
costFunction = zeros (1, NUMBERDECISIONVARIABLES);
for linkId = 1:NUMBERLINKS
for flowId=1:NUMBERFLOWS
columnA = columnMatrixA (linkId , flowId , NUMBERFLOWS);
costFunction (columnA) = linkCostPerGbps (linkId);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LINPROG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
options=optimset('display','off');
[flowRoutingVector,cost,exitFlagLP] = linprog (costFunction',A,b,[],[],[],[],[],options);
% . exitFlagLP: The value returned by the linprog solver
% 1 LINPROG converged to a solution X. (Means the flows can be carried)
% 0 Maximum number of iterations reached.
% -2 No feasible point found.
% -3 Problem is unbounded.
% -4 NaN value encountered during execution of algorithm.
% -5 Both primal and dual problems are infeasible.
% -7 Search direction became too small; no further progress can be
% made.
switch exitFlagLP
case 1
exitFlag = 0;
case 0
exitFlag = 1;
case -2
exitFlag = 2;
flowTable =[];
flowRoutingMatrix=[];
return
case -3
exitFlag = 2;
flowTable =[];
flowRoutingMatrix=[];
return
case -4
exitFlag=2;
flowTable =[];
flowRoutingMatrix=[];
return
case -5
exitFlag=2;
flowTable =[];
flowRoutingMatrix=[];
return
case -7
exitFlag=2;
flowTable =[];
flowRoutingMatrix=[];
return
end
for cont=1:length(flowRoutingVector)
if (flowRoutingVector(cont)<1E-6)
flowRoutingVector(cont) = 0;
end
end
flowRoutingMatrix = zeros (NUMBERFLOWS , NUMBERLINKS);
for flowId=1:NUMBERFLOWS
for linkId=1:NUMBERLINKS
flowRoutingMatrix(flowId,linkId) = flowRoutingVector (columnMatrixA (linkId , flowId , NUMBERFLOWS));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Receives a flowId (1,2,...) and a linkId (1,2,...) and return the index
% of the column in the constraints matrix for the decision variable
% associated to that flow and link
function columnID = columnMatrixA (linkId , flowId , NUMBERFLOWS)
columnID = flowId + NUMBERFLOWS * (linkId-1);