001/******************************************************************************* 002 * Copyright (c) 2016 Pablo Pavon Mariņo. 003 * All rights reserved. This program and the accompanying materials 004 * are made available under the terms of the GNU Lesser Public License v2.1 005 * which accompanies this distribution, and is available at 006 * http://www.gnu.org/licenses/lgpl.html 007 ******************************************************************************/ 008 009 010 011 012package com.net2plan.examples.ocnbook.offline; 013 014import java.io.File; 015import java.io.PrintWriter; 016import java.util.ArrayList; 017import java.util.Collection; 018import java.util.Collections; 019import java.util.HashMap; 020import java.util.HashSet; 021import java.util.LinkedList; 022import java.util.List; 023import java.util.Map; 024import java.util.Map.Entry; 025import java.util.Random; 026import java.util.Set; 027 028import cern.colt.matrix.tdouble.DoubleFactory1D; 029import cern.colt.matrix.tdouble.DoubleFactory2D; 030import cern.colt.matrix.tdouble.DoubleMatrix1D; 031import cern.colt.matrix.tdouble.DoubleMatrix2D; 032import cern.jet.math.tdouble.DoubleFunctions; 033 034import com.jom.DoubleMatrixND; 035import com.jom.OptimizationProblem; 036import com.net2plan.interfaces.networkDesign.Demand; 037import com.net2plan.interfaces.networkDesign.IAlgorithm; 038import com.net2plan.interfaces.networkDesign.Link; 039import com.net2plan.interfaces.networkDesign.Net2PlanException; 040import com.net2plan.interfaces.networkDesign.NetPlan; 041import com.net2plan.interfaces.networkDesign.NetworkElement; 042import com.net2plan.interfaces.networkDesign.Node; 043import com.net2plan.interfaces.networkDesign.ProtectionSegment; 044import com.net2plan.interfaces.networkDesign.Route; 045import com.net2plan.interfaces.networkDesign.SharedRiskGroup; 046import com.net2plan.utils.Constants.RoutingType; 047import com.net2plan.utils.Constants.SearchType; 048import com.net2plan.utils.DoubleUtils; 049import com.net2plan.utils.InputParameter; 050import com.net2plan.utils.IntUtils; 051import com.net2plan.utils.Pair; 052import com.net2plan.utils.Quadruple; 053import com.net2plan.utils.Triple; 054 055/** 056 * This algorithm is devoted to solve the several network planning problems in an optical WDM network (fiber placement, RWA, under different recovery schemes), appearing in the case study in the book section mentioned below. 057 * @net2plan.description 058 * @net2plan.keywords WDM, Topology assignment (TA), Flow assignment (FA), GRASP, JOM 059 * @net2plan.ocnbooksections Section 12.10 060 * @net2plan.inputParameters 061 * @author Pablo Pavon-Marino 062 */ 063public class Offline_tcfa_wdmPhysicalDesign_graspAndILP implements IAlgorithm 064{ 065 /* Stat */ 066 private double stat_totalCost , stat_costinks , stat_costCirc , stat_costNodes; 067 private int stat_numLinks , stat_numCirc , stat_numNodesDeg2 , stat_numNodesDeg3; 068 private int stat_numItGRASP , stat_numLSIterationsReducingCost , stat_numSolverCalls , stat_numSolverCallsOk; 069 private double stat_totalTimeSecs , stat_totalTimeCreateProgramsSecs , stat_totalTimeSolverSecs; 070 071 /* Parameters */ 072 private DoubleMatrix1D linkCost_e; 073 private Random rng; 074 private NetPlan netPlan; 075 076 /* Control */ 077 private int Efm, N, nSRGs , D , R; 078 private DoubleMatrix1D numCircH_d; 079 private DoubleMatrix2D A_er , A_dr , A_rs , A_se , Aout_ne , Ain_ne, Abid_dd , Abid_rr , Abid_ee; 080 private Map<Demand,Demand> opposite_d; 081 private Map<Route,Route> opposite_r; 082 private int [][] routeList_d; // for each demand, the array of route ids assigned to it 083 private int [][] traversedLinks_r; // for each demand, the array of route ids assigned to it 084 085 086 private InputParameter algorithm_randomSeed = new InputParameter ("algorithm_randomSeed", (long) 1 , "Seed of the random number generator"); 087 private InputParameter algorithm_outputFileNameRoot = new InputParameter ("algorithm_outputFileNameRoot", "tcfaWDM_graspAndILP" , "Root of the file name to be used in the output files. If blank, no output"); 088 private InputParameter algorithm_maxExecutionTimeInSeconds = new InputParameter ("algorithm_maxExecutionTimeInSeconds", (double) 300 , "Algorithm maximum running time in seconds" , 0 , false , Double.MAX_VALUE , true); 089 private InputParameter tcfa_circuitCost = new InputParameter ("tcfa_circuitCost", (double) 3.0 , "Cost of each circuit" , 0 , true , Double.MAX_VALUE , true); 090 private InputParameter tcfa_circuitCapacity_Gbps = new InputParameter ("tcfa_circuitCapacity_Gbps", (double) 10.0 , "Capacity of an optical circuit" , 0 , true , Double.MAX_VALUE , true); 091 private InputParameter tcfa_srgType = new InputParameter ("tcfa_srgType" , "#select# perBidirectionalLinkBundle noFailure", "Determines how the SRGs are initialized. The design must be tolerant to single SRG failures"); 092 private InputParameter tcfa_srgMttfPer1000Km_hrs = new InputParameter ("tcfa_srgMttfPer1000Km_hrs" , (double) 4380 , "Mean Time To Fail (MTTF) of the SRGs defined is this value multiplied by the link length divided by 1000", 0 , false , Double.MAX_VALUE , true); 093 private InputParameter tcfa_srgMttr_hrs = new InputParameter ("tcfa_srgMttr_hrs" , (double) 24 , "Mean Time To Repair (MTTR) of the SRGs", 0 , true , Double.MAX_VALUE , true); 094 private InputParameter tcfa_recoveryType = new InputParameter ("tcfa_recoveryType" , "#select# 1+1 shared restoration", "Determines the type of network recovery mechanism in the network."); 095 private InputParameter tcfa_maxPathLengthInKm = new InputParameter ("tcfa_maxPathLengthInKm" , (double) 10000 , "Maximum length that a lightpath can have", 0 , true , Double.MAX_VALUE , true); 096 private InputParameter tcfa_maxPathNumberOfHops = new InputParameter ("tcfa_maxPathNumberOfHops" , (int) 4 , "Maximum number of hops that a lightpath can make", 1 , Integer.MAX_VALUE); 097 private InputParameter tcfa_maxNumberPathsPerDemand = new InputParameter ("tcfa_maxNumberPathsPerDemand" , (int) 200 , "Maximum number of paths in the candidate path list, for each demand", 1 , Integer.MAX_VALUE); 098 private InputParameter tcfa_linkCapacity_numCirc = new InputParameter ("tcfa_linkCapacity_numCirc" , (int) 80 , "Number of wavelengths available per fiber", 1 , Integer.MAX_VALUE); 099 private InputParameter tcfa_linkCostPerKm = new InputParameter ("tcfa_linkCostPerKm" , (double) 0.4 , "Cost of one km of optical fiber", 0 , true , Double.MAX_VALUE , true); 100 private InputParameter tcfa_costNodeDegree2 = new InputParameter ("tcfa_costNodeDegree2" , (double) 4.0 , "Cost of one OADM of degree lower or equal than two", 0 , true , Double.MAX_VALUE , true); 101 private InputParameter tcfa_costNodeDegreeHigherThan2 = new InputParameter ("tcfa_costNodeDegreeHigherThan2" , (double) 4.0 , "Cost of one OADM of degree higher than two", 0 , true , Double.MAX_VALUE , true); 102 private InputParameter tcfa_networkDiameterKmToNormalize = new InputParameter ("tcfa_networkDiameterKmToNormalize" , (double) 4500 , "The link lengths are normalized so the longest link has this value", 0 , true , Double.MAX_VALUE , true); 103 private InputParameter tcfa_maxNumIterations = new InputParameter ("tcfa_maxNumIterations", (int) 10000 , "Maximum number of iterations" , 1 , Integer.MAX_VALUE); 104 private InputParameter algorithm_solverName = new InputParameter ("solverName", "#select# glpk cplex ipopt", "The solver name to be used by JOM. GLPK and IPOPT are free, CPLEX commercial. GLPK and CPLEX solve linear problems w/w.o integer contraints. IPOPT is can solve nonlinear problems (if convex, returns global optimum), but cannot handle integer constraints"); 105 private InputParameter algorithm_solverLibraryName = new InputParameter ("solverLibraryName", "", "The solver library full or relative path, to be used by JOM. Leave blank to use JOM default."); 106 private InputParameter algorithm_maxSolverTimeInSeconds = new InputParameter ("maxSolverTimeInSeconds", (double) 2 , "Maximum time granted to the solver in to solve each problem instance. If this time expires, the solver returns the best solution found so far (if a feasible solution is found)"); 107 private InputParameter tcfa_bidirectionalLinks = new InputParameter ("tcfa_bidirectionalLinks", true , "If true, the topology of fibers deployed is constrained to be bidirectional (the same number of fibers in both directions)"); 108 109 @Override 110 public String executeAlgorithm(NetPlan netPlan, Map<String, String> algorithmParameters, Map<String, String> net2planParameters) 111 { 112 /* Initialize all InputParameter objects defined in this object (this uses Java reflection) */ 113 InputParameter.initializeAllInputParameterFieldsOfObject(this, algorithmParameters); 114 if (netPlan.getNumberOfNodes() == 0) throw new Net2PlanException ("The input design has no nodes"); 115 116 117 try{ 118 this.netPlan = netPlan; 119 /* Initializations */ 120 121 netPlan.removeAllLinks(); netPlan.setRoutingType(RoutingType.SOURCE_ROUTING); 122 this.N = netPlan.getNumberOfNodes(); 123 this.rng = new Random (algorithm_randomSeed.getLong ()); 124 125 /* Initialize links */ 126 for (Node n1 : netPlan.getNodes()) for (Node n2 : netPlan.getNodes()) if (n1 != n2) netPlan.addLink(n1, n2, 0.01, netPlan.getNodePairEuclideanDistance(n1, n2), 200000 , null); 127 final double linkLengthKmFromEuclideanDistanceFactor = tcfa_networkDiameterKmToNormalize.getDouble() / netPlan.getVectorLinkLengthInKm().getMaxLocation() [0]; 128 for (Link e : netPlan.getLinks()) e.setLengthInKm(e.getLengthInKm() * linkLengthKmFromEuclideanDistanceFactor); 129 this.Efm = netPlan.getNumberOfLinks(); if (Efm != N*(N-1)) throw new RuntimeException ("Bad, Efm" + Efm + ", N: " + N); 130 this.linkCost_e = DoubleFactory1D.dense.make(Efm); for (Link e : netPlan.getLinks ()) linkCost_e.set(e.getIndex () , tcfa_linkCostPerKm.getDouble() * e.getLengthInKm()); 131 132 /* If there is no traffic, create uniform traffic matrix */ 133 if (netPlan.getNumberOfDemands() == 0) 134 for (Node n1 : netPlan.getNodes()) for (Node n2 : netPlan.getNodes()) if (n1 != n2) netPlan.addDemand(n1, n2, tcfa_circuitCapacity_Gbps.getDouble(), null); 135 136 137 /* Convert the demand offered traffic in the upper multiple of circuit capacity */ 138 DoubleMatrix2D trafficMatrix = netPlan.getMatrixNode2NodeOfferedTraffic(); 139 netPlan.removeAllDemands(); 140 for (Node n1 : netPlan.getNodes()) for (Node n2 : netPlan.getNodes()) if (n1.getIndex () > n2.getIndex ()) 141 { 142 final double maxTrafficBidir = Math.max(trafficMatrix.get(n1.getIndex (), n2.getIndex ()), trafficMatrix.get(n2.getIndex (), n1.getIndex ())); 143 netPlan.addDemand(n1, n2, tcfa_circuitCapacity_Gbps.getDouble() * Math.ceil(maxTrafficBidir / tcfa_circuitCapacity_Gbps.getDouble()) , null); 144 netPlan.addDemand(n2, n1, tcfa_circuitCapacity_Gbps.getDouble() * Math.ceil(maxTrafficBidir / tcfa_circuitCapacity_Gbps.getDouble()) , null); 145 } 146 147 /* If 1+1 then decouple demands in one per channel */ 148 if (tcfa_recoveryType.getString ().equals("1+1")) 149 { 150 NetPlan npcopy = netPlan.copy(); 151 netPlan.removeAllDemands(); 152 for (Demand d : npcopy.getDemands()) 153 { 154 final int numChannels = (int) Math.round (d.getOfferedTraffic() / tcfa_circuitCapacity_Gbps.getDouble ()); 155 for (int cont = 0 ; cont < numChannels ; cont ++) netPlan.addDemand(netPlan.getNodeFromId (d.getIngressNode().getId ()), netPlan.getNodeFromId (d.getEgressNode().getId ()), tcfa_circuitCapacity_Gbps.getDouble () , d.getAttributes()); 156 } 157 } 158 159 /* Initialize routes and demands to make it bidirectional */ 160 Pair<Map<Demand,Demand>,Map<Route,Route>> pair = initializeNetPlanLinksBidirDemandsAndRoutes (netPlan); 161 this.opposite_d = pair.getFirst(); 162 this.opposite_r = pair.getSecond(); 163 164 /* Initialize demand info */ 165 this.D = netPlan.getNumberOfDemands(); 166 this.numCircH_d = netPlan.getVectorDemandOfferedTraffic().assign(DoubleFunctions.div(this.tcfa_circuitCapacity_Gbps.getDouble ())).assign(DoubleFunctions.rint); 167 System.out.println("Total number of circuits to establish: " + numCircH_d.zSum()); 168 169 //if (1==1) throw new RuntimeException ("Stop"); 170 this.R = netPlan.getNumberOfRoutes(); 171 //System.out.println ("numCircH_d: " + numCircH_d); 172 173 /* Initialize SRG information */ 174 initializeSRGs(netPlan); 175 this.nSRGs = netPlan.getNumberOfSRGs(); 176 this.routeList_d = new int [D][]; for (Demand d : netPlan.getDemands ()) routeList_d [d.getIndex ()] = getIndexes (d.getRoutes ()); 177 this.traversedLinks_r = new int [R][]; for (Route r : netPlan.getRoutes ()) traversedLinks_r [r.getIndex ()] = getIndexes (r.getSeqLinksRealPath()); 178 179 180 /* Initialize aux arrays for speed-up */ 181 this.A_dr = netPlan.getMatrixDemand2RouteAssignment(); 182 for (Demand d : netPlan.getDemands ()) if (d.getRoutes().isEmpty()) throw new Net2PlanException ("A demand has no assigned routes"); 183 this.A_er = netPlan.getMatrixLink2RouteAssignment(); 184 this.A_se = DoubleFactory2D.dense.make(1 + nSRGs , Efm , 1.0); // 1 if link OK, 0 if fails 185 for (int contSRG = 0 ; contSRG < nSRGs ; contSRG ++) 186 for (Link e : netPlan.getSRG (contSRG).getLinks()) 187 A_se.set(contSRG+1 , e.getIndex () , 0.0); 188 this.A_rs = DoubleFactory2D.dense.make(R , 1 + nSRGs , 1.0); // 1 if link OK, 0 if fails 189 for (Route r : netPlan.getRoutes ()) for (Link e : r.getSeqLinksRealPath()) for (SharedRiskGroup srg : e.getSRGs()) 190 A_rs.set(r.getIndex (),1+srg.getIndex () , 0.0); 191 192 /* Check if the problem may have a solution: demands have at least one path in any failure state */ 193 if (A_dr.zMult(A_rs, null).getMinLocation() [0] == 0) { System.out.println("A_dr.zMult(A_rs, null): " + A_dr.zMult(A_rs, null)); throw new Net2PlanException ("Some demands cannot be routed in some failure state. We need more paths! (e.g. extend the path reach)"); } 194 195 this.Aout_ne = netPlan.getMatrixNodeLinkOutgoingIncidence(); 196 this.Ain_ne = netPlan.getMatrixNodeLinkIncomingIncidence(); 197 this.Abid_dd = DoubleFactory2D.sparse.make(D,D); for (Entry<Demand,Demand> entry : opposite_d.entrySet()) { Abid_dd.set(entry.getKey().getIndex (), entry.getValue().getIndex (), 1.0); Abid_dd.set(entry.getValue().getIndex (), entry.getKey().getIndex (), 1.0); } 198 this.Abid_rr = DoubleFactory2D.sparse.make(R,R); for (Entry<Route,Route> entry : opposite_r.entrySet()) { Abid_rr.set(entry.getKey().getIndex (), entry.getValue().getIndex (), 1.0); Abid_rr.set(entry.getValue().getIndex (), entry.getKey().getIndex (), 1.0); } 199 this.Abid_ee = DoubleFactory2D.sparse.make(Efm,Efm); for (Link e : netPlan.getLinks ()) { Abid_ee.set(e.getIndex (), oppositeLink(e).getIndex (), 1.0); Abid_ee.set(oppositeLink(e).getIndex (), e.getIndex (), 1.0); } 200 201 final long algorithmInitialtime = System.nanoTime(); 202 final long algorithmEndtime = algorithmInitialtime + (long) (algorithm_maxExecutionTimeInSeconds.getDouble() * 1E9); 203 204 /* Up stage */ 205 ArrayList<Integer> shuffledDemands = new ArrayList<Integer> (D); for (int d = 0 ; d < D ; d ++) shuffledDemands.add(d); 206 Collections.shuffle(shuffledDemands , rng); 207 if (tcfa_recoveryType.getString ().equals("shared")) 208 { 209 DoubleMatrix1D best_xr = null; DoubleMatrix1D best_pe = null; double bestCost = Double.MAX_VALUE; 210 int iterationCounter = 0; 211 while ((System.nanoTime() < algorithmEndtime) && (iterationCounter < tcfa_maxNumIterations.getInt ())) 212 { 213 DoubleMatrix1D thisIt_xr = greedyAndLocalSearch_shared(algorithmEndtime , netPlan); 214 DoubleMatrix1D thisIt_pe = A_er.zMult(thisIt_xr, null); thisIt_pe.assign(DoubleFunctions.div(tcfa_linkCapacity_numCirc.getInt ())).assign(DoubleFunctions.ceil); 215 double thisItCost = computeCost_shared(thisIt_xr, thisIt_pe); 216 if (thisItCost < bestCost) { best_xr = thisIt_xr; best_pe = thisIt_pe; bestCost = thisItCost; } 217 iterationCounter ++; 218 } 219 stat_numItGRASP = iterationCounter; 220 for (Route r : netPlan.getRoutes ()) r.setCarriedTraffic(tcfa_circuitCapacity_Gbps.getDouble () * best_xr.get(r.getIndex ()) , tcfa_circuitCapacity_Gbps.getDouble () * best_xr.get(r.getIndex ())); 221 for (Link e : netPlan.getLinks ()) e.setCapacity(tcfa_linkCapacity_numCirc.getInt () * tcfa_circuitCapacity_Gbps.getDouble () * best_pe.get(e.getIndex())); 222 } else if (tcfa_recoveryType.getString ().equals("1+1")) 223 { 224 DoubleMatrix1D best_xr = null; DoubleMatrix1D best_x2r = null; DoubleMatrix1D best_pe = null; double bestCost = Double.MAX_VALUE; 225 int iterationCounter = 0; 226 while ((System.nanoTime() < algorithmEndtime) && (iterationCounter < tcfa_maxNumIterations.getInt ())) 227 { 228 Pair<DoubleMatrix1D,DoubleMatrix1D> thisItSolution = greedyAndLocalSearch_11(algorithmEndtime , netPlan); 229 DoubleMatrix1D thisIt_xr = thisItSolution.getFirst(); 230 DoubleMatrix1D thisIt_x2r = thisItSolution.getSecond(); 231 DoubleMatrix1D thisIt_pe = computePe_11(thisIt_xr, thisIt_x2r); 232 double thisItCost = computeCost_11(thisIt_xr, thisIt_x2r , thisIt_pe); 233 if (thisItCost < bestCost) { best_xr = thisIt_xr; best_x2r = thisIt_x2r; best_pe = thisIt_pe; bestCost = thisItCost; } 234 iterationCounter ++; 235 } 236 stat_numItGRASP = iterationCounter; 237 for (Demand d : netPlan.getDemands ()) 238 { 239 if (Math.abs(best_xr.viewSelection(routeList_d [d.getIndex ()]).zSum() - 1) > 1E-3) throw new RuntimeException ("Bad"); 240 if (Math.abs(best_x2r.viewSelection(routeList_d [d.getIndex ()]).zSum() - 1) > 1E-3) throw new RuntimeException ("Bad"); 241 Route primaryRoute = null; ProtectionSegment protSegment = null; 242 for (Route r : d.getRoutes ()) if (Math.abs(best_xr.get(r.getIndex ()) - 1) <= 1e-3) { primaryRoute = r; primaryRoute.setCarriedTraffic(tcfa_circuitCapacity_Gbps.getDouble () , tcfa_circuitCapacity_Gbps.getDouble ()); break; } 243 for (Route r : d.getRoutes ()) if (Math.abs(best_x2r.get(r.getIndex ()) - 1) <= 1e-3) { protSegment = netPlan.addProtectionSegment(r.getSeqLinksRealPath() , tcfa_circuitCapacity_Gbps.getDouble () , null); break; } 244 primaryRoute.addProtectionSegment(protSegment); 245 } 246 for (Link e : netPlan.getLinks ()) e.setCapacity(tcfa_linkCapacity_numCirc.getInt () * tcfa_circuitCapacity_Gbps.getDouble () * best_pe.get(e.getIndex())); 247 248 } else if (tcfa_recoveryType.getString ().equals("restoration")) 249 { 250 DoubleMatrix2D best_xrs = null; DoubleMatrix1D best_pe = null; double bestCost = Double.MAX_VALUE; 251 int iterationCounter = 0; 252 while ((System.nanoTime() < algorithmEndtime) && (iterationCounter < tcfa_maxNumIterations.getInt ())) 253 { 254 DoubleMatrix2D thisIt_xrs = greedyAndLocalSearch_restoration(algorithmEndtime , netPlan); 255 DoubleMatrix1D thisIt_pe = computePe_restoration (thisIt_xrs); 256 double thisItCost = computeCost_shared(thisIt_xrs.viewColumn(0), thisIt_pe); 257 if (thisItCost < bestCost) { best_xrs = thisIt_xrs; best_pe = thisIt_pe; bestCost = thisItCost; } 258 iterationCounter ++; 259 } 260 stat_numItGRASP = iterationCounter; 261 for (Route r : netPlan.getRoutes ()) r.setCarriedTraffic(tcfa_circuitCapacity_Gbps.getDouble () * best_xrs.get(r.getIndex (),0) , tcfa_circuitCapacity_Gbps.getDouble () * best_xrs.get(r.getIndex (),0)); 262 for (Link e : netPlan.getLinks ()) e.setCapacity(tcfa_linkCapacity_numCirc.getInt () * tcfa_circuitCapacity_Gbps.getDouble () * best_pe.get(e.getIndex())); 263 } else throw new RuntimeException ("Bad"); 264 265 /* Remove unused routes and links */ 266 for (Route r : new HashSet<Route>(netPlan.getRoutes())) if (r.getCarriedTraffic() == 0) r.remove(); 267 for (Link e : new HashSet<Link>(netPlan.getLinks())) if (e.getCapacity() == 0) e.remove(); 268 269 Quadruple<Double,Double,Double,Double> q = computeCost (netPlan); 270 stat_totalCost = q.getFirst(); 271 stat_costinks = q.getSecond(); 272 stat_costNodes = q.getThird(); 273 stat_costCirc = q.getFourth(); 274 stat_numLinks = (int) (netPlan.getVectorLinkCapacity().zSum () / (tcfa_linkCapacity_numCirc.getInt () * tcfa_circuitCapacity_Gbps.getDouble ())); 275 stat_numCirc = (int) Math.round(netPlan.getDemandTotalCarriedTraffic() / tcfa_circuitCapacity_Gbps.getDouble ()); if (tcfa_recoveryType.getString ().equals("1+1")) stat_numCirc *=2; 276 for (Node n : netPlan.getNodes()) if (n.getOutgoingLinks().size() > 2) stat_numNodesDeg3 ++; else stat_numNodesDeg2 ++; 277 stat_totalTimeSecs = (System.nanoTime() - algorithmInitialtime)*1e-9; 278 checkSolution (netPlan); 279 280 double averageLinkUtilization = netPlan.getVectorLinkUtilizationIncludingProtectionSegments().zSum() / stat_numLinks; 281 final String fileNameStem = algorithm_outputFileNameRoot.getString() + "_" + tcfa_srgType.getString() + "_c" + tcfa_circuitCapacity_Gbps.getDouble () + "_t" + algorithm_maxExecutionTimeInSeconds.getDouble() + "_r" + tcfa_recoveryType.getString (); 282 try 283 { 284 PrintWriter pw = new PrintWriter (new File (fileNameStem + "_allResults.txt")); 285 pw.println(stat_totalCost ); 286 pw.println( stat_costinks ); 287 pw.println( stat_costNodes ); 288 pw.println( stat_costCirc ); 289 pw.println( stat_numLinks ); 290 pw.println( stat_numCirc ); 291 pw.println( stat_numNodesDeg2 ); 292 pw.println( stat_numNodesDeg3 ); 293 pw.println( stat_numItGRASP ); 294 pw.println( stat_numLSIterationsReducingCost ); 295 pw.println( stat_numSolverCalls ); 296 pw.println( stat_numSolverCallsOk ); 297 pw.println( stat_totalTimeSecs ); 298 pw.println( stat_totalTimeCreateProgramsSecs ); 299 pw.println( stat_totalTimeSolverSecs ); 300 pw.println( averageLinkUtilization ); 301 pw.close (); 302 } catch (Exception e) { e.printStackTrace(); throw new RuntimeException ("Not possible to write in File " + fileNameStem + "_allResults.txt"); } 303 304 305 netPlan.setAttribute("stat_totalCost", ""+stat_totalCost); 306 netPlan.setAttribute("stat_costinks", ""+stat_costinks); 307 netPlan.setAttribute("stat_costNodes", ""+stat_costNodes); 308 netPlan.setAttribute("stat_costCirc", ""+stat_costCirc); 309 netPlan.setAttribute("stat_numLinks", ""+stat_numLinks); 310 netPlan.setAttribute("stat_numCirc", ""+stat_numCirc); 311 netPlan.setAttribute("stat_numNodesDeg2", ""+stat_numNodesDeg2); 312 netPlan.setAttribute("stat_numNodesDeg3", ""+stat_numNodesDeg3); 313 netPlan.setAttribute("stat_numItGRASP", ""+stat_numItGRASP); 314 netPlan.setAttribute("stat_numLSIterationsReducingCost", ""+stat_numLSIterationsReducingCost); 315 netPlan.setAttribute("stat_numSolverCalls", ""+stat_numSolverCalls); 316 netPlan.setAttribute("stat_numSolverCallsOk", ""+stat_numSolverCallsOk); 317 netPlan.setAttribute("stat_totalTimeSecs", ""+stat_totalTimeSecs); 318 netPlan.setAttribute("stat_totalTimeCreateProgramsSecs", ""+stat_totalTimeCreateProgramsSecs); 319 netPlan.setAttribute("stat_totalTimeSolverSecs", ""+stat_totalTimeSolverSecs); 320 321 netPlan.saveToFile(new File (fileNameStem + ".n2p")); 322 323// if (netPlan.getDemandTotalBlockedTraffic() > 1E-3) throw new RuntimeException ("Bad"); 324// if (netPlan.getLinksOversubscribed().size() > 0) throw new RuntimeException ("Bad"); 325 System.out.println("JOM sol: cost: " + stat_totalCost + ", num links: "+ stat_numLinks); 326 327 return "Ok! cost: " + stat_totalCost + ", num links: "+ stat_numLinks; 328 329 } catch (Exception e) { e.printStackTrace(); throw new RuntimeException ("BAD OUT"); } 330 331 } 332 333 @Override 334 public String getDescription() 335 { 336 return "This algorithm is devoted to solve the several network planning problems appearing in the case study in Section 12.10 of the book mentioned below. The context is the planning of a WDM optical network, given a set of node locations, and a set of end-to-end optical connections (lightpaths) to establish. The algorithm should decide on the fibers to deploy, and the routing of the lightpaths on the fibers, minimizing the total network cost. The cost model includes the fiber renting cost (paid to a dark fiber provider), OADM cost (separating between OADMs of degree two, and degree higher than two), and the cost of the transponders. The network must be tolerant to a set of user-defined failures (defined by SRGs). The algorithm permits choosing between two recovery types: 1+1, and shared protection. In the latter, the number of lightpaths between two nodes that survive to any of the failures is at least the number of optical connections requested for that node pair. For solving the problem, the algorithm uses a GRASP heuristic, to govern the iterative search that is performed using small ILPs solved with JOM in each iteration."; 337 } 338 339 @Override 340 public List<Triple<String, String, String>> getParameters() 341 { 342 /* Returns the parameter information for all the InputParameter objects defined in this object (uses Java reflection) */ 343 return InputParameter.getInformationAllInputParameterFieldsOfObject(this); 344 } 345 346 private Quadruple<Double,Double,Double,Double> computeCost (NetPlan np) 347 { 348 double costLinks = 0; 349 for (Link e : np.getLinks ()) 350 costLinks += linkCost_e.get(e.getIndex ()) * e.getCapacity() / (tcfa_linkCapacity_numCirc.getInt () * tcfa_circuitCapacity_Gbps.getDouble ()); 351 352 double costCircuits = 0; 353 for (Route r : np.getRoutes ()) 354 costCircuits += tcfa_circuitCost.getDouble () * (r.getCarriedTraffic() / tcfa_circuitCapacity_Gbps.getDouble ()); 355 for (ProtectionSegment p : np.getProtectionSegments()) 356 costCircuits += tcfa_circuitCost.getDouble () * (p.getReservedCapacityForProtection() / tcfa_circuitCapacity_Gbps.getDouble ()); 357 double costNodes = 0; 358 for (Node n : np.getNodes()) 359 { 360 final int nodeDegree = (int) Math.max(n.getOutgoingLinks().size(), n.getIncomingLinks().size()); 361 costNodes += (nodeDegree <= 2)? tcfa_costNodeDegree2.getDouble () : tcfa_costNodeDegreeHigherThan2.getDouble (); 362 } 363 final double totalCost = costLinks+ costCircuits + costNodes; 364 365 System.out.println("Total cost: " + totalCost + ", links %: " + (costLinks/totalCost) + ", circ %: " + (costCircuits/totalCost) + ", nodes %: " + (costNodes/totalCost)); 366 367 System.out.println("-- From np cost total :" + totalCost + ", links: " + costLinks + ", nodes: " + costNodes + ", circ: " + costCircuits); 368 369 return Quadruple.of(totalCost,costLinks,costNodes,costCircuits); 370 } 371 372 private double computeCost_shared (DoubleMatrix1D x_r , DoubleMatrix1D p_e) 373 { 374 double costLinks = p_e.zDotProduct(linkCost_e); 375 final double costCircuits = tcfa_circuitCost.getDouble () * x_r.zSum(); 376 double costNodes = 0; 377 DoubleMatrix1D degree_n = Aout_ne.zMult(p_e, null); 378 for (int n = 0; n < N ; n ++) 379 costNodes += (degree_n.get(n) <= 2)? tcfa_costNodeDegree2.getDouble () : tcfa_costNodeDegreeHigherThan2.getDouble (); 380 final double totalCost = costLinks+ costCircuits + costNodes; 381 //System.out.println("-- shared cost total :" + totalCost + ", links: " + costLinks + ", nodes: " + costNodes + ", circ: " + costCircuits); 382 return totalCost; 383 } 384 private double computeCost_11 (DoubleMatrix1D x_r , DoubleMatrix1D x2_r , DoubleMatrix1D p_e) 385 { 386 double costLinks = p_e.zDotProduct(linkCost_e); 387 final double costCircuits = tcfa_circuitCost.getDouble () * (x_r.zSum() + x2_r.zSum()); 388 double costNodes = 0; 389 DoubleMatrix1D degree_n = Aout_ne.zMult(p_e, null); 390 for (int n = 0; n < N ; n ++) 391 costNodes += (degree_n.get(n) <= 2)? tcfa_costNodeDegree2.getDouble () : tcfa_costNodeDegreeHigherThan2.getDouble (); 392 final double totalCost = costLinks+ costCircuits + costNodes; 393 //System.out.println("-- shared cost total :" + totalCost + ", links: " + costLinks + ", nodes: " + costNodes + ", circ: " + costCircuits); 394 return totalCost; 395 } 396 397 private Link oppositeLink (Link e) { return e.getNetPlan ().getNodePairLinks(e.getDestinationNode(), e.getOriginNode() , false).iterator().next(); } 398 399 private Pair<Map<Demand,Demand>,Map<Route,Route>> initializeNetPlanLinksBidirDemandsAndRoutes (NetPlan np) 400 { 401 /* Remove lower half demands from np */ 402 np.removeAllRoutes(); np.removeAllProtectionSegments(); 403 for (Node n1 : np.getNodes()) for (Node n2 : np.getNodes()) if (n1.getIndex () > n2.getIndex ()) for (Demand d : np.getNodePairDemands(n1, n2,false)) d.remove (); 404 np.addRoutesFromCandidatePathList(null , "K" , "" + tcfa_maxNumberPathsPerDemand.getInt () , "maxLengthInKm" , "" + tcfa_maxPathLengthInKm.getDouble() , "maxNumHops" , "" + tcfa_maxPathNumberOfHops.getInt ()); 405 406 /* Add symmetric demands and routes */ 407 Map<Demand,Demand> opposite_d = new HashMap<Demand,Demand> (); 408 Map<Route,Route> opposite_r = new HashMap<Route,Route> (); 409 for (Demand d : new HashSet<Demand> (np.getDemands())) 410 { 411 final Demand opDemand = np.addDemand(d.getEgressNode(), d.getIngressNode(), d.getOfferedTraffic(), null); 412 opposite_d.put(d,opDemand); 413 opposite_d.put(opDemand,d); 414 for (Route r : new HashSet<Route> (d.getRoutes ())) 415 { 416 final Route oppRoute = np.addRoute(opDemand, r.getCarriedTraffic(), r.getOccupiedCapacity() , oppositeSeqLinks (r.getSeqLinksRealPath()), null); 417 opposite_r.put(r,oppRoute); opposite_r.put(oppRoute,r); 418 } 419 } 420 return Pair.of(opposite_d,opposite_r); 421 } 422 423 private void initializeSRGs (NetPlan np) 424 { 425 np.removeAllSRGs(); 426 if (tcfa_srgType.getString ().equals("perBidirectionalLinkBundle")) 427 { 428 for (Node n1 : np.getNodes()) for (Node n2 : np.getNodes()) if (n1.getIndex () < n2.getIndex()) 429 { 430 final double linkLengthKm = np.getNodePairLinks(n1, n2,false).iterator().next().getLengthInKm(); 431 final SharedRiskGroup srg = np.addSRG(tcfa_srgMttfPer1000Km_hrs.getDouble () * linkLengthKm / 1000, tcfa_srgMttr_hrs.getDouble (), null); 432 for (Link e : np.getNodePairLinks(n1, n2,true)) srg.addLink(e); 433 } 434 } 435 else if (tcfa_srgType.getString ().equals("noFailure")) 436 { 437 } else throw new Net2PlanException ("Wrong SRG type"); 438 } 439 440 private List<Link> oppositeSeqLinks (List<Link> s) 441 { 442 LinkedList<Link> oppList = new LinkedList<Link> (); 443 for (Link e : s) oppList.addFirst(oppositeLink(e)); 444 return (List<Link>) oppList; 445 } 446 447 private void checkSolution (NetPlan np) 448 { 449 if (np.getDemandsBlocked().size() != 0) 450 { 451 System.out.println("np.getDemandsBlocked(): " + np.getDemandsBlocked()); 452 System.out.println("np.getDemandOfferedTrafficMap(): " + np.getVectorDemandOfferedTraffic()); 453 System.out.println("np.getDemandBlockedTrafficMap(): " + np.getVectorDemandBlockedTraffic()); 454 throw new RuntimeException ("Bad"); 455 } 456 if (np.getLinksOversubscribed().size() != 0) 457 { 458 System.out.println("np.getLinksOversubscribed(): " + np.getLinksOversubscribed()); 459 System.out.println("np.getLinkCapacityMap(): " + np.getVectorLinkCapacity()); 460 System.out.println("np.getLinkCarriedTrafficMap(): " + np.getVectorLinkTotalCarriedTraffic()); 461 throw new RuntimeException ("Bad"); 462 } 463 464 /* Check bidirectional topology */ 465 for (Node n1 : np.getNodes()) 466 for (Node n2 : np.getNodes()) 467 if (n1.getIndex () > n2.getIndex ()) 468 { 469 Set<Link> thisLinks = np.getNodePairLinks(n1, n2,false); 470 Set<Link> oppLinks = np.getNodePairLinks(n2, n1,false); 471 if (thisLinks.size() > 1) throw new RuntimeException ("Bad"); 472 if (oppLinks.size() > 1) throw new RuntimeException ("Bad"); 473 if (tcfa_bidirectionalLinks.getBoolean()) 474 { 475 if (thisLinks.size() != oppLinks.size()) throw new RuntimeException ("Bad"); 476 if (thisLinks.size() > 0) 477 if (thisLinks.iterator().next().getCapacity() != oppLinks.iterator().next().getCapacity() ) throw new RuntimeException ("Bad"); 478 } 479 } 480 481 /* Check bidirectional routes */ 482 for (Route r : np.getRoutes()) 483 { 484 final Route oppRoute = opposite_r.get(r); 485 if (!np.getRoutes().contains(oppRoute)) throw new RuntimeException ("Bad"); 486 if (r.getCarriedTraffic() != oppRoute.getCarriedTraffic()) throw new RuntimeException ("Bad"); 487 } 488 489 if (tcfa_recoveryType.getString ().equals("1+1")) 490 { // one prot segment, and link disjoint 491 for (Route r : np.getRoutes()) 492 { 493 if (r.getPotentialBackupProtectionSegments().size() != 1) throw new RuntimeException ("Bad"); 494 final ProtectionSegment protSegment = r.getPotentialBackupProtectionSegments().iterator().next(); 495 List<Link> seqLinks = new ArrayList<Link> (r.getSeqLinksRealPath()); 496 seqLinks.retainAll(protSegment.getSeqLinks()); 497 if (!seqLinks.isEmpty()) 498 { 499 System.out.println ("Route : "+ r + " links: " + r.getSeqLinksRealPath() + ", prot segment " + protSegment + " links: " + protSegment.getSeqLinks() + ", get route carried traffic: " + r.getCarriedTraffic() + ", prot segment reserved traffic: " + protSegment.getReservedCapacityForProtection()); 500 throw new RuntimeException ("Bad"); 501 } 502 } 503 } 504 505 } 506 507 private DoubleMatrix1D jomStep_shared (Set<Integer> demandsChange , DoubleMatrix1D maxPe , DoubleMatrix1D x_r) 508 { 509 final long initTimeBeforeCreatingProgram = System.nanoTime(); 510 if (demandsChange.isEmpty()) throw new RuntimeException ("Bad"); 511 final int [] dChange = IntUtils.toArray(demandsChange); 512 int [] rChange = new int [0]; for (int d : dChange) rChange = IntUtils.concatenate(rChange, routeList_d [d]); 513 final int Rrest = rChange.length; 514 x_r.viewSelection(rChange).assign(0); 515 DoubleMatrix2D A_dRestrRest = A_dr.viewSelection(dChange, rChange); 516 DoubleMatrix2D A_erRest = A_er.viewSelection(null, rChange); 517 DoubleMatrix2D A_rRests = A_rs.viewSelection(rChange,null); 518 DoubleMatrix1D occupiedWithoutChangingR_e = A_er.zMult (x_r,null); 519 DoubleMatrix2D Abid_rRestrRest = DoubleFactory2D.sparse.make(Rrest,Rrest); 520 for (int rRest = 0 ; rRest < Rrest ; rRest ++) 521 { 522 final int r = rChange [rRest]; 523 final int opp_r = opposite_r.get(netPlan.getRoute(r)).getIndex (); 524 final int [] oppRrest = IntUtils.find(rChange, opp_r, SearchType.ALL); 525 if (oppRrest.length != 1) throw new RuntimeException ("Bad"); 526 Abid_rRestrRest.set(rRest , oppRrest [0] , 1.0); Abid_rRestrRest.set(oppRrest [0], rRest , 1); 527 } 528 OptimizationProblem op = new OptimizationProblem(); 529 op.setInputParameter("R",R); 530 op.setInputParameter("U",tcfa_linkCapacity_numCirc.getInt ()); 531 op.setInputParameter("c_e", linkCost_e.toArray() , "row"); 532 op.setInputParameter("c", tcfa_circuitCost.getDouble ()); 533 op.setInputParameter("d3Cost", tcfa_costNodeDegreeHigherThan2.getDouble ()); 534 op.setInputParameter("onesS", DoubleUtils.ones(1+nSRGs) , "row"); 535 op.setInputParameter("onesE", DoubleUtils.ones(Efm) , "row"); 536 op.setInputParameter("A_dRestrRest", new DoubleMatrixND (A_dRestrRest)); 537 op.setInputParameter("h_dRest", numCircH_d.viewSelection(dChange).toArray(), "row"); 538 op.setInputParameter("A_rRests", new DoubleMatrixND (A_rRests)); 539 op.setInputParameter("A_erRest", new DoubleMatrixND (A_erRest)); 540 op.setInputParameter("Aout_ne", new DoubleMatrixND (Aout_ne)); 541 op.setInputParameter("Abid_rRestrRest", new DoubleMatrixND (Abid_rRestrRest)); 542 op.setInputParameter("occup_e", occupiedWithoutChangingR_e.toArray(), "row"); 543 544 op.addDecisionVariable("d3_n", true, new int[] { 1, N }, 0, 1); 545 op.addDecisionVariable("p_e", true, new int[] { 1, Efm }, DoubleUtils.zeros(Efm) , (maxPe == null)? DoubleUtils.constantArray(Efm, Double.MAX_VALUE) : maxPe.toArray()); 546 op.addDecisionVariable("x_rRest", true, new int[] { 1, Rrest}, 0 , Double.MAX_VALUE); /* 1 if there is a link from node i to node j, 0 otherwise */ 547 548 op.setObjectiveFunction("minimize", "sum(c_e .* p_e) + c * sum (x_rRest) + d3Cost * sum(d3_n)"); 549 op.addConstraint("A_dRestrRest * diag(x_rRest) * A_rRests >= h_dRest' * onesS"); /* the flow-conservation constraints (NxD constraints) */ 550 op.addConstraint("A_erRest * x_rRest' + occup_e' <= U * p_e'"); /* the capacity constraints (E constraints) */ 551 op.addConstraint("Aout_ne * p_e' <= 2 + R * d3_n'"); /* the capacity constraints (E constraints) */ 552 op.addConstraint("x_rRest == x_rRest * Abid_rRestrRest"); 553 if (tcfa_bidirectionalLinks.getBoolean()) { op.setInputParameter("Abid_ee", new DoubleMatrixND (Abid_ee)); op.addConstraint("p_e == p_e * Abid_ee"); } 554 555 stat_totalTimeCreateProgramsSecs += eTime (initTimeBeforeCreatingProgram); 556 557 double solverTimeAllowed = algorithm_maxSolverTimeInSeconds.getDouble(); 558 do 559 { 560 final long t = System.nanoTime(); 561 try 562 { 563 stat_numSolverCalls ++; op.solve(algorithm_solverName.getString() , "solverLibraryName", algorithm_solverLibraryName.getString () , "maxSolverTimeInSeconds" , solverTimeAllowed); 564 } catch (Exception e) { System.out.println ("-- EXTRA TIME FOR SOLVER, NOW: " + solverTimeAllowed); solverTimeAllowed += 2; stat_totalTimeSolverSecs += eTime(t); continue;} 565 stat_totalTimeSolverSecs += eTime(t); 566 if (solverTimeAllowed != algorithm_maxSolverTimeInSeconds.getDouble()) System.out.println ("-- EXTRA TIME FOR SOLVER, NOW: " + solverTimeAllowed); 567 solverTimeAllowed *= 2; // duplicate the time if not feasible solution is found 568 } while (!op.solutionIsFeasible()); 569 stat_numSolverCallsOk ++; 570 final double [] new_x_rRest = op.getPrimalSolution("x_rRest").to1DArray(); 571 x_r.viewSelection(rChange).assign(new_x_rRest); 572 return x_r; 573 } 574 575 private Pair<DoubleMatrix1D,DoubleMatrix1D> jomStep_11 (Set<Integer> demandsChange , DoubleMatrix1D maxPe , DoubleMatrix1D x_r , DoubleMatrix1D x2_r) 576 { 577// System.out.println ("demandsChange :" + demandsChange); 578 final long initTimeBeforeCreatingProgram = System.nanoTime(); 579 580 if (demandsChange.isEmpty()) throw new RuntimeException ("Bad"); 581 final int [] dChange = IntUtils.toArray(demandsChange); 582// System.out.println ("dChange :" + Arrays.toString(dChange)); 583 int [] rChange = new int [0]; for (int d : dChange) rChange = IntUtils.concatenate(rChange, routeList_d [d]); 584 final int Rrest = rChange.length; 585 final int Drest = dChange.length; 586 x_r.viewSelection(rChange).assign(0); 587 x2_r.viewSelection(rChange).assign(0); 588 DoubleMatrix2D A_dRestrRest = A_dr.viewSelection(dChange, rChange); 589 DoubleMatrix2D A_erRest = A_er.viewSelection(null, rChange); 590 DoubleMatrix2D A_rRests = A_rs.viewSelection(rChange,null); 591 DoubleMatrix1D occupiedWithoutChangingR_e = A_er.zMult (x_r,null); 592 DoubleMatrix1D backup_occupiedWithoutChangingR_e = A_er.zMult (x2_r,null); 593 for (int e = 0;e < Efm; e ++) occupiedWithoutChangingR_e.set(e, occupiedWithoutChangingR_e.get(e) + backup_occupiedWithoutChangingR_e.get(e)); 594 DoubleMatrix2D Abid_rRestrRest = DoubleFactory2D.sparse.make(Rrest,Rrest); 595 for (int rRest = 0 ; rRest < Rrest ; rRest ++) 596 { 597 final int r = rChange [rRest]; final int opp_r = opposite_r.get(netPlan.getRoute (r)).getIndex (); 598 final int [] oppRrest = IntUtils.find(rChange, opp_r, SearchType.ALL); if (oppRrest.length != 1) throw new RuntimeException ("Bad"); 599 Abid_rRestrRest.set(rRest , oppRrest [0] , 1.0); Abid_rRestrRest.set(oppRrest [0], rRest , 1); 600 } 601 602 OptimizationProblem op = new OptimizationProblem(); 603 op.setInputParameter("R",R); 604 op.setInputParameter("U",tcfa_linkCapacity_numCirc.getInt ()); 605 op.setInputParameter("c_e", linkCost_e.toArray() , "row"); 606 op.setInputParameter("c", tcfa_circuitCost.getDouble ()); 607 op.setInputParameter("d3Cost", tcfa_costNodeDegreeHigherThan2.getDouble ()); 608 op.setInputParameter("onesS", DoubleUtils.ones(1+nSRGs) , "row"); 609 op.setInputParameter("onesE", DoubleUtils.ones(Efm) , "row"); 610 op.setInputParameter("A_dRestrRest", new DoubleMatrixND (A_dRestrRest)); 611 op.setInputParameter("h_dRest", numCircH_d.viewSelection(dChange).toArray(), "row"); 612 op.setInputParameter("A_rRests", new DoubleMatrixND (A_rRests)); 613 op.setInputParameter("A_erRest", new DoubleMatrixND (A_erRest)); 614 op.setInputParameter("Aout_ne", new DoubleMatrixND (Aout_ne)); 615 op.setInputParameter("Abid_rRestrRest", new DoubleMatrixND (Abid_rRestrRest)); 616 op.setInputParameter("occup_e", occupiedWithoutChangingR_e.toArray(), "row"); 617 618 op.addDecisionVariable("d3_n", true, new int[] { 1, N }, 0, 1); 619 op.addDecisionVariable("p_e", true, new int[] { 1, Efm }, DoubleUtils.zeros(Efm) , (maxPe == null)? DoubleUtils.constantArray(Efm, Double.MAX_VALUE) : maxPe.toArray()); 620 op.addDecisionVariable("x_rRest", true, new int[] { 1, Rrest}, 0 , 1); 621 op.addDecisionVariable("x2_rRest", true, new int[] { 1, Rrest}, 0 , 1); 622 623 op.setObjectiveFunction("minimize", "sum(c_e .* p_e) + d3Cost * sum(d3_n)"); 624 op.addConstraint("A_dRestrRest * x_rRest' == h_dRest'"); /* the flow-conservation constraints (NxD constraints) */ 625 op.addConstraint("A_dRestrRest * x2_rRest' == h_dRest'"); /* the flow-conservation constraints (NxD constraints) */ 626 op.addConstraint("A_erRest * (x_rRest+x2_rRest)' + occup_e' <= U * p_e'"); /* the capacity constraints (E constraints) */ 627 op.addConstraint("Aout_ne * p_e' <= 2 + R*d3_n'"); /* the capacity constraints (E constraints) */ 628 op.addConstraint("A_dRestrRest * diag(x_rRest+x2_rRest) * A_erRest' <= 1"); 629 op.addConstraint("x_rRest == x_rRest * Abid_rRestrRest"); 630 op.addConstraint("x2_rRest == x2_rRest * Abid_rRestrRest"); 631 if (tcfa_bidirectionalLinks.getBoolean()) { op.setInputParameter("Abid_ee", new DoubleMatrixND (Abid_ee)); op.addConstraint("p_e == p_e * Abid_ee"); } 632 633// System.out.println("Abid_rr sum filas: " + Abid_rRestrRest.zMult(DoubleFactory1D.dense.make(Rrest,1), null)); 634// System.out.println("Abid_rr sum cols: " + Abid_rRestrRest.zMult(DoubleFactory1D.dense.make(Rrest,1), null , 1 , 0 , true)); 635// System.out.println("Abid_rr non zeros: " + Abid_rRestrRest.zSum()); 636 stat_totalTimeCreateProgramsSecs += eTime (initTimeBeforeCreatingProgram); 637 638 double solverTimeAllowed = algorithm_maxSolverTimeInSeconds.getDouble (); 639 do 640 { 641 final long t = System.nanoTime(); 642 try 643 { 644 stat_numSolverCalls ++; op.solve(algorithm_solverName.getString() , "solverLibraryName", algorithm_solverLibraryName.getString () , "maxSolverTimeInSeconds" , solverTimeAllowed); 645 } catch (Exception e) { System.out.println ("-- EXTRA TIME FOR SOLVER, NOW: " + solverTimeAllowed); solverTimeAllowed += 2; stat_totalTimeSolverSecs += eTime(t); continue; } 646 if (solverTimeAllowed != algorithm_maxSolverTimeInSeconds.getDouble ()) System.out.println ("-- EXTRA TIME FOR SOLVER, NOW: " + solverTimeAllowed); 647 stat_totalTimeSolverSecs += eTime(t); 648 solverTimeAllowed *= 2; // duplicate the time if not feasible solution is found 649 //System.out.println ("Time to create program: " + (1E-9*(t-initTimeBeforeCreatingProgram)) + ", solver time " + ((System.nanoTime()-t)*1E-9) + " s ; "); 650 } while (!op.solutionIsFeasible()); 651 stat_numSolverCallsOk ++; 652 final double[] new_x_rRest = op.getPrimalSolution("x_rRest").to1DArray(); 653 final double[] new_x2_rRest = op.getPrimalSolution("x2_rRest").to1DArray(); 654 655 x_r.viewSelection(rChange).assign(new_x_rRest); 656 x2_r.viewSelection(rChange).assign(new_x2_rRest); 657 return Pair.of(x_r,x2_r); 658 } 659 660 private DoubleMatrix2D jomStep_restoration (Set<Integer> demandsChange , DoubleMatrix1D maxPe , DoubleMatrix2D x_rs) 661 { 662 final long initTimeBeforeCreatingProgram = System.nanoTime(); 663 if (demandsChange.isEmpty()) throw new RuntimeException ("Bad"); 664 final int [] dChange = IntUtils.toArray(demandsChange); 665 int [] rChange = new int [0]; for (int d : dChange) rChange = IntUtils.concatenate(rChange, routeList_d [d]); 666 final int Rrest = rChange.length; 667 final int Drest = dChange.length; 668 x_rs.viewSelection(rChange,null).assign(0); 669// DoubleMatrix2D A_dRestrRest = A_dr.viewSelection(dChange, rChange); 670 DoubleMatrix2D A_dRestrRest = DoubleFactory2D.sparse.make(Drest,Rrest); 671 for (int rRest = 0 ; rRest < Rrest ; rRest ++) 672 { 673 final int r = rChange [rRest]; final int dRest [] = IntUtils.find(dChange, netPlan.getRoute (r).getDemand().getIndex () , SearchType.ALL); 674 if (dRest.length != 1) throw new RuntimeException ("Bad"); 675 A_dRestrRest.set(dRest [0], rRest , 1.0); 676 } 677 678 DoubleMatrix2D A_erRest = A_er.viewSelection(null, rChange); 679 DoubleMatrix2D A_rRests = A_rs.viewSelection(rChange,null); 680 DoubleMatrix2D occupiedWithoutChangingR_es = A_er.zMult (x_rs,null); 681 DoubleMatrix2D Abid_rRestrRest = DoubleFactory2D.sparse.make(Rrest,Rrest); 682 for (int rRest = 0 ; rRest < Rrest ; rRest ++) 683 { 684 final int r = rChange [rRest]; final int opp_r = opposite_r.get(netPlan.getRoute (r)).getIndex (); 685 final int [] oppRrest = IntUtils.find(rChange, opp_r, SearchType.ALL); 686 if (oppRrest.length != 1) throw new RuntimeException ("Bad"); 687 Abid_rRestrRest.set(rRest , oppRrest [0] , 1.0); Abid_rRestrRest.set(oppRrest [0], rRest , 1); 688 } 689 690 OptimizationProblem op = new OptimizationProblem(); 691 op.setInputParameter("R",R); 692 op.setInputParameter("U",tcfa_linkCapacity_numCirc.getInt ()); 693 op.setInputParameter("c_e", linkCost_e.toArray() , "row"); 694 op.setInputParameter("c", tcfa_circuitCost.getDouble ()); 695 op.setInputParameter("d3Cost", tcfa_costNodeDegreeHigherThan2.getDouble ()); 696 op.setInputParameter("onesS", DoubleUtils.ones(1+nSRGs) , "row"); 697 op.setInputParameter("onesE", DoubleUtils.ones(Efm) , "row"); 698 op.setInputParameter("A_dRestrRest", new DoubleMatrixND (A_dRestrRest)); 699 op.setInputParameter("h_dRest", numCircH_d.viewSelection(dChange).toArray(), "row"); 700 op.setInputParameter("A_rRests", new DoubleMatrixND (A_rRests)); 701 op.setInputParameter("A_erRest", new DoubleMatrixND (A_erRest)); 702 op.setInputParameter("Aout_ne", new DoubleMatrixND (Aout_ne)); 703 op.setInputParameter("Abid_rRestrRest", new DoubleMatrixND (Abid_rRestrRest)); 704 op.setInputParameter("occup_es", occupiedWithoutChangingR_es.toArray()); 705 706 op.addDecisionVariable("d3_n", true, new int[] { 1, N }, 0, 1); 707 op.addDecisionVariable("p_e", true, new int[] { 1, Efm }, DoubleUtils.zeros(Efm) , (maxPe == null)? DoubleUtils.constantArray(Efm, Double.MAX_VALUE) : maxPe.toArray()); 708 op.addDecisionVariable("x_rRests", true, new int[] { Rrest, 1+nSRGs}, 0 , Double.MAX_VALUE); /* 1 if there is a link from node i to node j, 0 otherwise */ 709 710 op.setObjectiveFunction("minimize", "sum(c_e .* p_e) + d3Cost * sum(d3_n)"); 711 op.addConstraint("A_dRestrRest * x_rRests(all,0) == h_dRest'"); /* the flow-conservation constraints (NxD constraints) */ 712 op.addConstraint("A_dRestrRest * (x_rRests .* A_rRests) >= h_dRest' * onesS"); /* the flow-conservation constraints (NxD constraints) */ 713 op.addConstraint("A_erRest * x_rRests + occup_es <= U * p_e' * onesS"); /* the capacity constraints (E constraints) */ 714 op.addConstraint("Aout_ne * p_e' <= 2 + R*d3_n'"); /* the capacity constraints (E constraints) */ 715 op.addConstraint("x_rRests == Abid_rRestrRest * x_rRests"); 716 op.addConstraint("x_rRests >= diag(x_rRests(all,0)) * A_rRests"); // if not affected by failure, do not change 717 if (tcfa_bidirectionalLinks.getBoolean()) { op.setInputParameter("Abid_ee", new DoubleMatrixND (Abid_ee)); op.addConstraint("p_e == p_e * Abid_ee"); } 718 719 stat_totalTimeCreateProgramsSecs += eTime (initTimeBeforeCreatingProgram); 720 721 double solverTimeAllowed = algorithm_maxSolverTimeInSeconds.getDouble (); 722 do 723 { 724 final long t = System.nanoTime(); 725 try { stat_numSolverCalls ++; op.solve(algorithm_solverName.getString() , "solverLibraryName", algorithm_solverLibraryName.getString () , "maxSolverTimeInSeconds" , solverTimeAllowed); } catch (Exception e) { System.out.println ("-- EXTRA TIME FOR SOLVER, NOW: " + solverTimeAllowed); solverTimeAllowed += 2; stat_totalTimeSolverSecs += eTime(t); continue;} 726 if (solverTimeAllowed != algorithm_maxSolverTimeInSeconds.getDouble ()) System.out.println ("-- EXTRA TIME FOR SOLVER, NOW: " + solverTimeAllowed); 727 stat_totalTimeSolverSecs += eTime(t); 728 solverTimeAllowed *= 2; // duplicate the time if not feasible solution is found 729 //System.out.println ("Time to create program: " + (1E-9*(t-initTimeBeforeCreatingProgram)) + ", solver time " + ((System.nanoTime()-t)*1E-9) + " s ; "); 730 } while (!op.solutionIsFeasible()); 731 stat_numSolverCallsOk ++; 732 final DoubleMatrix2D new_x_rRests = op.getPrimalSolution("x_rRests").view2D(); 733 final DoubleMatrix1D new_pe = DoubleFactory1D.dense.make(op.getPrimalSolution("p_e").to1DArray()); 734 x_rs.viewSelection(rChange, null).assign(new_x_rRests); 735 return x_rs; 736 } 737 738 private DoubleMatrix1D greedyAndLocalSearch_shared (long algorithmEndtime , NetPlan np) 739 { 740 final long initGreedy = System.nanoTime(); 741 ArrayList<Pair<Node,Node>> shuffledNodePairs = new ArrayList<Pair<Node,Node>> (N*(N-1)/2); 742 for (Node n1 : np.getNodes ()) for (Node n2 : np.getNodes ()) if (n1.getIndex () < n2.getIndex ()) shuffledNodePairs.add(Pair.of(n1, n2)); 743 Collections.shuffle(shuffledNodePairs , rng); 744 DoubleMatrix1D current_xr = DoubleFactory1D.dense.make(R); 745 for (Pair<Node,Node> nodePair : shuffledNodePairs) 746 { 747 Set<Integer> thisDemand = new HashSet<Integer> (); 748 for (Demand d : np.getNodePairDemands(nodePair.getFirst(), nodePair.getSecond() , true)) 749 thisDemand.add(d.getIndex()); 750// System.out.println("thsiDemand: " + thisDemand); 751 if (jomStep_shared(thisDemand , null , current_xr) == null) throw new RuntimeException ("Returned null"); 752 } 753 DoubleMatrix1D current_pe = A_er.zMult(current_xr, null); current_pe.assign(DoubleFunctions.div(tcfa_linkCapacity_numCirc.getInt ())).assign(DoubleFunctions.ceil); 754 double current_cost = computeCost_shared(current_xr, current_pe); 755 //DoubleMatrix1D best_xr = current_xr.copy(); DoubleMatrix1D best_pe = current_pe.copy(); double bestCost = ; 756 757 System.out.println("---- Greedy solution Time " + (1E-9*(System.nanoTime() - initGreedy)) +", Cost : "+ current_cost + ", num links: " + current_pe.zSum()); 758 759 while (System.nanoTime() < algorithmEndtime) 760 { 761 //DoubleMatrix1D bestNeighbor_xr = null; DoubleMatrix1D bestNeighbor_pe = null; double bestNeighborCost = Double.MAX_VALUE; Pair<Long,Long> bestNeighborNodePair = null; 762 Collections.shuffle(shuffledNodePairs ,rng); 763 final long initLSIteration = System.nanoTime(); 764 765 boolean improvedSolution = false; 766 for (Pair<Node,Node> nodePair : shuffledNodePairs) 767 { 768// System.out.println("Pair: " + nodePair + ". Tabu list: " + tabuList_n1n2 + ", IS TABU: " + tabuList_n1n2.contains(nodePair)); 769 Set<Integer> demandsToChange = new HashSet<Integer> (); 770 for (Demand d : np.getNodePairDemands(nodePair.getFirst(), nodePair.getSecond() , true)) 771 demandsToChange.add(d.getIndex ()); 772 DoubleMatrix1D thisNeighbor_xr = current_xr.copy(); // do not modify the original one 773 if (jomStep_shared(demandsToChange , null , thisNeighbor_xr) == null) throw new RuntimeException ("Returned null"); 774 DoubleMatrix1D thisNeighbor_pe = A_er.zMult(thisNeighbor_xr, null); thisNeighbor_pe.assign(DoubleFunctions.div(tcfa_linkCapacity_numCirc.getInt ())).assign(DoubleFunctions.ceil); 775 final double thisNeighborCost = computeCost_shared (thisNeighbor_xr , thisNeighbor_pe); 776 if (thisNeighborCost < current_cost) 777 { 778 current_xr = thisNeighbor_xr; current_pe = thisNeighbor_pe; current_cost = thisNeighborCost; improvedSolution = true; 779 System.out.println("NEW BEST COST: " + current_cost + ", num links: " + current_pe.zSum()); 780 break; 781 } 782 } 783 784 System.out.println("- Local search TIME: " + ((System.nanoTime()-initLSIteration)*1e-9) + ", current cost: " + current_cost + ", num links: " + current_pe.zSum() + ", IMPROVED: " + improvedSolution); 785 786 if (!improvedSolution) break; else stat_numLSIterationsReducingCost ++; // end grasp iteration 787 } 788 789 return current_xr; 790 } 791 792 private Pair<DoubleMatrix1D,DoubleMatrix1D> greedyAndLocalSearch_11 (long algorithmEndtime , NetPlan np) 793 { 794 final long initGreedy = System.nanoTime(); 795 ArrayList<Pair<Node,Node>> shuffledNodePairs = new ArrayList<Pair<Node,Node>> (N*(N-1)/2); 796 for (Node n1 : np.getNodes ()) for (Node n2 : np.getNodes ()) if (n1.getIndex () < n2.getIndex ()) shuffledNodePairs.add(Pair.of(n1, n2)); 797 Collections.shuffle(shuffledNodePairs , rng); 798 DoubleMatrix1D current_xr = DoubleFactory1D.dense.make(R); 799 DoubleMatrix1D current_x2r = DoubleFactory1D.dense.make(R); 800 for (Pair<Node,Node> nodePair : shuffledNodePairs) 801 { 802 Set<Integer> thisDemand = new HashSet<Integer> (); 803 for (Demand d : np.getNodePairDemands(nodePair.getFirst(), nodePair.getSecond() , true)) 804 thisDemand.add(d.getIndex()); 805 final long t = System.nanoTime(); if (jomStep_11(thisDemand , null , current_xr , current_x2r) == null) throw new RuntimeException ("Returned null"); 806 //System.out.println("Node pair " + nodePair + " time increasing it: " + (1E-9*(System.nanoTime()-t)) + ", num demands: " + thisDemand.size()); 807 } 808 DoubleMatrix1D current_pe = computePe_11 (current_xr , current_xr); 809 double current_cost = computeCost_11(current_xr, current_x2r , current_pe); 810 811 System.out.println("---- Greedy solution Time " + (1E-9*(System.nanoTime() - initGreedy)) +", Cost : "+ current_cost + ", num links: " + current_pe.zSum()); 812 813 while (System.nanoTime() < algorithmEndtime) 814 { 815 Collections.shuffle(shuffledNodePairs ,rng); 816 final long initLSIteration = System.nanoTime(); 817 818 boolean improvedSolution = false; 819 for (Pair<Node,Node> nodePair : shuffledNodePairs) 820 { 821 Set<Integer> demandsToChange = new HashSet<Integer> (); 822 for (Demand d : np.getNodePairDemands(nodePair.getFirst(), nodePair.getSecond() , true)) 823 demandsToChange.add(d.getIndex()); 824 DoubleMatrix1D thisNeighbor_xr = current_xr.copy(); // do not modify the original one 825 DoubleMatrix1D thisNeighbor_x2r = current_x2r.copy(); // do not modify the original one 826 if (jomStep_11(demandsToChange , null , thisNeighbor_xr , thisNeighbor_x2r) == null) throw new RuntimeException ("Returned null"); 827 DoubleMatrix1D thisNeighbor_pe = computePe_11(thisNeighbor_xr, thisNeighbor_x2r); 828 final double thisNeighborCost = computeCost_11 (thisNeighbor_xr , thisNeighbor_x2r , thisNeighbor_pe); 829 if (thisNeighborCost < current_cost) 830 { 831 current_xr = thisNeighbor_xr; current_x2r = thisNeighbor_x2r; current_pe = thisNeighbor_pe; current_cost = thisNeighborCost; improvedSolution = true; 832 System.out.println("NEW BEST COST: " + current_cost + ", num links: " + current_pe.zSum()); 833 break; 834 } 835 } 836 837 System.out.println("- Local search TIME: " + ((System.nanoTime()-initLSIteration)*1e-9) + ", current cost: " + current_cost + ", num links: " + current_pe.zSum() + ", IMPROVED: " + improvedSolution); 838 839 if (!improvedSolution) break; else stat_numLSIterationsReducingCost ++; // end grasp iteration 840 } 841 842 return Pair.of(current_xr , current_x2r); 843 } 844 private DoubleMatrix2D greedyAndLocalSearch_restoration (long algorithmEndtime , NetPlan np) 845 { 846 final long initGreedy = System.nanoTime(); 847 ArrayList<Pair<Node,Node>> shuffledNodePairs = new ArrayList<Pair<Node,Node>> (N*(N-1)/2); 848 for (Node n1 : np.getNodes ()) for (Node n2 : np.getNodes ()) if (n1.getIndex () < n2.getIndex()) shuffledNodePairs.add(Pair.of(n1, n2)); 849 Collections.shuffle(shuffledNodePairs , rng); 850 DoubleMatrix2D current_xrs = DoubleFactory2D.dense.make(R,1+nSRGs); 851 for (Pair<Node,Node> nodePair : shuffledNodePairs) 852 { 853 Set<Integer> thisDemand = new HashSet<Integer> (); 854 for (Demand d : np.getNodePairDemands(nodePair.getFirst(), nodePair.getSecond(),true)) 855 thisDemand.add(d.getIndex ()); 856// System.out.println("thsiDemand: " + thisDemand); 857 if (jomStep_restoration(thisDemand , null , current_xrs) == null) throw new RuntimeException ("Returned null"); 858 } 859 DoubleMatrix1D current_pe = computePe_restoration(current_xrs); 860 double current_cost = computeCost_shared(current_xrs.viewColumn(0), current_pe); 861 862 System.out.println("---- Greedy solution Time " + (1E-9*(System.nanoTime() - initGreedy)) +", Cost : "+ current_cost + ", num links: " + current_pe.zSum() + ", SRGs: " + nSRGs); 863 864 while (System.nanoTime() < algorithmEndtime) 865 { 866 Collections.shuffle(shuffledNodePairs ,rng); 867 final long initLSIteration = System.nanoTime(); 868 869 boolean improvedSolution = false; 870 for (Pair<Node,Node> nodePair : shuffledNodePairs) 871 { 872 Set<Integer> demandsToChange = new HashSet<Integer> (); 873 for (Demand d : np.getNodePairDemands(nodePair.getFirst(), nodePair.getSecond() , true)) 874 demandsToChange.add(d.getIndex()); 875 DoubleMatrix2D thisNeighbor_xrs = current_xrs.copy(); // do not modify the original one 876 if (jomStep_restoration(demandsToChange , null , thisNeighbor_xrs) == null) throw new RuntimeException ("Returned null"); 877 DoubleMatrix1D thisNeighbor_pe = computePe_restoration(thisNeighbor_xrs); 878 final double thisNeighborCost = computeCost_shared (thisNeighbor_xrs.viewColumn(0) , thisNeighbor_pe); 879 if (thisNeighborCost < current_cost) 880 { 881 current_xrs = thisNeighbor_xrs; current_pe = thisNeighbor_pe; current_cost = thisNeighborCost; improvedSolution = true; 882 System.out.println("NEW BEST COST: " + current_cost + ", num links: " + current_pe.zSum()); 883 break; 884 } 885 } 886 887 System.out.println("- Local search TIME: " + ((System.nanoTime()-initLSIteration)*1e-9) + ", current cost: " + current_cost + ", num links: " + current_pe.zSum() + ", IMPROVED: " + improvedSolution); 888 889 if (!improvedSolution) break; else stat_numLSIterationsReducingCost ++; // end grasp iteration 890 } 891 892 return current_xrs; 893 } 894 895 private DoubleMatrix1D computePe_11 (DoubleMatrix1D x_r , DoubleMatrix1D x2_r) 896 { 897 DoubleMatrix1D p_e = A_er.zMult(x_r, null); 898 DoubleMatrix1D backup_occup_e = A_er.zMult(x2_r, null); 899 for (int e = 0 ; e < Efm ; e ++) p_e.set(e , p_e.get(e) + backup_occup_e.get(e)); 900 p_e.assign(DoubleFunctions.div(tcfa_linkCapacity_numCirc.getInt ())).assign(DoubleFunctions.ceil); 901 return p_e; 902 } 903 private DoubleMatrix1D computePe_restoration (DoubleMatrix2D x_rs) 904 { 905 DoubleMatrix2D p_es = A_er.zMult(x_rs, null); DoubleMatrix1D p_e = DoubleFactory1D.dense.make(Efm); 906 for (int e = 0; e < Efm ; e ++) p_e.set(e , p_es.viewRow(e).getMaxLocation() [0]); // the maximum per state s 907 p_e.assign(DoubleFunctions.div(tcfa_linkCapacity_numCirc.getInt ())).assign(DoubleFunctions.ceil); 908 return p_e; 909 } 910 911 private double eTime (long init) { return 1E-9*(System.nanoTime() - init); } 912 913 private static int [] getIndexes (Collection<? extends NetworkElement> col) { final int [] res = new int [col.size ()]; int cont = 0; for (NetworkElement el : col) res [cont++] = el.getIndex (); return res; } 914 915 916}