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}