No BSD License  

Highlights from
MatPlanWDM v0.5

image thumbnail

MatPlanWDM v0.5

by

 

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] = ...
%     MILP(traff_trafficMatrix, phys, parametersString)
% 
%>> 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.
%
%  . phys: Phys Structure. More information about netState in section 
%    "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)
% 
%   . netState: NetState Structure. More information about netState in
%     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] = ...
    MILP(traff_trafficMatrix, phys, parametersString)

    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] = ...
        processAlgorithmParameters (parametersString, cellOfDefaultParameters);   
    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
        incomingLinks = phys.incomingLinksToNode {i};
        for m=incomingLinks
            initial_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,1);
            end_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,2);
            pijmw_constraintMatrix (i,:,initial_mwPosition:end_mwPosition) = 1;
        end

        outgoingLinks = phys.outgoingLinksFromNode {i};
        for m=outgoingLinks
            initial_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,1);
            end_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,2);
            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 >=

        incomingLinks=phys.incomingLinksToNode {node}; % row vector of incoming links to the node    
        for m=incomingLinks
            initial_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,1);
            end_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,2);
            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 >=

        outgoingLinks = phys.outgoingLinksFromNode {node};
        for m=outgoingLinks
            initial_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,1);
            end_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,2);
            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
    for mw=1:NUMBERLINKWAVELENGTHS
            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,
            outgoingLinks = phys.outgoingLinksFromNode {i};
            for j = 1:N,
                if (i == j)
                    continue;
                end
                pijmw_constraintMatrix = zeros (N , N , NUMBERLINKWAVELENGTHS);
                bValue = 1;
                constraintType = -1;     % -1 <= ; 0 == ; +1 >=

                for m = outgoingLinks
                    initial_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,1);
                    end_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,2);
                    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
                    initial_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,1);
                    end_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,2);
                    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

            incomingLinks = phys.incomingLinksToNode {lpConservationNode};
            outgoingLinks = phys.outgoingLinksFromNode {lpConservationNode};
            nodeInputOrOutputLinks = union (incomingLinks,outgoingLinks);
            Wmax_thisNode = max(phys.numberWavelengthPerFiber (nodeInputOrOutputLinks));

            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

            incomingLinks = phys.incomingLinksToNode {lpConservationNode};
            outgoingLinks = phys.outgoingLinksFromNode {lpConservationNode};
            nodeInputOrOutputLinks = union (incomingLinks,outgoingLinks);

            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

                    for m = incomingLinks
                        initial_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,1);
                        end_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,2);
                        pijmw_constraintMatrix (sourceNodeOfLpAggregate , endNodeOfLpAggregate , initial_mwPosition:end_mwPosition) = 1;
                    end

                    for m = outgoingLinks
                        initial_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,1);
                        end_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,2);

                        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;

            outgoingLinks = phys.outgoingLinksFromNode{sourceNodeOfLpAggregate};
            for m=outgoingLinks
                initial_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,1);
                end_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,2);
                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
        for m=phys.outgoingLinksFromNode {i}
            initial_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,1);
            end_mwPosition = MWINITIALANDENDPOSITIONOFLINK (m,2);
            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(' Definition Ready');
    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] = ...
        callSolver_CPLEX (functionToOptimize,constraintsMatrix,bVector,N,NUMBERLINKWAVELENGTHS,NUMBERINGRESSEGRESSNODEPAIRS);

    tSolve1 = clock;
    timeOptimization=etime(tSolve1, tDef1);
    [hOpt, minOpt, secOpt]=sec2time(timeOptimization);
    fprintf('\n Optimization Ready');
    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

Contact us