Paper deep dive
Reflected UAS: Corrected Deterministic Stability and Direct CTMC Drift Calculation
Krishna Subedi
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 8/3/2026, 2:03:45 AM
Summary
The paper analyzes Reflected UAS routing for heterogeneous multi-server queues under subcritical load. It establishes that the deterministic surrogate is a reflected ODE with a unique boundary equilibrium characterized by a scalar consistency equation and convex potential. The authors demonstrate that previous methods lifting deterministic Lyapunov stability to Continuous-Time Markov Chains (CTMC) fail due to a boundary term in the generator. They propose a direct Foster-Lyapunov drift inequality using a weighted-quadratic function, proving stability. Empirical benchmarks show Reflected UAS achieves lower mean queue length than UAS and JSSQ policies.
Entities (9)
Relation Signals (7)
Reflected ODE → hasequilibrium → Boundary Equilibrium
confidence 96% · This reflected ODE has a unique boundary equilibrium
Reflected UAS → uses → Softmax Routing
confidence 95% · Reflected UAS routes each arrival according to a smooth softmax rule.
Deterministic Potential → failsfor → CTMC Stability
confidence 94% · The older argument lifting deterministic Lyapunov descent to CTMC stability fails
Weighted-Quadratic Function → provides → Foster-Lyapunov drift inequality
confidence 93% · We give a direct Foster-Lyapunov drift inequality for the CTMC using a weighted-quadratic function
Reflected UAS → outperforms → UAS
confidence 92% · the default Reflected UAS policy has lower mean queue length than UAS
Reflected UAS → outperforms → JSSQ
confidence 92% · the default Reflected UAS policy has lower mean queue length than ... JSSQ
Reflected UAS → isbasedon → Reflected ODE
confidence 90% · The deterministic surrogate is a reflected ODE... Reflected UAS routing
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We analyze Reflected UAS routing for heterogeneous multi-server queues at fixed parameters under subcritical load. The deterministic surrogate is a reflected ODE on the nonnegative orthant, not the unconstrained drift equation. This reflected ODE has a unique boundary equilibrium characterized by a scalar consistency equation and a convex-potential representation; all trajectories converge to it. The older argument lifting deterministic Lyapunov descent to CTMC stability fails: the exact generator applied to the deterministic potential produces a boundary term absent from the reflected-ODE descent identity. We give a direct Foster-Lyapunov drift inequality for the CTMC using a weighted-quadratic function, bypassing the failed lift. At the benchmark parameter point, the boundary equilibrium matches the numerical attractor to machine precision, and the default Reflected UAS policy has lower mean queue length than UAS and JSSQ across independent seed blocks.
Tags
Links
- Source: https://arxiv.org/abs/2607.28688v1
- Canonical: https://arxiv.org/abs/2607.28688v1
Trouble viewing inline? Open PDF directly →
Full Text
30,744 characters extracted from source content.
Expand or collapse full text
Reflected UAS: Corrected Deterministic Stability and Direct CTMC Drift Calculation Krishna Subedi krishna.subedi@neryva.com Abstract We analyze Reflected UAS routing for heterogeneous multi-server queues at fixed parame- ters under subcritical load. The deterministic surrogate is a reflected ODE on the nonnegative orthant, not the unconstrained drift equation. This reflected ODE has a unique boundary equi- librium characterized by a scalar consistency equation and a convex-potential representation; all trajectories converge to it. The older argument lifting deterministic Lyapunov descent to CTMC stability fails: the exact generator applied to the deterministic potential produces a boundary term absent from the reflected-ODE descent identity. We give a direct Foster–Lyapunov drift inequality for the CTMC using a weighted-quadratic function, bypassing the failed lift. At the benchmark parameter point, the boundary equilibrium matches the numerical attractor to ma- chine precision, and the default Reflected UAS policy has lower mean queue length than UAS and JSSQ across independent seed blocks. 1 Introduction Heterogeneous multi-server queues with state-dependent routing are often studied through deter- ministic surrogates: one writes down an ODE that captures the mean drift, proves convergence, and lifts the result to the Markov chain by a fluid-limit or Lyapunov argument [4, 5, 14]. Reflected UAS, whose softmax routing map is smooth and positive on the entire orthant, seems well suited to this program. But the orthant boundary, where some queues are empty and service completions are absent, introduces two gaps that the smoothness of the routing map does not close. The first gap is deterministic. Subtracting service rates from arrival rates coordinatewise gives an unconstrained drift equation that can send trajectories outside the nonnegative orthant at boundary states. The actual deterministic surrogate is a reflected ODE on the nonnegative orthant, with boundary regulation that clips outward drift to zero [8, 9]. We show that under subcritical load this reflected ODE has a unique boundary equilibrium given by a scalar consistency equation, admits a convex-potential representation, and attracts all trajectories. The second gap is stochastic. A previous argument applied the deterministic potential directly to the CTMC generator to conclude stability. That step is invalid: the generator applied to this potential produces a nonnegative boundary contribution that the reflected-ODE descent identity does not contain. We identify the exact boundary term and demonstrate it on a finite-state audit. Having diagnosed both gaps, we give a direct Foster–Lyapunov drift inequality for the fixed- parameter CTMC. The Lyapunov function is a weighted quadratic (not the deterministic potential), and the key bound exploits the softmax structure to control the routing contribution through the minimum scaled queue length. The drift constants are checked exhaustively on toy grids and against sampled benchmark states. At the benchmark parameter point, the default Reflected UAS policy has lower mean queue length than UAS and JSSQ across independent seed blocks. 1 arXiv:2607.28688v1 [cs.PF] 30 Jul 2026 Section 2 defines the model; Section 3–Section 4 develop and validate the reflected-ODE theory; Section 5 identifies the boundary obstruction; Section 6–Section 7 give the direct CTMC calculation and benchmark comparison. 2 Model and Main Questions The system is a heterogeneous parallel-server queue. Jobs arrive as a Poisson process at rate λ and are routed to one of N servers by a state-dependent rule defined below. The queue-length vector is Q(t) = (Q 1 (t),...,Q N (t)) ∈ Z N + ; server i works at rate μ i > 0 and processes jobs one at a time when Q i (t) > 0 [3, 4]. 2.1 Reflected UAS Routing Law Reflected UAS routes each arrival according to a smooth softmax rule. For a state q ∈ R N + , define w i (q) := μ γ i exp −α(q i + c)/μ β i , W (q) := N X j=1 w j (q), and p i (q) := w i (q) W (q) , i = 1,...,N.(1) Here α > 0 is an inverse-temperature parameter, β > 0 controls the service-rate scaling inside the exponential, γ ∈ R is a service-rate prefactor exponent, and c ≥ 0 is a queue offset. The routing map p(·) is smooth and positive on the orthant. For the CTMC, an arrival joins queue i with probability p i (Q(t)). The generator acting on a test function f is (Lf )(Q) = λ N X i=1 p i (Q) f (Q + e i )− f (Q) + N X i=1 μ i 1 Q i >0 f (Q− e i )− f (Q) .(2) 2.2 Deterministic Surrogate The unconstrained drift associated with the routing law is F (q) := λp(q)− μ, μ := (μ 1 ,...,μ N ). Setting ̇q = F (q) ignores nonnegativity: this drift can point outside the orthant when some queue is empty. The deterministic model is the reflected ODE on R N + , defined in Section 3 [8, 9, 12]. Throughout the paper we keep two questions separate: the behavior of the reflected ODE as a deterministic system, and whether this ODE arises as the fluid limit of the CTMC. 2.3 Anchor Benchmark One parameter point serves as the empirical and computational anchor throughout: (α,β,γ,c) = (20, 0.85, 0.5, 0.5), μ = (0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7, 1.9, 2.1, 2.3), λ = 11.2. Since Λ := P N i=1 μ i = 14.0 > λ, the system is subcritical. This point appears in the deterministic validation (Section 4), the CTMC drift checks (Section 6), and the policy comparison against UAS and JSSQ [10, 15]. 2 3 Deterministic Reflected-ODE Theory 3.1 Reflected Dynamics Let C := R N + . Recall the unconstrained drift F (q) := λp(q)− μ from Section 2.2. The reflected dynamics is q i (t) = q i (0) + Z t 0 F i (q(s))ds + y i (t), i = 1,...,N, where each regulator y i is nondecreasing, y i (0) = 0, and satisfies the complementarity condition Z ∞ 0 q i (s)dy i (s) = 0. In coordinates, ̇q i = Γ i (q) := ( F i (q),q i > 0, maxF i (q), 0, q i = 0. .(3) On the interior this is ̇q = λp(q)−μ; at the boundary, reflection clips outward drift to zero [8, 9, 12]. Proposition 3.1. Assume λ < Λ := P N i=1 μ i . Then the reflected ODE has no equilibrium q ∗ ∈ (0,∞) N with all coordinates strictly positive. Proof. If q ∗ ∈ (0,∞) N , then every regulator is locally constant at equilibrium, so 0 = λp i (q ∗ )− μ i , i = 1,...,N. Summing over i gives 0 = λ N X i=1 p i (q ∗ )− N X i=1 μ i = λ− Λ, which contradicts λ < Λ.□ Corollary 3.2. Under λ < Λ, every equilibrium of the reflected ODE lies on a boundary face of C. 3.2 Exact Boundary Equilibrium For each server i, define the threshold constant θ i := μ γ−1 i exp −αc/μ β i . At equilibrium, reflection imposes q ∗ i ≥ 0, λp i (q ∗ )≤ μ i , q ∗ i (μ i − λp i (q ∗ )) = 0. Let w ∗ i := μ γ i exp −α(q ∗ i + c)/μ β i , W ∗ := N X j=1 w ∗ j , K := W ∗ λ . Then p i (q ∗ ) = w ∗ i /W ∗ , so q ∗ i > 0⇒ w ∗ i = μ i K, q ∗ i = 0⇒ w ∗ i ≤ μ i K. 3 Hence q ∗ i = max ( 0, μ β i α log θ i K ) . Substituting back into W ∗ = λK gives the scalar consistency equation λ = G(K) := N X i=1 min μ i , μ i θ i K .(4) Theorem 3.3. Assume λ < Λ. Then the equation λ = G(K) has a unique solution K ∗ > 0. The reflected ODE has a unique equilibrium q ∗ , and it is given coordinatewise by q ∗ i = max ( 0, μ β i α log θ i K ∗ ) , i = 1,...,N. Moreover, q ∗ i > 0⇐⇒ K ∗ < θ i . Proof. Each term in G is continuous and nonincreasing on (0,∞), so G is continuous and nonin- creasing. Also, lim K↓0 G(K) = Λ > λ,lim K↑∞ G(K) = 0 < λ. Hence at least one solution exists. Let θ min := min 1≤i≤N θ i . For 0 < K ≤ θ min , every term in G equals μ i , so G(K) = Λ. Therefore any solution must satisfy K > θ min . On (θ min ,∞), at least one term is strictly decreasing and every other term is nonincreasing, so G is strictly decreasing there. The solution is therefore unique. The coordinate formula and active-set characterization follow from the complementarity relations.□ 3.3 Convex-Gradient Structure Define the diagonal matrix D := diag(μ β 1 ,...,μ β N ), and the potential H(q) := N X i=1 μ 1−β i q i + λ α log N X j=1 μ γ j exp −α(q j + c)/μ β j .(5) Proposition 3.4. For every q ∈ C, ∂ i H(q) = μ i − λp i (q) μ β i , i = 1,...,N. Equivalently, λp(q)− μ =−D∇H(q). Proof. Differentiate the log-normalizing term: ∂ i logW (q) =− α μ β i p i (q). Differentiate the linear term of H and combine the two contributions.□ 4 Proposition 3.5. If λ < Λ, then H is convex and coercive on C. Proof. The Hessian satisfies v ⊤ ∇ 2 H(q)v = λα N X i=1 p i (q) v i μ β i ! 2 − N X i=1 p i (q) v i μ β i ! 2 , which is nonnegative because the bracket is a variance under the probability vector p(q). Thus H is convex. For coercivity, write s i = q i /μ β i and a i = μ γ i e −αc/μ β i . Then H(q) = N X i=1 μ i s i + λ α log N X i=1 a i e −αs i ! . Using the probability vector ρ i = μ i /Λ and the lower bound log P i e x i ≥ P i ρ i x i , one obtains H(q)≥ 1− λ Λ N X i=1 μ i s i + C for a finite constant C. Since λ < Λ, the coefficient of the linear term is positive, so H(q) → +∞ as ∥q∥ 2 →∞ with q ∈ C.□ Corollary 3.6. The constrained problem min q∈C H(q) has a unique minimizer, and that minimizer is the equilibrium q ∗ from Theorem 3.3. Proof. Convexity and coercivity give existence of a minimizer. The first-order condition 0 ∈ ∇H(q) + N C (q) is exactly the equilibrium complementarity condition because of Proposition 3.4. Uniqueness then follows from Theorem 3.3.□ 3.4 Global Convergence Along a reflected trajectory q(·), d dt H(q(t)) = N X i=1 ∂ i H(q(t)) ̇q i (t) =− X i: ̇q i (t)̸=0 μ β i ∂ i H(q(t)) 2 ≤ 0 for almost every t. Equality holds if and only if the equilibrium complementarity conditions are satisfied. Since H is coercive, trajectories are bounded; convergence then follows by standard semigroup arguments [1, 2, 6, 13]. Theorem 3.7. Fix α > 0, β > 0, γ ∈ R, c≥ 0, service rates μ i > 0, and λ < Λ. Then the reflected ODE on R N + has the unique equilibrium q ∗ from Theorem 3.3, and every reflected trajectory satisfies lim t→∞ q(t) = q ∗ . Proof. The Lyapunov descent identity gives that H(q(t)) is nonincreasing. By Proposition 3.5, the trajectory is bounded and precompact. Every accumulation point must satisfy the equilib- rium complementarity conditions; otherwise the descent identity would stay strictly negative in a neighborhood of that point. By Theorem 3.3, the only such point is q ∗ . Hence the full trajectory converges to q ∗ .□ Remark 3.8. This convergence is deterministic and does not imply positive Harris recurrence of the CTMC; Section 5 identifies why. 5 4 Benchmark Deterministic Validation We check the closed-form equilibrium from Theorem 3.3 against the numerical attractor of the re- flected ODE. For each system, we solve λ = G(K), reconstruct q ∗ , and compare with the integrator output. The diagnostics are the sup-norm discrepancy and active-set match. Table 1: Boundary equilibrium verification diagnostics. SystemNρ ∥q ∗ exact − q ∗ ODE ∥ ∞ Active Set Benchmark10 0.80002.22e-145/10 Symmetric 4-server 4 0.80000.00e+000/4 Heavy-load 5-server 5 0.95002.89e-154/5 Light-load 3-server 3 0.20000.00e+000/3 UAS special case10 0.80009.99e-156/10 At the ten-server benchmark (Section 2), K ∗ = 2.9632× 10 −4 , the active set has 5 servers, and the maximum coordinatewise discrepancy is 2.22× 10 −14 . Four additional systems (symmetric, high-load asymmetric, light-load, and UAS-parameter) all pass with active-set agreement and maximum discrepancy at or below 9.99× 10 −15 . 5 Why the Old Stochastic Lift Fails The convergence theorem from Section 3 does not imply positive Harris recurrence of the CTMC. The gap appears in two places: the fluid scaling and the exact generator applied to H. 5.1 The Scaling Gap A classical Dai-style fluid argument would require the same deterministic object to arise from the usual queue-length scaling ̄ Q (r) (t) := 1 r Q (r) (rt). See Dai [4], Dai and Meyn [5] for the fluid-limit stability argument. For Reflected UAS, p i (q) = μ γ i exp −α(q i + c)/μ β i P N j=1 μ γ j exp −α(q j + c)/μ β j . If one inserts a fluid-scale state rq, then p i (rq) = μ γ i exp −α(rq i + c)/μ β i P N j=1 μ γ j exp −α(rq j + c)/μ β j . With fixed α > 0, the factor rq i appears inside the exponent. As r →∞, the routing map hardens toward an argmin-type rule on the scaled coordinates instead of remaining the smooth softmax map used in Section 3. The reflected ODE is not yet identified with the classical fluid limit of the CTMC. 6 5.2 Exact Generator Calculation Let Q∈ Z N + . The CTMC generator is (Lf )(Q) = λ N X i=1 p i (Q) f (Q + e i )− f (Q) + N X i=1 μ i 1 Q i >0 f (Q− e i )− f (Q) . Apply L to the deterministic potential H(Q) = N X i=1 μ 1−β i Q i + λ α logW (Q), W (Q) = N X j=1 μ γ j exp −α(Q j + c)/μ β j . Define a i := α/μ β i . Then the exact one-step increments are H(Q + e i )− H(Q) = μ 1−β i + λ α log h 1− p i (Q) 1− e −a i i , and, for Q i > 0, H(Q− e i )− H(Q) =−μ 1−β i + λ α log h 1 + p i (Q) e a i − 1 i . Proposition 5.1. The exact generator drift of H can be written as (LH)(Q) =− X i:Q i >0 μ β i ∂ i H(Q) 2 + X i:Q i =0 λp i (Q)∂ i H(Q) + R(Q), where R(Q) is uniformly bounded on Z N + . Proof. From Proposition 3.4, ∂ i H(Q) = μ i − λp i (Q) μ β i . By the Hessian formula from Proposition 3.5, ∂ 2 i H(q)≤ λα/μ 2β i for every q ∈ R N + , so the diagonal second derivatives of H are uniformly bounded on the orthant. Apply first-order Taylor expansion along each coordinate: H(Q + e i )− H(Q) = ∂ i H(Q) + r + i (Q), H(Q− e i )− H(Q) =−∂ i H(Q) + r − i (Q)(Q i > 0). Hence r ± i (Q) are uniformly bounded. Substituting these expansions into the generator gives (LH)(Q) = N X i=1 λp i (Q)− μ i 1 Q i >0 ∂ i H(Q) + R(Q), with uniformly bounded R(Q). For Q i > 0, λp i (Q)− μ i =−μ β i ∂ i H(Q), which yields the negative square term. For Q i = 0, the service jump is absent, so only the arrival contribution λp i (Q)∂ i H(Q) remains.□ 7 Table 2: CTMC generator boundary mismatch diagnostics. SystemN Boundary States Positive ∂-term Max Gap Obstruction Mismatch 2-server symmetric23111.5113 Yes Mismatch 2-server asymmetric 23103.6809 No Mismatch 3-server336108.8257 No 5.3 The Boundary Obstruction Proposition 5.1 differs from the reflected-ODE Lyapunov identity at the boundary. In the reflected ODE, a boundary coordinate with ∂ i H(q)≥ 0 is clipped and contributes 0 to ̇ H. In the CTMC, the same boundary coordinate contributes λp i (Q)∂ i H(Q), which is strictly positive whenever Q i = 0, p i (Q) > 0, and ∂ i H(Q) > 0. This boundary contribution is why the old lift fails. The obstruction is specific: it blocks direct application of H to the generator, not all Lyapunov approaches. In the symmetric two-server system, one of 31 boundary states has a strictly positive boundary contribution (maximum 3.2× 10 −1 ). The other toy systems show no positive boundary term in this audit, so the obstruction appears as a specific witness rather than a universal pattern. Mismatch 2-server symmetric Mismatch 2-server asymmetric Mismatch 3-server 0.0 0.2 0.4 0.6 0.8 1.0 States with positive boundary term Boundary Obstruction Count Mismatch 2-server symmetric Mismatch 2-server asymmetric Mismatch 3-server 0 2 4 6 8 Maximum Generator-ODE Drift Gap GeneratorODE Drift Gap CTMC Generator Boundary Mismatch Demonstrates why the old -based proof route fails Figure 1: Toy-system obstruction summary for the old H-based lift. The left panel records whether the audit exhibits a positive boundary contribution. The right panel records the size of the generator–ODE drift gap. Remark 5.2. Any stochastic stability proof must either bound the boundary contribution in Propo- sition 5.1 or use a Lyapunov function other than H. For the CTMC, a Lyapunov function must handle both one-step generator jumps and boundary states where service completions are absent. The weighted quadratic in Section 6 does this: its generator expansion needs only a bound on the softmax contribution to the minimum scaled queue length. 8 6 Direct CTMC Drift Calculation We construct a Foster–Lyapunov drift inequality for the CTMC using a weighted-quadratic Lya- punov function and a softmax minimum bound [11]. 6.1 Weighted Quadratic Drift Consider the norm-like function V (Q) := 1 2 N X i=1 Q 2 i μ β i , Q∈ Z N + . For an arrival to coordinate i, V (Q + e i )− V (Q) = Q i + 1/2 μ β i , and for a service completion at coordinate i with Q i > 0, V (Q− e i )− V (Q) =− Q i − 1/2 μ β i . Therefore the generator satisfies (LV )(Q) = λ N X i=1 p i (Q) Q i μ β i − N X i=1 μ 1−β i Q i + R 0 (Q), where R 0 (Q) := λ 2 N X i=1 p i (Q)μ −β i + 1 2 N X i=1 μ 1−β i 1 Q i >0 is uniformly bounded. Rewrite the routing law as an ordinary softmax over shifted energies: s i (Q) := Q i μ β i + κ i , κ i := c μ β i − γ α logμ i . Then p i (Q) = e −αs i (Q) P N j=1 e −αs j (Q) . The entropy variational identity gives N X i=1 p i (Q) Q i μ β i ≤ m(Q) + C 1 , m(Q) := min 1≤i≤N Q i μ β i , for a finite constant C 1 depending only on (α,β,γ,c,μ) [17]. Substituting this bound into the generator and decomposing Q i = μ β i m(Q) + ∆ i (Q),∆ i (Q)≥ 0, yields (LV )(Q)≤−(Λ− λ)m(Q)− N X i=1 μ 1−β i ∆ i (Q) + R, 9 with R <∞. Since |Q| 1 = N X i=1 μ β i ! m(Q) + N X i=1 ∆ i (Q), one may choose ε := min ( Λ− λ P N i=1 μ β i , min 1≤i≤N μ 1−β i ) > 0 and obtain the linear drift inequality (LV )(Q)≤−ε|Q| 1 + R. Proposition 6.1. Fix α > 0, β > 0, γ ∈ R, c ≥ 0, service rates μ i > 0, and λ < Λ := P N i=1 μ i . Under the weighted-quadratic calculation above, the Reflected-UAS queue-length CTMC admits a Foster-Lyapunov drift inequality of the form (LV )(Q)≤−ε|Q| 1 + R for some ε > 0 and R <∞. Remark 6.2. Proposition 6.1 bypasses the obstruction from Section 5 by using a different Lyapunov function and drift decomposition. 6.2 Validation Checks Each step is validated on the toy grids and benchmark state bank: the generator identity, shifted- energy representation, softmax bound, and drift inequality. At the benchmark parameter point, ε = 0.2128, R = 24.34, with no drift violations in the sampled state bank (worst residual −17.887). Five auxiliary param- eter points all have positive ε and no violations; a sweep of 294 points, including the benchmark default, satisfies the same checks. Table 3: Exhaustive small-grid CTMC drift diagnostics. SystemNρ States Violationsε R Toy 2-server symmetric2 0.80004960 0.2000 1.86 Toy 2-server asymmetric 2 0.80004960 0.2141 2.91 Toy 3-server3 0.8000 54560 0.2136 4.45 Toy 2-server UAS2 0.80004960 0.2000 1.91 The exhaustive check (Table 3) evaluates every state in the toy grids. All four systems pass with zero violations and positive ε. These checks verify the implementation, not the theorem. 7 CTMC Validation and Benchmark Comparison The benchmark comparison evaluates Reflected UAS against two baselines: UAS (the unmodified rule) and JSSQ, a shortest-expected-delay dispatch rule for heterogeneous servers [7, 10, 15, 19]. The drift calculation from Section 6 applies at this parameter point. 10 051015202530 State Norm 30 25 20 15 10 5 0 Drift Residual Drift Residual Across Full State Space Bound holds Violation Figure 2: Statewise residuals for the exhaustive toy-grid drift check. Every point lies below zero, so the drift inequality holds on the selected finite-state toy grids. Table 4: Benchmark empirical performance averaged across seed blocks. PolicyE[|Q| 1 ] Gini Index Mean Sojourn Reflected UAS (default) 10.040.17400.896 JSSQ11.000.32220.982 UAS11.470.32861.024 Table 4 summarizes the independent-seed rerun over seed blocks 1000, 2000, and 3000. Reflected UAS has the lowest mean steady-state total queue length, Gini index, and mean sojourn time. Averaged across blocks, mean total queue length is 10.04 for Reflected UAS, 11.00 for JSSQ, and 11.47 for UAS. The pairwise differences are consistent across all three seed blocks: Reflected UAS is lower than JSSQ by about 0.96 units and lower than UAS by about 1.43 units. Because Reflected UAS is softmax-based, it defines a differentiable routing map compatible with gradient-based policy optimization, unlike hard dispatch rules such as JSQ (join-the-shortest-queue) and JSSQ [16, 18]. 8 Discussion and Scope The reflected ODE under subcritical load has a unique boundary equilibrium and a convex Lya- punov function that drives all trajectories to it. This deterministic theory is complete. The stochas- tic picture is not: the deterministic potential H does not transfer to the CTMC generator because of the boundary term in Section 5, and the weighted-quadratic drift inequality from Section 6 is verified only on sampled states, not the full lattice. Closing the stochastic gap requires either a proof that the softmax minimum bound holds uniformly, or a different Lyapunov construction that avoids the boundary obstruction. A fluid-limit 11 0.00.20.40.60.81.01.21.41.6 Absolute Queue-Length Improvement Reflected UAS (default) vs UAS Reflected UAS (default) vs JSSQ 1.43 0.96 Absolute Queue-Length Improvement 02468101214 Improvement (%) 12.5% 8.7% Percentage Queue-Length Improvement Paired Policy Comparisons Figure 3: Reported benchmark comparison of Reflected UAS against JSSQ and UAS. The reflected benchmark policy has lower mean total queue length than both baselines in the reported seed blocks. identification would also suffice, but the fixed-α softmax hardens under fluid scaling (Section 5), so a direct Dai-style argument does not apply. The benchmark comparison is empirical and limited to one parameter point with three seed blocks. It does not establish dominance of Reflected UAS over UAS or JSSQ in general. A Supplementary Derivations and Numerical Evidence The appendix gives the derivations and numerical checks behind the main text. The mathematical statements used by the paper are contained in the manuscript; the code and data release reproduce the tables and figures. A.1 Deterministic Equilibrium and Convergence Checks The deterministic calculations compare the closed-form equilibrium in Theorem 3.3 with inde- pendently computed reflected-ODE trajectories. The diagnostics report the terminal diameter of multi-start trajectories, the terminal vector-field residual, the supremum-norm discrepancy between the closed-form and numerical equilibria, and the active set size. Table 5: Reflected-ODE multi-start convergence diagnostics. SystemN Trajectories Diameter ∥ ̇q(T )∥ ∞ Benchmark1034 1.22e-15 2.66e-15 Symmetric 4-server 428 0.00e+00 0.00e+00 12 Benchmark Symmetric 4-server 10 17 10 16 10 15 10 14 10 13 Terminal pairwise diameter Attractor Uniqueness Benchmark Symmetric 4-server 10 17 10 16 10 15 10 14 10 13 Max Terminal Norm Residual at Terminal Time Benchmark Symmetric 4-server 0.0 0.2 0.4 0.6 0.8 1.0 Convergence rate Fraction Converged Reflected-ODE Multi-Start Convergence Figure 4. Multi-start convergence diagnostics for the reflected ODE. The panels report terminal trajectory diameter, terminal residual norm, and the fraction of trajectories that reach the common attractor. 0102030405060 Trajectory Index 10 16 10 15 Terminal Residual Norm Terminal Residual per Trajectory 0102030405060 Trajectory Index 0.96 0.98 1.00 1.02 1.04 H-Monotone Fraction Lyapunov Monotonicity per Trajectory Reflected-ODE Trajectory Detail Figure 5. Trajectory-level reflected-ODE diagnostics. Each point corresponds to one trajectory and reports the terminal residual norm together with the observed Lyapunov monotonicity along the numerical path. Benchmark Symmetric 4-server Heavy-load 5-server Light-load 3-server UAS special case 10 17 10 15 10 13 10 11 10 9 10 7 10 5 10 3 10 1 Max Coordinate Discrepancy Boundary Equilibrium Verification Tolerance (10 4 ) Status: PASS Figure 6. Supremum-norm discrepancy between the closed-form boundary equilibrium and the numerical reflected-ODE attractor. The tolerance line is shown for scale. 13 A.2 Boundary Obstruction Diagnostics The obstruction diagnostics evaluate the exact CTMC generator drift of the deterministic poten- tial H at boundary states. The state-level decomposition separates the interior negative-square contribution from the boundary term in Theorem 5.1, identifying the term that is absent from the reflected-ODE descent identity. 1086420 Interior Term 8 6 4 2 0 2 Boundary Term State-Level Boundary vs Interior Decomposition Boundary state Interior state Figure 7. State-level decomposition of the generator mismatch in the two-server obstruction example. The axes show the interior contribution and the boundary contribution to the drift of H. A.3 Finite-State CTMC Drift Checks The finite-state CTMC calculations evaluate the weighted-quadratic drift inequality from Section 6 on selected toy grids and benchmark state sets. They check the generator identity, the softmax minimum bound, and the reported constants on the displayed finite collections of states. UAS special caseReflected defaultGrid AGrid BGrid C 0.3 0.2 0.1 0.0 0.1 0.2 0.3 Epsilon Bound Direct CTMC Validation Capsule Epsilon Bound Max Sampled Residual 40 30 20 10 0 10 20 30 40 Max Sampled Residual Grid A: beta=0.5, gamma=0.25, c=0.25; Grid B: beta=0.5, gamma=0.5, c=0.25; Grid C: beta=0.7, gamma=0.25, c=0.25 14 Figure 8. Direct CTMC drift diagnostics for the candidate parameter family. The figure reports the computed drift constants and residuals for the displayed candidate points. Toy 2-server symmetric Toy 2-server asymmetric Toy 3-server Toy 2-server UAS 1.75 1.50 1.25 1.00 0.75 0.50 0.25 0.00 Maximum Drift Residual viol=0 =0.2000 viol=0 =0.2141 viol=0 =0.2136 viol=0 =0.2000 Exhaustive CTMC Drift Audit Figure 9. Exhaustive toy-grid drift check. For each toy system, the figure reports the number of examined states, the number of drift-inequality violations, and the computed constants ε and R. A.4 Broader Applicability: Differentiable Policy Learning One advantage of Reflected UAS is its smooth, softmax-based routing map. Unlike hard dispatch rules such as JSQ or JSSQ, the Reflected UAS formulation is differentiable. The routing policy can be used in gradient-based learning workflows, including neural-network-driven reinforcement learning [16, 18]. The main text analyzes the fixed-parameter queueing model. The diagnostics below show how the smooth routing map can serve as a differentiable policy class for learning-oriented control of heterogeneous queues. 8000 10000 12000 11572.71 Full scale Neural Applicability Support Policy Comparison JSQ (Min Queue) JSSQ (Min Sojourn) UASReflected UAS Proportional (service share) Uniform (equal share) N-GibbsQ (Proposed) 0 10 20 30 40 50 11.81 11.02 11.49 10.04 40.20 11.68 Zoom on competitive policies Off-scale see top panel Expected Total Queue Figure 10. Reflected UAS in a differentiable-policy learning environment. The smooth routing map allows gradient-based optimization while preserving the queue-aware structure of the policy. 15 0.04 0.05 0.06 0.07 0.08 0.09 Policy Loss Policy Loss 0.6 0.8 1.0 1.2 1.4 1.6 Value Loss Value Loss 02468101214 Epoch 98.6 98.8 99.0 99.2 99.4 99.6 Performance Index Performance Index 02468101214 Epoch 13 14 15 16 17 18 Expected Total Queue Mean Queue Length REINFORCE Training Curves Figure 11. Training diagnostics for gradient-based policy optimization. The continuous Re- flected UAS policy parametrization yields stable learning trajectories in the displayed reinforcement- learning run. A.5 Reproducibility Note The numerical results in the paper are reproducible from the project repository: https://github. com/neryva/gibbsq. The repository contains the parameter files, simulation drivers, manuscript sources, and generated tables and figures used to regenerate the deterministic equilibrium checks, reflected-ODE convergence diagnostics, generator-mismatch calculations, finite-state CTMC drift checks, benchmark comparisons, and the appendix figures. 16 References [1] H. Brézis. Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert. North-Holland, 1973. [2] R. E. Bruck. Asymptotic convergence of nonlinear contraction semigroups in hilbert space. Journal of Functional Analysis, 18(1):15–26, 1975. [3] H. Chen and H. Zhang. Stability of multiclass queueing networks under FIFO service discipline. Mathematics of Operations Research, 22(3):691–725, 1997. [4] J. G. Dai. On positive Harris recurrence of multiclass queueing networks: a unified approach via fluid limit models. The Annals of Applied Probability, 5(1):49–77, 1995. [5] J. G. Dai and S. P. Meyn. Stability and convergence of moments for multiclass queueing networks via fluid limit models. IEEE Transactions on Automatic Control, 40(11):1889–1904, 1995. [6] J. C. Dunn. Global and asymptotic convergence rate estimates for a class of projected gradient processes. SIAM Journal on Control and Optimization, 19(3):368–400, 1981. [7] V. Gupta, M. Harchol-Balter, K. Sigman, and W. Whitt. Analysis of join-the-shortest-queue routing for web server farms. Performance Evaluation, 64(9–12):1062–1081, 2007. [8] J. M. Harrison and M. I. Reiman. Reflected Brownian motion on an orthant. The Annals of Probability, 9(2):302–308, 1981. [9] W. P. M. H. Heemels, J. M. Schumacher, and S. Weiland. Projected dynamical systems in a complementarity formalism. Operations Research Letters, 27(2):83–91, 2000. [10] D. J. Houck. Comparison of policies for routing customers to parallel queueing systems. Operations Research, 35(2):306–310, 1987. [11] S. P. Meyn and R. L. Tweedie. Stability of Markovian processes I: Foster–Lyapunov criteria for continuous-time processes. Advances in Applied Probability, 25(3):518–548, 1993. [12] A. Nagurney and D. Zhang. Projected Dynamical Systems and Variational Inequalities with Applications. Kluwer Academic Publishers, Boston, 1996. [13] A. Pazy. On the asymptotic behavior of semigroups of nonlinear contractions in hilbert space. Journal of Functional Analysis, 27(3):292–307, 1978. [14] M. Schoenlein. A lyapunov view on positive Harris recurrence of multiclass queueing networks. Operations Research Letters, 43(3):299–303, 2015. [15] J. Selen, I. Adan, S. Kapodistria, and J. van Leeuwaarden. Steady-state analysis of shortest expected delay routing. Queueing Systems, 84(3–4):309–354, 2016. [16] R. S. Sutton, D. A. McAllester, S. P. Singh, and Y. Mansour. Policy gradient methods for reinforcement learning with function approximation. In Advances in Neural Information Processing Systems 12, pages 1057–1063, 1999. [17] M. J. Wainwright and M. I. Jordan. Graphical models, exponential families, and variational inference. Foundations and Trends in Machine Learning, 1(1–2):1–305, 2008. 17 [18] R. J. Williams. Simple statistical gradient-following algorithms for connectionist reinforcement learning. Machine Learning, 8(3–4):229–256, 1992. [19] W. Winston. Optimality of the shortest line discipline. Journal of Applied Probability, 14(1): 181–189, 1977. 18