Mathematical Foundations

Only methods that appear in the user stories are included.

Contents

• Notation and terms

1. Break-even calculation
2. Pulsed arrival model
3. Probability distributions and parameter fitting
4. Dwell time and boarding capacity
5. Queueing at the transfer corridor
6. Delay propagation and recovery
7. Prediction error model
8. Hold duration search and optimization
9. Common random numbers
10. Confidence intervals and replication count
11. Sensitivity Analysis
12. Loop gain

---

Notation and terms

Notation

Symbol	Read as	Meaning
$\Sigma$	"sum of"	add up a list of values
$\prod$	"product of"	multiply a list of values together
$E[X]$	"expected value of X"	the average of $X$ over many repeats of the same situation
$\int$	"integral"	add up a changing quantity continuously over an interval, which is the area under its curve
$X \sim \text{Name}$	"X follows Name"	$X$ is a random value drawn from that distribution
$\hat{\theta}$	"theta hat"	an estimate calculated from data, as opposed to the true value
$\max(a, b)$, $\min(a, b)$	"max", "min"	the larger or the smaller of the two values
$\arg\max$	"argument of the maximum"	the input that gives the largest output, as opposed to the output itself
$f^{-1}$	"inverse of f"	the function that reverses $f$. It is not $1/f$
$\bar{X}$	"X bar"	the sample average of $X$
$Var(X)$	"variance of X"	how much the values of $X$ vary around their average
$Cov(X,Y)$	"covariance of X and Y"	how much $X$ and $Y$ vary together
$s$	"sample standard deviation"	estimated spread of observed results
$SE$	"standard error"	estimated uncertainty in a sample mean
$R$	"number of replications"	number of complete simulation runs
$\alpha$	"alpha"	significance level used when calculating a confidence interval
$\beta$	"beta"	probability of failing to detect a specified difference
$\delta$	"delta"	minimum meaningful difference the simulation should be able to detect
$\sigma_D$	"sigma D"	standard deviation of the paired differences
$G$	"loop gain"	The total effect of a disturbance after it travels around a feedback loop
$a_i$	"a sub i"	Attenuation factor for one step of disturbance propagation
$\rho(A)$	"spectral radius of A"	The largest absolute value of the eigenvalues of matrix $A$; used to assess whether repeated disturbances decay or grow
$\lambda_i$	"lambda sub i"	The $i$th eigenvalue of matrix $A$

Transit terms

Term	Meaning
Interchange	A station where passengers change between the REM and the Metro
Headway	The time between one train and the next on the same line
Dwell	The time a train spends stopped at a platform with its doors open. It is the normal stop, not the same as a hold
Hold	Extra time a train waits beyond its normal dwell so that more passengers can board
Onboard	Passengers who are already on the train
Alighting	Getting off a train
Terminus	The last station on a line, where the train turns around
Layover	Scheduled time at the terminus before the train starts its return trip
Slack	Spare time in the schedule that absorbs lateness. At the terminus this is the layover
Connection window	The time between the REM train arriving and the Metro's scheduled departure. This is the time passengers have to make the connection with no hold
Knock-on delay	Delay passed from one train to the trains behind it

Math terms

Term	Meaning
Distribution	How a quantity is spread across its possible values
Cumulative	A running total that only increases
Parameter	An adjustable number that controls the shape of a model
Probability density	The height of a distribution's curve. A higher curve at a point means values near that point are more likely
Likelihood	How probable the observed data would be if a given set of parameters were true
Quantile	The value below which a given fraction of the data falls. The median is the 50% quantile
Marginal	The extra amount from one more unit, for example the extra benefit from one more second of holding
Swept parameter	A parameter with no known value, so the simulation is run across a range of values
Optimization	Finding the best value or decision according to a defined objective
Discrete optimization	Optimization where only specific choices are allowed, such as 0, 10, 20, 30, 40, and 50 seconds
Objective	The quantity that the optimization tries to improve
Candidate	One possible value or strategy being evaluated
Rolling horizon	Repeatedly optimizing over a forward window and re-solving as the simulation moves forward
Variance	A measure of how much simulation results vary around their average
Covariance	A measure of how two quantities vary together
Standard deviation	A measure of the typical spread of values around their average
Standard error	An estimate of how much the sample mean would vary between repeated samples
Confidence interval	A range used to describe the uncertainty around an estimated average
Statistical power	The probability of detecting a specified difference when that difference is actually present
Replication	One complete run of the simulation
Paired difference	The difference between two strategies evaluated under the same random scenario
Random substream	A separate sequence of random values used for one source of randomness
Feedback loop	A chain of effects that eventually influences the part of the system where the chain started
Attenuation	Reduction in the size or effect of a disturbance as it propagates through the system
Loop gain	The total multiplication of a disturbance after it travels around a feedback loop
Spectral radius	The largest absolute value of the eigenvalues of a matrix; used to assess whether repeated effects decay or grow
Eigenvalue	A value describing how a particular direction or pattern is scaled by a matrix
Propagation matrix	A matrix describing how a disturbance at one part of the system affects another part
Stability	The tendency of disturbances to decay rather than grow as they propagate through the system

---

1. Break-even calculation

Concept: Cost-benefit comparison

Every holding decision is a trade. Waiting benefits the passengers who then make their connection, and costs everyone already on board. The decision is worthwhile only when the first outweighs the second.

Both sides are expressed in passenger-seconds, which makes them directly comparable: one person delayed sixty seconds is the same quantity as sixty people delayed one second.

An example that shows how the units work: If holding the train for 30 seconds lets 20 more people connect, and each of them would otherwise wait 5 minutes for the next train, the saving is 20 × 300 = 6,000 passenger-seconds. The 180 passengers already on board each lose 30 seconds, a cost of 180 × 30 = 5,400 passenger-seconds. The hold is worthwhile, but only by a small margin.

General form

A standard cost-benefit decision rule: act when the expected gain exceeds the expected loss.

$$\sum \text{benefit} > \sum \text{cost}$$

Our form

$$\text{Benefit}(h) = N(h) \cdot H \cdot w_{\text{wait}}$$

$$\text{Cost}(h) = L \cdot h \cdot w_{\text{onboard}} + \text{downstream cost}$$

$$\text{Hold if } \text{Benefit}(h) > \text{Cost}(h)$$

Symbol	Meaning
$h$	hold duration in seconds
$N(h)$	passengers who board because of the hold
$H$	headway they would otherwise wait
$L$	passengers already on board
$w$	weighting of waiting time against onboard time

The weighting $w$ is needed because waiting on a platform and sitting on a train are not experienced equally, even when they last the same amount of time. Waiting is generally worse, so the two are weighted differently. The ratio comes from the rider survey.

Why we use it

Holding a train helps some passengers and delays others. The break-even calculation puts both sides into one number, so TransitSync can decide whether a hold is worth it.

It uses passenger-seconds instead of headcount. For example, 30 passengers who each save 20 minutes gain 600 passenger-minutes, while 800 passengers who each lose 40 seconds lose about 533, so the small group wins.

Because the decision is one clear inequality, it can be:

• checked by hand
• reported with both the benefit and the cost shown
• negative, which is a valid result when the cost is larger

Where we use it

This is the core of the decision engine. It is implemented as a pure function with no dependency on simulator state, so it can be unit tested against values calculated by hand and explained independently in the final report.

---

2. Pulsed arrival model

Concept: Non-homogeneous point process

The observation

Transferring passengers do not arrive at a steady rate. A REM train arrives, and over the next few minutes a crowd emerges from the corridor: fast walkers first, then the bulk of the group, then a thinning tail.

Arrivals occur in waves tied to train events rather than as a continuous trickle.

General form

The counting process $N(t)$, is a running tally of how many events have occurred, and only ever increases, records how many events have occurred by time $t$. The rate at which that tally is increasing at a given moment is called the intensity, written $\lambda$.

Homogeneous — the rate stays constant, so the expected count grows in a straight line:

$$E[N(t)] = \lambda t$$

Non-homogeneous — the rate changes over time, so the expected count is the rate accumulated over the interval:

$$E[N(t)] = \int_0^t \lambda(u), du$$

$E[N(t)]$ is the expected count, meaning the average over many repeats. The integral adds up the intensity over every instant from 0 to $t$, and reduces to $\lambda t$ when the intensity is constant.

Transfer arrivals are the second kind. The rate jumps after each train arrival and falls to near zero in between.

Why the homogeneous case fails

Under a constant rate $\lambda$, the number of passengers caught by a hold of $h$ seconds is rate times duration:

$$N(h) = \lambda h$$

Substituting into the break-even condition from section 1:

$$\lambda h \cdot H \cdot w_{\text{wait}} > L \cdot h \cdot w_{\text{onboard}}$$

$$\lambda H \cdot w_{\text{wait}} > L \cdot w_{\text{onboard}}$$

The hold duration cancels out. The condition no longer depends on $h$, so there is no best duration to find: the model allows only an indefinite hold or none at all. Under a constant rate, each extra second of holding catches the same number of extra passengers, so the case for holding one more second never changes. There is no point at which holding stops being worthwhile and none at which it starts, so the model gives no way to choose a duration.

This is the key step. Under a constant arrival rate, each additional second of holding catches the same number of additional passengers, so the case for holding one more second never changes. There is no point at which holding stops being worthwhile, and no point at which it starts. The model gives no basis for choosing a duration, which is exactly what we are trying to determine.

Since $\lambda$ is small for a dispersed stream of arrivals while $L$ is comparatively large, the inequality is rarely satisfied. A constant-rate model therefore returns a negative result at every station, every time of day, and under every parameter setting.

Our form

Arrivals are driven by train events rather than by the clock. Each train releases a group of $M$ passengers whose walking times follow a distribution $F$.

$F(x)$ is a cumulative distribution: the fraction of the group that has finished walking by time $x$. If $F(100) = 0.4$, then 40% have reached the platform 100 seconds after the train arrived.

The passengers a hold rescues arrive after the connecting train was due to leave but before the hold ends:

$$N(h) = M \cdot \left[F(h + s) - F(s)\right]$$

Symbol	Meaning
$M$	size of the group released by one train
$F$	cumulative distribution of walking times, estimated from field measurement
$h$	hold duration under consideration
$s$	connection window: time between the arriving train and the connecting train's scheduled departure, in either direction (REM to Metro, or Metro to REM)

Example: $M = 60$, $s = 90$ s, $F(90) = 0.45$ and $F(120) = 0.70$. Without a hold, $60 \times 0.45 = 27$ passengers make it. Holding $h = 30$ s adds $60 \times (0.70 - 0.45) = 15$.

Where the values come from

$F$ is measured from the field sessions. $s$ is computed from the simulated schedule. $M$ is swept: transfer volumes are not published, so rather than guessing a value we run the simulation across the whole plausible range, at 40 transferring passengers, then 60, then 80, and report how the answer changes.

The field sessions exist to estimate $F$. At each interchange, arrivals onto the connecting platform are counted in 15-second blocks for six minutes after each train arrival. Repeated across roughly ten trains per session and across time bands, this produces the empirical shape of the walk, from which a continuous distribution is fitted by maximum likelihood (section 3).

Until those measurements are complete, the simulation runs on a placeholder distribution, labelled as assumed in the parameter provenance system. Replacing it is a configuration change rather than a code change.

Why this formulation admits an optimum

$F$ is a cumulative distribution, so it increases as the hold becomes longer and eventually approaches 1 as most of the transferring group has arrived.

Therefore, $N(h)$ increases with the hold duration, but the number of additional passengers rescued by each extra second of holding eventually decreases. Once most of the group has already reached the platform, additional holding provides little additional benefit.

At the same time, the cost of holding continues to increase because every passenger already on board experiences the additional delay.

This creates a trade-off: initially, additional holding can rescue many passengers, but eventually the marginal benefit becomes smaller than the marginal cost. The optimal hold duration is therefore the point where extending the hold further is no longer worthwhile.

The exact location of this point depends on the shape of the walking-time distribution $F$, which is why field measurement is important. A narrow distribution means passengers arrive within a relatively short window, so most potential connections are recovered quickly. A wider distribution with a longer tail means additional passengers continue arriving later, potentially making longer holds worthwhile.

The mean walking time alone is therefore not enough to characterize the decision. The distribution and especially its tail determine how many passengers can potentially be recovered by extending the hold.

Why we use it

How many extra passengers a hold catches depends on how arrivals are spread out after a train comes in. A constant rate cannot show this, so TransitSync models arrivals as a wave tied to each train.

This also connects the field work directly to the result, since the shape of $F$ decides where the best hold falls.

The model responds to conditions such as:

• the size of the transferring group
• the walking-time distribution at each station and direction
• the connection window

Where we use it

Passenger generation in the simulation engine. Transferring passengers are created as a batch on each train arrival event and spread over time according to their sampled walking times.

A regression test asserts that arrivals remain non-uniform. A later simplification to a constant rate would run without error and return plausible values while removing the effect under study.

---

3. Probability distributions and parameter fitting

Concept: Log-normal and gamma distributions; maximum likelihood estimation

The problem

Section 2 established that the model needs $F$, the distribution of walking times. The field sessions give us counts in 15-second blocks, which is essentially a bar chart of how many people arrived in each window.

A bar chart is not enough. The simulation needs to pick a walking time for each individual passenger, including values that fall between the blocks we measured, and values in the tail where we only saw one or two people. So we fit a smooth curve to the bars and draw from that instead.

The same applies to run times between stations. In every case the spread matters as much as the average, since the spread is what decides how many extra passengers a longer hold catches.

Choosing a distribution

Walking times have two properties that rule out most options.

There is a minimum. Nobody walks the corridor faster than some floor value. Times below that are impossible, not just unlikely.

The slow side has a longer tail than the fast side. Most people cluster near a typical time, but some take much longer: luggage, checking signage, waiting behind people on an escalator. There is no equivalent on the fast side, because walking speed has a limit.

A normal distribution has neither property. It is symmetric and has no floor, so it would allow negative walking times and would underestimate how many slow walkers there are. Since slow walkers are exactly the people a hold rescues, that error would make holding look less useful than it is.

The log-normal and gamma families both start at zero and have a longer tail on the right, which makes them the sensible candidates. We fit both and compare.

A log-normal curve is one where the logarithm of the values forms a bell curve: it starts at zero, peaks, and has a long right tail. A gamma curve is similar and differs mainly in how fast the tail dies away. Both are fitted and compared.

The two families are similar and differ mainly in the behaviour of the tail. It is not necessary to choose between them in advance, since fitting both and comparing the results requires little additional work.

What our field data actually gives us

Our observations are not exact walking times. They are counts within intervals.

If 12 passengers are counted during the block from 60 to 75 seconds, what we know is that those 12 arrived somewhere in that window. We do not know whether a given one took 61 seconds or 74. This is called interval data, and it is a consequence of how the measurement is done: counting people into blocks is practical on a platform, timing individuals is not.

For the initial fit, we represent each block by its midpoint. The 60 to 75 second block becomes 12 observations at 67.5 seconds.

This is an approximation. It lets us use standard fitting methods on data that was not collected as individual timings, at the cost of some precision within each block. The fitted distribution should therefore be read as an estimate of the underlying walking-time distribution, not as an exact record of what each passenger did.

The approximation is acceptable here because the blocks are narrow relative to the spread we are trying to capture. A 15-second block against a distribution spanning several minutes loses little. It would be a poor approximation if the blocks were wide or the distribution narrow.

General form

Once we have picked a family, we need its parameters. A log-normal distribution is described by two numbers, $\mu$ and $\sigma$. Different values give different curve shapes, and we want the pair that best matches what we observed.

Maximum likelihood estimation finds them by asking a simple question: for which parameter values would our observed data have been most likely to happen?

In practice this means we try different values for the distribution parameters and keep the ones that make our observed walking-time data most likely. It is a search over parameter values, scored by how well each candidate explains the data.

Two pieces of notation appear below:

• $\prod_i$ means multiply together, over every observation $i$. It is the multiplication equivalent of $\sum_i$.
• $\arg\max$ means the input that produces the maximum, as opposed to $\max$, which is the maximum value itself. If $g(3) = 10$ is the largest value $g$ takes, then $\max g = 10$ but $\arg\max g = 3$. Here we want the parameters, not the score, so we use $\arg\max$.

For a candidate set of parameters $\theta$, we multiply together how likely each individual observation was. That product is called the likelihood:

$$L(\theta) = \prod_i f(x_i ; \theta)$$

$$\hat{\theta} = \arg\max_\theta L(\theta)$$

The hat on $\hat{\theta}$ marks it as an estimate derived from data, rather than a true value we know.

In practice we take the logarithm first, which turns the product into a sum:

$$\hat{\theta} = \arg\max_\theta \sum_i \log f(x_i ; \theta)$$

This gives the same answer, because the logarithm is increasing: whatever maximises the product also maximises its logarithm. It avoids a practical problem, which is that multiplying hundreds of small probabilities together drives the result to zero in floating point, while adding their logarithms does not.

Our form

$$T_{\text{walk}} \sim \text{Lognormal}(\mu, \sigma^2)$$

The symbol $\sim$ means "follows": the walking time of a transferring passenger is a random value drawn from a log-normal curve.

$$\hat{\mu}, \hat{\sigma} = \arg\max \sum_i \log f(t_i ; \mu, \sigma)$$

Symbol	Meaning
$T_{\text{walk}}$	walking time for one transferring passenger
$t_i$	approximate walking times, taken as block midpoints
$\mu, \sigma$	parameters of the fitted curve

For a log-normal curve the median is $e^{\mu}$. As an example, with a median of 110 seconds, $\mu = \ln(110) \approx 4.7$. The real values are acquired from fitting the field data. A larger $\sigma$ gives a longer slow tail.

In practice a library runs the search, for example scipy.stats.lognorm.fit. SciPy names the parameters differently: its shape is $\sigma$ and its scale is $e^{\mu}$.

Symbol	Meaning
$T_{\text{walk}}$	walking time for one transferring passenger
$t_i$	approximate walking times, taken as block midpoints
$\mu, \sigma$	parameters of the fitted distribution

Checking the fit

Fitting always returns parameters, even when the family we chose was the wrong one. So the fit has to be checked rather than trusted.

Q-Q plot. Q-Q means quantile-quantile. A quantile is the value below which a given fraction of the data falls, so the median is the 50% quantile. Suppose 10% of passengers in the data had arrived by 72 s, while the fitted curve says its 10% quantile is 75 s. That is one point, (75, 72), with the curve's value on the horizontal axis and the data's on the vertical. Repeating at 20%, 30% and so on gives the full plot.

A good fit puts the points on the diagonal where the two values are equal. Points above it mean real walking times are longer than the curve predicts, and points below it mean shorter. Points that follow the diagonal through the middle but curve upward at the far right mean the slow tail is heavier than the curve allows. That last case matters most here, since slow walkers are the ones a hold rescues. Because the data are in 15-second blocks, the points form steps, which is expected.

Kolmogorov-Smirnov test. The cumulative curve of the data and the cumulative curve of the fit are drawn on the same axes, and the statistic is the largest vertical gap between them. If 70% of passengers had arrived by 120 s in the data and the fit says 64%, the gap there is 0.06. A small statistic means a close match. The test also gives a p-value, but with binned data and parameters estimated from the same data it is only approximate, so the size of the gap is the main result.

We fit both log-normal and gamma, compare these two checks, and record which one we chose and why.

Why we use it

Walking times differ from person to person, and that variation decides how many extra passengers a longer hold catches. A single average is not enough to compare a 30 second hold with a 60 second one.

Fitting a curve turns a few counts in 15-second blocks into something the simulation can draw from, including values between the blocks and in the slow tail. The same method is used for run times between stations.

Checking the fit matters because a curve that misses the slow tail would make holding look less useful than it is.

Where we use it

Fitted once from the field measurements and archived schedule data, then written into the configuration file. The simulation reads those parameters and samples from them. No fitting happens while the simulation runs.

That separation matters for reproducibility: a run from November and a run from March give the same answer because the parameters are a committed file rather than something recalculated from data that may have changed.

Implemented in: US-1.3 — Passenger generation and transfer arrivals, tasks Fit walking time distribution from field data and Implement walking time sampling. Also US-1.2, task Implement run time sampling.

---

4. Dwell time and boarding capacity

Concept: Throughput constraint with a fixed service window; linear regression with demand

The observation

The two networks handle doors differently, so they need different models.

The REM is automated. A train stops, the doors open for a set period, and they close on schedule regardless of how many people are still boarding. Dwell is effectively constant.

The Metro is not. The front-cab operator controls the doors manually, watching the platform and closing when boarding looks finished. Dwell therefore varies with demand, not through a mechanism but through a human decision.

On the REM, a hold costs exactly the hold duration and nothing more, but anyone who cannot get through the doors in the fixed window is left behind. On the Metro, a hold that brings a large group to the platform also lengthens the subsequent dwell, so the cost exceeds the hold itself.

It is also worth noting that the Metro operator is already making an informal version of the decision this project formalises: someone in the cab sees people approaching and waits a moment. Our system does not introduce holding to a network without it. It asks whether that judgment could be made better with information the operator does not have.

Dwell is the time a train spends at the platform with its doors open. It is distinct from a hold: dwell is the normal stop that occurs regardless, while a hold is additional time added on top of it.

General form

Where service is provided at a fixed rate over a fixed window, the quantity served is capped by their product:

$$\text{served} = \min(\text{demand},\ \text{rate} \times \text{window})$$

Where service time grows with demand, multiple linear regression applies:

$$y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \varepsilon$$

Our form

REM — fixed dwell, capacity-constrained boarding

$$d = d_{\text{scheduled}}$$

$$B_{\text{actual}} = \min\left(B_{\text{waiting}},\ C_{\text{door}} \cdot d_{\text{scheduled}} \cdot n_{\text{doors}}\right)$$

Metro — demand-dependent dwell

$$d = d_{\min} + \beta_b B + \beta_a A$$

Symbol	Meaning
$d_{\text{scheduled}}$	fixed dwell on the REM, an operator setting
$d_{\min}$	minimum Metro dwell with no passenger activity
$B, A$	passengers boarding and alighting
$B_{\text{waiting}}$	passengers attempting to board
$C_{\text{door}}$	boarding throughput per door, passengers per second
$n_{\text{doors}}$	doors available on the platform side
$\beta_b, \beta_a$	per-passenger boarding and alighting coefficients

In both cases, passengers who cannot board are returned to the platform queue and wait for the next train.

Why the two models differ

Modelling the two networks identically would misrepresent one of them. A fixed-dwell model applied to the Metro would remove the cost of concentrating passengers; a demand-dependent model applied to the REM would invent a cost its automation does not incur.

That the two operators behave differently is itself relevant to the project. Coordinating networks run on different operating principles is part of what makes the cross-operator case distinct from coordination within a single agency.

Where the numbers come from

$d_{\text{scheduled}}$ is an operator setting, published or directly observable.

$C_{\text{door}}$ is measured in the field by counting passengers through a door over a known interval, using the same method as the corridor throughput in section 5.

$\beta_b$ and $\beta_a$ are taken from TCRP Report 13, Rail Transit Capacity, which fits them by regression against observed rail dwell times. Montréal-specific coefficients are not published, so these are treated as published rather than measured, and are swept in the sensitivity analysis. Field sessions include a small set of Metro dwell timings at peak and off-peak as a plausibility check: if observed dwell does not vary between the two, the demand term is unnecessary and the Metro can be modelled with a fixed dwell as well.

Train capacity figures for the Azur and MR-73 rolling stock are published by the manufacturer and operator.

Why we use it

Holding would look cheaper than it is if every train had a fixed stop and anyone could board. This section adds the limits that prevent that.

The two networks are modelled differently because they are operated differently:

• REM: the doors close on schedule, so a hold does not lengthen the stop, but a large group may not fit through the doors in time and some are left behind.
• Metro: the operator controls the doors, so a hold that brings a large group also lengthens the stop and delays the trains behind it.

Passengers left behind feed directly into the missed-connection results.

Where we use it

Train boarding in the simulation engine. The boarding constraint is applied at every stop, and dwell is recalculated per stop on the Metro side.

---

5. Queueing at the transfer corridor

Concept: Deterministic queueing; cumulative arrival and departure curves

The observation

The transfer corridor only admits so many people per minute. When a large group arrives at once, the excess has to wait rather than passing through instantly. This is why a crowded transfer takes longer than an empty one even though the walking distance is identical.

General form

Queueing can be analysed by plotting two running totals against time: $A(t)$, how many people have arrived at the bottleneck so far, and $D(t)$, how many have passed through.

$$Q(t) = A(t) - D(t)$$

$$\text{People waiting} = \text{people who have arrived} - \text{people who have already gotten through}$$

The vertical gap between the two curves is the queue length. The horizontal gap is how long a given person waited. Little's law gives the related steady-state result, that average queue length equals arrival rate times average wait:

$$L = \lambda W$$

We do not use this directly, since we track individuals rather than averages, but it is the standard relationship and worth recognizing.

Our form

Because we simulate individuals rather than averages, we do not need a statistical queueing model. We track the two curves directly, with the departure curve limited by the corridor's capacity:

$$D(t) \text{ has slope} \leq C$$

$$Q(t) = A(t) - D(t)$$

$$\text{Delay}(k) = D^{-1}(k) - A^{-1}(k)$$

Symbol	Meaning
$A(t)$	cumulative arrivals at the bottleneck
$D(t)$	cumulative departures through it
$C$	corridor capacity in passengers per minute
$k$	an individual passenger

The last line needs unpacking. $A(t)$ answers "how many people have arrived by time $t$". Reading it backwards, $A^{-1}(k)$ answers "at what time had $k$ people arrived", which is the moment passenger number $k$ showed up. Likewise $D^{-1}(k)$ is the moment passenger $k$ got through the bottleneck.

Note this is the inverse function, not one divided by $A$.

So the difference between the two is how long that passenger spent waiting. In code this is simply the gap between two timestamps on the same passenger record, which is why the simulation does not need the curves at all. They are the way to reason about it on paper.

Why we use it

A crowded transfer takes longer than an empty one even though the walking distance is the same, because the corridor can only pass so many people per minute.

This matters most at high transfer volume, which is when holding looks most attractive. Modelling the queue shows what happens as volume passes the corridor's capacity, instead of assuming everyone walks through freely.

Each passenger's timestamps are enough to produce the queue, so the only new number needed is the corridor capacity $C$, measured on site.

Where we use it

Passenger movement between platforms. The capacity $C$ is measured on site at each interchange, as described in the field data collection plan.

---

6. Delay propagation and recovery

Concept: Recurrence relation with a maximum operator; saturating absorption

The observation

Two mechanisms determine how a delay spreads along a line.

The first is the minimum safe separation between consecutive trains. A train cannot close the gap to the one ahead of it, so if that train is running late, the following train must slow down as well. This is how a single delayed train becomes several.

The second is terminal recovery. The schedule allocates a layover at the end of each line before a train begins its return trip. A train arriving slightly late can still depart on time, since the layover absorbs the delay. A train arriving well behind schedule cannot, and carries the remaining lateness into its next trip.

General form

A recurrence relation defines each value in terms of the preceding one. Where a lower bound applies, a maximum operator appears:

$$x_n = \max\left(f(x_{n-1}),\ c\right)$$

This is the same idea as recursion: each value is computed from the one before it, working along the sequence. The maximum operator adds a floor, so that no value can fall below a set limit.

Saturating absorption describes a quantity absorbed up to a fixed limit, with only the excess passing through:

$$y = \max(0,\ x - k)$$

Our form

Applying the separation constraint to successive trains:

$$\text{arrival}i = \max\left(\text{natural}i,\ \text{arrival}{i-1} + h{\min}\right)$$

$$\text{departure}_i = \max\left(\text{arrival}i + \text{dwell}i,\ \text{departure}{i-1} + h{\min}\right)$$

Applying absorption at the terminus:

$$\text{carried} = \max(0,\ \text{delay} - \text{slack})$$

Symbol	Meaning
$h_{\min}$	minimum safe separation between consecutive trains
$\text{natural}_i$	when train $i$ would arrive with no train ahead of it
$\text{slack}$	scheduled layover at the terminus

Why these expressions carry so much

The first pair constitutes the entire knock-on delay model. Each train takes the later of two values: when it would have arrived unobstructed, or the earliest arrival permitted by the train ahead. Applied recursively down the line, this propagates a single delay to every following train with no additional mechanism required.

The absorption expression determines whether delay accumulates over a service day. If holds consume less than the available layover, lateness dissipates at each terminus. If they consume more, the excess carries forward and the line runs progressively later as the day proceeds. The effect is only observable over several hours, which is why full-day simulation forms part of the evaluation.

The separation constraint is enforced by the simulator rather than by any control strategy, so no strategy can bypass it.

Why we use it

A hold does not only affect the held train. The train behind cannot close the gap to it, so it slows as well, and that delay can pass along the line.

These two rules let TransitSync measure:

• how much knock-on delay a hold causes
• how long the line takes to recover
• whether lateness builds up over a full service day

The separation rule is enforced by the simulator, not by any strategy, so a strategy cannot look better by ignoring it.

Where we use it

Train movement in the simulation engine.

---

7. Prediction error model

Concept: Residuals; error variance that depends on a covariate (how wrong the controller's predictions are)

The observation

Trains do not run exactly to schedule. Run times, earlier delays and crowding all vary from trip to trip, so the actual arrival time is only known once the train arrives. Before that, a controller has the schedule or a live estimate, and both can be wrong.

The controller therefore works from a prediction, and predictions are wrong by varying amounts. A forecast made 30 seconds before an arrival is usually close. One made 5 minutes before is much less reliable.

If the controller were given the true arrival time, it would make better decisions than any real system could. The results would then describe an unachievable ceiling instead of something that could be deployed. The simulation therefore gives the controller a prediction with realistic error and never the true value.

General form

A residual is the gap between what happened and what was predicted:

$$e_i = y_i - \hat{y}_i$$

Here $y_i$ is the actual value and $\hat{y}_i$ is the prediction. The hat marks a value calculated from a model, as in section 3.

The spread of the residuals is described by their variance, $\sigma^2$. Its square root, $\sigma$, is the standard deviation. It is in the same units as the errors (here seconds) and is easier to picture, since it is roughly the typical size of an error. For a bell curve, about two thirds of errors fall within one standard deviation of the average.

If the spread is the same in every situation, the errors are homoscedastic. If it depends on something else, they are heteroscedastic, and the variance is written as a function of that quantity, here called $z$:

$$\text{homoscedastic: } \text{Var}(e) = \sigma^2$$

$$\text{heteroscedastic: } \text{Var}(e \mid z) = \sigma^2(z)$$

The symbol $\mid$ means "given", so $\text{Var}(e \mid z)$ is the variance of $e$ once the value of $z$ is known.

Our form

The quantity the error depends on is the prediction horizon $\tau$: how many seconds ahead of the arrival the prediction was made. Here the error is defined as predicted minus actual, the opposite sign to the general form, so that a positive $e$ means the prediction was later than reality:

$$e = t_{\text{predicted}} - t_{\text{actual}}$$

Errors observed in real data are modelled as:

$$e \mid \tau \sim \text{Normal}\left(\mu(\tau),\ \sigma^2(\tau)\right)$$

Normal is the bell curve with average $\mu(\tau)$ and variance $\sigma^2(\tau)$, and $\sim$ means "follows", as in section 3.

The controller is given the true time plus a scaled error drawn from that distribution:

$$\text{estimate} = t_{\text{actual}} + \kappa \cdot e$$

Symbol	Meaning
$\tau$	prediction horizon: how far ahead of the arrival the prediction is made
$\mu(\tau), \sigma(\tau)$	average error and spread of the error at that horizon
$e$	one error drawn from the distribution above
$\kappa$	error multiplier
$t_{\text{actual}}$	the true arrival time, known to the simulation but hidden from the controller

Worked example: a REM train will actually arrive in 120 seconds. Suppose, with illustrative values, that the error spread is $\sigma = 20$ s at a horizon of 120 s and $\sigma = 5$ s at a horizon of 30 s, with $\mu$ close to 0. With $\kappa = 1$, the controller's estimate is within about 20 s of the truth roughly two thirds of the time, so it will usually think the train arrives between 100 and 140 seconds from now. With $\kappa = 0.5$ the typical error halves to 10 s. With $\kappa = 0$ the estimate is exactly 120 s.

The same approach is applied to the estimated number of transferring passengers and to the onboard load. Each has its own multiplier, so the effect of each kind of error can be seen separately.

Where the values come from

$\mu(\tau)$ and $\sigma(\tau)$ are estimated by comparing what was predicted with what actually happened, grouped by how far ahead each prediction was made. This needs an archive that contains both. If the archived data only records what happened and not what was predicted, the errors cannot be measured, and the values are taken from published research on transit arrival prediction instead. No real-time data is published for the Metro, so errors for Metro arrivals are assumed in the same way. Errors in passenger counts are also assumed, since no archive records them. Each case is recorded in the parameter provenance system as measured, published or assumed.

$\kappa$ is swept. Like $M$ in section 2, it has no single correct value, so the simulation is run across a range of values. $\kappa = 1$ reproduces the error level seen in the real data, $\kappa < 1$ simulates better prediction than exists today, and $\kappa > 1$ simulates worse.

The sweep gives the project's accuracy threshold: how large the error can get before coordination stops being worthwhile. Better train position data is what CBTC provides (the new signalling system being installed on the Metro's blue line), so this turns "better data would help" into a number.

Why we use it

A real controller never knows exactly when a train will arrive. Giving the simulated controller a prediction with realistic error means the results describe what a deployable system could achieve.

Sweeping the error multiplier $\kappa$ also answers a question that is otherwise only argued: how accurate do predictions need to be before coordination is worthwhile?

Arrival times, passenger counts and onboard loads each have their own multiplier, so the project can see which improvement matters more.

Perfect-information mode ($\kappa = 0$) gives the upper bound, so the cost of imperfect information can be shown.

Where we use it

The observation layer between the simulation and the control strategy. Every value the controller receives passes through it, and the controller cannot reach the true values behind it.

Perfect-information mode sets $\kappa = 0$. It is run once and reported explicitly as an upper bound, not as something a real system could achieve.

---

8. Hold duration search and optimization

Concept

Discrete optimization and model predictive control.

The best hold duration is not necessarily obvious. A longer hold may allow more passengers to connect, but it also delays passengers who are already onboard and may create knock-on delays for trains behind.

Instead of choosing one fixed hold duration, TransitSync evaluates several possible durations and selects the one that gives the best overall result.

Discrete optimization

The simulation tests candidate hold durations in fixed increments up to the maximum allowed hold.

For example, if the maximum hold is 60 seconds and the increment is 10 seconds, the candidates are:

$$ 0, 10, 20, 30, 40, 50, 60 $$

For each candidate, the simulation calculates the objective value.

The best candidate is:

$$ h^* = \arg\max_{h \in X} J(h) $$

where:

• $h$ = candidate hold duration
• $X$ = set of allowed hold durations
• $J(h)$ = objective value for hold duration $h$
• $h^*$ = selected hold duration

The objective can combine the benefit to connecting passengers with the cost to onboard passengers and downstream trains.

Our form

For each candidate hold duration:

1. Estimate how many additional passengers make the connection.
2. Calculate passenger time saved by those additional connections.
3. Calculate the time cost to passengers already onboard.
4. Calculate any downstream or knock-on delay.
5. Combine these values into the objective.
6. Select the candidate with the highest objective value.

For example:

$$ J(h) = \text{Benefit}(h) - \text{Cost}(h) $$

The selected hold is:

$$ h^*=\arg\max_h J(h) $$

The simulator enforces a maximum hold, so an impractical duration cannot be applied.

Rolling-horizon optimization

A hold that is best at one station can delay a train that passengers need further along, so the best hold depends on the holds chosen at other stations.

TransitSync can use a rolling horizon.

At each decision point, the system looks ahead over a fixed time window and evaluates the possible decisions for several upcoming trains and stations.

It then applies the current decision and advances the simulation.

When the next decision point is reached, the forward window moves ahead and the optimization is performed again using the latest information.

This is similar to model predictive control (MPC): the system repeatedly predicts what will happen over a future window, chooses the best current action, then re-optimizes as new information becomes available.

Why we use it

Discrete optimization allows TransitSync to compare practical hold durations instead of assuming that one fixed value is always best.

The rolling horizon allows the decision to adapt to changing conditions such as:

• transfer passenger volume
• train headways
• onboard passenger load
• existing delays
• downstream congestion
• the connection window

Where it is used

Rolling-horizon optimization is used in the full-day simulation when decisions need to be updated as the simulated day progresses.

---

9. Common random numbers

Concept

Variance reduction, paired difference testing, and independent random substreams.

The simulation contains random events such as passenger walking times, transfer group sizes, run times between stations, and disruptions. If two strategies are tested using completely different random scenarios, some of the difference in their results may come from randomness rather than from the strategies themselves.

Common random numbers solve this by giving different strategies the same underlying random scenario.

For example, when comparing:

• a 30-second hold
• no hold

both simulations can use the same passenger arrivals, walking times, train conditions, and disruptions. The only intended difference is the control strategy.

General form

For two strategies $A$ and $B$, let:

$$ D = Y_A - Y_B $$

The variance of their difference is:

$$ \text{Var}(D) = \text{Var}(Y_A) + \text{Var}(Y_B) - 2,\text{Cov}(Y_A, Y_B) $$

When the two outputs are positively correlated, the covariance term reduces the variance of the difference.

This makes the comparison more precise without requiring as many simulation replications.

Our form

Each replication $r$ has its own master seed. Both strategies in that replication use the same seed, so they see the same random scenario, while different replications use different seeds and so cover different scenarios. The difference in each replication is:

$$ D_r = Y_{\text{hold},r} - Y_{\text{no-hold},r} $$

where:

• $r$ = replication number
• $Y_{\text{hold},r}$ = result from the hold strategy in replication $r$
• $Y_{\text{no-hold},r}$ = result from the no-hold strategy in replication $r$
• $D_r$ = difference between the two strategies under the same random scenario

The average difference across $R$ replications is:

$$ \bar{D} = \frac{1}{R}\sum_{r=1}^{R} D_r $$

Independent random substreams

The simulation can separate random-number generation into independent substreams for different sources of randomness, such as:

• passenger walking times
• transfer group sizes
• run times between stations
• disruptions
• onboard passenger loads
• prediction errors (section 7)

This prevents changing one random input from unintentionally shifting all the other random draws.

Prediction errors especially need their own stream, because different strategies make different numbers of decisions. If the errors shared a stream with the walking times, a strategy that decides more often would shift every walking time after it, and the two runs would no longer have the same passengers.

For example, changing the walking-time distribution should not cause a completely different sequence of disruption events.

Why we use it

Common random numbers make strategy comparisons more precise because both strategies experience the same underlying conditions.

This is especially useful when comparing:

• hold vs. no hold
• different hold durations
• different passenger-volume assumptions
• alternative control strategies

The goal is not to make simulation variance disappear. The goal is to reduce the random variation in the difference between strategies, making it easier to determine whether an observed improvement is caused by the strategy itself.

Where it is used

Common random numbers are used in the simulation and decision engine when comparing alternative strategies.

---

10. Confidence intervals and replication count

Concept

Interval estimation and statistical power analysis.

A single simulation run is not enough to determine how well a strategy performs because the result depends on random events.

For example, one simulation day might contain unusually short walking times and few disruptions, while another might contain longer walking times and several delays.

TransitSync therefore runs the simulation multiple times and reports:

• the average result
• the amount of variation between runs
• a confidence interval around the average
• the number of replications used

General form

Suppose a simulation produces $R$ results:

$$ Y_1,Y_2,\ldots,Y_R $$

The sample mean is:

$$ \bar{Y}=\frac{1}{R}\sum_{r=1}^{R}Y_r $$

The sample standard deviation is:

$$ s=\sqrt{\frac{1}{R-1}\sum_{r=1}^{R}(Y_r-\bar{Y})^2} $$

The standard error of the mean is:

$$ SE=\frac{s}{\sqrt{R}} $$

A confidence interval for the mean is:

$$ \bar{Y}\pm t_{\alpha/2,R-1}SE $$

where:

• $R$ = number of simulation replications
• $\bar{Y}$ = average simulation result
• $s$ = sample standard deviation
• $SE$ = standard error
• $t_{\alpha/2,R-1}$ = critical value from the $t$ distribution
• $\alpha$ = significance level

A 95% confidence interval uses $\alpha=0.05$.

Comparing two strategies

For TransitSync, the most important comparison is often between two strategies under the same random scenario.

Using the paired differences from Section 9:

$$ D_r=Y_{\text{hold},r}-Y_{\text{no-hold},r} $$

The average difference is:

$$ \bar{D}=\frac{1}{R}\sum_{r=1}^{R}D_r $$

The confidence interval is:

$$ \bar{D}\pm t_{\alpha/2,R-1} \frac{s_D}{\sqrt{R}} $$

where $s_D$ is the sample standard deviation of the paired differences.

If the confidence interval for the paired difference does not include zero, there is evidence that the two strategies have different average outcomes.

Choosing the number of replications

More replications generally produce a more precise estimate.

Because:

$$ SE=\frac{s}{\sqrt{R}} $$

the uncertainty decreases approximately with $1/\sqrt{R}$.

This means that increasing the number of replications improves precision, but with diminishing returns.

TransitSync can first run a pilot set of replications to estimate the variability of the results. That estimate can then be used to determine whether more replications are needed.

Statistical power

When comparing two strategies, the simulation should have enough replications to detect a difference that is meaningful for the project.

A simplified planning formula for paired differences is:

$$ R\approx \left( \frac{(z_{1-\alpha/2}+z_{1-\beta})\sigma_D} {\delta} \right)^2 $$

where:

• $R$ = required number of replications
• $\sigma_D$ = estimated standard deviation of the paired differences
• $\delta$ = minimum difference that the project wants to detect
• $\alpha$ = significance level
• $\beta$ = probability of failing to detect the specified difference
• $1-\beta$ = statistical power

The value of $\delta$ should represent a difference that is practically meaningful for TransitSync, rather than simply a very small numerical difference.

Why we use it

Confidence intervals show how uncertain the simulation results are.

Replication planning helps ensure that TransitSync runs enough scenarios to distinguish a meaningful strategy improvement from normal random variation.

This is especially important when two strategies produce results that are close to each other.

Where it is used

Confidence intervals and replication analysis are used in the results and analysis layer when comparing:

• hold vs. no hold
• different hold durations
• alternative control strategies
• different passenger-volume assumptions
• different operating conditions

The final results should record the number of replications, average outcome, standard deviation, confidence interval, and random-seed configuration so that the analysis can be reproduced.

---

11. Sensitivity analysis

Concept

One-at-a-time sensitivity analysis and finite-difference approximation.

Several inputs to the simulation are uncertain, such as transfer volume, walking-time spread, and prediction error. Some of them change the result a lot and some barely change it.

Sensitivity analysis finds out which is which. It changes one input at a time across the range of values it could realistically take, keeps every other input at its default, and records how much the result moves.

General form

The sensitivity of a result $Y$ to an input $x_i$ is how much $Y$ changes when $x_i$ changes and everything else stays the same. In calculus this is the partial derivative:

$$ \frac{\partial Y}{\partial x_i} $$

"Partial" means all the other inputs are held fixed.

The simulation is not a formula that can be differentiated, so the change is estimated instead. Move the input by a small amount $\Delta$ and see how far the result moves. This is called a finite difference.

$$ \frac{\partial Y}{\partial x_i} \approx \frac{Y(x_i + \Delta) - Y(x_i)}{\Delta} $$

Our form

TransitSync uses a range-based sensitivity measure rather than estimating the derivative directly. It compares the result at the low and high ends of a parameter's realistic range:

$$ S_i = \left| Y(x_i^{\text{high}}) - Y(x_i^{\text{low}}) \right| $$

where:

• $Y$ = the result being measured, by default the average net benefit of holding
• $x_i^{\text{low}}, x_i^{\text{high}}$ = the lowest and highest values that parameter $i$ could realistically take
• $S_i$ = how much the result changes because of parameter $i$

For each parameter:

1. Set every parameter to its default.
2. Run the simulation with the parameter at its low value.
3. Run it again with the parameter at its high value.
4. Calculate $S_i$.

The parameters are then sorted by $S_i$ and drawn as horizontal bars from longest to shortest. This is called a tornado chart, because of its shape.

For example, if changing transfer volume from 40 to 120 passengers changes the result by 3,400 passenger-seconds, while changing corridor capacity from 20 to 40 passengers per minute changes it by only 100, then transfer volume matters far more than corridor capacity.

Why we use it

Not every uncertain input can be measured, so the analysis shows where the uncertainty matters:

• A parameter with a large $S_i$ that is labelled assumed is the one to measure next.
• A parameter with a small $S_i$ can be fixed at its default, so the answer does not depend on how well it is known. It has little effect over the tested range, so it may not need further measurement or a detailed sweep.

It also keeps the full sweep manageable. The number of settings in a sweep is the number of values for each parameter multiplied together. With 10 parameters at 5 values each, that is almost 10 million settings.

The one-at-a-time analysis needs far fewer. Each parameter is tested at two values, a low and a high, so 10 parameters need 10 × 2 = 20 settings. One more setting, with every parameter at its default, gives the reference line on the chart. That makes 20 + 1 = 21 settings.

Each setting is still run several times (section 10). With 20 runs per setting, the 21 settings come to 420 runs, compared with almost 200 million for the full sweep.

Every run uses the same seeds (section 9), so a difference between the low and high results comes from the parameter and not from random luck.

One limit is that changing one parameter at a time cannot detect cases where a parameter only matters when another has a certain value. The break-even surface (section 2) covers the main one, transfer volume against headway.

Where it is used

Sensitivity analysis is used in the results and analysis layer. It is run early, before the full sweep, because its result decides which parameters the sweep needs to include.

---

12. Loop gain

Concept: Feedback system stability; spectral radius

The observation

A delay or disruption at one station can affect the next station, which can affect another train, which can eventually feed back into the original part of the network.

This creates a feedback loop. If each step of the loop reduces the disturbance, the effect eventually dies out. If the loop passes enough of the disturbance back into the system, the effect can persist or grow.

Loop gain measures this effect.

General form

For a simple loop with per-hop attenuation factors $a_1,a_2,\ldots,a_n$, the total loop gain is the product:

$$ G = \prod_{i=1}^{n} a_i $$

The interpretation is:

• $G < 1$: the disturbance decays each time it travels around the loop.
• $G = 1$: the disturbance neither grows nor decays.
• $G > 1$: the disturbance grows as it circulates through the loop.

For a network with multiple interacting loops, the equivalent measure is the spectral radius of the system matrix $A$:

$$ \rho(A)=\max_i |\lambda_i| $$

where $\lambda_i$ are the eigenvalues of $A$.

The stability condition is:

$$ \rho(A)<1 $$

A spectral radius below 1 means that repeated propagation of a disturbance becomes smaller over time.

Our form

If holds at different stations could trigger one another, the influence of each on the next would be written as a propagation factor. The numbers below are made up to show the arithmetic.

For example, suppose a disturbance passes through three stages with attenuation factors:

$$ a_1=0.8,\qquad a_2=0.6,\qquad a_3=0.5 $$

The loop gain is:

$$ G=0.8(0.6)(0.5)=0.24 $$

Only 24% of the original disturbance remains after one complete cycle. Repeating the cycle makes the disturbance progressively smaller.

If instead:

$$ G=1.2 $$

the disturbance becomes 20% larger after each cycle, indicating an unstable feedback loop.

For multiple interacting station loops, the propagation factors form a matrix:

$$ A= \begin{bmatrix} a_{11} & a_{12} & \cdots\ a_{21} & a_{22} & \cdots\ \vdots & \vdots & \ddots \end{bmatrix} $$

The spectral radius $\rho(A)$ provides a single measure of whether disturbances tend to decay or amplify across the network.

Why we would use it

This provides a mathematical way to distinguish recovering systems from systems where delays can reinforce one another.

A system with:

$$ \rho(A)<1 $$

tends to absorb disturbances over repeated propagation, while:

$$ \rho(A)>1 $$

indicates that disturbances can amplify without additional recovery mechanisms.

It is included because the project's scope can grow. Once holds at one station can trigger holds at another, a hold that looks harmless locally could circulate through several stations, and loop gain is the check for whether that dies out or grows.

Where we would use it

Loop gain is not calculated in the current scope. A feedback loop through the controller needs one hold to be able to trigger another, and the controller only protects REM to Metro connections. A Metro line held at McGill arrives late at Berri, but no controller operates there, so no second hold is triggered and the chain ends.

Loops would appear if Metro-to-Metro coordination were added, which is listed as future work, and loop gain would be calculated before turning it on. Until then, the one-hop constraint in US-2.4 limits how far a hold can spread, and delay building up over a day is covered by the layover rule in section 6.