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
012import java.io.File;
013import java.util.HashMap;
014import java.util.List;
015import java.util.Map;
016import java.util.Random;
017
018import cern.colt.matrix.tdouble.DoubleFactory1D;
019import cern.colt.matrix.tdouble.DoubleFactory2D;
020import cern.colt.matrix.tdouble.DoubleMatrix1D;
021import cern.colt.matrix.tdouble.DoubleMatrix2D;
022import cern.jet.math.tdouble.DoubleFunctions;
023
024import com.net2plan.examples.ocnbook.offline.Offline_ca_wirelessTransmissionPower;
025import com.net2plan.interfaces.networkDesign.Link;
026import com.net2plan.interfaces.networkDesign.Net2PlanException;
027import com.net2plan.interfaces.networkDesign.NetPlan;
028import com.net2plan.interfaces.simulation.IEventProcessor;
029import com.net2plan.interfaces.simulation.SimEvent;
030import com.net2plan.libraries.NetworkPerformanceMetrics;
031import com.net2plan.libraries.WirelessUtils;
032import com.net2plan.utils.GradientProjectionUtils;
033import com.net2plan.utils.InputParameter;
034import com.net2plan.utils.Pair;
035import com.net2plan.utils.TimeTrace;
036import com.net2plan.utils.Triple;
037
038/** 
039 * This module implements a distributed primal-gradient based algorithm for adjusting the transmission power of the links in a wireless network subject to interferences, to maximize the network utility enforcing a fair allocation of the resources.
040 *
041 * Ths event processor is adapted to permit observing the algorithm performances under user-defined conditions, 
042 * including asynchronous distributed executions, where signaling can be affected by losses and/or delays, and/or measurement errors. 
043 * The time evolution of different metrics can be stored in output files, for later processing. 
044 * As an example, see the <a href="../../../../../../graphGeneratorFiles/fig_sec9_6_powerAssignmentPrimal.m">{@code fig_sec9_6_powerAssignmentPrimal.m}</a> MATLAB file used for generating the graph/s of the case study in the 
045 * <a href="http://eu.wiley.com/WileyCDA/WileyTitle/productCd-1119013356.html">book</a> using this algorithm.
046 * 
047 * To simulate a network with this module, use the {@code Online_evGen_doNothing} generator.
048 * 
049 * @net2plan.keywords Transmission power optimization, Wireless, Distributed algorithm, Primal gradient algorithm, Capacity assignment (CA)
050 * @net2plan.ocnbooksections Section 9.6
051 * @net2plan.inputParameters 
052 * @author Pablo Pavon-Marino
053 */
054public class Online_evProc_powerAssignmentPrimal extends IEventProcessor 
055{
056        private InputParameter signaling_isSynchronous = new InputParameter ("signaling_isSynchronous", false , "true if all the distributed agents involved wake up synchronously to send the signaling messages");
057        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);
058        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);
059        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);
060        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);
061        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);
062        private InputParameter update_isSynchronous = new InputParameter ("update_isSynchronous", false , "true if all the distributed agents involved wake up synchronousely to update its state");
063        private InputParameter update_averageInterUpdateTime = new InputParameter ("update_averageInterUpdateTime", 1.0 , "Average time between two updates of an agent" , 0 , false , Double.MAX_VALUE , true);
064        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);
065        
066        private InputParameter gradient_maxGradientAbsoluteNoise = new InputParameter ("gradient_maxGradientAbsoluteNoise", 0.01 , "Max value of the added noise to the gradient coordinate in absolute values" , 0 , true , Double.MAX_VALUE , true);
067        private InputParameter gradient_gammaStep = new InputParameter ("gradient_gammaStep", 10.0 , "Gamma step in the gradient algorithm" , 0 , false , Double.MAX_VALUE , true);
068        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);
069        private InputParameter gradient_maxGradientCoordinateChange = new InputParameter ("gradient_maxGradientCoordinateChange", 1000.0 , "Maximum change in an iteration of a gradient coordinate" , 0 , false , Double.MAX_VALUE , true);
070
071        private InputParameter simulation_randomSeed = new InputParameter ("simulation_randomSeed", (long) 1 , "Seed of the random number generator");
072        private InputParameter simulation_outFileNameRoot = new InputParameter ("simulation_outFileNameRoot", "powerTransmissionAsignmentPrimal" , "Root of the file name to be used in the output files. If blank, no output");
073        private InputParameter simulation_maxNumberOfUpdateIntervals = new InputParameter ("simulation_maxNumberOfUpdateIntervals", 100.0 , "Maximum number of update intervals in average per agent" , 0 , false , Double.MAX_VALUE , true);
074        
075        private InputParameter control_pathLossExponent = new InputParameter ("control_pathLossExponent", 3.0 , "Exponent in the model for propagation losses" , 0 , true , Double.MAX_VALUE , true);
076        private InputParameter control_worseRiseOverThermal_nu = new InputParameter ("control_worseRiseOverThermal_nu", 10.0 , "Worse case ratio between interference power and thermal power at any receiver. Used to set the common thermal power in the receivers" , 0 , true , Double.MAX_VALUE , true);
077        private InputParameter control_interferenceAttenuationFactor_nu = new InputParameter ("control_interferenceAttenuationFactor_nu", 1.0e6 , "The interference power received in natural units is divided by this to reduce its effect" , 1 , true , Double.MAX_VALUE , true);
078        private InputParameter control_maxTransmissionPower_logu = new InputParameter ("control_maxTransmissionPower_logu", 3.0 , "The maximum link transmission power in logarithmic units (e.g. dBm)");
079        private InputParameter control_minTransmissionPower_logu = new InputParameter ("control_minTransmissionPower_logu", 0.0 , "The minimum link transmission power in logarithmic units (e.g. dBm)");
080        private InputParameter control_fairnessFactor = new InputParameter ("control_fairnessFactor", 2.0 , "Fairness factor in utility function of link capacities" , 0 , true , Double.MAX_VALUE , true);
081
082        private static final int SIGNALING_WAKEUPTOSENDMESSAGE = 300;
083        private static final int SIGNALING_RECEIVEDMESSAGE = 301;
084        private static final int UPDATE_WAKEUPTOUPDATE = 302;
085
086        private Random rng;
087        private int E , N;
088        private DoubleMatrix2D mac_g_nu_ee;
089        private DoubleMatrix1D mac_transmissionPower_logu_e; 
090        private DoubleMatrix1D mac_previousTransmissionPower_logu_e; 
091        private double mac_receptionThermalNoise_nu;
092        
093        private DoubleMatrix2D control_mostUpdatedMe2ValueKnownByLink1_e1e2; 
094                        
095        private TimeTrace stat_traceOf_u_e;
096        private TimeTrace stat_traceOf_p_e;
097        private TimeTrace stat_traceOf_objFunction;
098
099        private NetPlan currentNetPlan , copyInitialNetPlan;
100        
101        @Override
102        public String getDescription()
103        {
104                return "This module implements a distributed primal-gradient based algorithm for adjusting the transmission power of the links in a wireless network subject to interferences, to maximize the network utility enforcing a fair allocation of the resources.";
105        }
106
107        @Override
108        public List<Triple<String, String, String>> getParameters()
109        {
110                /* Returns the parameter information for all the InputParameter objects defined in this object (uses Java reflection) */
111                return InputParameter.getInformationAllInputParameterFieldsOfObject(this);
112        }
113
114        @Override
115        public void initialize(NetPlan currentNp, Map<String, String> algorithmParameters, Map<String, String> simulationParameters, Map<String, String> net2planParameters)
116        {
117                /* Initialize all InputParameter objects defined in this object (this uses Java reflection) */
118                InputParameter.initializeAllInputParameterFieldsOfObject(this, algorithmParameters);
119
120                this.currentNetPlan = currentNp;
121                this.copyInitialNetPlan = currentNp.copy ();
122                this.E = currentNp.getNumberOfLinks ();
123                this.N = currentNp.getNumberOfNodes ();
124                if (E == 0) throw new Net2PlanException ("The input design should have links");
125                
126                if (currentNp.getNumberOfLayers() != 1) throw new Net2PlanException ("This algorithm works in single layer networks");
127
128                this.rng = new Random (simulation_randomSeed.getLong ());
129
130                /* Initialize the gains between links, normalizing them so that the maximum gain is one */
131                this.mac_g_nu_ee = WirelessUtils.computeInterferenceMatrixNaturalUnits (currentNetPlan.getLinks () , control_interferenceAttenuationFactor_nu.getDouble() , control_pathLossExponent.getDouble());
132                //System.out.println("NOT normalized Gnu_ee: " + Gnu_ee);
133                final double maxGainValue = mac_g_nu_ee.getMaxLocation() [0];
134                mac_g_nu_ee.assign (DoubleFunctions.div(maxGainValue));
135                //System.out.println("normalized mac_g_nu_ee: " + mac_g_nu_ee);
136                
137                /* Initialize the thermal noise at the receivers, to have a worse case ROT (rise over thermal) */
138                double worseInterferenceReceivedAtMaxPower_nu = WirelessUtils.computeWorseReceiverInterferencePower_nu (control_maxTransmissionPower_logu.getDouble()  , mac_g_nu_ee);
139
140                /* Adjust the thermal noise in the receivers so that we have a given ROT */
141                this.mac_receptionThermalNoise_nu = worseInterferenceReceivedAtMaxPower_nu / control_worseRiseOverThermal_nu.getDouble();
142
143                /* Initialize the transmission power variables */
144                this.mac_transmissionPower_logu_e = DoubleFactory1D.dense.make (E , control_minTransmissionPower_logu.getDouble());
145                this.mac_previousTransmissionPower_logu_e = DoubleFactory1D.dense.make (E , control_minTransmissionPower_logu.getDouble());
146                        
147                /* Update the netplan object with the resulting capacities */
148                for (Link e : currentNp.getLinks())
149                        e.setCapacity(Math.log(computeSINR_e (e)));
150
151                this.control_mostUpdatedMe2ValueKnownByLink1_e1e2 = DoubleFactory2D.dense.make (E,E);
152                for (Link e1 : currentNp.getLinks())
153                {
154                        for (Link otherLink : currentNp.getLinks())
155                                if (e1 != otherLink)
156                                        control_mostUpdatedMe2ValueKnownByLink1_e1e2.set(e1.getIndex () , otherLink.getIndex(), computeMeFactor_e (otherLink));
157                }
158                
159                /* Initially all links receive a "wake up to transmit" event, aligned at time zero or y asynchr => randomly chosen */
160                for (Link e : currentNp.getLinks())
161                {
162                        final double signalingTime = (signaling_isSynchronous.getBoolean())? signaling_averageInterMessageTime.getDouble() : Math.max(0 , signaling_averageInterMessageTime.getDouble() + signaling_maxFluctuationInterMessageTime.getDouble() * (rng.nextDouble() - 0.5));
163                        this.scheduleEvent(new SimEvent (signalingTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_WAKEUPTOSENDMESSAGE , e));
164                        final double updateTime = (update_isSynchronous.getBoolean())? update_averageInterUpdateTime.getDouble() : Math.max(0 , update_averageInterUpdateTime.getDouble() + update_maxFluctuationInterUpdateTime.getDouble() * (rng.nextDouble() - 0.5));
165                        this.scheduleEvent(new SimEvent (updateTime , SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_WAKEUPTOUPDATE , e));
166                }
167
168                /* Intialize the traces */
169                this.stat_traceOf_u_e = new TimeTrace ();
170                this.stat_traceOf_p_e = new TimeTrace ();
171                this.stat_traceOf_objFunction = new TimeTrace ();
172                this.stat_traceOf_u_e.add(0.0, this.currentNetPlan.getVectorLinkCapacity());
173                this.stat_traceOf_p_e.add(0.0, this.mac_transmissionPower_logu_e.copy ());
174                this.stat_traceOf_objFunction.add(0.0 , NetworkPerformanceMetrics.alphaUtility(currentNetPlan.getVectorLinkCapacity() , control_fairnessFactor.getDouble()));
175
176                /* */
177//              System.out.println("Initialize control_mostUpdatedMe2ValueKnownByLink1_e1e2: " + this.control_mostUpdatedMe2ValueKnownByLink1_e1e2);
178//              System.out.println("mac_g_nu_ee: " + this.mac_g_nu_ee);
179        }
180
181        @Override
182        public void processEvent(NetPlan currentNetPlan, SimEvent event)
183        {
184                final double t = event.getEventTime();
185                switch (event.getEventType())
186                {
187                case SIGNALING_RECEIVEDMESSAGE: // a link receives the signaling informatio
188                {
189                        final Pair<Link,Pair<Link,Double>> signalInfo = (Pair<Link,Pair<Link,Double>>) event.getEventObject();
190                        final Link eMe = signalInfo.getFirst();
191                        final Link otherLink = signalInfo.getSecond().getFirst();
192                        final double otherLink_me = signalInfo.getSecond().getSecond();
193                        control_mostUpdatedMe2ValueKnownByLink1_e1e2.set (eMe.getIndex () , otherLink.getIndex () , otherLink_me);
194                        break;
195                }
196                
197                case SIGNALING_WAKEUPTOSENDMESSAGE: // A link sends signaling information
198                {
199                        final Link eMe = (Link) event.getEventObject();
200
201                        Pair<Link,Double> infoToSignal = Pair.of(eMe , this.computeMeFactor_e(eMe));
202                        if (rng.nextDouble() >= this.signaling_signalingLossProbability.getDouble()) // the signaling may be lost => lost to all nodes
203                                for (Link otherLink : currentNetPlan.getLinks ())
204                                        if (otherLink != eMe)
205                                        {
206                                                final double signalingReceptionTime = t + Math.max(0 , signaling_averageDelay.getDouble() + signaling_maxFluctuationInDelay.getDouble() * (rng.nextDouble() - 0.5));
207                                                this.scheduleEvent(new SimEvent (signalingReceptionTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_RECEIVEDMESSAGE , Pair.of(otherLink , infoToSignal)));
208                                        }
209
210                        /* Re-schedule when to wake up again */
211                        final double signalingTime = signaling_isSynchronous.getBoolean()? t + signaling_averageInterMessageTime.getDouble() : Math.max(t , t + signaling_averageInterMessageTime.getDouble() + signaling_maxFluctuationInterMessageTime.getDouble() * (rng.nextDouble() - 0.5));
212                        this.scheduleEvent(new SimEvent (signalingTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_WAKEUPTOSENDMESSAGE , eMe));
213                        break;
214                }
215
216                case UPDATE_WAKEUPTOUPDATE: // a link updates its power 
217                {
218                        final Link eMe = (Link) event.getEventObject(); 
219                        
220                        final double currentTransmissionPower_logu = this.mac_transmissionPower_logu_e.get(eMe.getIndex ());
221                        double gradientThisLink = computeGradient (eMe) + 2*gradient_maxGradientAbsoluteNoise.getDouble()*(rng.nextDouble()-0.5);
222                        double nextTransmissionPower_logu = currentTransmissionPower_logu + this.gradient_gammaStep.getDouble() * gradientThisLink;
223                        
224                        /* heavy ball */
225                        nextTransmissionPower_logu += this.gradient_heavyBallBetaParameter.getDouble() * (currentTransmissionPower_logu - mac_previousTransmissionPower_logu_e.get(eMe.getIndex ()));
226                        /* projection */
227                        nextTransmissionPower_logu = GradientProjectionUtils.euclideanProjection_boxLike(nextTransmissionPower_logu, this.control_minTransmissionPower_logu.getDouble(), this.control_maxTransmissionPower_logu.getDouble());
228
229                        if (gradient_maxGradientCoordinateChange.getDouble() > 0)
230                                nextTransmissionPower_logu = GradientProjectionUtils.scaleDown_maxAbsoluteCoordinateChange (currentTransmissionPower_logu , nextTransmissionPower_logu ,  gradient_maxGradientCoordinateChange.getDouble());
231                        
232                        this.mac_previousTransmissionPower_logu_e.set(eMe.getIndex (), this.mac_transmissionPower_logu_e.get(eMe.getIndex ()));
233                        this.mac_transmissionPower_logu_e.set(eMe.getIndex (), nextTransmissionPower_logu);
234
235                        //System.out.println("Link " + eIdMe + ": mac_transmissionPower_logu_e values Before: " + currentTransmissionPower_logu + ", after: " + nextTransmissionPower_logu);
236
237                        /* Send event next recomputing time */
238                        final double updateTime = update_isSynchronous.getBoolean()? t + update_averageInterUpdateTime.getDouble() : Math.max(t , t + update_averageInterUpdateTime.getDouble() + update_maxFluctuationInterUpdateTime.getDouble() * (rng.nextDouble() - 0.5));
239                        this.scheduleEvent(new SimEvent (updateTime , SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_WAKEUPTOUPDATE,  eMe));
240
241                        /* Update in currentNetPlan the capacity values */
242                        for (Link e : this.currentNetPlan.getLinks())
243                                e.setCapacity(Math.log(computeSINR_e (e)));
244
245                        this.stat_traceOf_u_e.add(t, this.currentNetPlan.getVectorLinkCapacity());
246                        this.stat_traceOf_p_e.add(t, this.mac_transmissionPower_logu_e.copy ());
247                        this.stat_traceOf_objFunction.add(t , NetworkPerformanceMetrics.alphaUtility(currentNetPlan.getVectorLinkCapacity() , control_fairnessFactor.getDouble()));
248
249                        if (t > this.simulation_maxNumberOfUpdateIntervals.getDouble() * this.update_averageInterUpdateTime.getDouble()) { this.endSimulation (); }
250                        
251                        break;
252                }
253
254                default: throw new RuntimeException ("Unexpected received event");
255                }
256                
257                
258        }
259
260        public String finish (StringBuilder st , double simTime)
261        {
262                if (simulation_outFileNameRoot.getString().equals("")) return null;
263                stat_traceOf_u_e.printToFile(new File (simulation_outFileNameRoot.getString() + "_ue.txt"));
264                stat_traceOf_objFunction.printToFile(new File (simulation_outFileNameRoot.getString() + "_objFunc.txt"));
265                stat_traceOf_p_e.printToFile(new File (simulation_outFileNameRoot.getString() + "_pe.txt"));
266                /* compute optimum solution */
267                Map<String,String> param = new HashMap<String,String> ();
268                param.put("solverName", "ipopt");
269                param.put("solverLibraryName", "");
270                param.put("maxSolverTimeInSeconds", "-1");
271                param.put("alphaFairnessFactor", "" + this.control_fairnessFactor.getDouble());
272                param.put("pathLossExponent", "" + this.control_pathLossExponent.getDouble());
273                param.put("worseRiseOverThermal_nu", "" + this.control_worseRiseOverThermal_nu.getDouble());
274                param.put("interferenceAttenuationFactor_nu", "" + this.control_interferenceAttenuationFactor_nu.getDouble());
275                param.put("maxTransmissionPower_logu", "" + this.control_maxTransmissionPower_logu.getDouble());
276                param.put("minTransmissionPower_logu", "" + this.control_minTransmissionPower_logu.getDouble());
277                new Offline_ca_wirelessTransmissionPower ().executeAlgorithm(copyInitialNetPlan , param , null);
278                final double optimumNetUtilityJOM = NetworkPerformanceMetrics.alphaUtility(copyInitialNetPlan.getVectorLinkCapacity() , control_fairnessFactor.getDouble());
279                DoubleMatrix1D optimum_ue = copyInitialNetPlan.getVectorLinkCapacity();
280                DoubleMatrix1D optimum_pe = DoubleFactory1D.dense.make (E);
281                for (Link e : copyInitialNetPlan.getLinks ()) optimum_pe.set (e.getIndex () , Double.parseDouble (e.getAttribute("p_e")));
282                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_ue.txt"), optimum_ue);
283                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_pe.txt"), optimum_pe);
284                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_objFunc.txt"), optimumNetUtilityJOM);
285                return null;
286        }
287        
288
289        private double computeGradient (Link thisLink)
290        {
291                final double u_e = thisLink.getCapacity();
292                final DoubleMatrix1D infoIKnow_me = control_mostUpdatedMe2ValueKnownByLink1_e1e2.viewRow(thisLink.getIndex ());
293                
294                double gradient = Math.pow(u_e, -this.control_fairnessFactor.getDouble());
295                double accumFactor = 0;
296                for (Link epp : this.currentNetPlan.getLinks())
297                        if (epp != thisLink)
298                                accumFactor += this.mac_g_nu_ee.get(thisLink.getIndex (),epp.getIndex ()) * infoIKnow_me.get(epp.getIndex ());
299                accumFactor *= Math.exp(this.mac_transmissionPower_logu_e.get(thisLink.getIndex ()));
300                //System.out.println("Gradient: positive factor: " + gradient + ", negative factor: " + accumFactor);
301                gradient -= accumFactor;
302                return gradient;
303        }
304
305        private double computeMeFactor_e (Link e)
306        {
307                final double u_e = e.getCapacity();
308                final double snr_e = Math.exp(u_e);
309                return  Math.pow(u_e,-this.control_fairnessFactor.getDouble()) * snr_e / (Math.exp(this.mac_transmissionPower_logu_e.get(e.getIndex ())) * this.mac_g_nu_ee.get(e.getIndex (),e.getIndex ()));
310        }
311        
312        private double computeSINR_e (Link e)
313        {
314                final double receivedPower_nu = Math.exp(this.mac_transmissionPower_logu_e.get(e.getIndex ())) * this.mac_g_nu_ee.get(e.getIndex (),e.getIndex ());
315                double interferencePower_nu = this.mac_receptionThermalNoise_nu;
316                for (Link eInt : this.currentNetPlan.getLinks ())
317                        if (eInt != e) interferencePower_nu += Math.exp(this.mac_transmissionPower_logu_e.get(eInt.getIndex ())) * this.mac_g_nu_ee.get(eInt.getIndex (),e.getIndex ());
318//              System.out.println ("SINR link " + e + ": " + receivedPower_nu / interferencePower_nu);
319//              System.out.println ("receiver power link " + e + ": " + receivedPower_nu + ", total interf power: " + interferencePower_nu + "thermal noise: " + mac_receptionThermalNoise_nu);
320                return receivedPower_nu / interferencePower_nu;
321        }
322
323}