The OR Formalization Drive

Why many-server queues in heavy traffic

Telephone exchanges, call centers, hospital wards and cloud server pools are queues with a large number $s$ of parallel servers. For such systems two classical approximations pull in opposite directions. Holding $s$ fixed and letting the traffic intensity $\rho \to 1$ gives the conventional heavy-traffic limit, in which almost every customer waits. Letting $s$ grow with $\rho$ fixed below one gives a limit in which almost nobody waits. Real systems sit in between: a positive fraction of customers wait, and waits are short. Halfin and Whitt (Operations Research 29 (1981) 567–588) identified the scaling that produces this middle regime and the limit laws it yields. The scaling is now called the Halfin–Whitt or quality-and-efficiency-driven (QED) regime, and it underlies square-root staffing: run $n \approx R + \beta\sqrt R$ servers for an offered load of $R$ Erlangs. Borst, Mandelbaum and Reiman (Operations Research 52 (2004) 17–34) built staffing rules for large call centers on it.

Timeline. Erlang (1917) gave the delay probability of the M/M/s queue. Iglehart (1965) proved diffusion limits for M/M/s queues with the number of servers growing and the traffic intensity fixed. Halfin and Whitt (1981) proved that the delay probability has a limit strictly between 0 and 1 exactly when $(1-\rho_n)\sqrt n$ converges to a positive constant (their Proposition 1), derived the limit of the stationary queue length (Theorem 1), and extended the process-level limit to GI/M/s queues. Later work extended the regime to many-server queues with general service times and abandonment.

Setting

An M/M/s queue has Poisson arrivals at rate $\lambda > 0$, $s \ge 1$ servers, and exponential service times with rate $\mu > 0$. Its traffic intensity is $\rho = \lambda/(s\mu)$. The number $Q(t)$ of customers in the system (waiting or in service) is a birth–death process on $\{0,1,2,\dots\}$ with birth rate $\lambda$ and death rate $\min(k,s)\mu$ in state $k$. When $\rho < 1$, $Q(t)$ converges in distribution to $Q(\infty)$, whose law $p_k = P(Q(\infty) = k)$ is the unique probability distribution solving the balance equations of this chain. The probability of delay $\alpha = P(Q(\infty) \ge s)$ is the Erlang-C formula.

The mission studies a sequence of such queues: queue $n \ge 1$ has $s_n = n$ servers, the fixed service rate $\mu$, and arrival rate $\lambda_n$ with $\rho_n = \lambda_n/(n\mu) < 1$ and $\lambda_n \to \infty$. Its stationary queue length is $Q_n(\infty)$, and the scaled queue length is

$$X_n = \frac{Q_n(\infty) - n}{\sqrt n}.$$

Throughout, $\Phi$ and $\varphi$ denote the standard normal distribution function and density, and for $\beta > 0$

$$\alpha(\beta) = \big[1 + \sqrt{2\pi}\,\beta\,\Phi(\beta)\,e^{\beta^2/2}\big]^{-1} = \frac{\varphi(\beta)}{\varphi(\beta) + \beta\Phi(\beta)} \in (0,1).$$

The standing assumption of the paper's Section 2 is the heavy-traffic condition

$$\lim_{n\to\infty} (1-\rho_n)\sqrt n = \beta, \qquad \beta > 0. \tag{2.2}$$

Formalization targets

Goal: Theorem 1

Under (2.2), with $\alpha = \alpha(\beta)$,

$$X_n \Rightarrow X, \qquad P(X \ge 0) = \alpha,\quad P(X > x \mid X \ge 0) = e^{-x\beta}\ (x \ge 0),\quad P(X \le x \mid X \le 0) = \frac{\Phi(\beta + x)}{\Phi(\beta)}\ (x \le 0).$$

The limit is exponential with rate $\beta$ above zero and a normal law truncated at zero below it, with masses $\alpha$ and $1-\alpha$.

Milestones, in the order the paper proves them

1. The stationary law (1.1)–(1.3) of one M/M/s queue, and the identification (1.2) of $P(Q(\infty) \ge s)$ with the Erlang-C formula.
2. Lemma 1: recursions expressing the partial moments $\sum_{k<s} k^m p_k$ and $\sum_{k \ge s} k^m p_k$ through lower ones and $\alpha$.
3. Proposition 1: $P(Q_n(\infty) \ge n) \to \alpha \in (0,1)$ if and only if (2.2) holds, and then $\alpha = \alpha(\beta)$.
4. Proposition 2: for $\delta > 0$ and $(n - \delta_n)/\sqrt n \to \delta$ with $\delta_n \le n$,

$$P(Q_n(\infty) \le \delta_n \mid Q_n(\infty) \le n) \to \frac{\Phi(\beta-\delta)}{\Phi(\beta)},\qquad \sqrt n\,P(Q_n(\infty) = [\delta_n] \mid Q_n(\infty) \le n) \to \frac{\varphi(\beta-\delta)}{\Phi(\beta)},$$

and for $(\delta_n - n)/\sqrt n \to \delta$ with $\delta_n \ge n$,

$$P(Q_n(\infty) \ge \delta_n \mid Q_n(\infty) \ge n) \to e^{-\delta\beta},\qquad \sqrt n\,P(Q_n(\infty) = [\delta_n] \mid Q_n(\infty) \ge n) \to \beta e^{-\delta\beta}.$$

5. Corollary 1: the first four moments and the variance of $X_n$ converge to explicit functions of $\alpha$ and $\beta$, for example $E X_n \to -\beta + \alpha/\beta$.

Significance

Theorem 1 is the stationary half of the Halfin–Whitt regime. It turns the Erlang-C formula, a ratio of sums with $n$ terms that is opaque for large $n$, into a two-parameter description of congestion: $\alpha$ is the fraction of customers who wait, $\beta$ is the scaled safety margin of servers, and the conditional laws above give the distribution of the number waiting and of the number of idle servers. Proposition 1 and Theorem 1 are what square-root staffing rules evaluate; Corollary 1 supplies the mean and variance used in the paper's numerical approximations.

All results of this mission are proved in the paper, by direct calculation from (1.1)–(1.3), Stirling's formula and the central limit theorem for Poisson variables. None has a machine-checked proof. Proposition 1 and the two Section 1 facts are posed on the platform as open theorems from the Gross et al. queueing textbook series; this mission adds the conditional and local limits, the weak-convergence theorem and the moment results on the same definitions, so that a closed development would give a fully checked derivation of the Halfin–Whitt limit from the balance equations.

Difficulty

The stationary law is explicit, so nothing here needs a stochastic process. The difficulty is analytic and uniform: each statement is a limit of ratios of sums of $n$ or infinitely many terms, $(n\rho_n)^k/k!$ below $n$ and a geometric tail above $n$, in which $\rho_n \to 1$ and $n \to \infty$ at linked rates. The lower half needs a central limit theorem for Poisson laws whose parameter $n\rho_n$ moves with $n$ and is evaluated at a moving point; the local limits (2.10) and (2.12) need Stirling's formula with explicit control of $\log \rho_n$ to second order. Weak convergence then needs tightness, or a direct argument from the conditional limits, and the moment limits need uniform integrability, which the fixed-$s$ moment recursions do not give by themselves. Holding $\rho$ fixed, or taking limits in $n$ and $\rho$ one after the other, gives degenerate answers ($\alpha = 0$ or $\alpha = 1$); the two limits must be taken together.

Formalization scope

The Lean development lives in the namespace HalfinWhitt81.Stationary. The law of $Q_n(\infty)$ is a sequence p n : ℕ → ℝ satisfying IsSteadyState (fun _ => lam n) (mmcDeath μ n) (p n) of the referenced definition QueueingFundamentals.BirthDeath.Balance (nonnegative, summing to one, solving the balance equations of the M/M/n chain); the Erlang-C formula and $\alpha(\beta)$ are erlangC and halfinWhittAlpha of QueueingFundamentals.BirthDeath.Erlang, and $\Phi$, $\varphi$ are Mathlib's cdf (gaussianReal 0 1) and gaussianPDFReal 0 1. These, and the three open theorems for (1.1)–(1.3), (1.2) and Proposition 1, come from the Gross et al. (2008) series and are referenced, not restated. The related Borst–Mandelbaum–Reiman items (DimCallCenters.*.lemma_4_1, for the continuous extension of Erlang-C) state a different result and are not used.

Conventions committed to:

• hypotheses on queue $n$ (rates, steady state, $\delta_n \le n$ or $\delta_n \ge n$) are imposed for $n \ge 1$, and only limits are asserted; $\lambda_n \to \infty$ is kept as a hypothesis although it follows from (2.2);
• the standing assumption "(2.1) or, equivalently, (2.2)" is the hypothesis (2.2) with $\beta > 0$, and $\alpha$ is written as $\alpha(\beta)$, equal to the limit in (2.1) by Proposition 1;
• probabilities of $Q_n(\infty)$ are sums of p n (probLE, probGE, probEqInt); a conditional probability $P(A \mid B)$ with $A \subseteq B$ is the ratio $P(A)/P(B)$; $[x]$ is Int.floor, with $P(Q = m) = 0$ for $m < 0$;
• $X_n$ uses $n^{-1/2}$: the page's (2.13) prints $n^{1/2}$, a misprint; (2.12) prints the limit $\beta e^{-\beta}$, a misprint for $\beta e^{-\delta\beta}$, which is what is formalized;
• "$X_n \Rightarrow X$" is stated as: there is a Borel probability measure $\nu$ on $\mathbb R$ with the three properties of Theorem 1 such that $\sum_k p_n(k)\,g((k-n)/\sqrt n) \to \int g\,d\nu$ for every bounded continuous $g$; the exponential clause is for $x \ge 0$ and the normal clause for $x \le 0$;
• infinite sums whose convergence is not implied by the hypotheses (the moment series of Lemma 1 and Corollary 1) have their convergence in the conclusion.

The goal cannot be met trivially: the law of $Q_n(\infty)$ is the M/M/n steady state, not a free sequence; the limit measure is asserted to exist together with its three properties and the convergence, so no vacuous "for every $\nu$" reading is possible; and the statement $P(X_n \ge 0) \to \alpha$ alone is Proposition 1, already posed.

A complete development needs: the closed form of the M/M/n steady state, Poisson central limit and local limit theorems with a moving parameter, Stirling's formula with error terms, a criterion for weak convergence from convergence of distribution functions, and uniform integrability for the moments. The Poisson limit theorems and the weak-convergence criterion are reusable beyond this mission. Proofs of the referenced open theorems, of individual milestones, and alternative routes to Theorem 1 (moment convergence, or the process-level limit) are all welcome.

Selected references

• S. Halfin and W. Whitt, Heavy-Traffic Limits for Queues with Many Exponential Servers, Operations Research 29(3) (1981) 567–588. https://doi.org/10.1287/opre.29.3.567
• R. B. Cooper, Introduction to Queueing Theory, Macmillan, 1972 (the source cited for (1.1)–(1.3)).
• D. L. Iglehart, Limiting Diffusion Approximations for the Many Server Queue and the Repairman Problem, Journal of Applied Probability 2(2) (1965) 429–441.
• S. Borst, A. Mandelbaum and M. I. Reiman, Dimensioning Large Call Centers, Operations Research 52(1) (2004) 17–34. https://doi.org/10.1287/opre.1030.0081
• D. Gross, J. F. Shortle, J. M. Thompson and C. M. Harris, Fundamentals of Queueing Theory, 4th ed., Wiley, 2008 (§2.4: the Erlang-C formula and the Halfin–Whitt limit).

Motivation

The contact process is the basic model of an infection spreading on a lattice: each infected site recovers at a constant rate, and each healthy site becomes infected at a rate that depends on how many of its neighbours are infected. It was introduced by T. E. Harris in Contact Interactions on a Lattice (Ann. Probab. 2, 1974), and it became one of the central examples of interacting particle systems, alongside the voter model and the stochastic Ising model (Liggett, Interacting Particle Systems, 1985, Ch. VI). Its main question is whether a finite infection dies out with probability one or survives forever with positive probability, and how that answer depends on the rates.

Harris's paper gives the first answers. It shows that, under natural conditions on the rates, the probability of survival is monotone, subadditive and submodular as a function of the initial set, and it uses these properties to prove a condition for certain extinction: Theorem 7.1. For the linear contact process (recovery rate $1$, infection rate $\lambda$ per infected neighbour) the condition is $(2d-1)\lambda<1$. Holley obtained Theorem 7.1 by different methods at about the same time (footnote 4, p. 981); Holley and Liggett (Ann. Probab. 3, 1975) later developed the duality theory of such processes.

Setting

Let $Z_d$ be the $d$-dimensional integer lattice, $d\ge1$. Sites $x,y$ are neighbours, $x\sim y$, if $\sum_i|x^i-y^i|=1$; $N_x$ is the set of the $2d$ neighbours of $x$. A configuration is a finite set $\xi\subset Z_d$ of infected sites; $\Xi_0$ is the set of all of them, $\xi(N_x)=|\xi\cap N_x|$, and $x\sim\xi$ means $x\notin\xi$ and $x$ has a neighbour in $\xi$.

The parameters are a recovery rate $\mu\ge0$ and infection rates $\lambda_0=0,\lambda_1,\dots,\lambda_{2d}\ge0$. The contact process $\{\xi_t\}$ started from a finite set is the continuous-time Markov chain on $\Xi_0$ that leaves $\xi$ at rate

$$q_\xi=\mu|\xi|+\sum_{x\sim\xi}\lambda_{\xi(N_x)}\qquad(4.4)$$

and jumps to $\xi\setminus x$ ($x\in\xi$) at rate $\mu$ and to $\xi\cup x$ ($x\sim\xi$) at rate $\lambda_{\xi(N_x)}$ (4.5). The empty set is absorbing. Write $P_t(\xi,\eta)$ for its transition function, $r(\xi,\eta)=q_{\xi\eta}/q_\xi$ for the transition matrix of its imbedded jump chain, and

$$p_t(\xi)=P_\xi\{\xi_t\neq\varnothing\},\qquad p_\infty(\xi)=\lim_{t\to\infty}p_t(\xi),\qquad m_t(\xi)=\mathscr E_\xi|\xi_t| .$$

The process is increasing, subadditive or submodular if, for all $\xi,\eta$ and $t\ge0$, respectively $p_t(\xi)\le p_t(\xi\cup\eta)$, $p_t(\xi\cup\eta)\le p_t(\xi)+p_t(\eta)$, or $p_t(\xi\cup\eta)+p_t(\xi\cap\eta)\le p_t(\xi)+p_t(\eta)$ ((5.2)–(5.4)).

Formalization targets

Goal: Theorem 7.1 (p. 981)

If

$$(2d-1)\lambda_k<k\mu,\qquad k=1,\dots,2d,$$

then

$$p_\infty(\xi)=0\qquad\text{for every finite }\xi\subset Z_d .$$

Milestones, in the order of the paper

1. Theorem 4.3 (p. 975): from a finite set the process has finitely many jumps on every $[0,t]$; equivalently $\sum_\eta P_t(\xi,\eta)=1$.
2. (5.1) (p. 976): $p_t(\xi)$ is non-increasing in $t$, so $p_\infty(\xi)$ is a limit.
3. Theorem 5.6(a) (p. 977): if $\lambda_k\uparrow$, the process is increasing.
4. Theorem 5.6(b) (p. 977): if moreover $\lambda_k/k\downarrow$, it is subadditive.
5. Lemma 5.8 (p. 978): if $\lambda_k\le\min_{j\ge k}\lambda'_j$, then $p_t\le p'_t$ and $m_t\le m'_t$.
6. Theorem 6.2 (p. 979): if $\{\lambda_i\}$ is concave and non-decreasing with $\lambda_0=0$, the process is submodular.
7. Change of time scale (p. 981): multiplying all rates by $c>0$ replaces $P_t$ by $P_{ct}$ and leaves $p_\infty$ unchanged.
8. First-step display (p. 981): for a subadditive process, a singleton $\{x\}$, and $n\ge1$, $p_\infty(\{x\})=\sum_\eta r^{(n)}(\{x\},\eta)p_\infty(\eta)\le p_\infty(\{x\})\,\mathscr E_x|\xi_{(n)}|$.
9. (7.2) (p. 981): for $\mu=1$, $\lambda_k=k\lambda$: $\pi_1=\frac{2d\lambda}{1+2d\lambda}\pi_2$, where $\pi_1,\pi_2$ are $p_\infty$ of a site and of a pair of neighbours.
10. (7.3)–(7.4) (p. 981): $\pi_2-\pi_1=(1+(2d-1)\lambda)\sum'r(\xi,\xi\cup x)(p_\infty(\xi\cup x)-\pi_2)\le(\pi_2-\pi_1)(2d-1)\lambda$.
11. Conclusion of the proof (p. 981): for $\mu=1$, $\lambda_k=k\lambda$, $\pi_1>0$ implies $(2d-1)\lambda\ge1$.

Significance

Theorem 7.1 gives a computable sufficient condition for extinction that holds in every dimension and for every non-linear infection law satisfying it. For the linear process it gives the lower bound $\lambda_c\ge1/(2d-1)$ on the critical infection rate; the paper's §9 shows conversely that the process survives for large $\lambda_k/\mu$, so the critical rate is non-trivial. The monotonicity, subadditivity and submodularity theorems of §§5–6 are the comparison tools on which much of the later theory of attractive spin systems was built, and Lemma 5.8 is the prototype of comparison between processes with ordered rates.

The results are proved in the paper, and the proof of Theorem 7.1 is half a page once §§5–6 are available. No machine-checked version of the contact process, of its comparison theorems, or of Theorem 7.1 is known to exist. A formalization produces a continuous-time Markov chain with unbounded rates and an absorbing state on a countable state space, its minimal transition function, a non-explosion criterion, first-step analysis for absorption probabilities, and coupling arguments for monotone comparison. The companion mission Contact Interactions on a Lattice II treats Theorem 7.6, the integrability of the expected size, on the same model.

Difficulty

The obvious argument is first-step analysis on singletons and pairs, (7.2) and (7.3). It fails as it stands because (7.3) involves $p_\infty$ of every three-point set $\xi\cup x$, and no finite system of first-step equations closes; some inequality between survival probabilities of nested sets is needed, and the milestones record the one the paper uses. Submodularity, subadditivity and monotonicity are not consequences of the rates by direct computation: the paper proves each by building several contact processes on one probability space (Tables 1–4), which in a formal development means constructing couplings of chains with unbounded rates and checking that each marginal has the right law. The reduction from general $\lambda_k$ to the linear case also needs the comparison Lemma 5.8, another coupling. Finally the rates grow linearly in $|\xi|$, so uniformization does not apply and non-explosion (Theorem 4.3) must be proved before transition probabilities are known to sum to one.

Formalization scope

The process is the countable-state chain of §4, not the Feller process on all subsets of $Z_d$ of §2: sites are Fin d → ℤ, configurations are Finsets, and $P_t$ is the minimal transition function of the rates (4.4)–(4.5), defined as the increasing limit of Feller's backward recursion. There is no measure-theoretic process object; $P_\xi$ and $\mathscr E_\xi$ enter only through $P_t$. Consequently Theorem 4.1 ($\Xi_0$ is stochastically closed) is built in, and every statement is made for finite initial sets, including those the page makes for arbitrary $\xi\in\Xi$ ((5.2)–(5.4), Lemma 5.8); for infinite $\xi$ the page itself notes $p_t(\xi)=1$.

The committed conventions:

• $P_t$, $p_t$, $p_\infty$ and $m_t$ take values in $[0,\infty]$; $p_t(\xi)=1-P_t(\xi,\varnothing)$ and $p_\infty(\xi)=\inf_{t\ge0}p_t(\xi)$, so $p_\infty=0$ means absorption at $\varnothing$ and cannot be produced by loss of mass to explosion.
• The rates are real; every theorem assumes $d\ge1$, $\mu\ge0$, $\lambda_0=0$ and $\lambda_k\ge0$ for $k\le2d$. Conditions on $\{\lambda_k\}$ range over $k=1,\dots,2d$ (or $0,\dots,2d$), never over all natural numbers.
• $2d-1$ is the real number $2d-1$.
• In the linear case the scalar rate is written $l$ ($\lambda$ is reserved in Lean), $\mu=1$ and $\lambda_k=kl$.

Trivializing formalizations are ruled out: the survival probability is not taken from an arbitrary semigroup or a uniformization (the rates are unbounded); $p_\infty$ is not a real infimum or limUnder, whose junk value $0$ would make the goal free; $\lambda$ is indexed by the number of infected, not vacant, neighbours, and $\lambda_0=0$ forbids spontaneous infection; and (7.2) is stated for every site and every neighbour, not for one fixed pair.

A complete development needs: non-explosion for chains with linearly growing rates; the Chapman–Kolmogorov and strong Markov properties of the minimal chain at jump times; lattice symmetries (translations and coordinate reflections) acting on the chain; and couplings of several contact processes realised as one chain on a product type space. The chain infrastructure and the coupling lemmas are reusable beyond this mission, in particular by its companion. Contributions to any milestone are welcome, including proofs of the milestones under the paper's hypotheses by methods other than Harris's.

Selected references

• T. E. Harris, Contact interactions on a lattice, The Annals of Probability 2(6), 969–988, 1974. https://doi.org/10.1214/aop/1176996493
• R. Holley and T. M. Liggett, Ergodic theorems for weakly interacting infinite systems and the voter model, The Annals of Probability 3(4), 643–663, 1975. https://doi.org/10.1214/aop/1176996306
• T. M. Liggett, Interacting Particle Systems, Grundlehren 276, Springer, 1985. https://doi.org/10.1007/978-1-4613-8542-4
• W. Feller, On the integro-differential equations of purely discontinuous Markoff processes, Transactions of the AMS 48, 488–515, 1940. https://doi.org/10.1090/S0002-9947-1940-0002697-3
• J. R. Norris, Markov Chains, Cambridge University Press, 1997, §2.8. https://doi.org/10.1017/CBO9780511810633

Motivation

The newsvendor problem is the basic model of stocking under uncertain demand: an order quantity is fixed before a random demand is observed, and every unit short or left over is penalized. Its solution is a quantile of the demand distribution. In practice that distribution is unknown and only a sample of past demands is available. The standard remedy, sample average approximation (SAA), replaces the expectation by the empirical average over the sample and orders the corresponding sample quantile. The practical question is how many observations make the SAA order nearly as good as the true optimum, without assuming anything about the demand distribution.

Levi, Roundy and Shmoys (Math. Oper. Res. 32(4), 2007) gave the first distribution-free answer, using Hoeffding's inequality. Levi, Perakis and Uichanco (Oper. Res. 63(6), 2015) improved it with Bernstein's inequality. The improvement matters most when the service level is high, which is the typical case in inventory practice. This mission formalizes that improved bound, Theorem 2 of the 2015 paper.

Setting

A retailer orders $q\in\mathbb R$ units before a real demand $D$ with law $\mu$ and cdf $F(q)=\Pr(D\le q)$ is realized. Unmet demand costs $b>0$ per unit (underage cost) and unsold stock costs $h>0$ per unit (overage cost). The expected cost is

$$C(q)=\mathbb E\big[b(D-q)^+ + h(q-D)^+\big],$$

which is finite when $\mathbb E|D|<\infty$. The critical quantile is $q^*=\inf\{q:F(q)\ge b/(b+h)\}$, and it minimizes $C$.

An order $q$ is $\epsilon$-optimal when its relative regret $(C(q)-C(q^*))/C(q^*)$ is at most $\epsilon$, i.e. $C(q)\le(1+\epsilon)C(q^*)$. The set of such orders is $S_\epsilon$. The function $C$ is convex with one-sided derivatives

$$\partial_+C(q)=-b+(b+h)F(q),\qquad \partial_-C(q)=-b+(b+h)\Pr(D<q).$$

The LRS interval is

$$S^{LRS}_\epsilon=\Big\{q:\partial_-C(q)\le\tfrac\epsilon3\min(b,h)\ \text{and}\ \partial_+C(q)\ge-\tfrac\epsilon3\min(b,h)\Big\}.$$

Given an i.i.d. sample $D^1,\dots,D^N$ from $\mu$, the empirical cdf is $\hat F_N(q)=\frac1N\sum_k\mathbf 1[D^k\le q]$, and the SAA solution is the sample quantile

$$\hat Q_N=\inf\{q:\hat F_N(q)\ge b/(b+h)\}.$$

It is a random variable through the sample.

Formalization targets

Goal: Theorem 2 (improved LRS bound)

For every demand law with $\mathbb E|D|<\infty$, every $N\ge1$ and every $0<\epsilon\le1$,

$$\Pr\big(\hat Q_N\notin S_\epsilon\big)\le 2\exp\Big(-\frac{N\epsilon^2}{18+8\epsilon}\cdot\frac{\min\{b,h\}}{b+h}\Big).$$

Milestones, in the order of the paper's proof

1. (§2, p. 7) $C$ is convex, and its one-sided derivatives are the two formulas above.
2. (Theorem EC.1) Bernstein's inequality for i.i.d. bounded variables: $\Pr\big(\frac1N\sum_iX^i-\mathbb E X^1\ge t\big)\le\exp\big(-Nt^2/(2\sigma^2+2tc/3)\big)$.
3. (Proposition EC.1) For every $\gamma>0$, $\Pr\big(\partial_-C(\hat Q_N)\le\gamma\text{ and }\partial_+C(\hat Q_N)\ge-\gamma\big)\ge1-2\exp\big(-3N\gamma^2/(6bh+8\gamma(b+h))\big)$.
4. (Display (3)) For $0<\epsilon\le1$, $S^{LRS}_\epsilon\subseteq S_\epsilon$.

Significance

The earlier LRS bound has the exponent $-\frac29N\epsilon^2\big(\min\{b,h\}/(b+h)\big)^2$. Theorem 2 replaces the square by the first power. When the critical ratio $b/(b+h)$ approaches $1$ (high service levels), $\min\{b,h\}/(b+h)$ is small, and the required sample size drops accordingly. Raising the service level from 95% to 99% multiplies the sample size the LRS bound requires by 25, but the sample size Theorem 2 requires only by 5 (p. 9). The bound holds for every demand distribution, with no density, continuity or support assumption.

The milestones are reusable on their own. One is the one-sided derivative formula for the expected newsvendor cost. Another is Bernstein's inequality for i.i.d. bounded variables, which Mathlib does not have; it has Hoeffding's inequality. The third is a distribution-free concentration statement for sample quantiles. The result is proved in the paper. To our knowledge none of these statements is formalized; this mission produces machine-checked versions.

Difficulty

A Hoeffding-type argument yields only the squared dependence on $\min\{b,h\}/(b+h)$. The improvement needs the variance $F(1-F)$ of the indicator $\mathbf 1[D\le q]$, which is small near extreme quantiles. Bernstein's inequality applies at a fixed point $q$. The sample quantile, however, is random, and $F$ can jump. The event that $\hat Q_N$ falls left of the target quantile therefore has to be reduced to deviations of $\hat F_N$ at deterministic points, with care at atoms of $D$, where $\partial_+C$ and $\partial_-C$ differ. The final step from the derivatives to the relative regret requires convexity of $C$ and a lower bound on $C(q^*)$ that holds for every distribution.

Formalization scope

• Model. The demand law is a probability measure μ on ℝ, with no sign restriction (the page's $q\ge0$ plays no role here). $F$ is Mathlib's cdf μ, and $\Pr(D<q)$ is μ (Set.Iio q). The realized cost is the published InventoryControl.newsboyLoss h b q d, and $C$ is a Bochner integral. Every statement that evaluates $C$ assumes Integrable id μ ($\mathbb E|D|<\infty$). This is the standing assumption under which $C$ is an expectation. It also rules out the trivializing formalization in which a non-integrable cost integrates to $0$ and every order is $\epsilon$-optimal.
• Quantiles. $q^*$ and $\hat Q_N$ are infima (sInf) of sets that are nonempty and bounded below under the hypotheses, never argmins or choice functions. The sample is x : Fin N → ℝ under the product measure Measure.pi (fun _ => μ), with $N\ge1$; the indices are 0-based.
• Probabilities. "With probability at least $1-p$" is stated as an upper bound $p$ on the outer measure of the failure set.
• $\epsilon$-optimality is the multiplicative form $C(q)\le(1+\epsilon)C(q^*)$, which avoids dividing by $C(q^*)$.
• Pinned hypotheses.
  • The goal and milestone 4 are stated for $0<\epsilon\le1$. The paper says "for any $\epsilon>0$", but both statements fail for large $\epsilon$. Take $b=h=1$, $N=1$ and $D\in\{0,M\}$ with $\Pr(D=M)=0.005$. Then $q^*=0$, and the order $M$ has relative regret $198$, yet at $\epsilon=150$ the bound promises a failure probability of about $1.9\cdot10^{-4}$.
  • In Theorem EC.1, the centred bound $|X^1-\mathbb EX^1|\le c$ is added to the printed $|X^1|\le c$. It holds for the indicators to which the paper applies the inequality.

A complete development needs the following:

• the one-sided derivatives of a convex integral functional;
• Bernstein's inequality for product measures;
• the quantile-to-cdf reduction for the empirical and the true cdf;
• the convexity estimate behind (3).

The Bernstein inequality and the sample-quantile concentration are useful beyond this mission. Contributions to any milestone are welcome.

Selected references

• R. Levi, G. Perakis, J. Uichanco, The Data-Driven Newsvendor Problem: New Bounds and Insights, Operations Research 63(6), 2015. https://doi.org/10.1287/opre.2015.1422. Statements are cited from the authors' accepted manuscript (MIT DSpace), body pp. 7–9 and e-companion pp. ec5–ec6.
• R. Levi, R. O. Roundy, D. B. Shmoys, Provably Near-Optimal Sampling-Based Policies for Stochastic Inventory Control Models, Mathematics of Operations Research 32(4), 2007. https://doi.org/10.1287/moor.1070.0272
• S. N. Bernstein, Theory of Probability, Moscow, 1927.
• P. H. Zipkin, Foundations of Inventory Management, McGraw-Hill, 2000.

Motivation

Comparison theorems for stochastic processes answer a practical question: if one system starts "lower" and moves "upward" less aggressively than another, does it stay below the other one for all time? Queueing, reliability and inventory models use such statements to order performance measures of two systems without computing either distribution. Results of this kind for real-valued Markov chains go back to Kalmykov (1962) and Daley (1968); O'Brien (1975) proved a comparison theorem for real sequences with general (non-Markov) dependence on the past.

Kamae, Krengel and O'Brien (1977) placed these results on a common foundation: an arbitrary partially ordered Polish space, where the state can be a vector, a path, a configuration or a measure. Their tool is a characterization of the stochastic order through monotone couplings, which goes back to Strassen (1965).

Timeline.

• 1962, Kalmykov: comparison of real Markov chains with stochastically monotone kernels.
• 1965, Strassen: existence of probability measures with given marginals; a coupling on a closed set $K \subseteq E \times E$ exists iff the marginals satisfy the matching inequalities.
• 1968, Daley: stochastically monotone Markov chains on $\mathbb R$.
• 1975, O'Brien: comparison theorem for real random sequences with history-dependent kernels.
• 1977, Kamae–Krengel–O'Brien: the order on a partially ordered Polish space, Theorem 1 (six characterizations), and the comparison theorem for sequences in products of such spaces (Theorem 2).

Setting

Let $E$ be a complete separable metric space with a closed partial order $\le$ (the set $\{(x,y) : x \le y\}$ is closed in $E \times E$) and its Borel $\sigma$-algebra. Write $\mathcal M(E)$ for the probability measures on $E$. A set $A \subseteq E$ is increasing if $x \in A$ and $x \le y$ imply $y \in A$; a function $f$ is increasing if $x \le y$ implies $f(x) \le f(y)$. For $P_1, P_2 \in \mathcal M(E)$, $P_1$ is stochastically smaller than $P_2$, written $P_1 \prec P_2$, if

$$\int f\,dP_1 \le \int f\,dP_2 \quad\text{for every bounded measurable increasing } f : E \to \mathbb R .$$

A stochastic kernel $k$ from $E_1$ to $E_2$ assigns to every $x \in E_1$ a probability measure $k(x,\cdot)$ on $E_2$, measurably in $x$. For $P_1 \in \mathcal M(E_1)$, $P_1 * k$ is the measure on $E_1 \times E_2$ with $(P_1 * k)(A_1 \times A_2) = \int_{A_1} k(x, A_2)\,P_1(dx)$, and $P_1^{k}$ is its second marginal. A kernel $k$ on $E \times E$ is upward if $k(x,\cdot)$ is concentrated on $\{y : y \ge x\}$ for every $x$.

For partially ordered Polish spaces $E_1, E_2, \dots$, the products $E^{n} = E_1 \times \cdots \times E_n$ and $E^\infty = \prod_i E_i$ carry the product topology and the coordinatewise order, and are again partially ordered Polish spaces. Given $P_1 \in \mathcal M(E_1)$ and kernels $p_n$ from $E^{n-1}$ to $E_n$ ($n \ge 2$), the measure $P_1 * p_2 * \cdots * p_n$ on $E^{n}$ is the law of $(X_1, \dots, X_n)$ when $X_1 \sim P_1$ and $X_n$ is drawn from $p_n(X_1, \dots, X_{n-1}, \cdot)$; its projective limit is the law of the whole sequence.

Formalization targets

Goal: Theorem 2 (the discrete-time comparison theorem)

Let $P_1, Q_1 \in \mathcal M(E_1)$ and let $p_n, q_n$ be stochastic kernels from $E^{n-1}$ to $E_n$, $n \ge 2$. If $P_1 \prec Q_1$ and

$$p_n(x^{n-1},\cdot) \prec q_n(y^{n-1},\cdot) \quad\text{whenever } x^{n-1} \le y^{n-1},$$

then there are random sequences $(X_n)$, $(Y_n)$ on one probability space, with initial laws $P_1$, $Q_1$ and conditional laws $p_n$, $q_n$, such that

$$P(X_i \le Y_i,\ i = 1, 2, \dots) = 1 .$$

Milestones

1. Theorem 1: for $P_1, P_2 \in \mathcal M(E)$, $P_1 \prec P_2$ is equivalent to each of: a coupling supported on $\{x \le y\}$; a representation $f(Z) \sim P_1$, $g(Z) \sim P_2$ with $f \le g$ and $Z$ real; random variables $X_1 \le X_2$ a.s. with laws $P_1, P_2$; $P_2 = P_1^{k}$ for an upward kernel $k$; and $P_1(B) \le P_2(B)$ for every closed increasing $B$.
2. Proposition 1: under the hypotheses of Theorem 2, $P_1 * p_2 * \cdots * p_n \prec Q_1 * q_2 * \cdots * q_n$ for every $n$.
3. Proposition 2: on $E^\infty$, if all finite-dimensional marginals satisfy $P^{(i)} \prec Q^{(i)}$, then $P \prec Q$.
4. Proposition 3: $\prec$ is preserved under weak convergence.
5. Proposition 4: $P_1 \prec P_2 \prec \cdots$ iff there are random elements $X_1 \le X_2 \le \cdots$ a.s. with $X_i \sim P_i$, iff there are $Z$ real and $f_1 \le f_2 \le \cdots$ with $f_i(Z) \sim P_i$.
6. Corollary 1 (i)–(iii), consequences of Theorem 2 when $E_1 = E_2 = \cdots$: for an increasing set $A$, $P(X_n \in A) \le P(Y_n \in A)$; first entrance times into $A$ satisfy $P(N_x < n) \le P(N_y < n)$; and $E f(X_n) \le E f(Y_n)$ for nondecreasing $f$ whenever the expectations exist.
7. Theorem 3: in a partially ordered Polish space with a compatible vector structure, kernels ordered through the increments they produce yield processes $(S_n)$, $(T_n)$ with $S_1 \le T_1$ and $S_{n+1} - S_n \le T_{n+1} - T_n$ for all $n$, almost surely.

Significance

The result. Theorem 2 turns an inequality between one-step transition laws into a pathwise inequality between whole trajectories. Every functional that is increasing in the path, such as hitting times of increasing sets, maxima, occupation counts or cumulative costs, is then ordered between the two processes, as Corollary 1 illustrates. Because the state space is any partially ordered Polish space, the theorem covers vector-valued queue lengths, networks, and processes whose state is itself a sequence, not just real chains. Theorem 1 is the basic tool for working with the order on such spaces: it converts between integrals, couplings, kernels and closed increasing sets.

Formalizing it. All results are proved in the paper (1977); none is machine-checked. The platform has the special case of Theorem 1 (i) $\Leftrightarrow$ (iv) for $E = \mathbb R^n$ as an open statement (PalmQueueing.Ordering.strassen_st); the general partially ordered Polish case, and the sequence comparison, are new. A formal development yields a reusable library for the stochastic order on general ordered spaces and its interaction with Mathlib's Ionescu-Tulcea construction of process laws (Kernel.trajMeasure).

Difficulty

The obvious argument for Theorem 2 couples the processes step by step: couple $X_1 \le Y_1$, then, given the two histories, couple $X_2 \le Y_2$, and so on. Each step needs a monotone coupling of $p_n(x^{n-1},\cdot)$ and $q_n(y^{n-1},\cdot)$ chosen measurably in the pair of histories; the existence of one coupling for each fixed pair (Theorem 1) does not by itself give a kernel.

Theorem 1's central implication, from the integral inequality to a coupling supported on the closed set $\{x \le y\}$, is Strassen's theorem. On $\mathbb R$ it follows from quantile functions; on a general ordered space there is no such formula. Passing from finite horizons to the infinite sequence is a further step, since an ordering of every finite-dimensional marginal must be turned into an ordering of the infinite-dimensional laws.

Formalization scope

• A partially ordered Polish space is [TopologicalSpace E] [PolishSpace E] [MeasurableSpace E] [BorelSpace E] [PartialOrder E] [OrderClosedTopology E]. This is the paper's standing assumption (Sec. 1, p. 899) and is carried by every theorem, also where a statement says only "Polish space". All sets and functions the paper quantifies over are measurable, also by the standing assumption.
• StochLE P₁ P₂ is defined through bounded measurable monotone functions, exactly as on p. 899. It is not defined through increasing sets, closed increasing sets or couplings, so none of the milestones holds by definition.
• Sequences are 0-based: the paper's $E_i$, $X_i$ are Lean's E (i - 1), X (i - 1), and the paper's kernel $p_n$ ($n \ge 2$) is p (n - 2) : Kernel (Π i : Iic (n - 2), E i) (E (n - 1)), the signature of Mathlib's Ionescu-Tulcea API.
• "Random sequences with initial law $P_1$ and conditional laws $p_n$" is encoded through their joint law, which these data determine: the law of $(X_n)$ is Kernel.trajMeasure P₁ p. Corollary 1 is stated about these laws directly. The goal therefore cannot be met by coupling only one-dimensional marginals or finite prefixes; it asserts the joint laws of the full sequences and one almost-sure event for all indices.
• Hypothesis (4) is required pointwise for all ordered pairs of histories; no stochastic monotonicity of the kernels is assumed.
• Corollary 1 (iii) takes expectations in the extended sense, $E f = E f^+ - E f^-$ in EReal, and "the expectation exists" means the two parts are not both infinite; expectations equal to $\pm\infty$ are covered.
• Theorem 1 (iii) and Proposition 4 (iii) leave the law of the real variable $Z$ free, as the paper does.
• Theorem 3 (p. 904) prints the conditional increment law of $T$ as $q_{n+1}(T_1, \dots, T_n; A + S_n)$; this is a misprint for $A + T_n$, and the Lean states the corrected form. Its compatible vector structure is [AddCommGroup E] [Module ℝ E] [ContinuousAdd E] [ContinuousSMul ℝ E] plus the hypothesis that translates of increasing sets are increasing.
• Not included: Sections 4–5 (continuous time, which needs the Skorohod space).

Welcome contributions: proofs of the milestones, in particular Strassen's coupling theorem on partially ordered Polish spaces, a measurable-selection lemma for monotone couplings, and general lemmas relating StochLE to increasing sets.

The source is the published version (The Annals of Probability), and printed page = PDF page + 898 throughout.

Selected references

• T. Kamae, U. Krengel, G. L. O'Brien, Stochastic Inequalities on Partially Ordered Spaces, The Annals of Probability 5(6), 1977, 899–912. https://doi.org/10.1214/aop/1176995659
• V. Strassen, The existence of probability measures with given marginals, The Annals of Mathematical Statistics 36(2), 1965, 423–439. https://doi.org/10.1214/aoms/1177700153
• G. L. O'Brien, The comparison method for stochastic processes, The Annals of Probability 3, 1975, 80–88. https://projecteuclid.org/journals/annals-of-probability/volume-3/issue-1
• D. J. Daley, Stochastically monotone Markov chains, Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete 10, 1968, 307–317 (as cited in the paper).
• G. I. Kalmykov, On the partial ordering of one-dimensional Markov processes, Theory of Probability and its Applications 7, 1962, 456–459 (as cited in the paper).

Motivation

A flexible manufacturing system is a set of machines through which parts of several types travel along fixed routes; each machine serves several buffers and must pay a set-up time whenever it switches from one buffer to another. Real-time scheduling decides, as the system evolves, which buffer each machine works on. Perkins and Kumar (IEEE Trans. Automat. Control 34, 1989) introduced simple distributed policies for this problem, of which the most natural is the clearing policy: a machine keeps working on a buffer until it is empty, and only then switches. They proved that every clear-a-fraction policy, a subclass of clearing policies, keeps every buffer bounded on acyclic systems whenever each machine has spare capacity, and left open whether clear-a-fraction policies stabilize all systems in which material flows around cycles.

Kumar and Seidman (IEEE Trans. Automat. Control 35(3), 1990, doi:10.1109/9.50339) answered no. Their Example 1 is a single part type that visits two machines in the order 1, 2, 2, 1. Every machine has spare capacity, yet under the clearing policy the buffer levels grow without bound, and they do so even when all set-up times are zero, so the instability comes from machines starving each other rather than from time lost to set-ups. Until then, instability had been suspected to require positive set-up times. Shortly afterwards Lu and Kumar exhibited instability of a static buffer-priority rule in a re-entrant network (IEEE Trans. Automat. Control 36, 1991); together these examples started the study of stability of multiclass queueing networks.

Setting

A manufacturing system has part types $p$ arriving at rates $d_p > 0$. Parts of type $p$ follow a route of length $n_p$: their $i$-th operation is at machine $\mu_{p,i}$, and they wait for it in buffer $b_{p,i}$, where each part needs processing time $\tau_{p,i} > 0$. Machine $m$ serves the buffers $B_m = \{b : \mu_b = m\}$, and switching from $b$ to $b'$ costs set-up time $\delta_{b,b'} \ge 0$.

Flows are continuous (fluid). The level of buffer $b$ at time $t \ge 0$ is $x_b(t) = x_b(0) + u_b(t) - y_b(t) \ge 0$, where $y_b(t)$ is its cumulative output and $u_b(t)$ its cumulative input: $d_p t$ for the first buffer of a route, and the output of the preceding buffer otherwise. Each machine works in runs: run $k$ is a set-up phase of length $\delta_{\beta_{k-1},\beta_k}$ followed by a processing phase on buffer $\beta_k$, during which the buffer is drained at rate $1/\tau_b$ while it is nonempty and passed through at its inflow rate when it is empty. The system is stable if $\sup_{0 \le t < \infty} x_b(t) < \infty$ for every buffer.

A clearing policy (Definition 1) is one in which a machine processing $b$ continues until the first time that $b$ is empty and some other buffer of the same machine is nonempty, and then commences a set-up for one of the nonempty buffers.

Example 1. One part type arrives at rate $d = 1$ and visits machine 1, machine 2, machine 2 and machine 1; its buffers are $1, 2, 3, 4$, so $B_1 = \{1, 4\}$ and $B_2 = \{2, 3\}$. Processing times are $\tau_1, \dots, \tau_4 > 0$, and $\delta_k$ is the time to set up to buffer $k$. The parameters satisfy the critical condition and the capacity condition

$$\tau_2 + \tau_4 > 1, \qquad \tau_1 + \tau_4 < 1, \qquad \tau_2 + \tau_3 < 1. \tag{3–5}$$

The initial state is $x(0) = (\xi, 0, 0, 0)$, with machine 1 set up for buffer 4 and machine 2 set up for buffer 3. Write

$$\lambda = \frac{\tau_4}{1-\tau_2} > 1, \quad \alpha = \frac{(\tau_4+1)(\delta_1+\delta_2)}{1-\tau_2} + \delta_3(\tau_4+1) + \delta_4, \quad \beta = \frac{\tau_4(\delta_1+\delta_2)}{1-\tau_2} + \tau_4\delta_3 + \delta_4.$$

Formalization targets

Goal: Example 1, both cases

Assume (3)–(5).

1. If $\delta_1, \dots, \delta_4 > 0$, there is $\xi_0$ such that for every $\xi \ge \xi_0$ ($\xi > 0$) a clearing trajectory from $(\xi, 0, 0, 0)$ exists, and every such trajectory has

$$\sup_{0 \le t < \infty} x_1(t) = +\infty.$$

2. If $\delta_1 = \dots = \delta_4 = 0$, the same holds for every $\xi > 0$.

Milestone: the Case 1 cycle map

For $\xi$ large enough, every clearing trajectory from $(\xi, 0, 0, 0)$ reaches, at $T_1 = (\lambda + \tau_2/(1-\tau_2))\xi + \alpha$,

$$x(T_1) = (\lambda\xi + \beta, 0, 0, 0),$$

with machines 1 and 2 again set up for buffers 4 and 3.

Milestone: the Case 2 magnification

With zero set-up times and any $\xi > 0$, every clearing trajectory reaches, at $t_5 = (\tau_2+\tau_4)\xi/(1-\tau_2)$,

$$x(t_5) = (\lambda\xi, 0, 0, 0),$$

with machines 1 and 2 again set up for buffers 4 and 3.

Significance

The example shows that the condition $\rho_m < 1$ on every machine, which is necessary for stability and sufficient for the existence of some stabilizing policy, does not make natural distributed policies stable once material flows around a cycle. The throughput of the line falls to $1/(\tau_2+\tau_4) < 1$ part per unit time although each machine could handle the demand. This motivates the paper's two positive results: sufficient conditions under which clear-a-fraction policies are stable (Theorem 1), and a supervisory mechanism that stabilizes any policy (Theorem 2), which are the subjects of the other missions of this series. The example is also an early instance of the phenomenon later studied as instability of multiclass fluid networks under work-conserving policies.

The paper's argument is a stage-by-stage computation of piecewise linear trajectories. No machine-checked version of it exists. A formal proof has to make precise what the paper leaves to the reader: that the clearing rule determines the trajectory, that the stage formulas are what that trajectory does, and that the cycle can be restarted. The formal model of runs, set-ups and the clearing rule built here is the same as in the other missions of the series.

Difficulty

The arithmetic of each cycle is routine once the trajectory is known. The difficulty is in the universal quantifier: the claim covers every clearing trajectory, and the clearing rule is defined implicitly, through "the first time thereafter" at which a buffer is empty and another one is nonempty. At several switching instants the buffer a machine switches to is empty and only starts to fill at that instant, and with zero set-up times a machine may begin a run at an instant where the switching condition already holds. Showing that each switch happens exactly when the paper says, and that the fluid levels then follow the printed formulas (including the reduced rate of a machine working on an empty buffer), is a uniqueness argument for a hybrid system, not a simulation. The existence half asks for the converse: an explicit trajectory, defined for all time, with infinitely many runs whose start times tend to infinity.

Formalization scope

• Time is real ($t \ge 0$); flows are fluid; there are no transport delays or assembly.
• The system is a general structure (part types Fin P, machines Fin M, buffers ⟨p, i⟩ with i : Fin (n p), paper index $i$ = Lean index $i+1$), instantiated as Example 1 with $d = 1$, route $(1,2,2,1)$, and $\delta_{b,b'} = \delta_{b'}$ for $b \ne b'$; staying on a buffer costs nothing.
• A trajectory is a schedule of runs per machine (possibly finitely many, the last lasting forever), with only finitely many run starts in any bounded interval. Processing obeys a rate cap ($y_b$ grows at most at rate $1/\tau_b$, and only while machine $\mu_b$ is in a processing phase of $b$) and runs at full rate while the buffer is nonempty.
• In Definition 1, a target buffer counts as "nonempty" when it is demanding: positive level, or inflow starting at that instant. The no-early-exit condition is imposed on the open processing interval. Under the literal reading (positive level) or a closed interval, the paper's own trajectories are not clearing, and the goal would hold vacuously; the existence clause in the goal rules out that trivialization.
• "Set up for buffer $b$ at time $T$" means the run in force on $(s_k, s_{k+1}]$.
• Unboundedness is stated for buffer 1: for every $C$ there is $t \ge 0$ with $x_1(t) > C$; "$\xi$ large enough" is $\exists \xi_0, \forall \xi \ge \xi_0$.

Useful contributions include lemmas about fluid trajectories that do not depend on the example (continuity of levels, the pass-through rate on an empty buffer, restarting a trajectory at a run boundary), which also serve the other missions of the series.

Selected references

• P. R. Kumar and T. I. Seidman, Dynamic instabilities and stabilization methods in distributed real-time scheduling of manufacturing systems, IEEE Trans. Automat. Control 35(3), 289–298, 1990. https://doi.org/10.1109/9.50339
• J. R. Perkins and P. R. Kumar, Stable, distributed, real-time scheduling of flexible manufacturing/assembly/disassembly systems, IEEE Trans. Automat. Control 34, 139–148, 1989 (reference [18] of the paper).
• S. H. Lu and P. R. Kumar, Distributed scheduling based on due dates and buffer priorities, IEEE Trans. Automat. Control 36, 1991.

Motivation

The newsvendor problem is the basic single-period inventory model: a decision maker orders $q$ units before a random demand $D$ is observed, pays $b>0$ per unit of unmet demand and $h>0$ per unit left over, and minimizes the expected cost $C(q)=E[b(D-q)^+ + h(q-D)^+]$. When the distribution of $D$ is known, the optimal order is the $b/(b+h)$ quantile of $D$. In practice the distribution is unknown and only a sample of past demands is available; the standard data-driven order is the sample average approximation (SAA), the $b/(b+h)$ quantile of the empirical distribution.

R. Levi, G. Perakis and J. Uichanco, The Data-Driven Newsvendor Problem: New Bounds and Insights (Operations Research 63(6), 2015), ask how many samples guarantee that the SAA order is $\epsilon$-optimal with high probability. Levi, Roundy and Shmoys (2007) gave a distribution-free bound. Levi, Perakis and Uichanco show that, asymptotically in $\epsilon$, the right measure of difficulty is a single scalar of the demand law, the weighted mean spread at the critical quantile, and that for the class of log-concave demand distributions (normal, uniform, exponential, logistic, Laplace and many others used in inventory theory) this scalar is bounded below by $\min(b,h)/(b+h)$. That uniform bound turns a distribution-dependent sample-size guarantee into one that holds for every log-concave demand, and is markedly tighter than the distribution-free one.

This mission formalizes the bound (Proposition 2) and the chain of lemmas the paper uses to prove it. The source is the authors' accepted manuscript of the paper (MIT DSpace), cited by its own page numbers.

Setting

Fix costs $b>0$, $h>0$ and the critical ratio $\beta=b/(b+h)\in(0,1)$. A demand law is given by a probability density $f:\mathbb R\to\mathbb R$: measurable, nonnegative, with $\int f=1$. Its cdf is $F(t)=\int_{(-\infty,t]}f$, and its $\beta$ quantile is

$$q^*=\inf\{q : F(q)\ge \beta\}.$$

The absolute mean spread (AMS, Definition 1) at $t$ is the gap between the conditional means above and below $t$,

$$\Delta(t)=E(D\mid D\ge t)-E(D\mid D\le t)=\frac{\int_{[t,\infty)}x f(x)\,dx}{1-F(t)}-\frac{\int_{(-\infty,t]}x f(x)\,dx}{F(t)},$$

and the weighted mean spread (WMS, Definition 2) is $\Delta(q^*)f(q^*)$.

A density is log-concave (Definition 3) if $\log f$ is concave on its support. Equivalently, $f(x)^a f(y)^{1-a}\le f(ax+(1-a)y)$ for all $x,y$ and $a\in[0,1]$; this form allows $f$ to vanish outside an interval.

For a log-concave $f$ and a point $t$ with $f(t)>0$, the number $\gamma_1$ is a supergradient of $\log f$ at $t$, written $\gamma_1\in\partial\log f(t)$, if $\log f(x)\le\log f(t)+\gamma_1(x-t)$ for all $x$ with $f(x)>0$. The paper splits the class $\mathbb L$ of log-concave densities into the subclasses $\mathbb L_{q^*,\gamma_0,\gamma_1}$ of densities with quantile $q^*$, $f(q^*)=\gamma_0$ and $\gamma_1\in\partial\log f(q^*)$, and minimizes the AMS over each subclass (problem (11)). The minimizer is the truncated exponential density

$$\tilde f(x)=\gamma_0e^{\gamma_1(x-q^*)}\ \text{ on }[\underline x,\overline x],\qquad \underline x=q^*+\tfrac1{\gamma_1}\log\!\big(1-\tfrac{\gamma_1}{\gamma_0}\beta\big),\quad \overline x=q^*+\tfrac1{\gamma_1}\log\!\big(1+\tfrac{\gamma_1}{\gamma_0}(1-\beta)\big),$$

and zero elsewhere (display (12)), with AMS $z^*_{q^*,\gamma_0,\gamma_1}$ in closed form.

The Lean development uses the names IsPdf, cdfOf, quantileOf, ams, IsLogSupergradient, tildeF and zStar for these objects, in the namespace DataDrivenNV.WMS.

Formalization targets

Goal: Proposition 2

For every log-concave density $f$,

$$\Delta(q^*)\,f(q^*)\ \ge\ \frac{\min(b,h)}{b+h}.$$

The statement has no further hypothesis: no moment condition, no continuity, no restriction on $\gamma_0,\gamma_1$.

Milestones, in the order of the paper's proof

1. Lemma 1: at the quantile, $-\frac{b+h}{h}\le\frac{\gamma_1}{\gamma_0}\le\frac{b+h}{b}$.
2. Lemma 2: $f(x)\le\gamma_0e^{\gamma_1(x-t)}$ for every $x$.
3. Lemma 3 (Domination Lemma): a density dominated on the support of another, with the same mass left of $t$, has the larger AMS at $t$.
4. Proposition 1: $\tilde f$ belongs to $\mathbb L_{q,\gamma_0,\gamma_1}$ and minimizes the AMS there.
5. The closed form: $\Delta_{\tilde f}(q)=z^*_{q,\gamma_0,\gamma_1}=\frac{\gamma_0}{\gamma_1^2}\big[(\frac{b+h}{h}+\frac{\gamma_1}{\gamma_0})\log(1+\frac{\gamma_1}{\gamma_0}\frac{h}{b+h})+(\frac{b+h}{b}-\frac{\gamma_1}{\gamma_0})\log(1-\frac{\gamma_1}{\gamma_0}\frac{b}{b+h})\big]$.
6. Lemma 4: three elementary inequalities in $\beta\in(0,1)$ and $\eta\in(-\frac1{1-\beta},\frac1\beta)$.

Significance

The proposition is what makes the paper's log-concave sample-size bound (Theorem 4) parameter-free: the probability that the biased SAA order is $\epsilon$-optimal is, asymptotically, at least $1-2\exp(-\frac14N\epsilon\,\Delta(q^*)f(q^*))$ (Theorem 3), and Proposition 2 replaces the unknown $\Delta(q^*)f(q^*)$ by $\min(b,h)/(b+h)$. A manager who only knows that demand is log-concave can then size a sample without estimating any parameter of the law.

The result is proved in the paper; nothing about it is open. As far as is known it has no machine-checked proof. Formalizing it produces a checked lower bound for an inventory-theory quantity together with reusable facts about log-concave densities on the line: exponential envelopes from a supergradient, the quantile-and-slope constraints of Lemma 1, and the comparison of conditional means behind the Domination Lemma.

Difficulty

The definitions are elementary, but the proof passes through an infinite-dimensional optimization problem over a class of densities. The obvious first step, "the AMS is minimized by the most concentrated density", has no direct meaning without fixing the density value and the slope of $\log f$ at the quantile; after fixing them, one needs the envelope of Lemma 2, a stochastic comparison (Lemma 3) and an explicit integral of a truncated exponential. The Domination Lemma as printed is false (a density whose support has a gap to the right of $t$ is a counterexample), so it is formalized under the hypothesis that the dominating density is positive exactly on an interval around $t$, which is how Proposition 1 uses it. The case $\gamma_1=0$ (uniform minimizer) and the two boundary values of $\gamma_1/\gamma_0$ (one-sided exponential minimizers) are not covered by the closed form (12) and must be handled separately in a proof of the goal. Finally, Lemma 1 rests on the monotonicity of the failure rate $f/(1-F)$ and the reversed hazard rate $f/F$ of log-concave laws, which must be proved from log-concavity.

Formalization scope

Densities are functions ℝ → ℝ with IsPdf f (measurable, nonnegative everywhere, integrable, total integral 1); the law of $D$ is never introduced separately. Log-concavity is the published definition ConvexOptimization.LogConcaveOn Set.univ f (power form, zeros allowed). Concavity of Real.log ∘ f is not used, because Real.log 0 = 0 would treat $\log f$ as $0$ off the support. The quantile is an sInf; for $\beta\in(0,1)$ its defining set is nonempty and bounded below. The AMS is the difference of two ratios of integrals. At the quantile both denominators are positive, and the goal and Proposition 1 assume no integrability of $xf(x)$, since log-concave densities have exponential tails. The value $f(q^*)$ is taken pointwise; $q^*$ lies in the interior of the support, where a log-concave density is continuous.

Conventions and disclosed departures from the page:

• "$\gamma_1\in\partial\log f$, the set of all subgradients" is read as the superdifferential of the concave $\log f$ on $\{f>0\}$.
• Lemma 3 carries three added hypotheses: $f_2$ vanishes outside some $[l,u]\ni t$ and is positive on $(l,u)$; $0<F_1(t)<1$; and $xf_1(x)$, $xf_2(x)$ are integrable. The first repairs the printed statement; the other two are the conditions under which Definition 1 makes sense for general densities.
• Proposition 1 and the closed form of $z^*$ are stated for $\gamma_1\neq0$ and $\gamma_1/\gamma_0$ strictly inside the interval of Lemma 1, where (12) is a finite interval.
• The goal has no such restriction.

Assumption 1 of the paper (monotonicity of $f$ beyond $q^*$) and the continuity assumption of §3 are hypotheses of Theorems 3–4 only and are not used here. A formalization in which $f(q^*)$ or $\Delta(q^*)$ takes a junk value (a non-integrable $xf(x)$, a zero denominator) would make the goal trivially true or false; the hypotheses above exclude this, and no statement assumes the conclusion of another.

Contributions welcome: proofs of the milestones in any order; general lemmas on log-concave densities on $\mathbb R$ (exponential tails, integrability of moments, monotone hazard rates, continuity in the interior of the support); and the boundary cases of problem (11).

Selected references

• R. Levi, G. Perakis, J. Uichanco, The Data-Driven Newsvendor Problem: New Bounds and Insights, Operations Research 63(6):1294–1306, 2015. Authors' accepted manuscript, MIT DSpace. https://doi.org/10.1287/opre.2015.1422
• R. Levi, R. O. Roundy, D. B. Shmoys, Provably Near-Optimal Sampling-Based Policies for Stochastic Inventory Control Models, Mathematics of Operations Research 32(4):821–839, 2007. https://doi.org/10.1287/moor.1070.0272