# MatPlanWDM v0.5

### Pablo Pavon Mariño (view profile)

29 Jan 2007 (Updated )

Educational network planning tool for the RWA problem in WDM networks (MILP and heuristic based)

MILP.m
```% MILP
%
%>> Usage: [exitMsg exitFlag netState valueFunction timeDefinition timeOptimization] = ...
%
%>> Abstract: This function and all its auxiliary functions in the
%   directory MILP solve the Virtual Topology Design Problem formulated as
%   Mixed-Integer Linear Programming (MILP) Problem. It uses the TOMLAB/CPLEX
%   software to solve the mathematical programming problem. The inputs to the
%   problem include the traffic matrix, the physical topology and the MILP
%   parameters. The objective is to design a virtual topology by optimizing
%   a certain cost function subject to a set of constraints. The values of
%   the output variables of the MILP formulation will yield the OPTIMAL
%   VIRTUAL TOPOLOGY.
%
% 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.
%
%    "Structure of phys variable" from Help.
%
%   parametersString: String of parameters defined by the user. The syntax
%    of this string consist of:
%
%    parameter1Name = parameter1Value (, parameter2Name = parameter2Value, ...)
%
%    The number of white spaces is not checked.
%
%    The default parameters are defines so:
%
%    {['Parameter Name'],   ['Default Value'], ['Data type'],['ValidRange']}
%
%    {['allowWavelengthConversion'],  ['true(1)'], ['boolean'],     [];...
%    ['allowTrafficLosses'],        ['false(1)'], ['boolean'],     []; ...
%    ['cost_GbpsLoss'],                  ['100'], ['numeric'], ['(0,inf)']; ...
%    ['cost_GbpsElectronicallySwitched'], ['10'], ['numeric'], ['[0,inf)']; ...
%    ['cost_opticalTxPlusRx'],             ['0'], ['numeric'], ['[0,inf)']; ...
%    ['cost_virtualPerUsedWavelength'],    ['0'], ['numeric'], ['[0,inf]']; ...
%    ['maximumAllowedUtilization'],        ['1'], ['numeric'], ['[0,1]']; ...
%    ['allowMultiLp'],               ['true(1)'], ['boolean'],     [];...
%    ['applyMaximumDistancePerLp'], ['false(1)'], ['boolean'],     [];...
%    ['maximumDistancePerLp'],           ['800'], ['numeric'], ['(0,inf)'];...
%    ['applyMaximumPhysicalHopsPerLp'], ['false(1)'], ['boolean'], [];...
%    ['maximumPhysicalHopsPerLp'],         ['6'], ['integer'],     ['[1,inf)'];
%    };
%
%
% o Out:
%  . exitMsg: Exit Message.
%           - If the virtual topology design was successful:
%               exitMsg='OK!!!'(exitFlag == 0);
%           - If is not possible to find a feasible solution:
%               exitMsg = 'No feasible solution found' (exitFlag == 1)
%           - Otherwise, exitMsg indicates other sitiation (exitFlag >= 2).
%             Exit Message according to the exit flag from the CPLEX solver.
%
%  . exitFlag:
%          - 0, if it is possible to find a optimal solution
%          - 1, if it is not possible to find a feasible solution
%          - 2, otherwise (a subpotimal solution)
%
%     section "Structure of netState variable" from Help.
%
%   . valueFunction(1x1): Objective function value at optimum. It is a
%     real value.
%
%    timeDefinition(1x1): Elapsed time in hours, minutes and seconds for
%     the mathematical definition of the MILP problem. It is the time
%     needed by MATLAB so as to build the matrices and vectors which
%     implement the problem constraints and the objective function.
%
%    timeOptimization(1x1): Elapsed time in hours, minutes and seconds for
%     the optimization of the MILP problem. It is the time needed by
%     TOMLAB/CPLEX so as to solve the MILP problem.
%
function [exitMsg exitFlag netState valueFunction timeDefinition timeOptimization] = ...

netState = struct ('flowRoutingMatrix' , [] ,  'lightpathRoutingMatrix' , [] , 'lightpathTable' , [] , 'flowTable' , [] , 'numberOfOccupiedTWCs' , [],'numberOfOccupiedTxs' , [],'numberOfOccupiedRxs' , []);

%We define the default parameter values
cellOfDefaultParameters = {['allowWavelengthConversion'],  ['true(1)'], ['boolean'],     [];...
['allowTrafficLosses'],        ['false(1)'], ['boolean'],     []; ...
['cost_GbpsLoss'],                  ['100'], ['numeric'], ['(0,inf)']; ...
['cost_GbpsElectronicallySwitched'], ['10'], ['numeric'], ['[0,inf)']; ...
['cost_opticalTxPlusRx'],             ['0'], ['numeric'], ['[0,inf)']; ...
['cost_virtualPerUsedWavelength'],    ['0'], ['numeric'], ['[0,inf]']; ...
['maximumAllowedUtilization'],        ['1'], ['numeric'], ['[0,1]']; ...
['allowMultiLp'],               ['true(1)'], ['boolean'],     [];...
['applyMaximumDistancePerLp'], ['false(1)'], ['boolean'],     [];...
['maximumDistancePerLp'],           ['800'], ['numeric'], ['(0,inf)'];...
['applyMaximumPhysicalHopsPerLp'], ['false(1)'], ['boolean'], [];...
['maximumPhysicalHopsPerLp'],         ['6'], ['integer'],     ['[1,inf)'];
};

%We process the string of the algorithm parameters
[exitFlagOfProcessParam exitMsgOfProcessParam userParam] = ...
if (exitFlagOfProcessParam ~= 0), error(['>>> Bad Algorithm Parameters:', exitMsgOfProcessParam]); end

%We Initalize the MILP variables
initializationMILPVariablesScript;

constraintsMatrix = sparse(0,NUMBERVARIABLES);
bVector = zeros (0,0);

fprintf('\n\n-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->');
fprintf('\nDefining the Constraints and the Objective Function %d.');

fprintf('\n');
fprintf('-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->\n\n');
fprintf('\n... ');

tDef0 = clock;

%********************1. We define the PROBLEM CONSTRAINTS *****************
%**************************************************************************

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% CONSTRAINTs 0: NOT ALLOWED VALUES FOR DECISION VARIABLES AND
%%%%% NO LOOPS IN ORIGIN AND DESTINATION NODES FOR LIGHTPATHS AND FLOWS

bValue = 0;
constraintType = 0; % -1 <= ; 0 == ; +1 >=
pijmw_constraintMatrix = zeros (N , N , NUMBERLINKWAVELENGTHS);
fijsd_constraintMatrix = zeros (N , N , N , N);  % (i,j,s,d)
for i=1:N
pijmw_constraintMatrix (i,:,initial_mwPosition:end_mwPosition) = 1;
end

pijmw_constraintMatrix (:,i,initial_mwPosition:end_mwPosition) = 1;
end
pijmw_constraintMatrix (i,i,:) = 1;
fijsd_constraintMatrix (i,i,:,:) = 1;
fijsd_constraintMatrix (:,:,i,i) = 1;
fijsd_constraintMatrix (:,i,i,:) = 1;
fijsd_constraintMatrix (i,:,:,i) = 1;
end
[constraintsMatrix , bVector] = addConstraint (constraintsMatrix , bVector , pijmw_constraintMatrix , fijsd_constraintMatrix , bValue , constraintType , NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD);

fprintf('\nConstraints Type 0 Defined\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% LIGHTPATH ROUTING AND WAVELENGTH CONSTRAINTS %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% CONSTRAINTs 1: VIRTUAL INDEGREE
for node = 1:N
pijmw_constraintMatrix = zeros (N , N , NUMBERLINKWAVELENGTHS);
bValue = phys.numberRxPerNode (node);
constraintType = -1; % -1 <= ; 0 == ; +1 >=

pijmw_constraintMatrix(:, node , initial_mwPosition:end_mwPosition) = 1;
end

[constraintsMatrix , bVector] = addConstraint (constraintsMatrix , bVector , pijmw_constraintMatrix , [] , bValue , constraintType , NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD);
end

fprintf('Constraints Type 1 Defined\n');

%%%%% CONSTRAINTs 2: VIRTUAL OUTDEGREE
for node = 1:N
pijmw_constraintMatrix = zeros (N , N , NUMBERLINKWAVELENGTHS);
bValue = phys.numberTxPerNode (node);
constraintType = -1;     % -1 <= ; 0 == ; +1 >=

pijmw_constraintMatrix (node , : , initial_mwPosition:end_mwPosition) = 1;
end

[constraintsMatrix , bVector] = addConstraint (constraintsMatrix , bVector , pijmw_constraintMatrix , [] , bValue , constraintType , NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD);
end

fprintf('Constraints Type 2 Defined\n');

%%%%% CONSTRAINTs 3: NO WAVELENGTH IN A LINK CARRIES MORE THAN ONE LP
pijmw_constraintMatrix = zeros (N , N , NUMBERLINKWAVELENGTHS);
bValue = 1;
constraintType = -1;     % -1 <= ; 0 == ; +1 >=
pijmw_constraintMatrix (:,:,mw) = 1;

[constraintsMatrix , bVector] = addConstraint (constraintsMatrix , bVector , pijmw_constraintMatrix , [] , bValue , constraintType , NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD);
end

fprintf('Constraints Type 3 Defined\n');

%%%%% CONSTRAINTs 4: CONSTRAINT OF MAX 1 LIGHPTATH BETWEEN EACH I,J PAIR
if not(userParam.allowMultiLp),
for i = 1:N,
for j = 1:N,
if (i == j)
continue;
end
pijmw_constraintMatrix = zeros (N , N , NUMBERLINKWAVELENGTHS);
bValue = 1;
constraintType = -1;     % -1 <= ; 0 == ; +1 >=

pijmw_constraintMatrix (i , j , initial_mwPosition:end_mwPosition) = 1;
end

[constraintsMatrix , bVector] = addConstraint (constraintsMatrix , bVector , pijmw_constraintMatrix , [] , bValue , constraintType , NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD);
end
end
end

fprintf('Constraints Type 4 Defined\n');

%%%%% CONSTRAINT 5: CONSTRAINT OF MAX NUMBER OF PHYSICAL HOPS PER
if not(userParam.allowMultiLp) & (userParam.applyMaximumPhysicalHopsPerLp),
for i = 1:N,
for j = 1:N,
if (i == j)
continue;
end
pijmw_constraintMatrix = zeros (N , N , NUMBERLINKWAVELENGTHS);
bValue = userParam.maximumPhysicalHopsPerLp;
constraintType = -1;     % -1 <= ; 0 == ; +1 >=

pijmw_constraintMatrix (i , j , :) = 1;

[constraintsMatrix , bVector] = addConstraint (constraintsMatrix , bVector , pijmw_constraintMatrix , [] , bValue , constraintType , NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD);
end
end
end

fprintf('Constraints Type 5 Defined\n');

%%%%% CONSTRAINT 6: CONSTRAINT OF MAX DISTANCE IN KM OF PHYSICAL HOPS
%%%%% OF LP
if not(userParam.allowMultiLp) & (userParam.applyMaximumDistancePerLp),
for i = 1:N,
for j = 1:N,
if (i == j)
continue;
end
pijmw_constraintMatrix = zeros (N , N , NUMBERLINKWAVELENGTHS);
bValue = userParam.maximumDistancePerLp;
constraintType = -1;     % -1 <= ; 0 == ; +1 >=

for m = 1:M
pijmw_constraintMatrix (i , j , initial_mwPosition:end_mwPosition) = phys.linkLengthInKm (m);
end

[constraintsMatrix , bVector] = addConstraint (constraintsMatrix , bVector , pijmw_constraintMatrix , [] , bValue , constraintType , NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD);
end
end
end

fprintf('Constraints Type 6 Defined\n');

%%%%% CONSTRAINT 7: ROUTING OF LIGHTPATHS WITH/WITHOUT WAVELENGTH CONTINUITY
%%%%% CONSTRAINT
bValue = 0;
if (not(userParam.allowWavelengthConversion))
for lpConservationNode = 1:N

for sourceNodeOfLpAggregate = 1:N
for endNodeOfLpAggregate = 1:N
if (sourceNodeOfLpAggregate == endNodeOfLpAggregate)
continue;
end
if (sourceNodeOfLpAggregate == lpConservationNode)
constraintType = -1;
elseif (endNodeOfLpAggregate == lpConservationNode)
constraintType = +1;
else
constraintType = 0;
end

for w = 1:Wmax_thisNode % the number of wavelengths of the fibre with more number of wavelengths
pijmw_constraintMatrix = zeros (N , N , NUMBERLINKWAVELENGTHS);

pijmw_constraintMatrix (sourceNodeOfLpAggregate , endNodeOfLpAggregate , MWPOSITIONSOFTHENODE_INCOMING{lpConservationNode,w}) = 1;
pijmw_constraintMatrix (sourceNodeOfLpAggregate , endNodeOfLpAggregate , MWPOSITIONSOFTHENODE_OUTGOING{lpConservationNode,w}) = -1;

[constraintsMatrix , bVector] = addConstraint (constraintsMatrix , bVector , pijmw_constraintMatrix , [] , bValue , constraintType , NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD);
end
end
end
end
else
bValue = 0;
for lpConservationNode = 1:N

for sourceNodeOfLpAggregate = 1:N
for endNodeOfLpAggregate = 1:N
if (sourceNodeOfLpAggregate == endNodeOfLpAggregate)
continue;
end

pijmw_constraintMatrix = zeros (N , N , NUMBERLINKWAVELENGTHS);

if (sourceNodeOfLpAggregate == lpConservationNode)
constraintType = -1;
elseif (endNodeOfLpAggregate == lpConservationNode)
constraintType = +1;
else
constraintType = 0;
end

pijmw_constraintMatrix (sourceNodeOfLpAggregate , endNodeOfLpAggregate , initial_mwPosition:end_mwPosition) = 1;
end

pijmw_constraintMatrix (sourceNodeOfLpAggregate , endNodeOfLpAggregate , initial_mwPosition:end_mwPosition) = -1;
end

[constraintsMatrix , bVector] = addConstraint (constraintsMatrix , bVector , pijmw_constraintMatrix , [] , bValue , constraintType , NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD);
end
end
end
end

fprintf('Constraints Type 7 Defined\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% TRAFFIC FLOW ROUTING CONSTRAINTS %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% CONSTRAINT 8: LP CAPACITY CONSTRAINTS
for sourceNodeOfLpAggregate = 1:N
for endNodeOfLpAggregate = 1:N
if (sourceNodeOfLpAggregate == endNodeOfLpAggregate)
continue;
end

pijmw_constraintMatrix = zeros (N , N , NUMBERLINKWAVELENGTHS);
fijsd_constraintMatrix = zeros (N , N , N , N);    % (i,j,s,d)
bValue = 0;
constraintType = -1;     % -1 <= ; 0 == ; +1 >=

fijsd_constraintMatrix (sourceNodeOfLpAggregate,endNodeOfLpAggregate,:,:)= 1;

pijmw_constraintMatrix (sourceNodeOfLpAggregate , endNodeOfLpAggregate , initial_mwPosition:end_mwPosition) = -userParam.maximumAllowedUtilization*phys.lightpathCapacity;
end

[constraintsMatrix , bVector] = addConstraint (constraintsMatrix , bVector , pijmw_constraintMatrix , fijsd_constraintMatrix , bValue , constraintType , NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD);
end
end

fprintf('Constraints Type 8 Defined\n');

%%%%% CONSTRAINT 9: CONSERVATION FLOW CONSTRAINTS
for conservationFlowNode = 1:N
for ingressNodeOfFlow = 1:N
for egressNodeOfFlow = 1:N
if (ingressNodeOfFlow == egressNodeOfFlow)
continue;
end

fijsd_constraintMatrix = zeros (N , N , N , N);    % (i,j,s,d)

if (ingressNodeOfFlow == conservationFlowNode)
if (userParam.allowTrafficLosses)
constraintType = -1;
else
constraintType = 0;
end
bValue = traff_trafficMatrix (ingressNodeOfFlow,egressNodeOfFlow);
elseif (egressNodeOfFlow == conservationFlowNode)
if (userParam.allowTrafficLosses)
constraintType = +1;
else
constraintType = 0;
end
bValue = - traff_trafficMatrix (ingressNodeOfFlow,egressNodeOfFlow);
else
constraintType = 0;
bValue = 0;
end

fijsd_constraintMatrix (conservationFlowNode,:,ingressNodeOfFlow,egressNodeOfFlow)= 1;
fijsd_constraintMatrix (:,conservationFlowNode,ingressNodeOfFlow,egressNodeOfFlow)= -1;
fijsd_constraintMatrix (conservationFlowNode,conservationFlowNode,ingressNodeOfFlow,egressNodeOfFlow)= 0;

[constraintsMatrix , bVector] = addConstraint (constraintsMatrix , bVector , [] , fijsd_constraintMatrix , bValue , constraintType , NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD);
end
end
end

fprintf('Constraints Type 9 Defined\n');

%*****************2. We define the OBJECTIVE FUNCTION**********************
%**************************************************************************

%% first, the contribution of electronic costs
fijsd_objectiveMatrix = userParam.cost_GbpsElectronicallySwitched * ones (N , N , N , N);

%% now, the contribution of virtual costs
pijmw_objectiveMatrix = userParam.cost_virtualPerUsedWavelength * ones (N , N , NUMBERLINKWAVELENGTHS);

%% now, the contribution of optical costs
for i=1:N
pijmw_objectiveMatrix (i,:,initial_mwPosition:end_mwPosition) = ...
pijmw_objectiveMatrix (i,:,initial_mwPosition:end_mwPosition) + userParam.cost_opticalTxPlusRx;
end
end

%% now the contribution of the cost of loosing traffic
if (userParam.allowTrafficLosses)
totalTrafficOffered = sum (sum (traff_trafficMatrix));
routedTraffic = 0;
for i_s=1:N
fijsd_objectiveMatrix (i_s,:,i_s,:) = fijsd_objectiveMatrix (i_s,:,i_s,:) - userParam.cost_GbpsLoss;
end
end

functionToOptimize = [array2vect(pijmw_objectiveMatrix) array2vect(fijsd_objectiveMatrix)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tDef1 = clock;
timeDefinition=etime(tDef1,tDef0);
[hDef, minDef, secDef]=sec2time(timeDefinition);
fprintf(['\n\n Elapsed time to define the problem: ',num2str(hDef),' hours ',num2str(minDef),' minutes ',num2str(secDef),' seconds \n']);

%-----------------------------********************3. TOMLAB CPLEX OPTMIZATION  ***********************---------------------------------------------------------------
%----------------------------------------------------------------------------------------------------------------------------------------------------

fprintf('\n\n-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->');
fprintf('\nTOMLAB / Cplex solving Problem %d.');
fprintf('\n');
fprintf('-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->-->\n\n');
fprintf('\n... ');

[exitMsg exitFlag x, slack, v, rc, valueFunction, ninf, sinf, basis, lpiter, ...
glnodes, confstat, iconfstat, sa] = ...

tSolve1 = clock;
timeOptimization=etime(tSolve1, tDef1);
[hOpt, minOpt, secOpt]=sec2time(timeOptimization);
fprintf(['\n\n Elapsed time to solve the problem: ',num2str(hOpt),' hours ',num2str(minOpt),' minutes ',num2str(secOpt),' seconds \n']);
fprintf('\n exitFlag = %d\n exitMsg = %s\n' ,exitFlag,exitMsg);

if exitFlag == 1,
netState.lightpathTable =[];
netState.lightpathRoutingMatrix=[];
netState.flowTable=[];
netState.flowRoutingMatrix=[];
valueFunction =[];
timeDefinition = [];
timeOptimization =[];
return;
end

%-----------------------********************4. We capture the solver CPLEX OUTPUTS****************************-----------------------------------------
%-----------------------------------------------------------------------------------------------------------------------------------------------------

% We capture into a programming variable the information about the
% Lightpath-Wavelength-Link Indicator contained in the solution vector x
% obtained by the solver.

pijmw = x (1:NUMBERLINKWAVELENGTHS * NUMBERINGRESSEGRESSNODEPAIRS);
fijsd = x (NUMBERLINKWAVELENGTHS * NUMBERINGRESSEGRESSNODEPAIRS + 1 :  NUMBERVARIABLES);

[pijmw fijsd] = roundSolutionVector (pijmw , fijsd);

if sum(pijmw)==0 | sum(fijsd)==0,
exitFlag = 1;
exitMsg = 'Null optimal integer solution found';
netState.lightpathTable =[];
netState.lightpathRoutingMatrix=[];
netState.flowTable=[];
netState.flowRoutingMatrix=[];
return
end

%We translate the variables which represent the configuration of the
%virtual topology from the format used in the MILP optimization into the
%standard format.
[netState.lightpathTable netState.lightpathRoutingMatrix] = trans_pijmw2LPRoutingMatrix(pijmw, phys);
[netState.flowTable netState.flowRoutingMatrix] = trans_fijsd2flowRoutingMatrix(fijsd, netState.lightpathTable , netState.lightpathRoutingMatrix , phys, traff_trafficMatrix);
netState.flowTable=[netState.flowTable -1*ones(size(netState.flowTable,1),3)];
netState.lightpathTable =[(1:size(netState.lightpathTable,1))' netState.lightpathTable];

numberOfNodes=length(phys.numberTxPerNode);
for nodei=1:numberOfNodes
netState.numberOfOccupiedTxs(nodei) = length(find(netState.lightpathTable(:,2)==nodei));
netState.numberOfOccupiedRxs(nodei) = length(find(netState.lightpathTable(:,3)==nodei));
end
numberOfLightpaths=size(netState.lightpathTable,1);
netState.numberOfOccupiedTWCs=zeros(numberOfNodes,1);
for lpID=1:numberOfLightpaths,
[sequenceOfNodeIDs, sequenceOfFiberIDs, sequenceOfWavelengthIDs] = lightpathPathComputation (lpID, phys, netState);
for i=1:length(sequenceOfFiberIDs)-1,
if  sequenceOfWavelengthIDs(i) ~= sequenceOfWavelengthIDs(i+1),
netState.numberOfOccupiedTWCs (sequenceOfNodeIDs(i+1)) = netState.numberOfOccupiedTWCs (sequenceOfNodeIDs(i+1)) + 1;
end
end
end

end

function [pijmw_result fijsd_result] = roundSolutionVector (pijmw , fijsd)
pijmw_result = pijmw;
fijsd_result = fijsd;

roundLimit = (max (fijsd)) / 1E10;
for i=1:length(fijsd),
if (fijsd(i) < roundLimit),
fijsd_result(i) = 0;
end
end

precision = 1E-1;
for i=1:length(pijmw),
if (abs (pijmw(i)-0) < precision),
pijmw_result(i) = 0;
elseif (abs (pijmw(i)-1) < precision),
pijmw_result(i) = 1;
else
disp ('CPLEX precision error in Pijmw variables');
error ('CPLEX precision error in Pijmw variables');
end
end
end
```