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