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 ******************************************************************************/
008package com.net2plan.examples.ocnbook.onlineSim;
009
010
011
012
013import java.io.File;
014import java.io.FileOutputStream;
015import java.io.PrintStream;
016import java.util.List;
017import java.util.Map;
018import java.util.Random;
019import java.util.Set;
020
021import cern.colt.matrix.tdouble.DoubleFactory1D;
022import cern.colt.matrix.tdouble.DoubleFactory2D;
023import cern.colt.matrix.tdouble.DoubleMatrix1D;
024import cern.colt.matrix.tdouble.DoubleMatrix2D;
025
026import com.jom.OptimizationProblem;
027import com.net2plan.interfaces.networkDesign.Demand;
028import com.net2plan.interfaces.networkDesign.Link;
029import com.net2plan.interfaces.networkDesign.Net2PlanException;
030import com.net2plan.interfaces.networkDesign.NetPlan;
031import com.net2plan.interfaces.networkDesign.Route;
032import com.net2plan.interfaces.simulation.IEventProcessor;
033import com.net2plan.interfaces.simulation.SimEvent;
034import com.net2plan.libraries.GraphUtils;
035import com.net2plan.libraries.NetworkPerformanceMetrics;
036import com.net2plan.utils.Constants.RoutingType;
037import com.net2plan.utils.GradientProjectionUtils;
038import com.net2plan.utils.InputParameter;
039import com.net2plan.utils.Pair;
040import com.net2plan.utils.Quadruple;
041import com.net2plan.utils.TimeTrace;
042import com.net2plan.utils.Triple;
043
044/** 
045 * This module implements a distributed primal-gradient based algorithm using a barrier function, for adapting the demand injected traffic (congestion control) in the network, to maximize the network utility enforcing a fair allocation of the resources.
046 *
047 * Ths event processor is adapted to permit observing the algorithm performances under user-defined conditions, 
048 * including asynchronous distributed executions, where signaling can be affected by losses and/or delays, and/or measurement errors. 
049 * The time evolution of different metrics can be stored in output files, for later processing. 
050 * As an example, see the <a href="../../../../../../graphGeneratorFiles/fig_sec9_4_congestionControlPrimal.m">{@code fig_sec9_4_congestionControlPrimal.m}</a> MATLAB file used for generating the graph/s of the case study in the 
051 * <a href="http://eu.wiley.com/WileyCDA/WileyTitle/productCd-1119013356.html">book</a> using this algorithm.
052 * 
053 * To simulate a network with this module, use the {@code Online_evGen_doNothing} generator.
054 * 
055 * @net2plan.keywords Bandwidth assignment (BA), Distributed algorithm, Primal gradient algorithm
056 * @net2plan.ocnbooksections Section 9.4
057 * @net2plan.inputParameters 
058 * @author Pablo Pavon-Marino
059 */
060public class Online_evProc_congestionControlPrimal extends IEventProcessor 
061{
062        private InputParameter signaling_isSynchronous = new InputParameter ("signaling_isSynchronous", false , "true if all the distributed agents involved wake up synchronously to send the signaling messages");
063        private InputParameter signaling_averageInterMessageTime = new InputParameter ("signaling_averageInterMessageTime", 1.0 , "Average time between two signaling messages sent by an agent" , 0 , false , Double.MAX_VALUE , true);
064        private InputParameter signaling_maxFluctuationInterMessageTime = new InputParameter ("signaling_maxFluctuationInterMessageTime", 0.5 , "Max fluctuation in time between two signaling messages sent by an agent" , 0 , true , Double.MAX_VALUE , true);
065        private InputParameter signaling_averageDelay = new InputParameter ("signaling_averageDelay", 0.0 , "Average time between signaling message transmission by an agent and its reception by other or others" , 0 , true , Double.MAX_VALUE , true);
066        private InputParameter signaling_maxFluctuationInDelay = new InputParameter ("signaling_maxFluctuationInDelay", 0.0 , "Max fluctuation in time in the signaling delay, in absolute time values. The signaling delays are sampled from a uniform distribution within the given interval" , 0 , true , Double.MAX_VALUE , true);
067        private InputParameter signaling_signalingLossProbability = new InputParameter ("signaling_signalingLossProbability", 0.05 , "Probability that a signaling message transmitted is lost (not received by other or others involved agents)" , 0 , true , Double.MAX_VALUE , true);
068        private InputParameter update_isSynchronous = new InputParameter ("update_isSynchronous", false , "true if all the distributed agents involved wake up synchronousely to update its state");
069        private InputParameter update_averageInterUpdateTime = new InputParameter ("update_averageInterUpdateTime", 1.0 , "Average time between two updates of an agent" , 0 , false , Double.MAX_VALUE , true);
070        private InputParameter update_maxFluctuationInterUpdateTime = new InputParameter ("update_maxFluctuationInterUpdateTime", 0.5 , "Max fluctuation in time in the update interval of an agent, in absolute time values. The update intervals are sampled from a uniform distribution within the given interval" , 0 , true , Double.MAX_VALUE , true);
071        private InputParameter gradient_diagonalScaling = new InputParameter ("gradient_diagonalScaling", false , "Whether diagonal scaling is applied or not");
072        private InputParameter gradient_maxGradientAbsoluteNoise = new InputParameter ("gradient_maxGradientAbsoluteNoise", 0.0 , "Max value of the added noise to the gradient coordinate in absolute values" , 0 , true , Double.MAX_VALUE , true);
073        private InputParameter gradient_gammaStep = new InputParameter ("gradient_gammaStep", 20.0 , "Gamma step in the gradient algorithm" , 0 , false , Double.MAX_VALUE , true);
074        private InputParameter gradient_heavyBallBetaParameter = new InputParameter ("gradient_heavyBallBetaParameter", 0.0 , "Beta parameter of heavy ball, between 0 and 1. Value 0 means no heavy ball" , 0 , true , 1.0 , true);
075        private InputParameter gradient_maxGradientCoordinateChange = new InputParameter ("gradient_maxGradientCoordinateChange", 1.0 , "Maximum change in an iteration of a gradient coordinate" , 0 , false , Double.MAX_VALUE , true);
076        private InputParameter gradient_penaltyMethod = new InputParameter ("gradient_penaltyMethod", "#select# barrier exterior-penalty" , "Use a barrier (interior penalty) or an exterior penalty for enforcing link capacity constraints");
077        private InputParameter gradient_exteriorPenaltyMuFactor = new InputParameter ("gradient_exteriorPenaltyMuFactor", 1000 , "Mu factor to apply in the exterior penalty method" , 0 , true , Double.MAX_VALUE , true);
078        private InputParameter gradient_interiorPenaltyEpsilonFactor = new InputParameter ("gradient_interiorPenaltyEpsilonFactor", 0.00001 , "Factor to multiply the delay penalization function" , 0 , true , Double.MAX_VALUE , true);
079        private InputParameter simulation_maxNumberOfUpdateIntervals = new InputParameter ("simulation_maxNumberOfUpdateIntervals", 800.0 , "Maximum number of update intervals in average per agent" , 0 , false , Double.MAX_VALUE , true);
080        private InputParameter simulation_randomSeed = new InputParameter ("simulation_randomSeed", (long) 1 , "Seed of the random number generator");
081        private InputParameter simulation_outFileNameRoot = new InputParameter ("simulation_outFileNameRoot", "congestionControlPrimal" , "Root of the file name to be used in the output files. If blank, no output");
082
083        private InputParameter control_minHd = new InputParameter ("control_minHd", 0.1 , "Minimum traffic assigned to each demand" , 0 , true , Double.MAX_VALUE , true);
084        private InputParameter control_maxHd = new InputParameter ("control_maxHd", 1.0E6 , "Maximum traffic assigned to each demand" , 0 , true , Double.MAX_VALUE , true);
085        private InputParameter control_fairnessFactor = new InputParameter ("control_fairnessFactor", 2.0 , "Fairness factor in utility function of congestion control" , 0 , true , Double.MAX_VALUE , true);
086        
087        private Random rng;
088        
089        private static final int SIGNALING_WAKEUPTOSENDMESSAGE = 300;
090        private static final int SIGNALING_RECEIVEDMESSAGE = 301;
091        private static final int UPDATE_WAKEUPTOUPDATE = 302;
092
093        private NetPlan currentNetPlan;
094        
095        private double control_epsilonOrMuFactor;
096        private DoubleMatrix1D control_priceFirstOrder_e;
097        private DoubleMatrix1D control_priceSecondOrder_e;
098        private DoubleMatrix1D control_previous_h_d;
099        private DoubleMatrix2D control_mostUpdatedLinkFirstOrderPriceKnownDemand_de; 
100        private DoubleMatrix2D control_mostUpdatedLinkSecondOrderPriceKnownDemand_de; 
101        private boolean control_isBarrierMethod;
102        private int D , E;
103        
104        private TimeTrace stat_traceOf_objFunction;
105        private TimeTrace stat_traceOf_hd;
106        private TimeTrace stat_traceOf_maxLinkTraffic;
107
108        @Override
109        public String getDescription()
110        {
111                return "This module implements a distributed primal-gradient based algorithm using a barrier function, for adapting the demand injected traffic (congestion control) in the network, to maximize the network utility enforcing a fair allocation of the resources.";
112        }
113
114        @Override
115        public List<Triple<String, String, String>> getParameters()
116        {
117                /* Returns the parameter information for all the InputParameter objects defined in this object (uses Java reflection) */
118                return InputParameter.getInformationAllInputParameterFieldsOfObject(this);
119        }
120
121        @Override
122        public void initialize(NetPlan currentNetPlan, Map<String, String> algorithmParameters, Map<String, String> simulationParameters, Map<String, String> net2planParameters)
123        {
124                /* Initialize all InputParameter objects defined in this object (this uses Java reflection) */
125                InputParameter.initializeAllInputParameterFieldsOfObject(this, algorithmParameters);
126
127                this.currentNetPlan = currentNetPlan;
128                if (currentNetPlan.getNumberOfLayers() != 1) throw new Net2PlanException ("This algorithm works in single layer networks");
129
130                /* Remove all routes, and create one with the shortest path in km for each demand */
131                currentNetPlan.removeAllUnicastRoutingInformation();
132                currentNetPlan.setRoutingType(RoutingType.SOURCE_ROUTING);
133                currentNetPlan.addRoutesFromCandidatePathList(currentNetPlan.getVectorLinkLengthInKm().toArray()  , "K" , "1");
134
135                this.control_isBarrierMethod = gradient_penaltyMethod.getString ().equals ("barrier");
136                this.control_epsilonOrMuFactor = control_isBarrierMethod? gradient_interiorPenaltyEpsilonFactor.getDouble() : gradient_exteriorPenaltyMuFactor.getDouble();
137                this.rng = new Random (simulation_randomSeed.getLong());
138                this.D = currentNetPlan.getNumberOfDemands ();
139                this.E = currentNetPlan.getNumberOfLinks ();
140                if ((E == 0) || (D == 0)) throw new Net2PlanException ("The input design should have links and demands");
141                
142                /* Initially all demands offer exactly hdMin */
143                control_previous_h_d = DoubleFactory1D.dense.make (D , control_minHd.getDouble());
144                currentNetPlan.setVectorDemandOfferedTraffic(control_previous_h_d);
145                for (Demand d : currentNetPlan.getDemands()) // carry the demand traffic
146                { 
147                        final Set<Route> routes = d.getRoutes(); 
148                        for (Route r : routes) r.setCarriedTraffic(d.getOfferedTraffic() / routes.size () , d.getOfferedTraffic() / routes.size ()); 
149                }
150                
151                /* Set the initial prices in the links */
152                this.control_priceFirstOrder_e = DoubleFactory1D.dense.make (E); for (Link e : this.currentNetPlan.getLinks ()) control_priceFirstOrder_e.set (e.getIndex () , computeFirstOrderPriceFromNetPlan(e));
153                this.control_priceSecondOrder_e = DoubleFactory1D.dense.make (E); for (Link e : this.currentNetPlan.getLinks ()) control_priceSecondOrder_e.set (e.getIndex () , computeSecondOrderPriceFromNetPlan(e));
154                
155                /* Initialize the information each demand knows of the prices of all the links */
156                this.control_mostUpdatedLinkFirstOrderPriceKnownDemand_de = DoubleFactory2D.dense.make (D,E);
157                this.control_mostUpdatedLinkSecondOrderPriceKnownDemand_de = DoubleFactory2D.dense.make (D,E);
158                for (Demand d : currentNetPlan.getDemands ())
159                {
160                        control_mostUpdatedLinkFirstOrderPriceKnownDemand_de.viewRow (d.getIndex ()).assign (control_priceFirstOrder_e);
161                        control_mostUpdatedLinkSecondOrderPriceKnownDemand_de.viewRow (d.getIndex ()).assign (control_priceSecondOrder_e);
162                }
163                
164                /* Initially all nodes receive a "wake up to transmit" event, aligned at time zero or y asynchr => randomly chosen */
165                for (Link e : currentNetPlan.getLinks())
166                {
167                        final double signalingTime = (signaling_isSynchronous.getBoolean())? signaling_averageInterMessageTime.getDouble() : Math.max(0 , signaling_averageInterMessageTime.getDouble() + signaling_maxFluctuationInterMessageTime.getDouble() * (rng.nextDouble() - 0.5));
168                        this.scheduleEvent(new SimEvent (signalingTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_WAKEUPTOSENDMESSAGE , e));
169                }
170                for (Demand d : currentNetPlan.getDemands())
171                {
172                        final double updateTime = (update_isSynchronous.getBoolean())? update_averageInterUpdateTime.getDouble() : Math.max(0 , update_averageInterUpdateTime.getDouble() + update_maxFluctuationInterUpdateTime.getDouble() * (rng.nextDouble() - 0.5));
173                        this.scheduleEvent(new SimEvent (updateTime , SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_WAKEUPTOUPDATE , d));
174                }
175
176                /* Intialize the traces */
177                this.stat_traceOf_hd = new TimeTrace ();
178                this.stat_traceOf_objFunction = new TimeTrace (); 
179                this.stat_traceOf_maxLinkTraffic = new TimeTrace ();
180                this.stat_traceOf_hd.add(0 , this.currentNetPlan.getVectorDemandOfferedTraffic());
181                this.stat_traceOf_objFunction.add(0 , NetworkPerformanceMetrics.alphaUtility(currentNetPlan.getVectorDemandOfferedTraffic() , control_fairnessFactor.getDouble()));
182                this.stat_traceOf_maxLinkTraffic.add(0.0, this.currentNetPlan.getVectorLinkTotalCarriedTraffic().getMaxLocation() [0]);
183        }
184
185        @Override
186        public void processEvent(NetPlan currentNetPlan, SimEvent event)
187        {
188                final double t = event.getEventTime();
189                switch (event.getEventType())
190                {
191                case SIGNALING_RECEIVEDMESSAGE: // A node receives from an out neighbor the q_nt for any destination
192                {
193                        final Quadruple<Demand,Link,Double,Double> signalInfo = (Quadruple<Demand,Link,Double,Double>) event.getEventObject();
194                        final Demand dMe = signalInfo.getFirst();
195                        final Link e = signalInfo.getSecond ();
196                        control_mostUpdatedLinkFirstOrderPriceKnownDemand_de.set (dMe.getIndex () , e.getIndex () , signalInfo.getThird ());
197                        control_mostUpdatedLinkSecondOrderPriceKnownDemand_de.set (dMe.getIndex () , e.getIndex () , signalInfo.getFourth ());
198                        break;
199                }
200                
201                case SIGNALING_WAKEUPTOSENDMESSAGE: // A node broadcasts signaling info to its 1 hop neighbors
202                {
203                        final Link eMe = (Link) event.getEventObject();
204
205                        /* Send the events of the signaling information messages to all the nodes */
206                        if (rng.nextDouble() >= this.signaling_signalingLossProbability.getDouble()) // the signaling may be lost => lost to all demands
207                                for (Route route : eMe.getTraversingRoutes())
208                                {
209                                        final Demand d = route.getDemand();
210                                        Quadruple<Demand,Link,Double,Double> infoToSignal = Quadruple.of(d , eMe ,  this.computeFirstOrderPriceFromNetPlan(eMe) , computeSecondOrderPriceFromNetPlan(eMe));
211                                        final double signalingReceptionTime = t + Math.max(0 , signaling_averageDelay.getDouble() + signaling_maxFluctuationInDelay.getDouble() * (rng.nextDouble() - 0.5));
212                                        this.scheduleEvent(new SimEvent (signalingReceptionTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_RECEIVEDMESSAGE , infoToSignal));
213                                }
214                        
215                        /* Re-schedule when to wake up again */
216                        final double signalingTime = signaling_isSynchronous.getBoolean()? t + signaling_averageInterMessageTime.getDouble() : Math.max(t , t + signaling_averageInterMessageTime.getDouble() + signaling_maxFluctuationInterMessageTime.getDouble() * (rng.nextDouble() - 0.5));
217                        this.scheduleEvent(new SimEvent (signalingTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_WAKEUPTOSENDMESSAGE , eMe));
218                        break;
219                }
220
221                case UPDATE_WAKEUPTOUPDATE: 
222                {
223                        final Demand dMe = (Demand) event.getEventObject();
224                        
225                        DoubleMatrix1D infoIKnow_priceFirstOrder_e = this.control_mostUpdatedLinkFirstOrderPriceKnownDemand_de.viewRow (dMe.getIndex ());
226                        DoubleMatrix1D infoIKnow_priceSecondOrder_e = this.control_mostUpdatedLinkSecondOrderPriceKnownDemand_de.viewRow (dMe.getIndex ());
227
228                        /* compute the demand price as weighted sum in the routes of route prices  */
229                        double demandWeightedSumLinkPrices = 0;
230                        double demandWeightedSumSecondDerivativeLinkPrices = 0;
231                        double demandCarriedTraffic = 0; 
232                        for (Route r : dMe.getRoutes ())
233                        {
234                                final double h_r = r.getCarriedTraffic();
235                                demandCarriedTraffic += h_r;
236                                for (Link e : r.getSeqLinksRealPath())
237                                {
238                                        demandWeightedSumLinkPrices += h_r * infoIKnow_priceFirstOrder_e.get(e.getIndex ());
239                                        demandWeightedSumSecondDerivativeLinkPrices += h_r * infoIKnow_priceSecondOrder_e.get(e.getIndex ());
240                                }
241                        }
242                        if (Math.abs(demandCarriedTraffic - dMe.getOfferedTraffic()) > 1E-3) throw new RuntimeException ("Not all the traffic is carried. demandCarriedTraffic: " + demandCarriedTraffic);
243                        demandWeightedSumLinkPrices /= demandCarriedTraffic;
244                        demandWeightedSumSecondDerivativeLinkPrices /= demandCarriedTraffic;
245                        
246                        /* compute the new h_d */
247                        final double old_hd = dMe.getOfferedTraffic();
248                        final double gradient_d = Math.pow(old_hd, -this.control_fairnessFactor.getDouble()) - this.control_epsilonOrMuFactor * demandWeightedSumLinkPrices +  2*gradient_maxGradientAbsoluteNoise.getDouble()*(rng.nextDouble()-0.5);
249                        double seconDerivativehd2ForDiagonalScaling = (gradient_diagonalScaling.getBoolean())? Math.abs(-this.control_fairnessFactor.getDouble() * Math.pow(old_hd, -this.control_fairnessFactor.getDouble()-1)- this.control_epsilonOrMuFactor * demandWeightedSumSecondDerivativeLinkPrices) : 1;
250                        if (!Double.isFinite(seconDerivativehd2ForDiagonalScaling)) seconDerivativehd2ForDiagonalScaling = 1; 
251                        
252                        double coordinateChange =  this.gradient_gammaStep.getDouble() / seconDerivativehd2ForDiagonalScaling * gradient_d + this.gradient_heavyBallBetaParameter.getDouble() * (old_hd - this.control_previous_h_d.get(dMe.getIndex ())  );
253                        
254                        if (gradient_maxGradientCoordinateChange.getDouble() > 0) 
255                                coordinateChange = Math.signum(coordinateChange) * Math.min(gradient_maxGradientCoordinateChange.getDouble() , Math.abs(coordinateChange));  
256
257                        final double new_hd = GradientProjectionUtils.euclideanProjection_boxLike (old_hd + coordinateChange , control_minHd.getDouble() , control_maxHd.getDouble ());
258
259                        control_previous_h_d.set (dMe.getIndex (), dMe.getOfferedTraffic());
260                        dMe.setOfferedTraffic(new_hd); 
261                        final Set<Route> routes = dMe.getRoutes(); // carry the demand traffic
262                        for (Route r : routes) r.setCarriedTraffic(dMe.getOfferedTraffic() / routes.size () , dMe.getOfferedTraffic() / routes.size ());
263
264                        final double updateTime = update_isSynchronous.getBoolean()? t + update_averageInterUpdateTime.getDouble() : Math.max(t , t + update_averageInterUpdateTime.getDouble() + update_maxFluctuationInterUpdateTime.getDouble() * (rng.nextDouble() - 0.5));
265                        this.scheduleEvent(new SimEvent (updateTime , SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_WAKEUPTOUPDATE,  dMe));
266
267                        this.stat_traceOf_hd.add(t, this.currentNetPlan.getVectorDemandOfferedTraffic());
268                        this.stat_traceOf_objFunction.add(t, NetworkPerformanceMetrics.alphaUtility(currentNetPlan.getVectorDemandOfferedTraffic() , control_fairnessFactor.getDouble()));
269                        this.stat_traceOf_maxLinkTraffic.add(t, this.currentNetPlan.getVectorLinkTotalCarriedTraffic().getMaxLocation() [0]);
270
271                        if (t > this.simulation_maxNumberOfUpdateIntervals.getDouble() * this.update_averageInterUpdateTime.getDouble()) { this.endSimulation (); }
272                        
273                        break;
274                }
275                        
276
277                default: throw new RuntimeException ("Unexpected received event");
278                }
279                
280                
281        }
282
283        public String finish (StringBuilder st , double simTime)
284        {
285                if (simulation_outFileNameRoot.getString().equals("")) return null;
286                stat_traceOf_hd.printToFile(new File (simulation_outFileNameRoot.getString() + "_hd.txt"));
287                stat_traceOf_objFunction.printToFile(new File (simulation_outFileNameRoot.getString() + "_objFunc.txt"));
288                stat_traceOf_maxLinkTraffic.printToFile(new File (simulation_outFileNameRoot.getString() + "_maxYe.txt"));
289                Pair<DoubleMatrix1D,Double> optPair = computeOptimumSolution ();
290                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_hd.txt"), optPair.getFirst());
291                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_objFunc.txt"), optPair.getSecond());
292                return null;
293        }
294        
295                        
296        private Pair<DoubleMatrix1D,Double> computeOptimumSolution ()
297        {
298                OptimizationProblem op = new OptimizationProblem();
299
300                /* Add the decision variables to the problem */
301                op.addDecisionVariable("h_d", false, new int[] {1, D}, control_minHd.getDouble() , control_maxHd.getDouble());
302
303                /* Set some input parameters */
304                op.setInputParameter("u_e", currentNetPlan.getVectorLinkCapacity() , "row");
305                op.setInputParameter("alpha", this.control_fairnessFactor.getDouble());
306                op.setInputParameter("R_de", currentNetPlan.getMatrixDemand2LinkAssignment());
307
308                /* Sets the objective function */
309                if (control_fairnessFactor.getDouble() == 1)
310                    op.setObjectiveFunction("maximize", "sum(ln(h_d))");
311                else if (control_fairnessFactor.getDouble() == 0)
312                    op.setObjectiveFunction("maximize", "sum(h_d)");
313                else
314                    op.setObjectiveFunction("maximize", "(1-alpha) * sum(h_d ^ (1-alpha))");
315
316                op.addConstraint("h_d * R_de <= u_e"); // the capacity constraints (E constraints)
317
318                /* Call the solver to solve the problem */
319                op.solve("ipopt");
320
321                /* If an optimal solution was not found, quit */
322                if (!op.solutionIsOptimal()) throw new Net2PlanException("An optimal solution was not found");
323
324                /* Retrieve the optimum solutions */
325                DoubleMatrix1D h_d = op.getPrimalSolution("h_d").view1D ();
326                return Pair.of(h_d ,NetworkPerformanceMetrics.alphaUtility(h_d , control_fairnessFactor.getDouble()));
327        }
328                        
329        
330        
331        /* Computes the price, and the price for diagonal scaling */
332        private double computeFirstOrderPriceFromNetPlan (Link e)
333        {
334                final double y_e = e.getCarriedTrafficIncludingProtectionSegments();
335                final double u_e = e.getCapacity();
336                if (u_e == 0) throw new RuntimeException ("Zero capacity in a link");
337                double price = 0;
338                if (control_isBarrierMethod)
339                {
340                        if (y_e / u_e < 1)
341                                price = u_e / Math.pow(u_e - y_e , 2);
342                        else
343                                price = Double.MAX_VALUE; 
344                }
345                else
346                        price = Math.max(0, y_e - u_e);
347                return price;
348        }
349
350        private double computeSecondOrderPriceFromNetPlan (Link e)
351        {
352                final double y_e = e.getCarriedTrafficIncludingProtectionSegments();
353                final double u_e = e.getCapacity();
354                if (u_e == 0) throw new RuntimeException ("Zero capacity in a link");
355                double priceDiagScaling = 0;
356                if (control_isBarrierMethod)
357                {
358                        if (y_e / u_e < 1)
359                                priceDiagScaling = 2 * u_e / Math.pow(u_e - y_e , 3);
360                        else
361                                priceDiagScaling = Double.MAX_VALUE; //2 * u_e / Math.pow(u_e - 0.999999 * u_e , 3);
362                }
363                else
364                        priceDiagScaling = (y_e <= u_e)? 0 : 2;
365                return priceDiagScaling;
366        }
367
368        
369}