On Computing the Mean Stationary Response Time in M/M/1 Fork-Join (n,k) Queues

Copyright (c) 2026 - Pierre M. Fiorini, Ph.D.

Busy-period grammar computation and paper-fragment generator.

This module is the executable companion to the paper

        "On Computing the Mean Stationary Response Time
        in M/M/1 Fork-Join $(n,k)$ Queues"

It provides two deliberately separate command-line workflows.

Copyright (c) 2026 - Pierre M. Fiorini, Ph.D.

===============================================================================
TWO EXECUTION OPTIONS
===============================================================================

Option 1: regular computation, no LaTeX generation
--------------------------------------------------

    python busy_period_grammar_fresh.py compute \
        --n 3 --k 2 --lambda-rate 0.4 --mu 1 1 1 \
        --tolerance 1e-10 --max-expansions 500000

This mode runs the analytic best-first frontier computation for one model and
prints the usual numerical diagnostics.  It does not write tables, figures,
PGFPlots code, or paper fragments.  Add --des to run an independent-replication
DES validation for the same model.  This is the normal one-case computational
workflow.

Option 2: paper generation with fresh values
--------------------------------------------

    python busy_period_grammar_fresh.py paper \
        --out generated_fresh \
        --des-reps-small 30 --des-warmup-small 5000 --des-measured-small 50000 \
        --des-reps-large 20 --des-warmup-large 3000 --des-measured-large 20000

This mode recomputes all analytic table entries from the best-first probability
mass frontier solver, runs fresh DES replications, and then writes every table,
TikZ diagram, and PGFPlots graph fragment used in the paper.  Unlike earlier
versions of this generator, paper mode does not populate validation tables from
frozen manuscript constants.  The displayed analytic values, DES means, DES
half-widths, absolute differences, pending masses, and expansion counts are all
written from the current execution.

For quick debugging, paper mode accepts --skip-des, --case-filter, and small DES
replication settings.  For manuscript-grade regeneration, use the DES settings
reported in the paper or larger values.

===============================================================================
MODEL IMPLEMENTED BY THE ANALYTIC FRONTIER SOLVER
===============================================================================

The analytic model is the infinite-buffer Markovian cancel-on-k fork-join queue.
Jobs arrive as a Poisson process of rate lambda.  Each job creates one task at
each of n FCFS M/M/1 branches.  Branch i has exponential service rate mu_i.  The
job departs when any k of its n copies complete service, and all unfinished
copies of that same job are cancelled immediately.

The state used by the analytic solver is a queue-length state, not a job-id
state.  This is valid because all arrivals are synchronous and all branches are
FCFS.  A horizontal queue depth represents one active job.  If a service
completion occurs at branch i and the branch length before completion is r, the
number of unfinished copies of that depth-r job is

    M_r(q) = #{j : q_j >= r}.                                      (eq:Mr)

After the mechanical completion, cancellation fires exactly when

    M_r(q) - 1 < n - k + 1.                                  (eq:canceltest)

When cancellation fires, all residual copies of that job are removed by

    q'_j = q_j - 1{q_j >= r}.                                  (eq:cancelvec)

When cancellation does not fire, only the completing branch is decremented.
These formulas are implemented by full_transitions() and grouped_transitions().

The full heterogeneous state is q=(q_1,...,q_n).  The total rate out of a
non-empty state is

    Sigma(q) = lambda + sum_{i:q_i>0} mu_i.                       (eq:sigma)

The active-job count is the number of occupied horizontal queue depths:

    N_sys(q) = max_i q_i.                                      (eq:Nsysfull)

For homogeneous or partially homogeneous systems, branches are partitioned into
rate classes.  For rate class g, c_{g,a} is the number of branches in group g
with queue length a.  The grouped state is

    x = (c_1,...,c_G),     sum_a c_{g,a} = n_g.               (eq:groupstate)

The grouped exit rate and active-job count are

    Sigma(x) = lambda + sum_g mu_g sum_{a>=1} c_{g,a},          (eq:groupsigma)
    N_sys(x) = max{a : sum_g c_{g,a} > 0}.                     (eq:groupsigma)

The grouped rank count is

    M_a(x) = sum_h sum_{b>=a} c_{h,b}.                          (eq:groupMr)

Grouping is exact symmetry lumping; it is not a mean-field approximation and it
is not a finite-buffer truncation.

===============================================================================
RENEWAL-REWARD EQUATIONS AND FRONTIER EXPANSION
===============================================================================

The empty state is regenerative.  Starting immediately after the initiating
arrival, let tau_0 be the first return time to the empty state.  The steady-state
mean job response time is computed by the renewal-reward identity

    E[T] = E[ integral_0^{tau_0} N_sys(X(t)) dt ] / E[A(tau_0)].   (eq:rr-intro)

The first-step residual reward equations are

    R(x) = N_sys(x)/Sigma(x) + sum_{y!=0} P(x,y) R(y),             (eq:Rfirst)
    J(x) = lambda/Sigma(x) + sum_{y!=0} P(x,y) J(y).              (eq:Jfirst)

The initiating arrival is counted separately.  Therefore the displayed response
ratio is

    E[T] = R(1) / (1 + J(1)).                                  (eq:grammarrr)

AnalyticalBestFirstFrontierSolver evaluates these equations forward over the
busy-period grammar.  It stores a sparse pending frontier measure m(x).  Popping
a state with mass m contributes

    m * N_sys(x)/Sigma(x)

to accumulated job-time reward, contributes

    m * lambda/Sigma(x)

to expected future admitted jobs, and pushes each nonterminal successor with its
embedded-chain probability.  Empty successors are absorbed.

The traversal order is best-first: the largest pending mass is expanded next.
This is an evaluation order, not a different stochastic model.  If the frontier
is fully drained with no pruning and no resampling, the result equals the exact
analytic regenerative ratio in exact arithmetic and is numerically exact up to
floating-point roundoff.

===============================================================================
ERROR-ACCOUNTING EQUATIONS
===============================================================================

After finitely many expansions, the current non-empty frontier F has residuals

    R_F = sum_{x in F} m(x) R(x),       J_F = sum_{x in F} m(x) J(x).  (eq:residuals)

With no pruning or resampling,

    T = (R_hat + R_F) / (J_hat + J_F).                         (eq:frontier_exact)

Subtracting the prefix estimate T_hat=R_hat/J_hat gives

    T - T_hat = (R_F J_hat - R_hat J_F)/(J_hat(J_hat+J_F)).     (eq:error_decomp)

Pending mass alone is not an error bound on an unbounded state space.  If U_R
and U_J are nonnegative reward envelopes satisfying the paper's inequalities,
then the computable residual interval is

    R_hat/(J_hat + Jbar_F) <= T <= (R_hat + Rbar_F)/J_hat.       (eq:interval)

When deterministic frontier resampling is used, the moved mass is

    D_move = 1/2 sum_x |F_tilde(x) - F(x)|,                       (eq:moved)

and bounded test functions satisfy the total-variation perturbation bound

    |sum_x f(x)F_tilde(x)-sum_x f(x)F(x)| <= osc_F(f) D_move.    (eq:tvbound)

The mass-conservation diagnostic is

    epsilon_mass = |m_absorbed + m_pending + m_pruned - 1|.       (eq:masserr)

===============================================================================
DES VALIDATION AND CONFIDENCE INTERVALS
===============================================================================

The script includes an event-level DES simulator with the same FCFS cancel-on-k
semantics as the analytic grammar.  The simulator uses job identifiers, branch
queues, per-branch in-service tokens, and a global event heap.  Tokens invalidate
stale service-completion events whenever an in-service residual copy is cancelled
by a job completion on another branch.

For R independent replications, after warm-up, replication r has mean

    Tbar_r = (1/M) sum_{j=1}^M T_{rj}.                           (eq:repmean)

The reported DES confidence interval is

    Tbar_DES +/- t_{0.975,R-1} s_R/sqrt(R),
    s_R^2 = (1/(R-1)) sum_r (Tbar_r - Tbar_DES)^2.               (eq:ci)

The DES interval is a simulation sampling interval only.  The analytic grammar
value is deterministic and is accompanied by frontier diagnostics rather than a
sampling interval.