View Javadoc
1   package org.opentrafficsim.road.gtu.generator.headway;
2   
3   import org.djunits.value.vdouble.scalar.Duration;
4   import org.djunits.value.vdouble.scalar.Time;
5   import org.opentrafficsim.base.parameters.ParameterException;
6   import org.opentrafficsim.core.distributions.Generator;
7   import org.opentrafficsim.core.distributions.ProbabilityException;
8   
9   import nl.tudelft.simulation.dsol.simulators.SimulatorInterface;
10  import nl.tudelft.simulation.jstats.distributions.DistNormal;
11  import nl.tudelft.simulation.jstats.streams.StreamInterface;
12  
13  /**
14   * Headway generation based on {@code Arrivals}.
15   * <p>
16   * Copyright (c) 2013-2018 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
17   * BSD-style license. See <a href="http://opentrafficsim.org/node/13">OpenTrafficSim License</a>.
18   * <p>
19   * @version $Revision$, $LastChangedDate$, by $Author$, initial version 13 dec. 2017 <br>
20   * @author <a href="http://www.tbm.tudelft.nl/averbraeck">Alexander Verbraeck</a>
21   * @author <a href="http://www.tudelft.nl/pknoppers">Peter Knoppers</a>
22   * @author <a href="http://www.transport.citg.tudelft.nl">Wouter Schakel</a>
23   */
24  public class ArrivalsHeadwayGenerator implements Generator<Duration>
25  {
26  
27      /** Arrivals. */
28      private final Arrivals arrivals;
29  
30      /** Simulator. */
31      private final SimulatorInterface.TimeDoubleUnit simulator;
32  
33      /** Random stream to draw headway. */
34      private final StreamInterface stream;
35  
36      /** Random headway generator. */
37      private final HeadwayDistribution distribution;
38  
39      /** First GTU. */
40      private boolean first = true;
41  
42      /**
43       * @param arrivals Arrivals; arrivals
44       * @param simulator SimulatorInterface.TimeDoubleUnit; simulator
45       * @param stream StreamInterface; random stream to draw headway
46       * @param distribution HeadwayDistribution; random headway distribution
47       */
48      public ArrivalsHeadwayGenerator(final Arrivals arrivals, final SimulatorInterface.TimeDoubleUnit simulator,
49              final StreamInterface stream, final HeadwayDistribution distribution)
50      {
51          this.arrivals = arrivals;
52          this.simulator = simulator;
53          this.stream = stream;
54          this.distribution = distribution;
55      }
56  
57      /**
58       * Returns a new headway {@code h} assuming that the previous vehicle arrived at the current time {@code t0}. The vehicle
59       * thus arrives at {@code t1 = t0 + h}. This method guarantees that no vehicle arrives during periods where demand is zero,
60       * while maintaining random headways based on average demand over a certain time period.<br>
61       * <br>
62       * The general method is to find {@code h} such that the integral of the demand pattern {@code D} from {@code t0} until
63       * {@code t1} equals {@code r}: &#931;{@code D(t0 > t1) = r}. One can think of {@code r} as being 1 and representing an
64       * additional vehicle to arrive. The headway {@code h} that results correlates directly to the mean demand between
65       * {@code t0} and {@code t1}.<br>
66       * <br>
67       * The value of {@code r} always has a mean of 1, but may vary between specific vehicle arrivals depending on the headway
68       * distribution. When assuming constant headways for any given demand level, {@code r} always equals 1. For exponentially
69       * distributed headways {@code r} may range anywhere between 0 and infinity.<br>
70       * <br>
71       * This usage of {@code r} guarantees that no vehicles arrive during periods with 0 demand. For example:
72       * <ul>
73       * <li>Suppose we have 0 demand between 300s and 400s.</li>
74       * <li>The previous vehicle was generated at 299s.</li>
75       * <li>The demand at 299s equals 1800veh/h (1 veh per 2s).</li>
76       * <li>For both constant and exponentially distributed headways, the expected next vehicle arrival based on this demand
77       * value alone would be 299 + 2 = 301s. This is within the 0-demand period and should not happen. It's also not
78       * theoretically sound, as the demand from 299s until 301s is not 1800veh/h on average.</li>
79       * <li>Using integration we find that the surface of demand from 299s until 300s equals 0.5 veh for stepwise demand, and
80       * 0.25 veh for linear demand. Consequently, the vehicle will not arrive until later slices integrate to an additional 0.5
81       * veh or 0.75 veh respectively. This additional surface under the demand curve is only found after 400s.</li>
82       * <li>In case the exponential headway distribution would have resulted in {@code r} &lt; 0.5 (stepwise demand) or 0.25
83       * (linear demand), a vehicle will simply arrive between 299s and 300s.</li>
84       * </ul>
85       * <br>
86       * @return Duration; new headway
87       * @throws ProbabilityException if the stored collection is empty
88       * @throws ParameterException in case of a parameter exception
89       */
90      @Override
91      public Duration draw() throws ProbabilityException, ParameterException
92      {
93          Time now = this.simulator.getSimulatorTime();
94          // initial slice times and frequencies
95          Time t1 = now;
96          double f1 = this.arrivals.getFrequency(t1, true).si;
97          Time t2 = this.arrivals.nextTimeSlice(t1);
98          if (t2 == null)
99          {
100             return null; // no new vehicle
101         }
102         double f2 = this.arrivals.getFrequency(t2, false).si;
103         // next vehicle's random factor
104         double rem = this.distribution.draw(this.stream);
105         if (this.first)
106         {
107             // first headway may be partially in the past, take a random factor
108             rem *= this.stream.nextDouble();
109             this.first = false;
110         }
111         // integrate until rem (by reducing it to 0.0, possibly in steps per slice)
112         while (rem > 0.0)
113         {
114             // extrapolate to find 'integration = rem' in this slice giving demand slope, this may beyond the slice length
115             double dt = t2.si - t1.si;
116             double t;
117             double slope = (f2 - f1) / dt;
118             if (Math.abs(slope) < 1e-12) // no slope
119             {
120                 if (f1 > 0.0)
121                 {
122                     t = rem / f1; // rem = t * f1, t = rem / f1
123                 }
124                 else
125                 {
126                     t = Double.POSITIVE_INFINITY; // no demand in this slice
127                 }
128             }
129             else
130             {
131                 // reverse of trapezoidal rule: rem = t * (f1 + (f1 + t * slope)) / 2
132                 double sqrt = 2 * slope * rem + f1 * f1;
133                 if (sqrt >= 0.0)
134                 {
135                     t = (-f1 + Math.sqrt(sqrt)) / slope;
136                 }
137                 else
138                 {
139                     t = Double.POSITIVE_INFINITY; // not sufficient demand in this slice, with negative slope
140                 }
141             }
142             if (t > dt)
143             {
144                 // next slice
145                 rem -= dt * (f1 + f2) / 2; // subtract integral of this slice using trapezoidal rule
146                 t1 = t2;
147                 t2 = this.arrivals.nextTimeSlice(t1);
148                 if (t2 == null)
149                 {
150                     return null; // no new vehicle
151                 }
152                 f1 = this.arrivals.getFrequency(t1, true).si; // we can't use f1 = f2 due to possible steps in demand
153                 f2 = this.arrivals.getFrequency(t2, false).si;
154             }
155             else
156             {
157                 // return resulting integration times
158                 return Duration.createSI(t1.si + t - now.si);
159             }
160         }
161         throw new RuntimeException("Exception while determining headway from Arrivals.");
162     }
163 
164     /**
165      * Headway distribution.
166      * <p>
167      * Copyright (c) 2013-2018 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved.
168      * <br>
169      * BSD-style license. See <a href="http://opentrafficsim.org/node/13">OpenTrafficSim License</a>.
170      * <p>
171      * @version $Revision$, $LastChangedDate$, by $Author$, initial version 5 dec. 2017 <br>
172      * @author <a href="http://www.tbm.tudelft.nl/averbraeck">Alexander Verbraeck</a>
173      * @author <a href="http://www.tudelft.nl/pknoppers">Peter Knoppers</a>
174      * @author <a href="http://www.transport.citg.tudelft.nl">Wouter Schakel</a>
175      */
176     public interface HeadwayDistribution
177     {
178 
179         /** Constant headway. */
180         HeadwayDistribution CONSTANT = new HeadwayDistribution()
181         {
182             @Override
183             public double draw(final StreamInterface randomStream)
184             {
185                 return 1.0;
186             }
187 
188             @Override
189             public String getName()
190             {
191                 return "CONSTANT";
192             }
193         };
194 
195         /** Exponential headway distribution. */
196         HeadwayDistribution EXPONENTIAL = new HeadwayDistribution()
197         {
198             @Override
199             public double draw(final StreamInterface randomStream)
200             {
201                 return -Math.log(randomStream.nextDouble());
202             }
203 
204             @Override
205             public String getName()
206             {
207                 return "EXPONENTIAL";
208             }
209         };
210 
211         /** Uniform headway distribution. */
212         HeadwayDistribution UNIFORM = new HeadwayDistribution()
213         {
214             @Override
215             public double draw(final StreamInterface randomStream)
216             {
217                 return 2.0 * randomStream.nextDouble();
218             }
219 
220             @Override
221             public String getName()
222             {
223                 return "UNIFORM";
224             }
225         };
226 
227         /** Triangular headway distribution. */
228         HeadwayDistribution TRIANGULAR = new HeadwayDistribution()
229         {
230             @Override
231             public double draw(final StreamInterface randomStream)
232             {
233                 double r = randomStream.nextDouble();
234                 if (r < .5)
235                 {
236                     return Math.sqrt(r * 2.0);
237                 }
238                 return 2.0 - Math.sqrt((1.0 - r) * 2.0);
239             }
240 
241             @Override
242             public String getName()
243             {
244                 return "TRIANGULAR";
245             }
246         };
247 
248         /** Triangular (left side, mean 2/3) and exponential (right side, mean 4/3) headway distribution. */
249         HeadwayDistribution TRI_EXP = new HeadwayDistribution()
250         {
251             @Override
252             public double draw(final StreamInterface randomStream)
253             {
254                 double r = randomStream.nextDouble();
255                 if (r < .5)
256                 {
257                     return Math.sqrt(r * 2.0); // left-hand side of triangular distribution, with mean 2/3
258                 }
259                 return 1.0 - Math.log((1.0 - r) * 2.0) / 3.0; // 1 + 1/3, where 1/3 is mean of exponential right-hand side
260                 // note: 50% with mean 2/3 and 50% with mean 1 + 1/3 gives a mean of 1
261             }
262 
263             @Override
264             public String getName()
265             {
266                 return "TRI_EXP";
267             }
268         };
269 
270         /** Log-normal headway distribution (variance = 1.0). */
271         HeadwayDistribution LOGNORMAL = new HeadwayDistribution()
272         {
273             /** Mu. */
274             private final double mu = Math.log(1.0 / Math.sqrt(2.0));
275 
276             /** Sigma. */
277             private final double sigma = Math.sqrt(Math.log(2.0));
278 
279             @Override
280             public double draw(final StreamInterface randomStream)
281             {
282                 return Math.exp(new DistNormal(randomStream, this.mu, this.sigma).draw());
283             }
284 
285             @Override
286             public String getName()
287             {
288                 return "LOGNORMAL";
289             }
290         };
291 
292         /**
293          * Draws a randomized headway factor. The average value returned is always 1.0. The returned value is applied on the
294          * demand pattern by (reversed) integration to derive actual headways.
295          * @param randomStream StreamInterface; random number stream
296          * @return randomized headway factor
297          */
298         double draw(StreamInterface randomStream);
299 
300         /**
301          * Returns the distribution name.
302          * @return distribution name
303          */
304         String getName();
305 
306     }
307 
308 }