Stochastic Systems

Motivation

Airlines sell seats on a single flight in several fare classes. Discount classes book first, so the carrier must decide how many seats to protect for later, higher-paying passengers. The classical answer, from Littlewood's rule for two classes to the nested optimality conditions of Brumelle and McGill (1993), computes optimal protection levels from a forecast of the demand distribution of each class. Those forecasts must be built from censored booking data, which is where most of the practical difficulty lies.

van Ryzin and McGill (2000) observed that the optimality conditions themselves can drive an adaptive rule. After each flight one only records whether a fill event occurred (whether demand reached each protection level) and moves the levels by a stochastic-approximation step. No forecast and no optimization is needed. Their Theorem 1 asserts that this rule converges almost surely to the optimal protection levels, with an explicit mean-square rate. It is one of the standard references for model-free revenue management. The convergence proof combines Robbins–Monro theory, the Robbins–Siegmund supermartingale lemma and an induction over fare classes.

Setting

There are $k+1$ fare classes with fares $f_1 > f_2 > \cdots > f_{k+1} > 0$ and discount ratios $r_{i+1} = f_{i+1}/f_1$. On a flight the class demands $X_1, \dots, X_{k+1}$ are independent, nonnegative and continuously distributed. Successive flights $X^1, X^2, \dots$ are independent with the same law.

A protection vector is $\theta = (\theta_1, \dots, \theta_k)$. The fill events are

$$A_i(\theta, X) = \{X_1 > \theta_1,\ X_1 + X_2 > \theta_2,\ \dots,\ X_1 + \cdots + X_i > \theta_i\}, \qquad i = 1, \dots, k.$$

A vector $\theta^*$ is characterised by the optimality condition (3), $P(A_i(\theta^*, X)) = r_{i+1}$ for $i = 1, \dots, k$. The adjustments are $H_i(\theta, X) = r_{i+1} - \mathbf 1(A_i(\theta, X))$, and their means are $h_i(\theta) = r_{i+1} - P(A_i(\theta, X))$. From an arbitrary $\theta^1$ the algorithm updates

$$\theta^{n+1} = \theta^n - \gamma_n H(\theta^n, X^n), \qquad \gamma_n = \frac{A}{n + B},\ A > 0,\ B \ge 0.$$

The interim protection levels $p_i(\theta) = \max\{\theta_j : 1 \le j \le i\}$ are what the airline actually uses to control bookings. The updates themselves monitor $A_i(\theta^n, X^n)$.

Formalization targets

Goal: Theorem 1 (p. 765)

Assume:

• A1: every $X_i$ has bounded support.
• A2 (windowed): for each $R > 0$ there is $\delta_R > 0$ with $(\theta_i - \theta^*_i)\,h_i(\theta_i, \theta^*_{i-1}, \dots, \theta^*_1) \ge \delta_R|\theta_i - \theta^*_i|^2$ whenever $|\theta_i - \theta^*_i| \le R$.
• A3: the distribution functions of the partial sums $X_1 + \cdots + X_i$ are Lipschitz.

Then for $i = 1, \dots, k$,

$$\theta^n_i \to \theta^*_i \ \text{a.s.}, \qquad p_i(\theta^*) = \theta^*_i,$$

and there are $\beta > 0$ and $C$ with

$$E|\theta^n_i - \theta^*_i|^2 \le C\,\gamma_n^{\beta/2^{i-1}}.$$

The exponent is left existential, as in the paper, so any improvement of the rate constant still proves the goal.

Milestones

The milestones follow the paper's own proof, in attack order:

• the base case $i = 1$ (p. 765);
• the Lipschitz bound on $h_{i+1}$ (p. 766);
• Lemma 3, boundedness of the iterates (p. 764);
• the almost-supermartingale inequality (12) (p. 766);
• Lemma 1 (p. 764);
• the a.s. summability of $\gamma_n|\theta^n_i - \theta^*_i|$ (p. 766);
• Lemma 2, Robbins–Siegmund (p. 764);
• the induction step for (9) (p. 766);
• Lemma 4, the rate recursion (p. 764).

Significance

Theorem 1 says that the nested optimality conditions of Brumelle and McGill can be reached by observing binary fill events alone. It also says that the limit is automatically ordered, $\theta^*_1 \le \cdots \le \theta^*_k$, so the interim levels coincide with the limit. The algorithm is distribution-free and needs neither demand forecasts nor uncensoring, which is why it is used in practice and as a baseline for later data-driven and learning-based revenue management.

The theorem is proved in the paper. As far as is known it has no machine-checked proof. Formalizing it requires a vector Robbins–Monro argument with coupled coordinates, the Robbins–Siegmund lemma, and an explicit polynomial rate. Two of the milestones are general probability facts reusable far beyond this mission: Lemma 1, and Robbins–Siegmund, which no platform item states yet. The formalization also corrects the record on two printed slips, described under Formalization scope.

Difficulty

Coordinate $i+1$ is not a Robbins–Monro process by itself. Its fill event $A_{i+1}(\theta, X)$ depends on all of $\theta_1, \dots, \theta_{i+1}$, so the drift of $\theta^n_{i+1}$ is perturbed by the errors of the lower coordinates. Applying a scalar convergence theorem coordinate by coordinate fails at this point. The perturbation has to be bounded via the Lipschitz assumption A3 and shown to be summable along the path, which uses the mean-square rate already established for class $i$. The rate and the almost-sure convergence must therefore be carried through the induction together. The exponent halves at each class. A second difficulty is that A2 is only usable on a bounded set, so boundedness of the iterates (Lemma 3) has to come first and must be uniform.

Formalization scope

Classes are indexed from $1$ in Lean ($\mathbb N$-indexed vectors whose index $0$ and indices above $k+1$ are never read). A flight's demand vector has the product law $\bigotimes_i \nu_i$, and the sample path has the countable product of copies of it. Lean's iterate … n is the paper's $\theta^{n+1}$, so the rate pairs it with $\gamma_{n+1}$. Continuity of a demand law is "every singleton is null". Expectations are Bochner integrals, and integrability is part of each rate conclusion. The fill event reuses the published NestedSeatAlloc.ProbCond.nestEvent.

Pinned and corrected statements:

• A2 as printed is unsatisfiable. It asks for one $\delta$ valid for all $\theta_i$. Since $|h_i| \le 1$, the inequality fails once $|\theta_i - \theta^*_i| > 1/\delta$, so the printed Theorem 1 is vacuously true. The mission uses the windowed A2 above. That is what the proof applies to bounded iterates, and what the paper's sufficient condition (a density bounded below near $\theta^*$) yields. A formalization with the printed A2, or with any other hypothesis no instance satisfies, is ruled out.
• The base case. The page claims the rate "for all $0 < \beta < 1$". That is false in general. The milestone states "for some $\beta > 0$", which is all Theorem 1 uses.
• Integrability. Lemma 2 assumes each $Z_n$ integrable, so that the conditional expectation is genuine. Lemma 3's "bounded (a.s.)" is read as a bound uniform in $n$.
• Standing assumptions. Nonnegative demands and positive fares are added where the argument uses them.

Not posed: §4's modifications (projection onto the capacity, randomized integer levels, booking lead times) and the simulations of §5.

Contributions welcome: proofs of the general lemmas (Lemma 1, Robbins–Siegmund) stated for arbitrary probability spaces, measurability facts for the recursion, and the elementary Lemma 4.

Selected references

• G. van Ryzin and J. McGill, Revenue Management Without Forecasting or Optimization: An Adaptive Algorithm for Determining Airline Seat Protection Levels, Management Science 46(6), 760–775, 2000. https://doi.org/10.1287/mnsc.46.6.760.11936
• S. L. Brumelle and J. I. McGill, Airline Seat Allocation with Multiple Nested Fare Classes, Operations Research 41(1), 127–137, 1993. https://doi.org/10.1287/opre.41.1.127
• H. Robbins and D. Siegmund, A convergence theorem for non negative almost supermartingales and some applications, in Optimizing Methods in Statistics, Academic Press, 233–257, 1971. https://doi.org/10.1016/B978-0-12-604550-5.50015-8
• A. Benveniste, M. Métivier and P. Priouret, Adaptive Algorithms and Stochastic Approximations, Springer, 1990. https://doi.org/10.1007/978-3-642-75894-2
• E. Lukacs, Stochastic Convergence, 2nd ed., Academic Press, 1975.

Motivation

A simulation metamodel is a cheap surrogate for an expensive stochastic simulation: after running the simulation at a few design settings, one predicts the mean response at settings never simulated. Kriging, the interpolation method of geostatistics and of the design and analysis of computer experiments (DACE; Sacks, Welch, Mitchell and Wynn 1989; Santner, Williams and Notz 2003), treats the unknown response surface as a realization of a random field and predicts by the best linear predictor. Kriging was developed for deterministic computer codes, where an observation is the exact response. A stochastic simulation returns the response plus sampling noise, and the noise may be correlated across settings when the experimenter uses common random numbers (CRN).

Ankenman, Nelson and Staum, Stochastic Kriging for Simulation Metamodeling (Proc. 2008 Winter Simulation Conference, pp. 362–370), extend kriging to this setting. The mission formalizes §2 of that paper: the form of the optimal linear predictor, called stochastic kriging, and the formula for its mean squared error, which shows how much the simulation noise costs. The source used throughout is the WSC 2008 proceedings paper, not the later Operations Research 58(2) 2010 article, whose numbering differs.

Setting

Fix a probability space $(\Omega,\mathcal F,P)$ and a set of design settings $\mathbf x$. On replication $j=1,2,\dots$ at setting $\mathbf x$ the simulation outputs (display (3), with constant trend)

$$\mathcal Y_j(\mathbf x)=\beta_0+\mathsf M(\mathbf x)+\varepsilon_j(\mathbf x).$$

The constant $\beta_0$ is the overall mean. $\mathsf M$ is a mean-zero random field, the extrinsic uncertainty imposed by the modeller. $\varepsilon_j(\mathbf x)$ is mean-zero sampling noise, the intrinsic uncertainty of the simulation; its variance may depend on $\mathbf x$, and noises at different settings may be correlated (CRN). All of these have finite second moments, and the field is uncorrelated with the noise.

An experiment design runs $n_i\ge1$ replications at each of $k$ design points $\mathbf x_1,\dots,\mathbf x_k$. The data are the sample means (4)

$$\bar{\mathcal Y}(\mathbf x_i)=\frac1{n_i}\sum_{j=1}^{n_i}\mathcal Y_j(\mathbf x_i),\qquad\bar{\mathcal Y}=(\bar{\mathcal Y}(\mathbf x_1),\dots,\bar{\mathcal Y}(\mathbf x_k))^\top,$$

and the target at a point $\mathbf x_0$ is the noise-free response $\mathsf Y(\mathbf x_0)=\beta_0+\mathsf M(\mathbf x_0)$.

Write $\Sigma_{\mathsf M}(\mathbf x,\mathbf x')=\mathrm{Cov}[\mathsf M(\mathbf x),\mathsf M(\mathbf x')]$. Then $\Sigma_{\mathsf M}$ is the $k\times k$ matrix of these covariances at the design points, $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$ is the vector $(\mathrm{Cov}[\mathsf M(\mathbf x_0),\mathsf M(\mathbf x_i)])_{i}$, and $\Sigma_\varepsilon$ is the $k\times k$ covariance matrix of the averaged noises $\sum_{j=1}^{n_i}\varepsilon_j(\mathbf x_i)/n_i$. A linear predictor (5) is $\lambda_0+\lambda^\top\bar{\mathcal Y}$ with arbitrary weights $(\lambda_0,\lambda)\in\mathbb R\times\mathbb R^k$; its mean squared error is $\mathrm E[(\lambda_0+\lambda^\top\bar{\mathcal Y}-\mathsf Y(\mathbf x_0))^2]$. The optimal MSE $\mathrm{MSE}^\star$ is the infimum of this quantity over all weights.

Formalization targets

Goal: display (7)

There is a rule $\Xi$, assigning a positive definite matrix $\Xi(\Sigma_{\mathsf M},\Sigma_\varepsilon)$ to every pair of positive definite matrices, such that in every model as above with $\Sigma_{\mathsf M}\succ0$ and $\Sigma_\varepsilon\succ0$, and at every $\mathbf x_0$,

$$\begin{aligned}\mathrm{MSE}^\star&=\Sigma_{\mathsf M}(\mathbf x_0,\mathbf x_0)-\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top[\Sigma_{\mathsf M}+\Sigma_\varepsilon]^{-1}\Sigma_{\mathsf M}(\mathbf x_0,\cdot)\\&=\Big[\Sigma_{\mathsf M}(\mathbf x_0,\mathbf x_0)-\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top\Sigma_{\mathsf M}^{-1}\Sigma_{\mathsf M}(\mathbf x_0,\cdot)\Big]+\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top\,\Xi\,\Sigma_{\mathsf M}(\mathbf x_0,\cdot).\end{aligned}$$

The bracket is the usual kriging MSE. The goal also states that $\mathrm{MSE}^\star$ strictly exceeds it whenever $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)\neq0$. The matrix $\Xi$ is existential and depends only on the two covariance matrices, not on $\mathbf x_0$.

Milestone: display (6)

If $\Sigma_{\mathsf M}+\Sigma_\varepsilon\succ0$, the stochastic kriging predictor

$$\widehat{\mathsf Y}(\mathbf x_0)=\beta_0+\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top[\Sigma_{\mathsf M}+\Sigma_\varepsilon]^{-1}(\bar{\mathcal Y}-\beta_0\mathbf 1_k)$$

is of the form (5), attains the smallest MSE among all predictors of that form, and is the only one that does.

Further items (not milestones)

• (9)–(10): the two-point example, with closed forms for the predictor and its MSE in terms of $\tau^2$, $r_{12}$, $r_0$, $\mathsf V$, $\rho$ and $n$.
• (14): the equicorrelated example, $\mathrm{MSE}^\star=\tau^2\big(1-k r_0^2/(1+(k-1)r+\gamma/n)\big)$.
• §3, p. 365: under the Gaussian Assumption 1, $\widehat{\mathsf Y}(\mathbf x_0)=\mathrm E[\mathsf Y(\mathbf x_0)\mid\bar{\mathcal Y}]$ almost surely.

Significance

Display (7) quantifies the price of simulation noise. With exact observations the best achievable error is the kriging MSE; with noisy observations it is larger by a positive definite form in the cross-covariance vector. Every point $\mathbf x_0$ that is correlated with the design therefore loses accuracy. The two examples make this concrete. In (10) the MSE increases with the intrinsic correlation $\rho$, which is why the paper concludes that common random numbers, a standard variance-reduction device for comparing systems, do not help prediction. Display (14) is the baseline against which the paper measures the cost of estimating the noise variance. The conditional-expectation statement shows that under Gaussian assumptions stochastic kriging is optimal among all predictors, not only linear ones.

The results are classical in substance: they are the best-linear-prediction formulas of second-order random-field theory, applied to data whose covariance is $\Sigma_{\mathsf M}+\Sigma_\varepsilon$. The paper states them without proof ("we can show"). No machine-checked version exists. Formalizing them yields a reusable, measure-theoretic best-linear-predictor theorem for random vectors defined from a model rather than postulated, together with the Loewner-order fact behind the positive definiteness of $\Xi$.

Difficulty

The algebra is short; the difficulty is in deriving the second-moment structure from the model, not assuming it. The covariance of $\bar{\mathcal Y}$ must be computed from the sample-mean definition, using bilinearity of covariance over finite sums of square-integrable variables, and shown to equal $\Sigma_{\mathsf M}+\Sigma_\varepsilon$. The cross-covariance with $\mathsf Y(\mathbf x_0)$ must be shown to equal $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$. The MSE of an arbitrary affine predictor then has to be written as a quadratic in $(\lambda_0,\lambda)$, and its infimum identified.

A first idea is to take "a random vector with mean $\beta_0\mathbf 1$ and covariance $\Sigma$" as the hypothesis. That proves a different, more abstract statement and leaves the model's content unproved.

For the second line of (7), the positive definiteness of $\Xi$ does not follow from a scalar argument. It is a matrix inequality, $(A+B)^{-1}\prec A^{-1}$ for $A,B\succ0$, and it fails without $\Sigma_\varepsilon\succ0$: when $\Sigma_\varepsilon=0$ the two lines coincide.

Formalization scope

• The model. Settings form an arbitrary type X; design points are x : Fin k → X, so the paper's index $i$ is i.val + 1. Replication $j\ge1$ is index j - 1 of a sequence ε : ℕ → X → Ω → ℝ.
• Moments and covariances. Field and noise are in $L^2$ and have mean zero. Covariances are Mathlib's covariance. The optimal MSE is an infimum over all of $\mathbb R\times\mathbb R^k$, with no unbiasedness constraint and no sign constraint.
• Hypotheses not printed in §2. Three are added and disclosed:
  1. uncorrelatedness of $\mathsf M$ and $\varepsilon$, the second-order content of Assumption 1's "independent of $\mathsf M$", without which the covariance of $\bar{\mathcal Y}$ is not $\Sigma_{\mathsf M}+\Sigma_\varepsilon$;
  2. square integrability;
  3. $\Sigma_{\mathsf M}\succ0$ and $\Sigma_\varepsilon\succ0$ in the goal. Lean's matrix inverse returns $0$ for a singular matrix, and $\Xi\succ0$ is false when $\Sigma_\varepsilon$ is singular.
• Generality kept. No Gaussianity, no independence across replications and no "no CRN" restriction is imposed on (6) and (7).
• Examples. $|r_{12}|<1$ in (10) and $r\le1$ in (14) are added, since both are correlations.
• Non-trivializing. A positive semidefinite $\Xi$, or a $\Xi$ chosen after $\mathbf x_0$, would make the goal trivial or strictly weaker; the statement fixes $\Xi$ as a function of $(\Sigma_{\mathsf M},\Sigma_\varepsilon)$ before the model and the point.

Infrastructure needed and reusable beyond this mission:

• covariance of averages of $L^2$ random variables;
• the MSE of an affine predictor as a quadratic form;
• minimization of a positive definite quadratic;
• the Loewner antitonicity of the inverse;
• for the conditional-expectation item, Gaussian conditioning via uncorrelated-hence-independent residuals.

Contributions to any of these are welcome.

Selected references

• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Proceedings of the 2008 Winter Simulation Conference, IEEE, pp. 362–370, 2008. http://www.informs-sim.org/wsc08papers/042.pdf
• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Operations Research 58(2), 371–382, 2010 (later journal version, not used here). https://doi.org/10.1287/opre.1090.0754
• J. Sacks, W. J. Welch, T. J. Mitchell, H. P. Wynn, Design and Analysis of Computer Experiments, Statistical Science 4(4), 409–423, 1989. https://doi.org/10.1214/ss/1177012413
• T. J. Santner, B. J. Williams, W. I. Notz, The Design and Analysis of Computer Experiments, Springer, 2003. https://doi.org/10.1007/978-1-4757-3799-8

Motivation

Stochastic simulation models of queues, supply chains and manufacturing systems are often too slow to run at every input setting of interest. A metamodel fitted to a moderate number of simulation runs then stands in for the simulation, for instance inside an optimization or a sensitivity analysis. Kriging, developed in geostatistics and in the design and analysis of deterministic computer experiments (DACE), predicts an unknown response surface by treating it as a realization of a Gaussian random field. Deterministic kriging interpolates its data exactly, which is wrong for stochastic simulation, where each run returns the response plus sampling noise.

Ankenman, Nelson and Staum's stochastic kriging separates the two sources of uncertainty: the extrinsic uncertainty of the unknown surface and the intrinsic uncertainty of the simulation output. Its optimal predictor requires the intrinsic noise variances, which are never known in practice and are replaced by sample variances computed from the same replications. This mission formalizes the paper's first key result: this plug-in step introduces no prediction bias.

The source is the Winter Simulation Conference 2008 proceedings version of the paper (pp. 362–370). The later Operations Research 58(2) (2010) article numbers its results differently and is not the reference here.

Setting

Fix a probability space $(\Omega,\mathcal F,P)$ and design variables $\mathbf x\in\mathbb R^d$. On replication $j=1,2,\dots$ at setting $\mathbf x$ the simulation returns

$$\mathcal Y_j(\mathbf x)=\beta_0+\mathsf M(\mathbf x)+\varepsilon_j(\mathbf x),$$

where $\beta_0\in\mathbb R$ is a constant trend, $\mathsf M$ is a mean-zero random field on $\mathbb R^d$ (the unknown surface) and $\varepsilon_j(\mathbf x)$ is the noise of replication $j$. The quantity to predict at any point $\mathbf x_0$, simulated or not, is the noise-free response $\mathsf Y(\mathbf x_0)=\beta_0+\mathsf M(\mathbf x_0)$.

An experiment design is a list of pairs $(\mathbf x_i,n_i)$, $i=1,\dots,k$, with distinct design points $\mathbf x_i$ and $n_i$ replications at $\mathbf x_i$. The data are summarized by the sample means $\bar{\mathcal Y}(\mathbf x_i)=\frac1{n_i}\sum_{j=1}^{n_i}\mathcal Y_j(\mathbf x_i)$, collected in the vector $\bar{\mathcal Y}\in\mathbb R^k$, and the sample variances

$$\mathcal S^2(\mathbf x_i)=\frac1{n_i-1}\sum_{j=1}^{n_i}\big(\mathcal Y_j(\mathbf x_i)-\bar{\mathcal Y}(\mathbf x_i)\big)^2 .$$

The extrinsic covariances are the $k\times k$ matrix $\Sigma_{\mathsf M}=(\mathrm{Cov}[\mathsf M(\mathbf x_h),\mathsf M(\mathbf x_i)])_{h,i}$ and the vector $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)=(\mathrm{Cov}[\mathsf M(\mathbf x_0),\mathsf M(\mathbf x_i)])_{i}$.

Assumption 1. $\mathsf M$ is a stationary Gaussian random field: all its finite-dimensional laws are multivariate normal with mean $0$, the covariance is $\tau^2R(\mathbf x-\mathbf x')$ with $\tau^2>0$ and $R(\mathbf 0)=1$, and the covariance matrix at finitely many distinct points is positive definite. At each design point the noises $\varepsilon_1(\mathbf x_i),\varepsilon_2(\mathbf x_i),\dots$ are i.i.d. $N(0,\mathsf V(\mathbf x_i))$ with $\mathsf V(\mathbf x_i)>0$; noises at different design points are independent (no common random numbers); and the whole noise family is independent of $\mathsf M$.

The intrinsic variance at a design point is estimated by $\widehat{\mathsf V}(\mathbf x_i)=\mathcal S^2(\mathbf x_i)$, giving $\widehat\Sigma_\varepsilon=\mathrm{Diag}\{\widehat{\mathsf V}(\mathbf x_1)/n_1,\dots,\widehat{\mathsf V}(\mathbf x_k)/n_k\}$ and the plug-in stochastic kriging predictor

$$\widehat{\widehat{\mathsf Y}}(\mathbf x_0)=\beta_0+\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top\big[\Sigma_{\mathsf M}+\widehat\Sigma_\varepsilon\big]^{-1}\big(\bar{\mathcal Y}-\beta_0\mathbf 1_k\big).\tag{13}$$

Formalization targets

Goal: Theorem 1 (p. 366)

Under Assumption 1, with $n_i\ge2$ replications at every design point, the prediction error is integrable and

$$\mathrm E\Big[\widehat{\widehat{\mathsf Y}}(\mathbf x_0)-\mathsf Y(\mathbf x_0)\Big]=0 .$$

Milestones (p. 365)

1. Under Assumption 1, $(\mathsf Y(\mathbf x_0),\bar{\mathcal Y}(\mathbf x_1),\dots,\bar{\mathcal Y}(\mathbf x_k))$ is multivariate normal.
2. Under Assumption 1, $\mathcal S^2(\mathbf x_i)$ has a scaled chi-squared distribution: $(n_i-1)\mathcal S^2(\mathbf x_i)/\mathsf V(\mathbf x_i)\sim\chi^2_{n_i-1}$.
3. Under Assumption 1, $\mathcal S^2(\mathbf x_i)$ is strongly consistent for $\mathsf V(\mathbf x_i)$: the sample variance of the first $m$ replications converges to $\mathsf V(\mathbf x_i)$ almost surely as $m\to\infty$.

Significance

Theorem 1 says that the practitioner's shortcut, estimating the noise variances from the replications and plugging them into the optimal predictor, keeps the predictor unbiased; the cost of not knowing the intrinsic variance is paid entirely in mean squared error. This is what makes the paper's subsequent comparison of the plug-in MSE with the known-variance MSE the right measure of the penalty for estimation.

None of these statements has a machine-checked proof, and Mathlib has neither the chi-squared law of the normal sample variance nor the independence of the normal sample mean and sample variance. A complete development produces these classical facts of normal sampling theory for the first time in Lean, together with a reusable model of noisy observations of a Gaussian random field.

Difficulty

The predictor is a nonlinear function of the data: its weight vector $[\Sigma_{\mathsf M}+\widehat\Sigma_\varepsilon]^{-1}\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$ is random and depends on the same replications as $\bar{\mathcal Y}$. Linearity of expectation therefore does not apply directly: the obvious argument, which treats the weights as constants, fails, and the dependence between the estimated variances, the sample means and the field has to be controlled exactly. The joint independence structure of Assumption 1 — across replications, across design points, and between noise and field — is load-bearing; a weaker pairwise or pointwise independence does not suffice. Integrability of the error must also be established, since the weights depend on random variances.

Formalization scope

Design points are indexed by Fin k (the paper's $i$ is i.val + 1) and lie in EuclideanSpace ℝ (Fin d); replication $j\ge1$ of the paper is index j - 1 of a sequence, so sums run over Finset.range. The field is M : EuclideanSpace ℝ (Fin d) → Ω → ℝ with Mathlib's IsGaussianProcess; covariances are ProbabilityTheory.covariance; stationarity is $\mathrm{Cov}[\mathsf M(\mathbf y),\mathsf M(\mathbf y')]=\tau^2R(\mathbf y-\mathbf y')$ with $R(\mathbf 0)=1$. The noise law is HasLaw … (gaussianReal 0 (V (x i))); independence is one iIndepFun over the family indexed by Fin k × ℕ plus one IndepFun between that whole family and the whole field. The variance function is ℝ≥0-valued with $\mathsf V(\mathbf x_i)>0$ assumed; replication counts satisfy $n_i\ge2$ wherever $\mathcal S^2$ appears and $n_i\ge1$ where only means appear. The chi-squared law $\chi^2_\nu$ is gammaMeasure (ν/2) (1/2).

The covariance quantities $\Sigma_{\mathsf M}$ and $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$ in (13) are the true ones; only $\Sigma_\varepsilon$ is estimated, and $\widehat\Sigma_\varepsilon$ is built from the same replications as $\bar{\mathcal Y}$. A formalization in which $\widehat\Sigma_\varepsilon$ is deterministic or independent of the data, in which the matrix $\Sigma_{\mathsf M}+\widehat\Sigma_\varepsilon$ may be singular (Lean's Matrix.inv returns $0$ there, reducing the claim to $\mathrm E[\mathsf M(\mathbf x_0)]=0$), or in which the expectation is asserted without integrability, would trivialize the goal and is ruled out: positive definiteness of $\Sigma_{\mathsf M}$ at distinct design points is part of Assumption 1, and integrability is part of the conclusion. $\mathbf x_0$ may coincide with a design point.

A complete development needs: the law of linear images of independent Gaussian vectors; independence of the sample mean and the residual vector for i.i.d. normal samples (Cochran's theorem); the chi-squared law of the normal sample variance; the strong law of large numbers for the sample variance; and a measurability and boundedness argument for the random weight vector. The normal-sampling results are reusable well beyond this mission and are welcome as separate contributions.

Selected references

• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Proceedings of the 2008 Winter Simulation Conference, IEEE, pp. 362–370, 2008. http://www.informs-sim.org/wsc08papers/042.pdf
• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Operations Research 58(2):371–382, 2010. https://doi.org/10.1287/opre.1090.0754
• T. J. Santner, B. J. Williams, W. I. Notz, The Design and Analysis of Computer Experiments, Springer, 2003. https://doi.org/10.1007/978-1-4757-3799-8

Motivation

A G/GI/N queue has $N$ identical servers, a general arrival process, a first-come-first-served waiting room and i.i.d. service times with a general distribution $F$ of mean $1$. In the Halfin–Whitt regime the number of servers grows while the traffic intensity $\rho^N$ approaches $1$ at rate $\sqrt N(1-\rho^N)\to\beta\in\mathbb R$; this is the standard asymptotic model for large call centers, where a positive fraction of customers wait but waiting times are short. Halfin and Whitt (Oper. Res. 1981) obtained a diffusion limit for exponential service times. Reed (arXiv:0912.2837, Ann. Appl. Probab. 2009) proved the corresponding limit for general service times: the centred and scaled queue length converges (Theorem 5.1) to the solution of a nonlinear stochastic convolution equation, (5.33) below.

Equation (5.33) is written in terms of an infinite-server limit, and it is not evident from its form that it reduces to the Halfin–Whitt diffusion when service is exponential. Section 5.5 of the paper gives a second representation, Corollary 5.2, in which the convolution against $F$ is replaced by a convolution against the renewal function of $F$. With exponential service the renewal function is $M(t)=t$, and the representation becomes the Halfin–Whitt diffusion. This mission formalizes Corollary 5.2 and the renewal-theoretic steps of its proof.

Timeline. Halfin–Whitt (1981): exponential service, diffusion limit. Reed (2009): general service with finite mean; Corollary 5.2 links the general limit back to 1981. Karlin–Taylor (1975) and Ross (1983): the renewal-type equation and the renewal equation that the proof cites.

Setting

Let $\mu$ be a probability measure on $[0,\infty)$ with mean $\int x\,d\mu(x)=1$ (the service-time law), $F(t)=\mu((-\infty,t])$ its distribution function and $G=1-F$ its tail. The equilibrium distribution is

$$F_e(x)=\int_0^x G(u)\,du,\qquad x\ge0.\tag{5.4}$$

Let $\mu^{*n}$ be the $n$-fold convolution of $\mu$ (the law of the $n$th renewal epoch $S_n$). The renewal measure is $dM=\sum_{n\ge1}\mu^{*n}$ and the renewal function is $M(t)=\sum_{n\ge1}F^{*n}(t)$, the expected number of renewals by time $t$.

A path $x:\mathbb R\to\mathbb R$ is càdlàg on $[0,\infty)$ if it is right continuous at every $t\ge0$ with finite left limits at every $t>0$; only its values on $[0,\infty)$ matter. Fix a càdlàg driving path $\zeta$ (in the paper $\tilde\zeta=\tilde M_Q+\tilde Q_I$, (5.40)) and $\beta\in\mathbb R$. Write $y^+=\max(y,0)$ and $y^-=\min(y,0)\le0$. The two equations are

$$q(t)=\zeta(t)-\beta F_e(t)+\int_0^t q^+(t-s)\,dF(s),\qquad t\ge0,\tag{5.33}$$ $$q(t)=\zeta(t)+\int_0^t\zeta(t-u)\,dM(u)-\beta t-\int_0^t q^-(t-u)\,dM(u),\qquad t\ge0.\tag{5.41}$$

All Stieltjes integrals are over the closed interval $[0,t]$.

Formalization targets

Goal: Corollary 5.2 (p. 29)

For every service law $\mu$, every càdlàg $\zeta$ and every $\beta$:

$$\forall q\ \text{càdlàg}:\quad q\ \text{solves (5.33)}\iff q\ \text{solves (5.41)},$$

and (5.41) has a càdlàg solution that is unique on $[0,\infty)$ among càdlàg paths.

Milestones, in proof order

1. (5.39): $M$ is finite, solves $M(t)=F(t)+\int_0^tM(t-u)\,dF(u)$, and is its unique locally bounded measurable solution.
2. (5.42)–(5.43): for locally bounded measurable $H$, $r=H+\int_0^tH(t-u)\,dM(u)$ is the unique locally bounded solution of $r(t)=H(t)+\int_0^tr(t-u)\,dF(u)$.
3. The identity $F_e(t)+\int_0^tF_e(t-s)\,dM(s)=t$ for $t\ge0$ (p. 30).
4. (5.44): every càdlàg solution of (5.33) satisfies (5.41) with $-\beta t$ replaced by $-\beta\bigl(F_e(t)+\int_0^tF_e(t-s)\,dM(s)\bigr)$.
5. (5.46) (further): for exponential service of rate $1$, $M(t)=t$ and (5.41) becomes $q(t)=\zeta(t)+\int_0^t\zeta(s)\,ds-\beta t-\int_0^tq^-(s)\,ds$.

Significance

The result. Corollary 5.2 expresses the many-server diffusion limit as an equation that is linear in the driving process $\zeta$ and whose only nonlinearity is the idle-server term $\int_0^tq^-(t-u)\,dM(u)$. This is how the paper verifies that its general-service limit agrees with Halfin and Whitt's for exponential service, and the paper's sequel interprets the integral term as the limiting idle time of the servers. Uniqueness in (5.41) also yields the paper's remark that $\tilde Q^N\Rightarrow\tilde Q_M$, the solution of (5.41).

The formalization. The result is proved in the paper, but only partly written: the printed proof shows that the solution of (5.33) satisfies (5.41), while the converse and the uniqueness are asserted. The mission states both. It also produces, independently of queueing, the renewal measure of a law on $[0,\infty)$ with possible atoms, its finiteness on bounded sets, the renewal equation and the solution of renewal-type equations, and the identity $F_e+F_e*dM=t$; none of these is in Mathlib or published on the platform. To our knowledge nothing in this mission has a machine-checked proof.

Difficulty

The forward direction is formal once (5.42)–(5.43) is available, but (5.42)–(5.43) itself needs the renewal measure to be finite on bounded intervals when $F$ may have an atom at $0$, and the uniqueness half needs an argument that does not assume the kernel has mass below one on any window starting at $0$. The converse direction, (5.41) $\Rightarrow$ (5.33), cannot be read off the printed proof, because (5.41) is not an equation of renewal type in $q$. The uniqueness of (5.41) is not a contraction argument either: the atom of $dM$ at $0$ has mass $F(0)/(1-F(0))$, which may exceed $1$, so a naive Picard estimate on $[0,\delta]$ fails.

Formalization scope

• Paths are functions $\mathbb R\to\mathbb R$; every equation is asserted for $t\ge0$ and every uniqueness is equality on $[0,\infty)$. Values at negative times are never read.
• The service law is a probability measure $\mu$ on $\mathbb R$ with $\mu((-\infty,0))=0$, integrable identity and mean $1$. The paper places "no additional restrictions on $F$ beyond a first moment".
• The renewal measure is defined as $\sum_{n\ge1}\mu^{*n}$ (Mathlib's additive convolution of measures), not through an i.i.d. sequence; the paper's "expected number of renewals" equals it by Tonelli.
• Stieltjes integrals are Lebesgue integrals over the closed $[0,t]$, so atoms at $0$ of $dF$ and $dM$ are included. $q^-$ is $\min(q,0)$, not Mathlib's negPart.
• Explicit readings: "unique strong solution" is existence and uniqueness among càdlàg paths, pathwise for each fixed $\zeta$. The paper's random statement follows by applying it to each sample path, which lies in $D[0,\infty)$ almost surely (p. 30). "Locally bounded" is bounded on every $[0,T]$, with measurability added for the functions in (5.39) and (5.42).
• Corrected slip: (5.42) prints $\int_0^tr(t)\,dF(t-u)$; the stated equation is $\int_0^tr(t-u)\,dF(u)$. Milestone (5.42)–(5.43) assumes $F$ carried by $[0,\infty)$ with $F(0)<1$, under which the cited result holds.
• Ruled out: the solution of (5.33) is not assumed to exist (it exists by the paper's Proposition 3.1 with $a=0$, $B=F$, which the solver must supply). $M$ is the renewal measure of the same $F$, never an arbitrary measure assumed to satisfy (5.39). Uniqueness is among càdlàg paths on $[0,\infty)$, so no non-measurable path makes an integral vanish.
• Not in the mission: Theorems 4.1 and 5.1 (weak convergence in the Skorohod space), whose limit this corollary rewrites, and the Brownian-motion claim after (5.46).
• Related platform items, credited but not the same object: BellmanDP.Inventory.renewal_equation_exists_unique (a renewal equation with a density kernel of mass below $1$); ChenWhitt93.Reflection.Basic (càdlàg paths on $[0,T]$). Reusable output: the renewal measure, (5.39) and (5.42)–(5.43) serve any renewal-theoretic development. Alternative proofs of the converse and of uniqueness are welcome.

Selected references

• J. Reed, The G/GI/N queue in the Halfin–Whitt regime, Ann. Appl. Probab. 19(6), 2009, 2211–2269. https://arxiv.org/abs/0912.2837 (v1), https://doi.org/10.1214/09-AAP609
• S. Halfin and W. Whitt, Heavy-traffic limits for queues with many exponential servers, Oper. Res. 29(3), 1981, 567–588. https://doi.org/10.1287/opre.29.3.567
• S. Karlin and H. M. Taylor, A First Course in Stochastic Processes, 2nd ed., Academic Press, 1975 (reference [11] of the paper).
• S. M. Ross, Stochastic Processes, Wiley, 1983 (reference [19] of the paper, Exercise 3.4).

Why the lost-sales state space matters

The lost-sales inventory system with a positive lead time is one of the basic models of inventory theory. Demand that cannot be met from stock is lost rather than backlogged, and orders arrive a fixed number $L$ of periods after they are placed. Unlike the backlog system, where the inventory position is a sufficient one-dimensional state and a base-stock policy is optimal, the lost-sales system needs the whole pipeline of outstanding orders as its state, and its optimal policy has no simple form. Computing optimal policies by dynamic programming therefore runs into the size of an $L$-dimensional state space.

Paul Zipkin's paper Old and New Methods for Lost-Sales Inventory Systems (Operations Research 56(5), 2008, doi:10.1287/opre.1070.0471) revisits this model, compares classical and new heuristics against optimal policies, and prepares the computation with a state-reduction result. Morton (1969, SIAM Review 11, Equation (83)) showed that under an optimal policy the state remains in a bounded set; Zipkin's §4 extends the argument to every policy in a class $Z(s)$ that contains the usual heuristics. This mission formalizes that section: the reduced state space $X(s)$, the policy class $Z(s)$, and Claim 1, which says that states outside $X(s)$ are transient.

Setting

Fix the lead time $L\ge 1$. The state is the $L$-vector $x=(x_0,\dots,x_{L-1})$: $x_0$ is the stock on hand after this period's delivery, and $x_l$ ($1\le l\le L-1$) is the order that arrives $l$ periods later, so $x_{L-1}$ is the most recent order. Given the state $x$, an order $z\ge 0$ and a demand $d\ge 0$, the next state is

$$x_+ = \bigl([x_0-d]^+ + x_1,\ x_2,\ \dots,\ x_{L-1},\ z\bigr),\qquad [a]^+=\max(a,0).$$

Demands in successive periods are independent with a common law $D$ on $[0,\infty)$.

The pipeline partial sums are $v_l=\sum_{m=l}^{L-1}x_m$ for $l=0,\dots,L-1$ and $v_L=0$; $v_0$ is the inventory position and $v_{L-1}=x_{L-1}$ the most recent order. A level vector is $s=(s_0,\dots,s_L)$ with $s_0\ge s_1\ge\dots\ge s_L\ge 0$. For a level vector $s$,

$$X(s)=\{x\ge 0:\ v_l\le s_l,\ l=0,\dots,L-1\},$$

and $Z(s)$ is the set of stationary policies $z$ (an order $z(x)\ge 0$ in each state $x\ge 0$) with $z(x)=0$ for $x\notin X(s)$ and $v_l+z(x)\le s_l$ for $l=0,\dots,L$ and $x\in X(s)$. The vector base-stock policy with parameter $s$ is

$$z_s(x)=\bigl[\min\{s_l-v_l,\ l=0,\dots,L\}\bigr]^+,$$

equation (5) of the paper with a general $s$ in place of Morton's $\bar s$. Under a stationary policy $z$, from an initial state $x$ and along a demand path $(d_t)$, the state process is $x_0=x$, $x_{t+1}=(x_t)_+$ with order $z(x_t)$ and demand $d_t$. In Lean these are next, v, IsLevelVector, Xs, Zs, vectorBaseStock and traj in the namespace LostSalesOldNew.StateReduction.

Formalization targets

Goal: Claim 1 (p. 1259)

For a level vector $s$, a policy $z\in Z(s)$, i.i.d. demands with law $D$ on $[0,\infty)$ and $P(d>0)>0$, and any initial state $x\ge 0$, almost surely

$$\exists\,T\ \forall\,t\ge T:\quad x_t\in X(s).$$

The paper's word is "transient"; the statement above says that every state outside $X(s)$, indeed the whole complement, is visited only finitely often.

Milestones

1. (5) and §4: for every level vector $s$, $z_s\in Z(s)$.
2. Invariance (proof of Claim 1): $x\in X(s)$, $z\in Z(s)$, $d\ge 0$ imply $x_+\in X(s)$.
3. Entrance (proof of Claim 1): along a demand path with $d_t\ge 0$ and $\sum_{t<n}d_t\to\infty$, the process started at any $x\ge 0$ visits $X(s)$.

A further item records the paper's remark that $X(s)$ and $Z(s)$ are increasing in $s$.

Significance

Claim 1 lets the dynamic program be solved on the compact set $X(\bar s)$ with policies restricted to $Z(\bar s)$, which is what makes the paper's exact computations (lead times up to four, §6) feasible. Since $Z(s)$ contains the vector base-stock policies and, for $s$ large enough, any policy that stops ordering in large states, the claim also says that every plausible policy eventually confines the state to a bounded set.

The result has a short proof in the paper and no machine-checked version. Formalizing it produces a reusable encoding of the lost-sales pipeline with general stationary policies, and makes explicit two points the paper leaves implicit: the probabilistic reading of "transient", and the demand condition under which the process must leave the complement of $X(s)$. Related platform material treats the backlog analogue (VeinottBaseStock.base_stock_level_absorbing, a base-stock level set is absorbing) and lost-sales dynamics under other conventions (CappedBaseStock_Model, from an empty initial state under capped base-stock policies with average cost); those formulations use another state convention and cost criterion and are not imported here.

Difficulty

The invariance step is algebra on partial sums. The content of the claim is the entrance step: outside $X(s)$ the policy orders nothing, so the process can only reach $X(s)$ through demand draining the stock and the pipeline. That requires the cumulative demand to diverge, which fails pathwise (demand identically zero leaves a state with $x_0>s_0$ fixed forever) and holds only almost surely, from independence and $P(d>0)>0$. The naive reading of "transient" as invariance alone, or as the pathwise statement with divergent demand assumed, misses this step.

Formalization scope

• States are Fin L → ℝ with $L\ge 1$ in the goal; theorems quantify over nonnegative states. The vector $(x,z)$ is Fin.snoc x z; $v$ takes a natural-number index and vanishes at $l\ge L$, so $v_L=0$.
• Level vectors are Fin (L+1) → ℝ, antitone with nonnegative last entry. $X(s)$ constrains $v_0,\dots,v_{L-1}$; $Z(s)$ constrains $v_l+z(x)$ for $l=0,\dots,L$.
• Policies are deterministic stationary functions of the state, with no measurability required. The zero condition and the sign condition of $Z(s)$ range over states $x\ge 0$ (the paper's state space); over all of $\mathbb R^L$ the vector base-stock policy would not be in $Z(s)$.
• Demand is a probability measure $D$ on $\mathbb R$ with $D((-\infty,0))=0$; the demand sequence is the coordinate process under Measure.infinitePi (fun _ : ℕ => D). The hypothesis $D((0,\infty))>0$ is added to the paper's assumptions; without it the claim is false.
• "Transient" is read as: almost surely the state lies in $X(s)$ from some period on.
• Ruled out: the goal is neither the invariance statement alone, which says nothing about states outside $X(s)$, nor the pathwise entrance statement with divergent cumulative demand assumed, which drops the probabilistic content.

A complete development needs the partial-sum algebra of the pipeline, the second Borel–Cantelli lemma or the strong law for i.i.d. nonnegative sequences under Measure.infinitePi, and induction along traj. The pipeline encoding is shared with the second mission of this series (lead-time monotonicity). Proofs of any milestone and of the goal are welcome.

Selected references

• P. Zipkin, Old and New Methods for Lost-Sales Inventory Systems, Operations Research 56(5):1256–1263, 2008. https://doi.org/10.1287/opre.1070.0471
• T. Morton, Bounds on the Solution of the Lagged Optimal Inventory Equation with No Demand Backlogging and Proportional Costs, SIAM Review 11:572–576, 1969 (as cited in Zipkin 2008).
• T. Morton, The Near-Myopic Nature of the Lagged-Proportional-Cost Inventory Problem with Lost Sales, Operations Research 19(7):1708–1716, 1971. https://doi.org/10.1287/opre.19.7.1708
• R. Levi, G. Janakiraman, M. Nagarajan, Provably Near-Optimal Balancing Policies for Stochastic Inventory Control Models with Lost Sales, Mathematics of Operations Research, 2008 (forthcoming at the time of Zipkin 2008).

Motivation

Queueing models often represent work arriving to and leaving a server by a real-valued input path. The resulting workload cannot fall below a barrier. Functional limit theorems describe what happens to such paths under scaling, but a limit for the input path is useful only when it can be carried through the workload transformation. Ward Whitt’s 1980 paper studies several path transformations that occur in stochastic-process limits and states conditions under which their limits can be computed in Skorohod topologies. Section 6 treats the running supremum and a reflecting barrier, including paths with a linear drift that changes with the scaling index. This mission formalizes the three-regime reflection result, Theorem 6.4, and the source results on its dependency path. Whitt (1980)

The result is already proved in the paper. Its formal Lean statement remains an open proof target here. The mission is therefore concerned with formalizing the known theorem, including its exact jump conditions and convergence modes, rather than proposing a new queueing claim. Whitt (1980), pp. 80–81

Setting

A càdlàg path on a real interval $T$ is right-continuous and has a limit from the left wherever $T$ approaches a time from the left. The collection of such real-valued paths is $D(T,\mathbb R)$. For a compact interval $[a,b]$, the $J_1$ distance between $x$ and $y$ allows an increasing homeomorphism $\lambda$ of time and measures the larger of its distance from the identity and the uniform distance between $x$ and $y\circ\lambda$. On a general interval, $J_1$ convergence is specified by restrictions to compact subintervals whose endpoints are continuity points of the limit or endpoints of $T$. The paper uses complete separable metric spaces for general path-space results; Section 6 specializes to real-valued paths and intervals with $0$ as a closed left endpoint. Whitt (1980), pp. 70, 80

For $x\in D(T,\mathbb R)$ and $t\in T$, define the running supremum $x^{\uparrow}(t)=\sup_{0\le s\le t}x(s)$ and running infimum $x^{\downarrow}(t)=\inf_{0\le s\le t}x(s)$. The reflecting-barrier map is $f(x)(t)=x(t)-x^{\downarrow}(t)$. This is the paper’s definition: it subtracts the full running infimum. The identity path is $e(t)=t$, and $\theta$ denotes the zero path. A path has no positive jumps if $x(t)\le x(t-)$ for each $t>0$ in $T$; no negative jumps reverses that inequality. Whitt (1980), pp. 80–81

Formalization targets

The goal is Theorem 6.4. Suppose $x_n-c_ne\to x$ in $J_1$ and $x(0)=0$. Its three conclusions are

$$\begin{aligned} c_n\to c &\implies f(x_n)\to f(x+ce) &&(J_1),\\ c_n\to+\infty &\implies f(x_n)-c_ne\to x &&(J_1),\\ c_n\to-\infty\ \text{and }x\text{ has no positive jumps} &\implies f(x_n)\to\theta &&(U). \end{aligned}$$

Here $U$ means uniform convergence on every compact subinterval. The goal keeps all three drift regimes together, as the printed theorem does. Its milestones are the finite-partition fact cited from Billingsley, the interval-splitting Lemma 2.2, continuity of addition at limits without shared discontinuities in Theorem 4.1, the running-supremum inequality in Theorem 6.1, the three-regime supremum theorem 6.2, and the paper’s explicit claim that $f$ is $J_1$-continuous. Whitt (1980), pp. 71, 79–81

Significance

The theorem specifies which output path appears after reflection as the linear centering changes. Finite drift leaves a reflected limit with drift included; large positive drift leaves the centered input limit; large negative drift drives the reflected path uniformly to zero when upward jumps of the limit are excluded. The distinction between $J_1$ and the stronger compact-uniform mode in the last case is part of the result. These statements make the paper’s barrier transformation usable when an input-process limit is already known. Whitt (1980), Theorem 6.4

A completed formal development would establish reusable path-space components: a precise $J_1$ convergence interface for real intervals, a theorem about addition with disjoint discontinuities, and continuity properties of running extrema and reflection. The present Lean items state these results and compile with proof placeholders. No machine-checked proof of these mission theorems is claimed here. The milestone list follows the claims Whitt states or cites, so later proofs can be attached to the same mathematical targets. Whitt (1980), §§2, 4, 6

Difficulty

A uniform estimate alone does not establish a $J_1$ limit: separate paths may converge using different changes of time, and addition can fail when the limiting paths jump together. Whitt’s Theorem 4.1 gives continuity precisely at pairs with disjoint discontinuity sets. In the reflection theorem, a linearly growing drift also changes which extrema matter, and the sign of a jump decides whether the large-drift conclusion can hold. The negative-drift conclusion asks for compact-uniform convergence, which is stronger than merely retaining $J_1$ convergence. These are real constraints on the statement; dropping a jump condition or replacing $U$ by $J_1$ changes the theorem. Whitt (1980), pp. 78–81

Formalization scope

Paths are total functions $\mathbb R\to S$, but only their values on the stated interval count. The Lean definitions require càdlàg membership of both prelimit and limit paths in each $J_1$ convergence claim. A time change on $[a,b]$ is strictly increasing, continuous, and onto; this represents the paper’s increasing homeomorphism. Uniform and $J_1$ distances take values in the extended nonnegative reals, so an unbounded supremum cannot silently become zero. Continuity points use continuity within $T$, and the discontinuity set likewise contains only points of $T$. Sequential convergence is used because the paper’s $J_1$ spaces are metrizable. These choices encode the paper’s path-space conventions rather than imposing conditions on values outside $T$. Whitt (1980), §§2, 4, 6

The compact form of Theorem 6.1 uses metric (2.1); the noncompact metric (2.2) is outside this mission’s definition layer. Theorem 4.1 retains the paper’s complete separable metric-space setting and metric compatibility of addition. Its measurability clause is omitted: this mission does not construct the Borel measurable space of $D(T,S)$, and a product measurable structure on all total functions would express a different claim. The source’s $f(x)=x-x^{\downarrow}$ is used on both prelimit and limit paths; it is not replaced by a clipped running infimum. Whitt (1980), pp. 70, 78–81

The definition layer, the addition theorem, the running-extremum results, and the barrier continuity statement are welcome proof targets. The theorem’s centered-path hypothesis remains a genuine $J_1$ convergence assertion, and the jump conditions apply to the limit path only. A vacuous convergence predicate or an independently assumed conclusion would not represent this mission’s goal.

Selected references

• Ward Whitt, Some Useful Functions for Functional Limit Theorems, Mathematics of Operations Research 5(1), 1980, pp. 67–85. DOI: 10.1287/moor.5.1.67

Why time reversal matters

Functional limit theorems describe the behavior of entire time-indexed paths under scaling. In queueing and other stochastic models, a path of departures, workload, or accumulated input may be related to another path by reversing a finite observation window. To carry a limit theorem through that change of viewpoint, one needs to know how reversal acts on the topology of the path space. Ward Whitt's 1980 paper studies composition, reflection, first passage, and time reversal as maps that transfer functional limit theorems to derived processes (Whitt 1980).

The time-reversal result is a precise metric statement. Ordinary reversal preserves the Skorohod $J_1$ distance exactly. A version shifted to begin at the group identity satisfies a bound with constant two. These facts give continuity and, when the random paths satisfy the hypotheses, make reversal available to the continuous mapping theorem. Whitt's introduction explains how deterministic path continuity connects to stochastic-process limits (Whitt 1980, pp. 67–69).

Paths, distances, and reversal

Let $S$ be a complete separable metric group, written additively, with identity $0$ and metric $m$. The metric is translation invariant: moving both points by the same group element leaves their distance unchanged. Let $D$ contain the càdlàg paths $x:[0,1]\to S$, meaning paths that are right continuous and possess left limits. Section 8 additionally requires $x(1)=x(1-)$, so the terminal value agrees with the value approached from earlier times. Let $D_\theta=\{x\in D:x(0)=0\}$. These are the spaces Whitt fixes for time reversal (Whitt 1980, p. 83).

The uniform distance is $\rho(x,y)=\sup_{0\le t\le1}m(x(t),y(t))$. A time change $\lambda\in\Lambda$ is an increasing homeomorphism of $[0,1]$ onto itself; $e(t)=t$ is the identity time change. The $J_1$ distance used in the paper is

$$d(x,y)=\inf_{\lambda\in\Lambda} \max\{\rho(\lambda,e),\rho(x,y\circ\lambda)\}.$$

Two paths may be close when their values can be matched after a small distortion of time. This is the paper's equation (2.1), rather than an unspecified equivalent $J_1$ metric (Whitt 1980, p. 70).

The reverse-time map reads the left limit at the reflected time:

$$R(x)(t)= \begin{cases} x((1-t)-),&0\le t<1,\\ x(0),&t=1, \end{cases} \qquad r(x)(t)=R(x)(t)-x(1).$$

The left limit is essential: the direct value $x(1-t)$ places jumps on the wrong side and generally does not give a right-continuous reverse path. The terminal convention gives $R(x)(0)=x(1)$, so $r(x)(0)=0$. Whitt also uses the reflected time change $(-r)(\lambda)(t)=1-\lambda(1-t)$ in the proof (Whitt 1980, pp. 83–84).

Formalization targets

The goal is Whitt's Theorem 8.1. For every $x,y\in D$,

$$d(Rx,Ry)=d(x,y),\qquad d(rx,ry)\le 2d(x,y).$$

The second inequality also applies when both paths lie in $D_\theta$, since $D_\theta\subseteq D$. The constant is the paper's printed constant, not a parameter to optimize in this mission (Whitt 1980, Theorem 8.1, p. 83).

The milestones record four claims from the theorem's proof: reversal preserves its path spaces and has inverse behavior on its natural domains; reflected time changes stay in $\Lambda$; reversal obeys exact and bounded uniform-distance relations; and reversal interacts with time changes through $(-r)(\lambda)$. Their descriptions reproduce the printed proof sentences. The proof calls the time-change cost a “minimum” once, although equation (2.1) and the displayed calculation use a maximum. The formal target follows equation (2.1) (Whitt 1980, pp. 70, 84).

What the result gives

An isometry preserves every $J_1$ distance, so ordinary reversal preserves convergent and Cauchy sequences in this metric. The bound for $r$ makes shifted reversal Lipschitz and therefore continuous. A functional limit theorem for input paths can then be translated into one for the corresponding reversed paths whenever the path-space hypotheses hold. The result is already proved in Whitt's paper; the open work here is a machine-checked development of its definitions, supporting claims, and metric theorem.

Formalizing the result also creates reusable interfaces for endpoint-constrained càdlàg paths, left-limit reversal, uniform path distance, and reflected time changes. Those interfaces may support later statements about reversed stochastic models. This mission does not assert weak convergence of probability laws or Borel measurability; those would require a separate measurable-space development for $D$.

Mathematical difficulty

The immediate idea $t\mapsto x(1-t)$ does not preserve càdlàg paths at jumps: a right-continuous jump becomes left-continuous. The correction uses left limits, but endpoint behavior then matters. Without $x(1)=x(1-)$, the shifted reverse path need not start at the identity. The $J_1$ assertion also concerns an infimum over all allowed time changes, so a uniform-distance statement alone does not establish the theorem. These are the points at which reflection of the time axis becomes a claim about the precise path-space topology (Whitt 1980, pp. 83–84).

Formalization scope

Lean represents paths as total functions $\mathbb R\to S$; all path conditions, metrics, and equality claims read only $[0,1]$. Càdlàg membership requires right continuity and left limits within the interval, with $x(1)=x(1-)$ included in $D$. Completeness and separability are explicit because the paper calls $S$ a complete separable metric space. The additive metric group is represented by an additive commutative group with an explicit translation-invariance hypothesis. This chooses the commutative reading exemplified by $S=\mathbb R^k$; the paper uses additive notation but does not separately state commutativity.

The paper's $\Lambda$ is encoded by strict increase, continuity, and surjectivity on $[0,1]$, which characterize increasing homeomorphisms there. The paper's $d$ is equation (2.1) on both sides of the goal. Lean uses extended nonnegative distances for $\rho$ and $d$, avoiding the default value a real supremum can assign to an empty or unbounded set. The claims do not become trivial by dropping $D$ membership: it is an explicit hypothesis, and the zero path over $\mathbb R$ satisfies it. A local, sorry-free check confirms that instance. Two easier statements are excluded on purpose: the goal is stated for the $J_1$ distance $d$, not the uniform distance $\rho$, and $R$ reads the left limit $x((1-t)-)$, not the value $x(1-t)$.

The proof's first sentence says the maps are one-to-one and onto. The map $r:D\to D_\theta$ forgets an additive offset and is not injective on all of $D$; the formal milestone expresses involution for $R$ on $D$ and for the restriction of $r$ to $D_\theta$. One identity printed for $r$ with $x\in D_\theta$ also holds for every $x\in D$; the formal milestone uses that broader domain and says so in its item description. Contributions are welcome for the left-limit and time-change infrastructure as well as the four listed statements and the goal.

Selected references

• Ward Whitt, Some Useful Functions for Functional Limit Theorems, Mathematics of Operations Research 5(1), 1980, pp. 67–85. DOI: 10.1287/moor.5.1.67.

Motivation

A service system with several servers can use a larger fraction of its capacity than a single-server system while still giving customers an acceptable wait. The design question is how much unused capacity is needed as the number of servers grows. Ward Whitt's study begins with a factory manager deciding whether to put four machines in one work area. It uses queueing models to connect the number of machines, their utilization, and the delays customers experience. The proposed utilization rule, $(1-\rho)\sqrt{s}=\gamma$, keeps the unused fraction $1-\rho$ proportional to $1/\sqrt{s}$ as the number of servers $s$ changes (Whitt 1992, pp. 708–710).

The paper distinguishes the chance that a customer waits at all from the length of the wait once a delay occurs. That distinction matters when a staffing choice must control a long wait rather than merely the fraction of customers who wait. The exact $M/M/s$ result in §3.1 gives a benchmark for the paper's later approximations to more general arrival and service processes (Whitt 1992, pp. 719–720).

Setting

An $M/M/s$ queue has Poisson arrivals at rate $\lambda$, $s$ identical servers, independent exponential service times, unlimited waiting room, and first-come first-served service. Time is measured so that each server's service rate is one. The server utilization is $\rho=\lambda/s$. We consider $s\geq1$ and $0<\lambda<s$, so the queue has a stationary distribution. Let $p_n$ be the stationary probability that $n$ customers are in the system. Its state process is a birth–death process: arrivals occur at rate $\lambda$ and departures from state $n$ occur at rate $\min(n,s)$.

Let $W$ be an arriving customer's steady-state time in the queue before service begins. An arrival finding fewer than $s$ customers starts service at once. An arrival finding $s+j$ customers waits for $j+1$ service completions while all servers are occupied. In the formal model, the law of $W$ has an atom at zero of mass $\sum_{n<s}p_n$ and, for each $j\geq0$, a gamma component with shape $j+1$, rate $s$, and weight $p_{s+j}$. The arrival-state weights use the Poisson-arrivals-see-time-averages property. This construction represents the FCFS waiting time independently of any claim about its eventual conditional distribution.

For an event $W>0$ of positive probability, $\mathcal L(W\mid W>0)$ denotes the conditional waiting-time law. The paper calls a customer delayed exactly when $W>0$. The scale parameter $\gamma$ is defined by the utilization equation $(1-\rho)\sqrt{s}=\gamma$; because $\rho<1$, this equation makes $\gamma$ positive.

Formalization targets

The supporting target is the §3.1 statement that a delayed customer's waiting time is exponential with rate $s(1-\rho)$, equivalently with mean $1/[s(1-\rho)]$:

$$\mathcal L(W\mid W>0)=\operatorname{Exp}\bigl(s(1-\rho)\bigr).$$

The goal is Proposition 3.1 under the utilization equation. It records both the entire conditional tail and its mean, for every $x\geq0$:

$$\mathbb P(W>x\mid W>0)=e^{-\gamma\sqrt{s}\,x}, \qquad \mathbb E[W\mid W>0]=\frac{1}{\gamma\sqrt{s}}.$$

The goal also asserts $\mathbb P(W>0)>0$, so these conditional quantities are defined. The scan of Proposition 3.1 prints the exponent without $x$ and the mean denominator with $\sqrt2$; the statement here uses the values determined by the preceding exponential-law sentence and by the paragraph following the proposition (Whitt 1992, pp. 719–720).

Significance

The result gives an exact delay distribution for the Markovian multi-server system. With $s$ servers and the utilization equation held fixed, the conditional tail at any positive threshold and the mean wait of a delayed customer both depend on $s$ through $\gamma\sqrt{s}$. Thus the same rule that organizes the probability-of-delay discussion has a definite implication for customers who actually wait. It is also the exact comparison point for the paper's later, explicitly approximate analysis of general $G/G/s$ queues (Whitt 1992, §3.2).

Formalizing this result supplies a reusable representation of the stationary FCFS waiting-time law as a mixture of Erlang distributions, together with a precise connection between birth–death state probabilities and delay. The related platform items QueueingFundamentals.BirthDeath.mmc_steady_state and QueueingFundamentals.BirthDeath.erlang_c_formula state the stationary state probabilities and the probability of delay, respectively; both are open and neither states this waiting-time law. The open QueueingFundamentals.BirthDeath.halfin_whitt concerns the heavy-traffic limit for the $M/M/c$ case, a different target. The current result is known mathematically from Whitt's 1992 paper; this mission asks for its machine-checked proof in the specified model.

Difficulty

The conditional distribution is not immediate from an individual service time. A delayed arrival can encounter any queue length $s+j$, so its waiting time is a gamma law whose shape varies with $j$. The stationary probabilities provide the weights of infinitely many such components. Showing that this whole mixture has one exponential law requires a relation among the stationary tail probabilities and control of the infinite sum. The same mixture must have a well-defined first moment for the conditional mean formula.

Formalization scope

Lean uses a natural number $s\geq1$, a real arrival rate $0<\lambda<s$, and a real sequence $p_n$ satisfying the published QueueingFundamentals.BirthDeath.IsSteadyState predicate with constant birth rate $\lambda$ and death rate $\min(n,s)$. That predicate includes nonnegativity, total mass one, and the global balance equations. Its published definition is imported as a reference. Service rate is exactly one, as fixed in §2.1 of the paper. The conditional waiting-time law is a measure on real time, with a point mass at zero and a countable gamma mixture. The positive-time restriction of this measure expresses conditioning in the milestone; real measure values and an integral express the tail and mean in the goal.

The formal hypotheses spell out positive arrival rate, at least one server, and subcritical utilization. These make the stationary law and the positive exponential rate meaningful. The goal states positive delay probability before dividing by it. The waiting-time law is built from stationary queue lengths and service completions, rather than being defined as exponential; proving that it becomes exponential is the mathematical task. Solvers may develop reusable measure-mixture identities, stationary-tail facts, and gamma or exponential distribution results. The approximation formulas of §3.2 and the paper's numerical tables are outside this mission.

Selected references

• Ward Whitt, Understanding the Efficiency of Multi-Server Service Systems, Management Science 38(5):708–723, 1992. DOI: 10.1287/mnsc.38.5.708.