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 011package com.net2plan.examples.ocnbook.offline; 012 013import java.io.File; 014import java.util.List; 015import java.util.Map; 016 017import cern.colt.matrix.tdouble.DoubleFactory1D; 018import cern.colt.matrix.tdouble.DoubleFactory2D; 019import cern.colt.matrix.tdouble.DoubleMatrix1D; 020import cern.colt.matrix.tdouble.DoubleMatrix2D; 021import cern.jet.math.tdouble.DoubleFunctions; 022import cern.jet.math.tdouble.DoublePlusMultSecond; 023 024import com.jom.OptimizationProblem; 025import com.net2plan.interfaces.networkDesign.Demand; 026import com.net2plan.interfaces.networkDesign.IAlgorithm; 027import com.net2plan.interfaces.networkDesign.Link; 028import com.net2plan.interfaces.networkDesign.Net2PlanException; 029import com.net2plan.interfaces.networkDesign.NetPlan; 030import com.net2plan.interfaces.networkDesign.Node; 031import com.net2plan.libraries.GraphUtils; 032import com.net2plan.utils.Constants.RoutingType; 033import com.net2plan.utils.InputParameter; 034import com.net2plan.utils.Quintuple; 035import com.net2plan.utils.TimeTrace; 036import com.net2plan.utils.Triple; 037 038/** 039 * Finds the routing and mocular capacities for a network that minimize the cost, using a dual decomposition approach 040 * 041 * The time evolution of different metrics can be stored in output files, for later processing. 042 * As an example, see the <a href="../../../../../../graphGeneratorFiles/fig_sec11_7_modularCapacitiesAndRouting_dualDecomp.m">{@code fig_sec11_7_modularCapacitiesAndRouting_dualDecomp.m}</a> MATLAB file used for generating the graph/s of the case study in the 043 * <a href="http://eu.wiley.com/WileyCDA/WileyTitle/productCd-1119013356.html">book</a> using this algorithm. 044 * @net2plan.description 045 * @net2plan.keywords Capacity assignment (CA), Flow assignment (FA), Modular capacities , Dual decomposition 046 * @net2plan.ocnbooksections Section 11.7 047 * @net2plan.inputParameters 048 * @author Pablo Pavon-Marino 049 */ 050public class Offline_cfa_modularCapacitiesAndRoutingDualDecomposition implements IAlgorithm 051{ 052 private double PRECISIONFACTOR; 053 private double numModulesUpperBound; 054 055 private int N , E; 056 057 private InputParameter moduleCapacity = new InputParameter ("moduleCapacity", 1.0 , "Capacity of one module" , 0 , false , Double.MAX_VALUE , true); 058 private InputParameter solverName = new InputParameter ("solverName", "#select# ipopt glpk cplex", "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"); 059 private InputParameter solverLibraryName = new InputParameter ("solverLibraryName", "" , "The solver library full or relative path, to be used by JOM. Leave blank to use JOM default."); 060 private InputParameter maxSolverTimeInSeconds = new InputParameter ("maxSolverTimeInSeconds", (double) -1 , "Maximum time granted to the solver to solve the problem. If this time expires, the solver returns the best solution found so far (if a feasible solution is found)"); 061 private InputParameter subgrad_gammaStep = new InputParameter ("subgrad_gammaStep", (double) 0.05 , "Gamma step in the algorithm iteration" , 0 , false , Double.MAX_VALUE , true); 062 private InputParameter subgrad_type = new InputParameter ("subgrad_type", "#select# 1-over-t constant 1-over-square-root-t", "Type of gradient algorithm. constant, 1-over-t, decreasing-square-t"); 063 private InputParameter simulation_outFileNameRoot = new InputParameter ("simulation_outFileNameRoot", "modularCapacitiesAndRoutingDualDecomp" , "Root of the file name to be used in the output files. If blank, no output"); 064 private InputParameter simulation_numIterations = new InputParameter ("simulation_numIterations", (int) 500 , "Number o iterations of the algorithm" , 1 , Integer.MAX_VALUE); 065 066 067 private DoubleMatrix2D tmColumnSumZero; 068 private DoubleMatrix2D A_ne; 069 070 private TimeTrace stat_bestObjFuncFeasible; 071 private TimeTrace stat_objFuncFeasible; 072 private TimeTrace stat_lowerBound; 073 private TimeTrace stat_bestLowerBound; 074 private TimeTrace stat_pie; 075 076 @Override 077 public String executeAlgorithm(NetPlan netPlan, Map<String, String> algorithmParameters, Map<String, String> net2planParameters) 078 { 079 /* Initialize all InputParameter objects defined in this object (this uses Java reflection) */ 080 InputParameter.initializeAllInputParameterFieldsOfObject(this, algorithmParameters); 081 082 if (netPlan.getNumberOfLayers() != 1) throw new Net2PlanException ("This algorithm works in single layer networks"); 083 084 /* Initialize some variables */ 085 this.N = netPlan.getNumberOfNodes(); 086 this.E = netPlan.getNumberOfLinks(); 087 this.PRECISIONFACTOR = Double.parseDouble(net2planParameters.get("precisionFactor")); 088 089 this.stat_objFuncFeasible = new TimeTrace (); 090 this.stat_bestObjFuncFeasible = new TimeTrace (); 091 this.stat_lowerBound = new TimeTrace (); 092 this.stat_bestLowerBound = new TimeTrace (); 093 this.stat_pie = new TimeTrace (); 094 095 netPlan.removeAllUnicastRoutingInformation(); 096 netPlan.setRoutingType (RoutingType.HOP_BY_HOP_ROUTING); 097 098 /* The diagonal in the traffic matrix contains minus the amount of traffic destined to that node */ 099 double [][] trafficMatrix = netPlan.getMatrixNode2NodeOfferedTraffic().toArray(); 100 for (int n1 = 0 ; n1 < N ; n1 ++) for (int n2 = 0 ; n2 < N ; n2 ++) if (n1 == n2) continue; else trafficMatrix [n2][n2] -= trafficMatrix [n1][n2]; 101 this.tmColumnSumZero= DoubleFactory2D.dense.make(trafficMatrix); 102 103 this.numModulesUpperBound = Math.ceil(netPlan.getVectorDemandOfferedTraffic().zSum() / moduleCapacity.getDouble()); 104 this.A_ne = netPlan.getMatrixNodeLinkIncidence(); 105 106 DoubleMatrix1D multipliers_eIP = DoubleFactory1D.dense.make(E, 1.0); 107 double highestDualCost = -Double.MAX_VALUE; 108 double lowestPrimalCost = Double.MAX_VALUE; 109 DoubleMatrix1D bestFeasible_ne = null; 110 DoubleMatrix2D bestRelaxedAndFeasible_xte = null; 111 for (int it = 1 ; it <= simulation_numIterations.getInt() ; it ++) 112 { 113 System.out.println("**** Iteration : " + it + ", mult: " + multipliers_eIP); 114 115 Quintuple<DoubleMatrix2D,DoubleMatrix1D,DoubleMatrix1D,Double,DoubleMatrix1D> q = solveSubproblems (multipliers_eIP); 116 final DoubleMatrix2D relaxedAndFeasible_xte = q.getFirst(); 117 final DoubleMatrix1D relaxed_ne = q.getSecond(); 118 final DoubleMatrix1D feasible_n_e = q.getThird(); 119 final double dualCost = q.getFourth(); 120 final DoubleMatrix1D gradient_e = q.getFifth (); 121 122 final double objFunc = feasible_n_e.zSum(); 123 124 highestDualCost = Math.max(highestDualCost , dualCost); 125 if (objFunc < lowestPrimalCost) 126 { 127 lowestPrimalCost = objFunc; 128 bestFeasible_ne = feasible_n_e.copy(); 129 bestRelaxedAndFeasible_xte = relaxedAndFeasible_xte.copy(); 130 } 131 132 if (highestDualCost > lowestPrimalCost + PRECISIONFACTOR) throw new RuntimeException ("Bad: highestDualCost: "+ highestDualCost + ", lowestPrimalCost: " + lowestPrimalCost); 133 134 System.out.println("* Feasible: Compute cost: " + objFunc); 135 136 this.stat_objFuncFeasible.add (it , new double [] { objFunc } ); 137 this.stat_bestObjFuncFeasible.add(it , new double [] { lowestPrimalCost } ); 138 this.stat_lowerBound.add(it , new double [] { dualCost } ); 139 this.stat_bestLowerBound.add(it , new double [] { highestDualCost } ); 140 this.stat_pie.add(it, multipliers_eIP.toArray() ); 141 142 double gamma = 0; 143 if (this.subgrad_type.getString ().equalsIgnoreCase("constant")) 144 gamma = this.subgrad_gammaStep.getDouble(); 145 else if (this.subgrad_type.getString ().equalsIgnoreCase("1-over-t")) 146 gamma = this.subgrad_gammaStep.getDouble() / it; 147 else if (this.subgrad_type.getString ().equalsIgnoreCase("1-over-square-root-t")) 148 gamma = this.subgrad_gammaStep.getDouble() / Math.sqrt(it); 149 else throw new Net2PlanException ("Unknown subgradient algorithm type"); 150 for (int e = 0 ; e < E ; e++) 151 multipliers_eIP.set(e , Math.max(0 , multipliers_eIP.get(e) + gamma * gradient_e.get(e) )); 152 } 153 154 saveNetPlan (netPlan , bestRelaxedAndFeasible_xte , bestFeasible_ne); 155 156 finish (netPlan); 157 158 return "Ok! Best solution. Cost: " + bestFeasible_ne.zSum(); 159 } 160 161 @Override 162 public String getDescription() 163 { 164 return "Given a network with a set of given nodes, and links, and a given known unicast offered traffic. This algorithm jointly computes (i) the routing of the traffic, and the capacity to assign to the links. The link capacities are constrained to be modular: an integer multiple of a link capacity quantum. Optimization target is minimizing the network cost, given by the number of capacity modules to install in the network. The problem in NP-hard. The algorithm implements a dual decomposition iterative, using a subgradient algorithm."; 165 } 166 167 @Override 168 public List<Triple<String, String, String>> getParameters() 169 { 170 /* Returns the parameter information for all the InputParameter objects defined in this object (uses Java reflection) */ 171 return InputParameter.getInformationAllInputParameterFieldsOfObject(this); 172 } 173 174 private Quintuple<DoubleMatrix2D,DoubleMatrix1D,DoubleMatrix1D,Double,DoubleMatrix1D> solveSubproblems (DoubleMatrix1D pi_e) 175 { 176 /* Upper layer */ 177 DoubleMatrix2D x_te = DoubleFactory2D.dense.make(N,E); 178 { 179 OptimizationProblem op = new OptimizationProblem(); 180 op.setInputParameter("pi_e", pi_e.toArray() , "column"); 181 op.setInputParameter("A_ne", A_ne.toArray()); 182 op.setInputParameter("TM", this.tmColumnSumZero.toArray()); 183 op.addDecisionVariable("x_te", false , new int[] { N , E }, 0, Double.MAX_VALUE); // dest-link formulation at IP layer 184 op.setObjectiveFunction("minimize", "sum (x_te * pi_e)"); 185 op.addConstraint("A_ne * (x_te') == TM"); // the flow-conservation constraints (NxN constraints) 186 /* Call the solver to solve the problem */ 187 op.solve(solverName.getString (), "solverLibraryName", solverLibraryName.getString () , "maxSolverTimeInSeconds" , maxSolverTimeInSeconds.getDouble ()); 188 189 /* If no solution is found, quit */ 190 if (op.feasibleSolutionDoesNotExist()) throw new Net2PlanException("The problem has no feasible solution"); 191 if (op.foundUnboundedSolution()) throw new Net2PlanException ("Found an unbounded solution"); 192 if (!op.solutionIsFeasible()) throw new Net2PlanException("A feasible solution was not found"); 193 x_te = op.getPrimalSolution("x_te").view2D(); 194 } 195 196 DoubleMatrix1D n_e = DoubleFactory1D.dense.make(E); 197 for (int e = 0 ; e < E ; e ++) 198 n_e.set(e , (1 - (moduleCapacity.getDouble() * pi_e.get(e)) >= 0)? 0 : numModulesUpperBound); 199 200 DoubleMatrix1D y_e = x_te.copy().zMult(DoubleFactory1D.dense.make(N, 1.0) , null , 1.0 , 0.0 , true); 201 DoubleMatrix1D gradient_e = y_e.copy().assign(n_e , DoublePlusMultSecond.plusDiv(-moduleCapacity.getDouble())); 202 203 /* Compute the dual cost */ 204 final double dualCost = x_te.copy().zMult(pi_e, null).zSum() - n_e.zDotProduct(pi_e) * moduleCapacity.getDouble(); 205 206 /* Compute a feasible solution from the relaxed x_te */ 207 DoubleMatrix1D feasible_n_e = y_e.copy().assign(DoubleFunctions.div(moduleCapacity.getDouble())).assign(DoubleFunctions.ceil); 208 209// System.out.println("carriedTraffic_e: " + carriedTraffic_e); 210// System.out.println("capacityLink_e: " + capacityLink_e); 211// System.out.println("barU_e: " + barU_e); 212 213 return Quintuple.of(x_te.copy () , n_e.copy() , feasible_n_e.copy() , dualCost , gradient_e); 214 } 215 216 private void finish (NetPlan np) 217 { 218 /* If no output file, return */ 219 if (simulation_outFileNameRoot.getString ().equals("")) return; 220 stat_objFuncFeasible.printToFile(new File (simulation_outFileNameRoot.getString () + "_objFunc.txt")); 221 stat_bestObjFuncFeasible.printToFile(new File (simulation_outFileNameRoot.getString () + "_bestObjFunc.txt")); 222 stat_lowerBound.printToFile(new File (simulation_outFileNameRoot.getString () + "_lowerBound.txt")); 223 stat_bestLowerBound.printToFile(new File (simulation_outFileNameRoot.getString () + "_bestLowerBound.txt")); 224 stat_pie.printToFile(new File (simulation_outFileNameRoot.getString () + "_pie.txt")); 225 226 227 double lowerBoundNumLps_1 = 0; 228 for (Node s : np.getNodes()) 229 { 230 double outTraffic = 0; for (Demand d : s.getOutgoingDemands()) outTraffic += d.getOfferedTraffic(); 231 lowerBoundNumLps_1 += Math.ceil(outTraffic / moduleCapacity.getDouble()); 232 } 233 double lowerBoundNumLps_2 = 0; 234 for (Node t : np.getNodes()) 235 { 236 double inTraffic = 0; for (Demand d : t.getIncomingDemands()) inTraffic += d.getOfferedTraffic(); 237 lowerBoundNumLps_2 += Math.ceil(inTraffic / moduleCapacity.getDouble()); 238 } 239 double lowerBoundNumLps_3 = 0; 240 for (Node t : np.getNodes()) 241 for (Node s : np.getNodes()) 242 { 243 if (s == t) continue; 244 final int numHopsSP = GraphUtils.getShortestPath(np.getNodes () , np.getLinks(), s, t, null).size(); 245 for (Demand d : np.getNodePairDemands(s,t,false)) lowerBoundNumLps_3 += numHopsSP * d.getOfferedTraffic(); 246 } 247 lowerBoundNumLps_3 = Math.ceil(lowerBoundNumLps_3/moduleCapacity.getDouble()); 248 249 TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString () + "_heuristicLBCost.txt") , new double [] { lowerBoundNumLps_1 , lowerBoundNumLps_2 , lowerBoundNumLps_3 } ); 250 251 } 252 253 void saveNetPlan (NetPlan netPlan , DoubleMatrix2D x_te , DoubleMatrix1D n_e) 254 { 255 /* Set the routing at the IP layer */ 256 netPlan.setRoutingType(RoutingType.HOP_BY_HOP_ROUTING); 257 netPlan.removeAllForwardingRules(); 258 netPlan.setRoutingFromDestinationLinkCarriedTraffic(x_te , true); 259 for (Link e : netPlan.getLinks()) 260 e.setCapacity(moduleCapacity.getDouble() * n_e.get(e.getIndex ())); 261 for (Demand d : netPlan.getDemandsBlocked()) 262 if (d.getBlockedTraffic() > PRECISIONFACTOR) throw new RuntimeException ("Bad"); 263 for (Link e : netPlan.getLinksOversubscribed()) 264 if (e.getOccupiedCapacityIncludingProtectionSegments() - e.getCapacity() > PRECISIONFACTOR) throw new RuntimeException ("Bad"); 265 } 266 267}