An ODE for an Overloaded X Model Involving a Stochastic Averaging Principle

We study an ordinary differential equation (ODE) arising as the many-server heavy-traffic fluid limit of a sequence of overloaded Markovian queueing models with two customer classes and two service pools. The system, known as the X model in the call-center literature, operates under the fixed-queue-ratio-with-thresholds (FQR-T) control, which we proposed in a recent paper as a way for one service system to help another in face of an unanticipated overload. Each pool serves only its own class until a threshold is exceeded; then one-way sharing is activated with all customer-server assignments then driving the two queues toward a fixed ratio. For large systems, that fixed ratio is achieved approximately. The ODE describes system performance during an overload. The control is driven by a queue-difference stochastic process, which operates in a faster time scale than the queueing processes themselves, thus achieving a time-dependent steady state instantaneously in the limit. As a result, for the ODE, the driving process is replaced by its long-run average behavior at each instant of time; i.e., the ODE involves a heavy-traffic averaging principle (AP).

occurs. We want the control to automatically detect the overload. The FQR-T control is designed to prevent sharing of customers (i.e., sending customers to be served at the other-class service pool) when sharing is not needed, and automatically activate sharing when the system becomes overloaded due to a sudden shift in the arrival rates.
This paper is the third in a series of five papers. First, In [15] we initiated study of this overload-control problem and proposed the FQR-T control; see [15] for a discussion of related literature. We used a heuristic stationary fluid approximation to derive the optimal control when a convex holding cost is charged to the two queues during the overload incident. Within that framework, we showed with simulations that FQR-T outperforms the best fixed allocation of servers, even when the new arrival rates are known. The stationary point of the fluid model was derived using a heuristic flow-balance argument, which equates the rate of flow into the system to the rate of flow out of the system, when the system is in steady state.
Second, in [16] we applied the heavy-traffic AP as an engineering principle in order to justify the ODE considered here to describe the transient fluid approximation of the X system under FQR-T after an overload has occurred. We observed that the FQR-T control is driven by a queue-difference stochastic process, which operates in a faster time scale than the queueing processes themselves, so that it should achieve a time-dependent steady state instantaneously in the MS-HT limit, i.e., as the scale (arrival rate and number of servers) increases; see §3.1. We argued heuristically that the ODE should arise as the limit of a properly-scaled sequence of overloaded X-model systems, provided that the driving process is replaced by its long-run average behavior at each instant of time. We performed simulation to justify the approximations.
The present paper and the next two provide mathematical support. The present paper establishes important properties of the ODE suggested in [16]. The fourth and fifth papers prove limits. In [17] we prove that the fluid approximation, as a deterministic function of time, arises as the MS-HT limit of a sequence of X models; i.e., we prove a functional weak law of large numbers (FWLLN). This FWLLN is based on the AP; see [4,8] for previous examples. In [18] we prove the corresponding functional central limit theorem (FCLT) that describes the stochastic fluctuations about the deterministic fluid path.
We prove convergence to the ODE in [17] by the standard two-step procedure, described in Ethier and Kurtz [5]: (i) establishing tightness and (ii) uniquely characterizing the limit process. The tightness argument follows familiar lines, but characterizing the limit process turns out to be challenging.
Indeed, characterizing the limit process depends on the results here. Thus, the present paper provides a crucial ingredient for the limits established in [17,18].
The AP makes the ODE unconventional. The AP creates a singularity region, causing the ODE not to be continuous in its full state space. Hence, classical results of ODE theory, such as those establishing existence, uniqueness and stability of solutions, cannot be applied directly. Moreover, existing algorithms for numerically solving ODE's cannot be applied directly either, since the solution to the ODE requires that the time-dependent steady state of the fast-time-scale process (FTSP) be computed at each instant. Nevertheless, we establish the existence of a unique solution to the ODE, show that there exists a unique stationary point; and show that the fluid process converges to its stationary point as time evolves. Moreover, we show that the convergence to stationarity is exponentially fast. The key is a careful analysis of the FTSP, which we represent as a quasi-birth-and-death (QBD) process. Finally, we provide a numerical algorithm for solving the ODE based on the matrix-geometric method [10].
Here is how the rest of this paper is organized: The next two sections provide background. In §2 we elaborate on the X queueing model and the FQR-T control; that primarily is a review of [15]. In §3 we provide a brief overview of the MS-HT scaling and a heuristic explanation of the AP. In §4 we introduce the ODE that we study in subsequent sections. In §5 we state out main result, establishing the existence of a unique solution. In §6 we establish properties of the FTSP, which depends on the state of the ODE, and whose steady-state distribution influences the evolution of the ODE. In §7 we define the state space of the ODE, and prove the main theorem about existence of a unique solution. In §8 we establish the existence of a unique stationary point and show that the fluid solution converges to that stationary point as time evolves. In §9 we prove that a solution converges to stationarity exponentially fast. In §10 we provide conditions for global state space collapse, i.e., for having the AP operate for all t ≥ 0. In §11 we develop an algorithm to numerically solve the ODE (given an initial condition), based on the theory developed in the previous sections. We conclude in §12 with one postponed proof.
Additional material appears in an appendix, available from the authors' web pages. There we analyze the system with an underloaded initial state, and show that the approximating fluid models then lead to our main ODE in finite time. We elaborate on the algorithm and give two more numerical examples. We also provide a few omitted proofs. Finally, we mention remaining open problems.

Preliminaries.
This section reviews the highlights of [15], starting with a definition of the original X queueing model, for which the ODE serves as an approximation.
2.1. The Original Queueing Model. The Markovian X model has two classes of customers, arriving according to independent Poisson processes with ratesλ 1 andλ 2 . There are two queues, one for each class, in which customers that are not routed to service immediately upon arrival wait to be served. Customers are served from each queue in order of arrival. Each class-i customer has limited patience, which is assumed to be exponentially distributed with rate θ i , i = 1, 2. If a customer does not enter service before he runs out of patience, then he abandons the queue. The abandonment keep the system stable for all arrival and service rates.
There are two service pools, with pool j having m j homogenous servers (or agents) working in parallel. This X model was introduced to study two large systems that are designed to operate independently under normal loads, but can help each other in face of unanticipated overloads. We assume that all servers are cross-trained, so that they can serve both classes. The service times depend on both the customer class i and the server type j, and are exponentially distributed; the mean service time for each class-i customer by each pool-j agent is 1/µ i,j . All service times, abandonment times and arrival processes are assumed to be mutually independent. The FQR-T control described below assigns customers to servers.
We assume that, at some unanticipated point of time, the arrival rates change, with at least one increasing. We further assume that the staffing cannot be changed (in the time scale under consideration) to respond to this unexpected change of arrival rates. Hence, the arrival processes change from Poisson with ratesλ 1 andλ 2 to Poisson processes with unknown (but fixed) rates λ 1 and λ 2 , whereλ i < m i /µ i,i , i = 1, 2 (normal loading), but λ i > µ i,i m i for at least one i (the unanticipated overload). Without loss of generality, we assume that pool 1 (and class-1) is the overloaded (or more overloaded) pool. The fluid model (ODE) is an approximation for the system performance after the overload has occurred, so that we start with the new arrival rate pair (λ 1 , λ 2 ).

The FQR-T Control for the Original Queueing
Model. The FQR-T control is based on two positive thresholds, k 1,2 and k 2,1 , and the two queueratio parameters, r 1,2 and r 2,1 . We define two queue-difference stochastic processesD 1,2 (t) ≡ Q 1 (t) − r 1,2 Q 2 (t) andD 2,1 ≡ r 2,1 Q 2 (t) − Q 1 (t). As long asD 1,2 (t) ≤ k 1,2 andD 2,1 (t) ≤ k 2,1 we consider the system to be normally loaded (i.e., not overloaded) so that no sharing is allowed. Hence, in that case, the two classes operate independently. Once one of these inequalities is violated, the system is considered to be overloaded, and sharing is initialized. For example, ifD 1,2 (t) > k 1,2 , then class 1 is judged to be overloaded and service-pool 2 is allowed to start helping queue 1. As soon as the first class-1 customer starts his service in pool 2, we drop the threshold k 1,2 , but keep the other threshold k 2,1 . Now, the sharing of customers is done as follows: If a type-2 server becomes available at time t, then it will take its next customer from the head of queue 1 ifD 1,2 (t) > 0. Otherwise, it will take its next customer from the head of queue 2. If at some time t after sharing has started queue 1 empties, orD 2,1 (t) = k 2,1 then the threshold k 1,2 is reinstated. The control works similarly if class 2 is overloaded, but with pool-1 servers helping queue 2, and with the threshold k 2,1 dropped once it is crossed.
In addition, we impose the condition of one-way sharing: we allow sharing in only one direction at any time. Thus, in the example above, where sharing is done with pool 2 helping class 1, we do not later allow pool 1 to help class 2 until there are no more pool-2 agents serving class-1 customers. Sharing is initiated with pool 1 helping class 2 whenD 2,1 (t) > k 2,1 and there are no pool-2 agents serving class-1 customers. And similarly in the other direction.
Let Q i (t) be the number of customers in the class-i queue at time t, and let Z i,j (t) be the number of class-i customers being served in pool j at time t, i, j = 1, 2. Let q i (t) and z i,j (t) be the fluid approximations of Q i (t) and Z i,j (t), respectively. With the assumptions on the X system and the FQR-T control, the six-dimensional stochastic process (Q i (t), Z i,j (t); i, j = 1, 2) describing the overloaded system becomes a continuous-time Markov chain (CTMC) (with stationary transition rates).
Once sharing is initialized, the control makes the overloaded X model operate as an overloaded N model, and keeps the two queues at approximately the target ratio, e.g., if queue 1 is being helped, then Q 1 (t) ≈ r 1,2 Q 2 (t). If sharing is done in the opposite direction, then r 2,1 Q 2 (t) ≈ Q 1 (t) for all t ≥ 0. That is substantiated by simulation experiments, some of which are reported in [15,16].
In addition to the thresholds k 1,2 and k 2,1 , discussed above, the model also includes shifting constants κ 1,2 and κ 2,1 . The shifting constants may be introduced after the threshold is dropped, because it may be dictated by the optimal ratio function in [15]. Let q * i and z * i,j , i, j = 1, 2 denote the fluid steady state values of q i (t) and z i,j (t). (We will show that a unique steady state, or stationary point, exists for the fluid approximation in §8 below.) If the optimal relation between the steady state fluid queues is q * 1 = r * 1,2 q * 2 +κ 1,2 for some κ 1,2 ∈ R, where r * 1,2 denotes the fluid optimal ratio, (assuming that pool 2 needs to help class 1), as is the case when the holding cost is separable and quadratic with non-zero constant and linear terms, then we use the shifted FQR-T control. Shifted FQR-T centers about κ 1,2 instead at about zero. For example, if class 1 is overloaded, then every server takes his new customer from the head of queue 1 ifD i,j (t) > κ 1,2 . Otherwise, it takes the new customer from the head of its own class queue. We call that control shifted FQR-T since it keeps the two queues at a fixed ratio, but shifted by the constant κ 1,2 . We can think of FQR-T as the special case of shifted FQR-T with κ 1,2 = 0. The beauty of the control is that, for large-scale service systems, FQR-T and shifted FQR-T tend to achieve their purpose; i.e., they keep the two queues approximately in fixed relation. In the stochastic system this means that the two-dimensional vector (Q 1 (t), Q 2 (t)) evolves approximately as a one-dimensional process. In the fluid model this approximation becomes exact; We no longer need to consider the three-dimensional process x(t) ≡ (q 1 (t), q 2 (t), z 1,2 (t)), since it is enough to consider z 1,2 (t) together with only one of the queues. The other queue is determined by the first via the statespace collapse (SSC) equation q 1 (t) = r i,j q 2 (t) + κ i,j , depending on which way the sharing is performed. In [17] SSC is shown to hold asymptotically in the MS-HT limit.
3. The Many-Server Heavy-Traffic Fluid Limit. In this section we briefly describe the convergence of the sequence of stochastic systems to the fluid limit, as established in [17]. Without loss of generality we assume that class 1 is overloaded, and receives help from service-pool 2. (Class 2 may also be overloaded, but less than class 1, so that pool 2 should be serving some class-1 customers.) 3.1. Many-Server Heavy-Traffic (MS-HT) Scaling. To develop the fluid limit in [17], we consider a sequence of X systems, indexed by n (denoted by superscript), with arrival rates and number of servers growing proportionally to n, i.e., with the service and abandonment rates held fixed. We then define the associated fluid-scaled stochastic processes For each system n, there are threshold k n 1,2 and k n 2,1 , scaled so that The first scaling by n is chosen to make the thresholds asymptotically negligible in MS-HT fluid scaling, so they detect overloads immediately when they occur. The second scaling by √ n is chosen to make the thresholds asymptotically infinite in MS-HT diffusion scaling, so that asymptotically the thresholds will not be exceeded under normal loading. It is significant that MS-HT scaling shows that we should be able to simultaneously satisfy both conflicting objectives in large systems.
There are also the shifting thresholds κ n i,j , arising from consideration of separable quadratic cost functions; see §2.2, but we do not specify their scale. If sharing is taking place, then at some time it was activated by sending the first class-1 customer to service pool 2. We thus need only consider κ n 1,2 and the weighted-difference processD n 1, Hence, we redefine the difference process. Let where κ ≡ κ 1,2 and r ≡ r * 1,2 . With the new definition in (3.4), we allow κ n to be of any order less than or equal to O(n); in particular, we assume that κ n /n → κ for 0 ≤ κ < ∞. There are two principle cases: κ = 0 and κ > 0. The first case produces FQR (after sharing has began); the second case produces shifted FQR (shifted by the constant κ n ).
With the new process D n in (3.4), we can apply the same FQR routing rule for both the FQR and shifted FQR cases: if D n (t) > 0, then every newly available agent (in either pool) takes his new customer from the head of the class-1 queue. If D n (t) ≤ 0, then every newly available agent takes his new customer from the head of his own queue.

3.2.
A Heuristic View of the AP. The AP is concerned with the system behavior when sharing is taking place; i.e., when some, but not all, of the pool 2 agents are serving class 1. That takes place when q 1 = rq 2 +κ. In that situation, it can be shown that the queue-difference process D n in (3.4) is an order O(1) process, without any spatial scaling, i.e., for each t, the sequence of unscaled random variables {D n (t) : n ≥ 1} turns out to be stochastically bounded (or tight) in R. That implies that D n operates in a time scale that is different from the other processes Q n i and Z n 1,2 , which are scaled by dividing by n in (3.2). With the MS-HT scaling in (3.1), in order for the two queues to change significantly (in a relative sense, which is captured by the scaling in (3.2)), there needs to be O(n) arrivals and departures from the queues. In contrast, the difference process D n can never go far from 0, because it has drift pointing towards 0 from both above and below. Thus, the difference process oscillates more and more rapidly about 0 as n increases. Thus, over short time intervals in which X n remains nearly unchanged for large n, the process D n moves rapidly in its state space, nearly achieving a local steady state. As n increases, the speed of the difference process increases, so that in the limit, it achieves a steady state instantaneously. That steady state is a local steady state, because it depends on x(t), the fluid limit x at time t.
To formalize this separation of time scales, we define a family of timeexpanded difference processes: for each n ≥ 1 and t ≥ 0, let Dividing s by n in (3.5) allows us to examine what is happening right after time t in the faster time scale. For each t, a different process D n t is defined. For every t ≥ 0 and s > 0, the time increment [t, t + s/n) becomes infinitesimal in the limit. A main result in [17] as n → ∞, where the limit D t ≡ {D t (s) : s ≥ 0} is the FTSP. For each n, the control depends on whether or not D n (t) > 0. In turn, the limiting ODE depends on the corresponding steady-state probability of the FTSP, 4. The ODE. We now specify the ODE, which is the main subject of this paper. We assume that class 1 is overloaded, even after receiving help from pool 2. Hence both pools are fully busy and some pool-2 agents are helping class 1, so that z 1,1 (t) = m 1 , z 2,1 (t) = 0 and z 1,2 (t) + z 2,2 (t) = m 2 . As a consequence, we only need consider z 1,2 among these four variables.
Recall that the shifting constant satisfies κ ≥ 0. We consider the restricted state space S ≡ [κ, ∞) × [0, ∞) × [0, m 2 ]. We thus avoid the transient region in which q 1 < rq 2 + κ with q 2 = 0, whereq 1 (t) > 0 andq 2 (t) < 0, but q 2 remains at 0 while q 1 increases to the shifting constant κ. The restricted state space, with q 1 ≥ κ is shown to be the space of the fluid limit of the system in [17]. We will also show in Theorem 5.1 below that the ODE cannot leave this restricted state space.
It is convenient to specify the conditions on the model parameters in terms of the steady-state formulas for the queues in isolation. For that purpose, let q a i be the length of fluid-queue i and let s a i be the amount of spare service capacity in service-pool i, in steady state, when there is no sharing, i = 1, 2. The quantities q a i and s a i are well known, since they are the steady state quantities of the fluid model for the Erlang-A model (M/M/m i + M ) with arrival-rate λ i , service-rate µ i,i and abandonment-rate θ i ; see Theorem 2.3 in [21], especially equation (2.19), and §5.1 in [15]. In particular, where (x) + ≡ max{x, 0}. It is easy to see that q a i s a i = 0, i = 1, 2. We thus make the following assumption, which is assumed to hold henceforth.
Assumption A.
(I) The model parameters satisfy We now explain these assumptions. Clearly, a sufficient condition for both pools to be overloaded is to have s a 1 = s a 2 = 0, i.e., to have no spare service capacity in either pool in their individual steady states. However, if s a 2 > 0, both pools can still be overloaded, provided that enough class-1 fluid is processed in pool 2. To have the solution be eventually in S, we require that θ 1 (q a 1 − κ) ≥ µ 1,2 s a 2 . This condition ensures that service pool 2 is also full of fluid when sharing is taking place, i.e., z 1,2 (t) + z 2,2 (t) = m 2 for all t ≥ 0 (assuming that pool 2 is full at time 0). To see why, note that when service-pool 2 has spare service capacity (s a 2 > 0), sharing will be activated if q a 1 > κ, because q a 2 = 0. Now, the maximum amount of class-1 fluid that pool 2 can process, while still processing all of the class-2 fluid (so that q 2 is kept at zero), is µ 1,2 s a 2 . hence, µ 1,2 s a 2 is the minimal amount of class-1 fluid that should flow to pool 2. On the other hand, θ 1 q a 1 = λ 1 − µ 1,1 m 1 is equal to the "extra" class-1 fluid input, i.e., all the class-1 fluid that pool 1 cannot process. Some of this "extra" class-1 fluid might abandon (if q 1 > 0). The minimal amount of class-1 fluid that abandons is θ 1 κ (but κ can be equal to zero).
We thus require that all the class-1 fluid, that is not served in pool 1, minus the minimal amount of class-1 fluid that abandons, is larger than µ 1,2 s a 2 . With this requirement, pool 2 is assured to be full, assuming that it is initialized full. (If pool 2 is not initialized full, then it will fill up after some finite time period; see the appendix.) Remark 5.1. (class 1 need not be more overloaded than class 2) In this paper we are interested in analyzing the ODE (4.2) as given. Hence, in Assumption A we do not assume that class 1 is more overloaded than class 2; i.e., we do not require that q a 1 − κ ≥ rq a 2 . This extra assumption is not needed for our results for the specified ODE. In contrast, this assumption is needed in order to show that the ODE holds as the fluid limit, with class 1 receiving help; see Assumption 1 in [17].
We exploit Assumption A to show that the boundaries of S in R 3 play no role.
We give the proof in §8.4 after the necessary tools have been developed. Our main result establishes the existence of a unique solution.
Theorem 5.2. (existence and uniqueness) For any w 0 ∈ S, there exists a unique function x : [0, ∞) → S such that, (i) for all t ≥ 0, there exist δ(t) > 0 such that x is right-differentiable at t, differentiable on (t, t + δ(t)) and satisfies the IVP (4.3) based on the ODE (4.1) over [t, t + δ(t)) with initial value x(t), and (ii) x is continuous and differentiable almost everywhere.
Theorem 5.2 has two parts: First, there is (i) establishing the local existence and uniqueness of a conventional differentiable solution on each interval [t, t + δ(t)), for which it suffices to consider a single t, e.g., t = 0. Second, there is (ii) justifying an overall continuous solution.
We prove Theorem 5.2 in the next two sections. The proof is tied to the characterization of π 1,2 in (4.2) and (3.7), and thus the FTSP D t . We need to determine conditions for the FTSP D t to be positive recurrent, so that the AP holds, and then calculate its steady-state distribution in order to find π 1,2 . Moreover, we need to establish topological properties of the function π 1,2 , such as continuity and differentiability. We discuss the FTSP D t next.
6. The Fast-Time-Scale Process. Recall that the FTSP D t is the limit of D n t without any scaling (see (3.6)), where D n t is the time-expanded difference process defined in (3.5) associated with the queue-difference stochastic process D n ≡ (Q n 1 − κ n ) − rQ n 2 in (3.4). Since there is no scaling of space, the state space for the FTSP D t is the countable lattice {±j ± kr : j, k ∈ Z} in R. To see this, first observe from (3.4) that D n has state space {±j ± kr − κ n : j, k ∈ Z}. Next, because of the subtraction in (3.5), D n t has state space {±j ± kr : j, k ∈ Z}. Finally, because of the convergence in (3.6), the FTSP D t has this same state space.
6.1. The Fast-Time-Scale CTMC. We fix a time t and assume that we are given the value x(t) ≡ (q 1 (t), q 2 (t), z 1,2 (t)). In order to simplify the analysis we assume that r is rational. That clearly is without any practical loss of generality. Specifically, we assume that r = j/k for some positive integers j and k without any common factors. We then multiply the process by k, so that all transitions can be expressed as ±j or ±k in the state space Z. In that case, the FTSP D t ≡ {D t (s) : s ≥ 0} becomes a CTMC.
Let λ (j) ) and µ (k) + (m, x(t)) be the transition rates of the FTSP D t for transitions of +j, +k, −j and −k, respectively, when D t (s) = m > 0. Similarly, we define the transitions when These rates are the limits of the rates of D n t as n → ∞ withX n (t) ⇒ x(t). First, for D t (s) = m ∈ (−∞, 0], the upward rates are corresponding, first, to a class-1 arrival and, second, to a departure from the class-2 queue, caused by a type-2 agent service completion (of either customer type) or by a class-2 customer abandonment. Similarly, the downward rates are corresponding, first, to a departure from the class-1 customer queue, caused by a class-1 agent service completion or by a class-1 customer abandonment, and, second, to a class-2 arrival. Next, for D t (s) = m ∈ (0, ∞), we have upward rates corresponding, first, to a class-1 arrival and, second, to a departure from the class-2 customer queue caused by a class-2 customer abandonment. The downward rates are corresponding, first, to a departure from the class-1 customer queue, caused by (i) a type-1 agent service completion, (ii) a type-2 agent service completion (of either customer type), or (iii) by a class-1 customer abandonment and, second, to a class-2 arrival.

6.2.
Representing the FTSP D t as a QBD. Further analysis is simplified by exploiting matrix geometric methods, as in [10]. In particular, we represent the integer-valued CTMC D t ≡ {D t (s) : s ≥ 0} just constructed as a homogeneous continuous-time quasi-birth-and-death (QBD) process, as in Definition 1.3.1 and §6.4 of [10]. In passing, note that the special case r = 1 is especially tractable, because then the QBD process reduces to an ordinary birth-and-death (BD) process.
To represent D t as a QBD process, we must re-order the states appropriately. We order the states so that the infinitesimal generator matrix Q can be written in block-tridiagonal form, as in Definition 1.3.1 and (6.19) of [10] (imitating the shape of a generator matrix of a BD process). In particular, we write where the four component submatrices B, A 0 , A 1 and A 2 are all 2m × 2m submatrices for m ≡ max {j, k}. In particular, These 2m × 2m matrices B, A 0 , A 1 and A 2 in turn can be written in block-triangular form composed of four m × m submatrices, i.e.,

. (All matrices are also functions of x(t).)
To achieve this representation, we need to re-order the states into levels. The main idea is to represent transitions of D t above and below 0 within common blocks. Let L(n) denote level n, n = 0, 1, 2, . . . We assign original states φ(n) to positive integers n according to the mapping: Then we order the states in levels as follows With this representation, the generator-matrix Q can be written in the form (6.5) above, where A 1 groups all the transitions within a level, A 0 groups the transitions from level L(n) to level L(n+1) and A 2 groups all transitions from level L(n) to level L(n − 1). Matrix B groups the transitions within the boundary level L(0), and is thus different than A 1 .
To illustrate, consider an example with r = 0.8, so that we can choose j = 4 and k = 5, yielding m = 5. The states are ordered in levels as follows Then the submatrices B µ , B λ , A + i and A − i , which form the block matrices B and A i , i = 0, 1, 2, have the form in (6.12) below, where (We solve a full numerical example with these matrices in §11.3.) Henceforth, we refer to D t interchangeably as the QBD or the FTSP.
6.3. Positive Recurrence. We show that positive recurrence depends only on the constant drift rates in the two regions, as one would expect. Let δ + and δ − be the drift in the positive and negative region, respectively; i.e., let Proof. We employ the theory in §7 of [10], modified for the continuoustime QBD. We first construct the aggregate matrices We then observe that the aggregate matrix A is reducible, so we need to consider the component matrices A + and A − , which both are irreducible CTMC infinitesimal generators in their own right. Let ν + and ν − be the unique stationary probability vectors of A + and A − , respectively, e.g., with ν + A + = 0 and ν + 1 = 1. The theory concludes that our QBD is positive recurrent if and only if In our application it is easy to see that both ν + and ν − are the uniform probability vector, attaching probability 1/m to each of the m states, from which the conclusion follows directly.
The alternative cases are simplified by the following relation: Hence there are only two cases in which the drift does not point inward: (i) δ + (x(t)) ≥ 0 and δ − (x(t)) > 0, (ii) δ − (x(t)) ≤ 0 and δ + (x(t)) < 0. In both cases the behavior is unambiguous: In case (i), clearly π 1,2 (x(t)) = 1; in case (ii), clearly π 1,2 (x(t)) = 0. 6.4. Computing π 1,2 . When the QBD is positive recurrent, the stationary vector of the QBD can be expressed as α ≡ {α n : n ≥ 0} ≡ {α n,j : n ≥ 0, 1 ≤ j ≤ m}, where α n ≡ (α + n , α − n ) for each n, with α + n and α − n both being 1 × m vectors. Then the desired probability π 1,2 can be expressed as (6.14) where 1 denotes a column vector with all entries 1 of the right dimension (here m×1), while 1 + represents a 2m×1 column vector, with m 1 ′ s followed by m 0 ′ s. By Theorem 6.4.1 and Lemma 6.4.3 of [10], the steady-state distribution has the matrix-geometric form where R is the 2m × 2m rate matrix, which is the minimal nonnegative solutions to the quadratic matrix equation A 0 + RA 1 + R 2 A 2 = 0, and can be found efficiently by existing algorithms, as in [10] (see §11 below). Since the matrices A 0 , A 1 and A 2 have the block-diagonal form in (6.6), so does R, with submatrices R + and R − . Since the spectral radius of the rate matrix R is strictly less than 1 (Corollary 6.2.4 of [10]), the sum of powers of R is finite, yielding Also, by Lemma 6.3.1 of [10], the boundary probability vector α 0 in (6.15) is the unique solution to the system Finally, given the above, and using (6.14), we see that the desired quantity π 1,2 can be represented as For further analysis, it is convenient to have alternative representations for π 1,2 (x). Let the vector 1 have the appropriate dimension in (6.19) below.
Proof. (a) When r = 1, the FTSP D t ≡ D t (x) evolves as an M/M/1 queue in each of the regions D t (s) > 0 and D t ≤ 0. Thus, we can look at the system at the successive times at which D t transitions from state 0 to state 1, and then again from state 1 to state 0. That construction produces an alternating renewal process of occupation times in each region, where these occupation times are distributed as the busy periods of the corresponding M/M/1 queues. Hence, π 1,2 (x) can be expressed as where T + (x) is the busy period of the M/M/1 queue in the upper region, while T − (x) is the busy period of the M/M/1 queue in the lower region. By the definition of A, these mean busy periods are finite in each region. In particular, are the constant drift rates up (away from the boundary) and down (toward the boundary) in the upper region in (6.3) and (6.4), depending on state x, while λ − (x) and µ − (x) are the constant drift rates down (away from the boundary) and up (toward the boundary) in the lower region in (6.1) and (6.2); e.g., (b) We first observe that we can reason as in the case r = 1, using a regenerative argument. We can let the regeneration times be successive transitions from one specific QBD state in level 0 with D t ≤ 0 to a specific state in level 1 where D t > 0. The intervals between successive transitions will be i.i.d. random variables with finite mean. Hence, we can represent π 1,2 (x) just as in (6.20), but where now T + (x) is the total occupation time in the upper region with D t (s) > 0 during a regeneration cycle, while T − (x) is the total occupation time in the lower region with D t (s) ≤ 0 during a regeneration cycle. Each of these occupation times can be broken up into first passage times. For example, T + (x) is the sum of first passage times from some state at level 0 to some other state in level 1 where D t (s) > 0. The regenerative cycle will end when the starting and ending states within levels 0 and 1 are the designated pair associated with the specified regeneration time. The successive pairs (i,j) of starting and ending states within the levels 0 and 1 evolve according to a positive-recurrent finite-state discrete-time Markov chain.
With the QBD representation, we can determine when the FTSP D t is positive recurrent, for a given x(t), using (6.10), and then numerically calculate π 1,2 . That allows us to numerically solve the ODE (4.1) in §11. We will also use the representations (6.17), (6.18), (6.19) and other QBD properties to deduce topological properties of π 1,2 .
7. Existence and Uniqueness of Solutions. This section is devoted to proving Theorem 5.2. For the local existence and uniqueness in Theorem 5.2 (i), we will show that the function Ψ in (4.2) is locally Lipschitz continuous in Theorem 7.1 below. That allows us to apply the classical Picard-Lindelöf theorem to deduce the desired existence and uniqueness of solutions to the IVP (4.3); see Theorem 2.2 of Teschl [19] or Theorem 3.1 in [9]. Afterwards, in §7.2 we establish the global properties in Theorem 5.2 (ii).
The boundary subset S b is a hyperplane in the state space S, and is therefore a closed subset. It is the subset of S in which SSC and the AP are taking place (in fluid scale). In S b the function π 1,2 can assume its full range of values, 0 ≤ π 1,2 (x) ≤ 1.
The region S + above the boundary is an open subset of S. For all x ∈ S + , π 1,2 (x) = 1. The region S − below the boundary is also an open subset of S. For all x ∈ S − , π 1,2 (x) = 0. It is important to keep in mind that, in order for S − to be a proper subspace of S, both service pools must be constantly full (in the fluid limit). Thus, if x ∈ S − , then z 1,1 = m 1 and z 1,2 + z 2,2 = m 2 (but q 1 and q 2 are allowed to be equal to zero).
It is immediate that the function Ψ in (4.2) is Lipschitz continuous on S + and S − , because π 1,2 (x) = 1 when x ∈ S + , and π 1,2 (x) = 0 when x ∈ S − , so that Ψ is linear in each region. However, Ψ is not linear on S b . To analyze Ψ on S b , we exploit properties of the QBD introduced in §6. We partition S b into three subsets, depending on the drift rates in (6.9). Let A be the set of all x ∈ S b for which the QBD is positive recurrent, as given in (6.10); i.e., let Let the other two subsets be By the relation (6.13), there are no other alternatives; i.e., From the continuity of the QBD drift-rates in (6.9), if follows that A is an open and connected subset of S b . Hence, A can be regarded as an open connected subset of R 2 Just as for the open subsets S + and S−, if the initial value is in A, then the ODE will remain within A over some initial subinterval. In contrast, the situation is more complicated in A + and A − .
There is potential movement out of the region S b only from the sets A + and A − . To understand what can happen, let d(x(t)) ≡ q 1 (t) − rq 2 (t) and d ′ (x(t)) ≡q 1 (t) − rq 2 (t), from (4.2). On A + and A − , the possibilities can be determined from the following lemma.
Proof. Substitute the appropriate values of π 1,2 (x(t)) into (4.2) and compute δ ± (x) from (6.1)-(6.4), recalling that r ≡ j/k. We next separate equality from strict inequality for the weak inequalities in Lemma 7.1. For that purpose, we decompose the sets A + and A − by letting Suppose that a solution exists for the ODE over a sufficiently small interval starting at Proof. We only treat the first two cases, because the reasoning for the last two is the same. If x(0) ∈ A + + , then d ′ (x(0)) > 0 by Lemma 7.1, which implies the result. If x(0) ∈ A + 0 , then d ′ (x(0)) = 0 by Lemma 7.1. To see why we cannot have x(0+) ∈ A − ∪ S − , note that then π 1,2 (x) would jump from 1 to 0, which would cause a jump in d ′ (x) because of Lemma 7.1 and the inequality in (6.13).
We now are ready to establish local Lipschitz continuity.
Note that all of S is covered by the five cases in Theorem 7.1. Also note that, in each case, when we conclude that Ψ is Lipschitz continuous on Ω 1 within Ω 2 , if a solution exists starting at some point in Ω 1 , then it will necessarily remain in Ω 2 for a short interval, by the reasoning above. We postpone the relatively long proof until §12.

Global Existence and Uniqueness.
This section is devoted to completing the proof of Theorem 5.2 by establishing global existence and uniqueness. We first observe that, in general, one overall differentiable solution to the ODE over [0, ∞) may not exist. From either S − or S + , the solution x can hit S b , i.e., one of the sets A, A − and A + , and have drifts that are inconsistent with the drifts in the new destination set. For example, in S + we necessarily have π 1,2 (x) = 1. However, in general there is nothing preventing ). The probability π 1,2 (x(t)) jumps instantaneously from 1 to some value strictly between 0 and 1 when A is hit. A numerical example is given in the appendix. To treat that case, We can start a new ODE at this hitting time of A.
We first show that the possible values of x are contained in a compact subset of S, provided that the initial values of the queue lengths are constrained. That is accomplished by proving that a solution to the IVP (4.3) is bounded. We use the notation: a ∨ b ≡ max{a, b}. In particular, the following upper bounds for the fluid queues hold: Proof. Since 0 ≤ z 1,2 ≤ m 2 and q i ≥ 0 in S, we only need to establish (7.5). To do so, it suffices to consider the bounding function describing the queue-length process of each queue in a modified system with no service processes, so that all the fluid output is due to abandonment, which produces a simple one-dimensional ODE for each queue; for the remaining details, see §D in the appendix.
Proof of Theorem 5.2 (ii). It follows from Theorem 5.2 (i) established above, and Theorems 7.1 and 7.2, that any solution x on [0, δ) can be extended to an interval [0, δ ′ ), δ ′ > δ (even δ ′ = ∞), with the solution {x(t) : t ∈ [0, δ ′ )} again being unique, provided that that the solution x makes no transitions from S − S b to A, causing a discontinuity in π 1,2 (x) and thus Ψ in (4.2). (See Theorem 3.3 in [9] and its proof for supporting details.) Moreover, the solution in S + or S − has a left limit at the time it hits A. The left limit exists because, by Theorem 7.2, the solution is bounded, and because the derivative in either S + or S − is bounded, by (4.2). At each such hitting time, a new ODE is constructed starting in A. That ensures the overall continuity of x. In general, there can be accumulation points of such hitting times of the set A from S−S b . However, any such accumulation point t must be in either A + or A − . That is so, because there then are sequences {t i n : n ≥ 1}, i = 1, 2 with x(t 1 n ) ∈ S − S b and x(t 2 n ) ∈ A for all n with t i n ↑ t and x(t i n ) → x(t) ∈ S b as n → ∞ for i = 1, 2. Finally, by Theorem 7.1, the function Ψ is locally Lipschitz continuous at each point in A + ∪ A − . Hence, the solution x must actually be differentiable at each of these accumulation times of hitting times. As a consequence, x is continuous and differentiable almost everywhere throughout [0, ∞).
In the proof of Theorem 5.2 (ii), just completed, we have also established the following result.  . If x(t) = x * for all t, then we say that the fluid solution is stationary, or in steady state.
8.1. Characterization of the Stationary Point. By definition, a stationary point x * ∈ S satisfies Ψ(x * ) = 0. From (4.2), we see that this gives a system of three equations with three unknowns, namely, q * 1 , q * 2 and z * 1,2 . The apparent fourth variable π * 1,2 ≡ π 1,2 (x * ) is a function of the other three variables and its value is determined by x * . In principle, the three equations in Ψ(x) = 0 can be solved directly to find all the roots of Ψ. However, π * 1,2 is a complicated function of x * having the complicated closed-form expression in (6.14) and (6.17).
Theorem 8.1 below states that, if there exists a stationary point for the fluid ODE (4.2), then this point is unique, and must have the specified form. The uniqueness of x * is proved by treating π * 1,2 as a fourth variable, and adding a fourth equation to the three equations Ψ(x) = 0. However, it does not prove that a stationary point exists. In general, the solution π * 1,2 we get from the system of four equations may not equal to π 1,2 (x * ), for the function π 1,2 defined in (3.7). The existence of a stationary point is proved in the next section.
The proof of existence is immediate from the proof of uniqueness when π 1,2 (x * ) is known in advance to be 0 or 1, with the value determined. That occurs everywhere except the region A; it occurs in the two regions S + and S − , but it also occurs in S b − A. Since the QBD is not positive recurrent in S b − A, it follows that π 1,2 (x * ) can only assume one of the values, 0 or 1, achieving the same value as in the neighboring region S + or S − . (We omit detailed demonstration.) But we will have to work harder in A.
Consider the following three equations with the three unknowns: z, q 1 (z) and q 2 (z). (here q 1 and q 2 are treated as functions of the variable z, not to be confused with the fluid solution which is a function of the time argument t.) Notice that q 1 (z) is decreasing with z, whereas q 2 (z) is increasing with z. Thus, there exists a unique solution to these three equations, which has z as in (8.1). We can recover x * from the solution to (8.4), and by doing so show that x * is unique and is always in one of the three regions S − , S + or S b (so that x * ∈ S). Let (q 1 (z), q 2 (z), z) be the unique solution to (8.4). First assume that z > m 2 , which implies that q 2 (z) > 0, and, by the third equation, q 1 (z) > κ. By replacing z with m 2 , q 1 (·) is increased and q 2 (·) is decreased (but is still positive), so that q 1 (m 2 ) − rq 2 (m 2 ) > κ (and, trivially, q 1 (m 2 ) > κ, q 2 (m 2 ) > 0). This implies that x * ≡ (q 1 (m 2 ), q 2 (m 2 ), m 2 ) ∈ S + and, if it is indeed a solution to Ψ(x) = 0, then x * is the unique stationary point for the ODE. Now assume that the unique solution to (8.4) has z < 0. By replacing z with 0 we have q 1 (0) < q 1 (z) and q 2 (0) > q 2 (z), which imply that q 1 (0) − rq 2 (0) < κ. Now, since q 1 (0) = q a 1 we have that q 1 (0) ≥ κ by Assumption A. This implies that q 1 (z) > κ, which further implies that rq 2 (z) = q 1 (z) − κ > 0, so that rq 2 (0) > rq 2 (z) > 0. Taking x * ≡ (q 1 (0), q 2 (0), 0), we see that x * ∈ S − , and if x * is indeed a solution to Ψ(x) = 0, then x * is the unique stationary point for the ODE.
Finally, assume that the solution x(z) ≡ (q 1 (z), q 2 (z), z) to (8.4) has 0 ≤ z ≤ m 2 . To conclude that x(z) is in S b we need to show that q(z), q 2 (z) ≥ 0, so that q * 1 = q 1 (z) and q * 2 = q 2 (z) are legitimate queue-length solutions. We now show that is the case under Assumption A.
An immediate consequence of the proof of Theorem 8.1 is that, in order to find the candidate stationary point x * , one has to solve the three equations in (8.4). The next corollary summarizes the values x * may take, depending on its region; the proof appears in the appendix.
1. If x * ∈ S b , then, for z defined in (8.1), If x * ∈ S + , as in (ii), then the system does not have enough service capacity to keep the weighted difference between the two queues at κ, even when all agents are working with class 1. In this case, the only output from queue 2 is due to abandonment, since no class-2 fluid is being served (in steady state). Queue 2 is then equivalent to the fluid approximation for an M/M/∞ system with service rate θ 2 and arrival rate λ 2 . On the other hand, queue 1 is equivalent to an overloaded inverted-V model: a system in which one class, having one queue, is served by two different service pools.
The next corollary gives necessary and sufficient conditions for x * to be in each region. It shows that the region of x * can be determined from rate considerations alone. We give the proof in the appendix.
x * ∈ A if and only if both inequalities are strict.
Remark 8.1. (most likely region in applications) It follows from Corollary 8.2 that, in applications, A is the most likely region for the stationary point when the system is overloaded, provided that the arrival rates are about 10 − 50% larger than planned during an overload incident. Typically, a much higher overload is needed in order for the stationary point to be in S + . As an example, consider the canonical example from [15]: There are 100 servers in each pool, serving their own class at rates µ 1,1 = µ 2,2 = 1. Type-2 servers serve class-1 customers at rate µ 1,2 = 0.8. Also, θ 1 = θ 2 = 0.3, r = 0.8 and κ = 0. Suppose that class 2 is not overloaded with λ 2 = 90. Then, for the stationary point to be in S + , we need to have λ 1 > µ 1,1 m 1 + µ 1,2 m 2 + θ 1 rλ 2 /θ 2 = 252, i.e., the class-1 arrival rate is 252% larger than the total service rate of pool 1. If λ 2 > 90, especially if pool 2 is also overloaded, then λ 1 needs to be even larger than that.

Existence of a Stationary Point.
We have just established uniqueness of the stationary point in S, and characterized it. In the process, we have also established existence in S − A. Now we will establish existence of the stationary point in A. First, we calculate the drift rates at x * ∈ A. Proof. Substitute x * in Corollary 8.1 (i) into (6.9), using (6.1)-(6.4).
We now are ready to prove existence.
Proof. We will prove that there must exist at least one stationary point. Given that result, by Theorem 8.1 and Corollary 8.1, there must be exactly one fixed point and that must be the x * given there. To establish existence, we will apply the Brouwer fixed point theorem. It concludes that a continuous function mapping a convex compact subset of Euclidean space R k into itself has at least one fixed point. We will let our domain be the set Choose η sufficiently small that x * ∈ C(η); that is easily ensured by Lemma 8.1. Since the rates in (6.1)-(6.4) and the drift in (6.9) are linear functions of x, we see that C(η) is a convex subset of A for each η > 0. Since the inequalities in (8.10) are weak, the set is closed. The intersection with B guarantees that the set C(η) is also bounded. Hence, C(η) is compact. By Theorem 5.2, for any x(0) ∈ C(η), there exists a unique solution to the ODE over [0, δ] for some positive δ. Hence, for any t with 0 < t < δ, the map from x(0) to x(t) is continuous; see §2.4 of [19]. Let x * L ≡ (q * 1,L , q * 2,L , z * 1,2,L ) and and z 1,2,t ≡ z 1,2,t ∨ z * 1,2,L ∧ z * 1,2,U , for i = 1, 2. We can choose η > 0 and ǫ > 0 sufficiently small so that, first, x * ∈ C(η) and, second, that x i,t ∈ C(η) for each x(0) ∈ C(η). Hence, the pair (C(η), φ t ) satisfies the conditions for the Brouwer fixed point theorem. Hence, there exists x(0) ∈ C(η) such that x(t) = x(0). Now let {t n : n ≥ 1} be a sequence of time points decreasing toward 0. We can apply the argument above to deduce that, for each n, there exists x n (0) in C(η) such that x n (t n ) = x n (0). However, from the ODE, we have the relation |x(t) − x(0) − Ψ(x(0))t| ≤ M t 2 for all sufficiently small t. Since {x n (0) : n ≥ 1} is bounded, there exists a convergent subsequence. Let x(0) be the limit of that convergent subsequence. For that limit, we necessarily have Ψ(x(0)) = 0. Hence, that x(0) must be a stationary point for the ODE. By Theorem 8.1, we must have x(0) = x * . 8.3. Global Asymptotic Stability. Having a unique stationary point does not imply that a fluid solution necessarily converges to that point as t → ∞. It does not even guarantee that a solution to the IVP (4.3) is asymptotically stable in the sense that, if x(0) − x * < ǫ, then x(t) → x * as t → ∞, no matter how small ǫ is. In fact, there is not even a guarantee that x(t) will remain in the ǫ-neighborhood of x * for all t ≥ 0. We will establish all of these properties in Theorem 8.3 below by showing that x * in §8.1 is globally asymptotically stable, as defined below: Definition 8.2. (global asymptotic stability) A point x * is said to be globally asymptotically stable if it is a stationary point and if, for any initial state x(0) and any ǫ > 0, there exists a time T ≡ T (x(0), ǫ) ≥ 0 such that x(t) − x * < ǫ for all t ≥ T .
Global asymptotic stability goes beyond simple convergence by also requiring that the limit be a stationary point. The proof of Theorem 8.3 relies on Lyapunov stability theory for deterministic dynamical systems; see Chapter 4 of Khalil [9]. Let E be an open and connected subset of R n containing the origin. We use standard vector notation to denote the inner product of vectors a, b ∈ R n , i.e., a·b = n i=1 a i b i . Definition 8.3. (Lie derivative) For a continuously differentiable function V : E → R, and a function Ψ : E → R n , the Lie derivative of V along Ψ is defined byV where ∇V ≡ ( ∂V ∂x 1 , . . . , ∂V ∂xn ) is the gradient of V .
Definition 8.4. (Lyapunov-function candidate) A continuously differentiable function V : E → R is a Lyapunov-function candidate if: In proving Theorem 8.3 we use the following theorem, which is Theorem 4.2 pg. 124 in [9]: Notice that, under the conditions of Theorem 8.4, the Lyapunov-function candidate V provides a form of monotonicity: We necessarily have V (0) = 0 and V (x(t)) strictly decreasing in t for x(t) = 0. To elaborate, we introduce the notion of a V -ball. We say that β V (α) is the α V -ball with center at x * and radius α if If x(t 0 ) ∈ β V (α) for some α ≥ 0 and t 0 ≥ 0, then x(t) ∈ β V (α) for all t ≥ t 0 .
Proof of Theorem 8.3. Theorem 8.4 applies directly only within one region, starting at a point in S + , S − , A, A − or A + . However, we will show that the same Lyapunov function V can be used in all regions, leading to global decrease of V as x * is being approached. Let x be the unique solution to (4.3). Let x * ≡ (q * 1 , q * 2 , z * 1,2 ) be the stationary point for the system (4.1). Without loss of generality, we perform a change of variables and define a new system whose unique stationary point is x = 0. To this end, let y = x − x * so thatẏ =ẋ = Ψ(x). Hence, Ψ(x) = Ψ(y + x * ) ≡ g(y) and we have that g(0) = Ψ(0 + x * ) = Ψ(x * ) = 0. That is, if x * is a stationary point for the original systemẋ = Ψ(x), then the stationary point for the new system,ẏ = g(y), is y * = 0. We distinguish between two cases: (i) µ 1,2 > µ 2,2 and (ii) µ 1,2 ≤ µ 2,2 .
By Theorem 8.4, y * = 0 is globally asymptotically stable for the modified system g(y). Hence, x * is globally asymptotically stable for the original system Ψ(x). That is, for every initial value x(0) we have that x(t) → x * .

Staying in S.
We also use the Lyapunov argument to prove Theorem 5.1, i.e., show that the solution to the ODE can never leave S.
We finally consider the remaining more challenging subcase:q 1 (t) < 0 andq 2 (t) < 0. We will show that this subcase is not possible. To that end, we further divide this case into three subcases: However, x(t) cannot be in A − , since then π 1,2 (x(t)) = 0, which implies that q 1 (t) is nondecreasing (plug π 1,2 (x(t)) = 0 and q 1 (t) = κ into the ODE (4.2)). Moreover, x(t) cannot be in A + , since then π 1,2 (x(t)) = 1, which implies that q 2 (t) is strictly increasing. Now assume the remaining possibility, x(t) ∈ A, and recall that Ψ is Lipschitz continuous in A, so that the Lyapunov argument holds over [t, t + η), for some η > 0. Specifically, the Lyapunov function V is monotone increasing in x(t), because x * > 0. (The inequality holds componentwise.) If µ 1,2 > µ 2,2 , then we take the Lyapunov function V 1 (x(t)) = q 1 (t) + q 2 (t). The monotonicity of V 1 at x(t) implies that at least one of the queues must be increasing, which contradicts the assumption that the derivative of both queues is negative at t. If µ 1,2 ≤ µ 2,2 , then we take the Lyapunov function V 2 (x(t)) = Cq 1 (t) + q 2 (t) + (C − 1)z 1,2 (t). We then choose C = 1 + ǫ with ǫ small enough, such thatV 2 (x(t)) < 0 (assuming the derivatives of both queues are strictly negative at t). Once again, this contradicts the positive monotonicity of V at x(t). This concludes the proof.

Exponential Stability.
Definition 9.1. (exponential stability) A stationary point x * is said to be exponentially stable if there exist two real constants ϑ, β > 0 such that for all t ≥ 0 and for all x(0), where · is a norm on R n .
Theorem 9.1. (exponential stability of the origin) Suppose that all the conditions of Theorem 8.4 are satisfied. In addition, assume that there exist positive constants K 1 , K 2 , K 3 and p such that Then the origin is exponentially stable, and We use the L 1 norm: Theorem 9.2. (exponential stability of x * ) Each x * in S is exponentially stable.
Proof. As in the proof of Theorem 8.3, Theorem 9.2 applies directly only within one region, starting at a point in S + , S − , A, A − or A + . However, again, the same Lyapunov function V can be used in all regions.
(ii) We use the Lyapunov function 10. Conditions for State-Space Collapse. In this section we give ways of verifying that x lies entirely in A, given that x(0) and x * are both in A. In the appendix we provide conditions for the solution to eventually reach A after an initial transient. The results here are intended to apply after this initial transient period has concluded.
For x * ∈ A, we will now show that there exist α > 0 and T ≡ T (α), such that global SSC can be inferred once x(T ) − x * < α. We exploit the drift rates at stationarity, defined by δ * + ≡ δ + (x * ) and δ * − ≡ δ − (x * ). It follows from the expressions in (6.9) that (10.5) δ Thus, if 0 < z * 1,2 < m 2 , then the positive recurrence condition (6.10) holds at the stationary point x * . (This agrees with (8.3) which has 0 < π * 1,2 < 1 if and only if 0 < z * 1,2 < m 2 .) In the next theorem we give explicit expressions for α. For reasonable rates, such as will hold in applications, α is quite large. In fact, in the numerical example considered in §11.3 we show that, typically in applications, α is so large, that we can infer that x lies entirely in A without even solving the IVP; i.e., x(0) ∈ β V (α).

In both cases, if there exists
Proof. We use β V (α), the α V -ball with center at x * and radius α, in (8.12). To find a proper α for the V -ball β V (α), we once again use the conditions (10.3) and (10.4). We first show how to find α for the case µ 2,2 = Bµ 1,2 for some B ≥ 1, i.e., when µ 1,2 ≤ µ 2,2 . Recall (proof of Theorem 8.3) that in this case, V 2 (x) = Cx 1 + x 2 + (C − 1)x 3 is a Lyapunov function for any C > B. Also, the Lyapunov function was defined for the modified system in which the origin was the stationary point.
Let x * = (q * 1 , q * 2 , z * 1,2 ) be the stationary point in A. First assume that, at some time T , V 2 (x(T )) = ǫ 1 , i.e., Cq 1 (T ) + q 2 (T ) + (C − 1)z 1,2 (T ) = ǫ 1 . If x(t) ∈ β V 2 (ǫ 1 ) for all t > T , then it must hold that To make sure δ + (x(t)) < 0, we use (10.3), reorganizing the terms. We need to have By (10.6), the above inequality holds if Plugging in the expressions for q * 1 , q * 2 and z * 1,2 , we see that we need to find an ǫ 1 > 0 such that We can take C as large as needed, so that the only term that matters on the LHS is rθ 2 ǫ 1 . Hence, we need to have Similarly, to make sure that δ − (x(t)) > 0, we use (10.4), reorganizing the terms. We need to have Using (10.6) again (with a different ǫ 2 ), we see that it suffices to show that Once again, plugging in the values of q * 1 , q * 2 and z * 1,2 , and taking C as large as needed, we can choose ǫ 2 > 0 such that Hence, we can take α as in (i). For the second case, when µ 1,2 > µ 2,2 , we use the Lyapunov function V 1 (x) = x 1 + x 2 . Using similar reasoning as above, we get Hence, in this case we can take α in (ii).

A Numerical Algorithm to Solve the IVP.
11.1. Computing π 1,2 (x) at a point x. The QBD structure in §6.2 allows us to use established efficient numerical algorithms from [10] to solve for the steady state of the QBD to compute π 1,2 (x), for any given x ≡ x(t) ∈ A. We start by computing the rate matrix R ≡ R(x). (To simplify notation, we drop the argument x, with the understanding that all matrices, are functions of x.) By Proposition 6.4.2 of [10], R is related to matrices G and U via In addition, the matrices G and R are the minimal nonnegative solutions to the quadratic matrix equations Hence, if can compute the matrix G, then the rate matrix R can be found via (11.1). Once R is known, we use (6.16) to compute α 0 . With α 0 and R in hand, π 1,2 (x) is easily computed via (6.17). It remains to compute the matrix G. We use the logarithmic reduction (LR) algorithm in §8.4 of [10], modified to the continuous case, as in §8.7 of [10]. The LR algorithm is quadratically convergent and is numerically well behaved. These two properties are important, because the matrix R(x) needs to be computed for many values of x when we numerically solve the IVP (4.3). From our experience with this algorithm, it takes fewer than ten iterations to achieve a 10 −6 precision (when calculating G).
The forward Euler algorithm is known to have an error proportional to the step size h, and to be relatively numerically unstable at times, but it was found to be adequate. It would be easy to apply more sophisticated algorithms, such as general linear methods, which have a smaller error, and can be more numerically stable. The only adjustment required is to replace the Euler step in (11.3) by the alternative method.
In the numerical example in §11.3 below we let the ratio be r = 0.8 = 4/5, so that all the matrices, used in the computations for π 1,2 , are of size 10×10. It took less than 10 seconds for the algorithm to terminate (using a relatively slow, 1 GB memory, laptop). The same example, but with r = 20/25, so that the matrices are now 50 × 50, the algorithm took less than a minute to terminate. Moreover, the answers to both trials were exactly the same, up to the 7th digit. In both cases, we performed 5000 Euler steps (each of size h = 0.01, so that the termination time is T = 50). It is easily seen that π 1,2 had to be calculated for over 4500 different points, starting at the time π 1,2 becomes positive (see Figure 2 in the following example).
The validity of the solution can be verified by comparing it to simulation results, as in the example below and others in [15,16]. There are two other ways to verify the validity: First, we can check that the solution converges to the stationary point x * , which can be computed explicitly using (8.2). Second, within A we can see that the two queues keep at the target ratio r, even though this relation between the two queues is not forced explicitly by the algorithm. 11.3. A Numerical Example. We now provide a numerical example of the algorithm for solving the ODE in (4.1). In addition, we added the sample paths of the stochastic processes Q n 1 and Z n 1,2 , after scaling as in (3.2), on top of the trajectories of the solution to their fluid counterparts q 1 and z 1,2 .
The model has the same target ratio r = 0.8 as in the example in §6.2 with component rate matrices in (6.12). We chose a large queueing system with scaling factor n = 1000, so that the stochastic fluctuations do not to hide the general structure of the simulated sample paths. We let the ODE model parameters be m 1 = m 2 = 1, λ 1 = 1.3, λ 2 = 0.9, µ 1,1 = µ 2,2 = 1, µ 1,2 = µ 2,1 = 0.8, θ 1 = θ 2 = 0.3 and κ = 0. The associated queueing model has the same parameters µ i,j and θ i , but the other parameters are multiplied by n. The plots are shown without dividing by n.
We ran the algorithm and the simulation for 50 time units. We used an Euler step of size h = 0.01, so we performed 5000 Euler iterations. In each Euler iteration we performed several iterations to calculate the matrix G in (11.1), which is used to calculate the instantaneous steady-state probability π 1,2 . Figures 1-4 show q 1 (t)/q 2 (t), π 1,2 (x(t), q 1 (t) and z 1,2 (t) as functions of time t for a system initialized empty. After a short period, the pools fill up. Then q 1 (t) starts to grow, and immediately then fluid (customers) starts flowing to pool 2, causing z 1,2 (t) to grow. Figures 1-4 show that, for practical purposes, steady state is achieved for t ∈ [10,20].
In Figure 1 we see that once S b is hit, the ratio between the queues is kept at the target ratio 0.8. This is an evidence for the validity of the numerical solution, and a strong demonstration of the AP. In Figure 2 we see that initially, while q 1 = 0, π 1,2 = 0. This lasts until z 2,2 (t) + z 1,2 (t) = m 2 , at which time the space S is hit, specifically S b ), and the averaging begins. Once S b is hit, π 1,2 becomes almost constant, even before the system reaches steady state. Thus the functions q 1 , q 2 and z 1,2 have exponential form, supporting the results of §9.
We got x(t n ) ≡ (q 1 (t n ), q 2 (t n ), z 1,2 (t n )) = (0.3639, 0.4550, 0.2385) and Before solving the ODE, we can apply Theorem 10.2 to conclude that the solution will remain in A after it first hits A.
12. Proof of Theorem 7.1. We have previously observed that the first two conclusions involving S + and S − are valid. We now prove the three conclusions involving A, A + and A − . We will use the fact that a function mapping a convex compact subset of R m to R n is Lipschitz on that domain if it has a bounded derivative. Since we can always work with balls in R m (which are convex with compact closure), that in turn implies that a function mapping an open subset of R m to R n is locally Lipschitz whenever it has a bounded derivative on each ball in the domain; e.g., see Lemma 3.2 of [9]. The three sets A, A + and A − are convex.
The key is what happens in A.
For understanding, it is helpful to first verify this theorem in the special case r = 1, where the QBD process reduces to a BD process. Thus we first give a proof for that special case.
Proof for the special case r = 1. We use the fact that the ODE remains within A if it starts in A, so we are regarding A as an open connected convex subset of R 2 . The key component of the function Ψ in A is π 1,2 . We exploit the explicit representations in (6.18) and (6.21). From (6.1)-(6.4), the partial derivatives of λ ± (x) and µ ± (x) with respect to the three components of x, i.e., q 1 , q 2 and z 1,2 , are constants. From (6.18) and (6.21), we see that the partial derivatives of π 1,2 (x) with respect to each of the three components of x exist, are finite and continuous. That takes care of A.
We next consider A − and A + ; the reasoning for these two cases is essentially the same, with (6.18) making it quite elementary. We see that π 1,2 (x) → 0 and these partial derivatives approach finite limits as x → x b ∈ A − for x ∈ A, while π 1,2 (x) → 1 and these partial derivatives approach finite limits as x → x b ∈ A + for x ∈ A. In both cases we have a conventional heavy-traffic limit: ρ ± (x) ↑ 1 as x → x b . Hence, the partial derivatives of π 1,2 (x) are continuous and bounded on S b . As a consequence, for any ǫ-ball in S − S − about x in A + , there exists a constant K such that |π 1,2 (x 1 ) − π 1,2 (x 2 )| ≤ K x 1 − x 2 3 for all x 1 and x 2 in the ǫ-ball, where · 3 is the maximum norm on R 3 . A similar statement applies to A − . Hence we have completed the proof for r = 1. In closing, note that we cannot conclude that π 1,2 (x) is even continuous on all of S, because for x ∈ A we may have a sequence {x n : n ≥ 1} with x n ∈ S + for all n (or x n ∈ S − for all n), with x n → x as n → ∞, π 1,2 (x n ) = 1 for all n (or = 0), while 0 < π 1,2 (x) < 1.
We now treat the general case.
Proof of Theorem 7.1 in the general case. We first consider A. As in the case r = 1, we use the fact that the ODE remains within A if it starts in A, so we are regarding A as an open connected convex subset of R 2 . We will look at π 1,2 , and thus the QBD, as a function of the variable x ∈ A, which is an element of R 3 . By the definition of the matrices A 0 , A 1 and A 2 in (6.6) (see also the example in §6.2), these matrices are twice differentiable with respect to any of their elements. By the definition of the rates in (6.1)-(6.4), which are the elements of the matrices A 0 , A 1 and A 2 , these matrix elements in turn have constant partial derivatives with respect to each of the three real components of x at each x ∈ A, i.e., with respect to q 1 , q 2 and z 1,2 . It follows from Theorem 2.3 in He [7] that the rate matrix R in (6.15), which is the minimal nonnegative solution to the quadratic matrix equation A 0 + RA 1 + R 2 A 2 = 0, is also twice differentiable with respect to the matrix elements of A 0 , A 1 and A 2 , and thus also with respect to the three real components of x at each x ∈ A.
It thus suffices to look at the derivatives with respect to one of the elements of the matrices A 0 , A 1 or A 2 . It follows from the normalizing expression in (6.16) and the differentiability of R, that α 0 is also differentiable. Hence, from (6.17), we see that π 1,2 is differentiable at each x ∈ A, with By differentiating (6.16), we have so that α ′ 0 is continuous. The continuity of R ′ and α ′ 0 with respect to one of the elements of the matrices A 0 , A 1 or A 2 implies that the derivative π ′ 1,2 with respect to one of the elements of the matrices A 0 , A 1 or A 2 is finite and continuous on A, which in turn implies that the partial derivatives with respect to the three real components of x at each x ∈ A are finite and continuous as well. Hence, Ψ is locally Lipschitz continuous on A, as claimed.
We next show that π 1,2 and thus Ψ are locally Lipschitz continuous in neighborhoods of points in A + within S − S − and of points in A − within S − S + . We will only consider A + , because the two cases are essentially the same. In both cases, the situation is complicated starting from (12.1) because the entries of α 0 (x) become negligible, while the entries of (I − R) −1 (x) explode as x → x b . However, the two different limits cancel their effect. We exploit (6.19). The representation in (6.19) is convenient because now All key asymptotics take place in R + .
Since the crucial asymptotics involves only R + , we see that we only need carefully consider one of the two regions, in this case the upper one. To obtain results about R + , from a process perspective, it suffices to replace the given QBD by a new QBD with the upper region and reflection at the lower boundary. The new QBD model involving only R + is equivalent to a relatively simple single-server queue. The net input is a linear combination of four Poisson processes, and so has stationary and independent increments. The queue length process in the revised model is an elementary M AP/M SP/1 queue, as in §4 of [1], which has as QBD representation with rate matrix R + .
For the asymptotics, the key quantities are the spectral radii of the matrices R + (x) and R − (x), say η + (x) and η − (x), and the way that these depend on the drifts δ + (x) and δ − (x) as x → x b . The spectral radius η + (x) is the unique root in the interval (0, 1) of the equation det[A + 0 (x) + A + 1 (x)η + A + 2 (x)η 2 ] = 0, and similarly for η − (x); see (39) on p. 241 of [12], the Appendix of [13] and §4 of [1]. We see that η In general, we can represent powers of the matrix R (and similarly for R + and R − ) asymptotically as (12.3) R n = vuη n + o(η n ) as n → ∞, where u and v are the left and right eigenvectors of the eigenvalue η, respectively, normalized so that u1 = 1 and uv = 1. Moreover, as η → 1, the matrix inverse (I − R) −1 is dominated by these terms. Hence, we can do a heavy-traffic expansion of η + (x) and the related quantities as x → x b ∈ A + with x ∈ A, as in [2]; see the Appendix of [13]. As x → x b , all quantities in (6.19) have finite continuous limits as We then obtain where h is a continuous function on A + . Hence, there exist constants K 1 and K 2 such that for all x sufficiently close to x b . Finally, we can apply the triangle inequality with (12.7) to obtain |π 1,2 ( Hence, π 1,2 (x) and thus Ψ are locally Lipschitz continuous on A + within S − S − . Hence the proof is complete.
show that, if indeed x * ∈ A ∪ S + , then T < ∞, where The transient period of the fluid system can be divided into two distinct periods: The first transient period, on the interval [0, T ), lasts until the fluid limit hits S b . The second transient period is the one starting at the hitting time T , and is described by the ODE (4.2). This period was analyzed in the previous sections. The first transient period is described by different ODE's, depending on the state of the system. These ODE's, for the initial condition in (B.1), are given in the proof of Theorem B.1 below.
We shall prove that T < ∞ under the extra assumption that at no time during [0, T ) is z 2,1 > 0. The assumption can be verified directly by solving the fluid model of the first transient period. We discuss this condition after the proof of Theorem B.1.
Proof. We start by developing the ODE to describe the system before hitting S. As before, we do not consider the original queueing model and prove convergence to the appropriate fluid limit, but instead we develop the ODE directly. We first consider the case in s a 2 > 0 (so that q a 2 = 0), i.e., class 2 experiences no overload by itself (before pool 2 starts serving class-1 fluid). First, there is an initial period in which the pools are being filled with fluid. It is easy to see that as long as neither pool is full, the pool-content functions z i,i (t) behave as the fluid approximations for the number in system at time t in an M/M/∞ queueing model with arrival rate λ i and service rate µ i,i , i = 1, 2; e.g., see [14] (where it assumed that λ = µ, so that λ/µ = 1). Therefore, the system evolution is described by the pair of ODE'ṡ z 1,1 (t) = λ 1 − µ 1,1 z 1,1 (t), z 1,1 (0) = ζ 1 z 2,2 (t) = λ 2 − µ 2,2 z 2,2 (t), z 2,2 (0) = ζ 2 , and the unique solution to each ODE is These ODE's describe the dynamics of the two classes until one of the pools is full, i.e., until the time 1. Hence, from the beginning (time 0), z 1,1 increases until time t ′ 1 ≥ t 1 at which z 1,1 = m 1 . Then q 1 increases, satisfying (B.5) with q 1 (t ′ 1 ) = 0. By the assumption on x * , and following Corollary 8.2, there exists a time T < ∞ such that q 1 (T ) − rq 2 (T ) = κ. This is because rq 2 (t) ≤ rq a 2 < q a 1 − κ for all t ≤ T . On the other hand, it follows trivially from the solution to (B.5), that q a 1 is the globally asymptotically stable point of (B.5). Hence, for every ǫ > 0, there exists t ǫ such that q 1 (t) > q a 1 − ǫ for all t ≥ t ǫ . (This is because, by the initial conditions, q 1 (t) ≤ q a 1 for all t). Thus, we can find ǫ > 0 such that The second scenario of the second case has pool 1 filled first at time t 1 , so that q 1 starts increasing according to (B.5). If q 1 reaches κ before q 2 starts increasing, then we have the same behavior as when s a 2 > 0. However, if at time t 2 in (B.4) q 2 > 0, then the two queues will continue increasing independently until time T . Once again, (B.6) can be shown to hold, so that T < ∞.
We can easily calculate the exact value of x(T ) and use it to calculate the QBD drift rates δ + (x(T )) and δ − (x(T )) to find whether the positiverecurrence condition (6.10) holds at T , so that x(T ) ∈ A.
Remark B.1. (sharing in the wrong direction) In Theorem B.1 we assumed that we never have z 2,1 > 0. The reason is that, if z 2,1 ever does become positive, then the fluid x never hits the region S. To see that this is so, suppose that for some time t 4 sharing is initialized, with class-2 fluid flowing to service pool 1. Then z 2,1 is increasing until a time t 5 at which q 1 (t 5 ) − rq 2 (t 5 ) = κ, and the AP begins to operate. At that time, z 2,1 will start decreasing according to the ODĖ z 2,1 (t) = −µ 2,1 z 2,1 (t), t ≥ t 5 , whose unique solution is Hence z 2,1 remains strictly positive for all t ≥ t 5 , and S is never hit. Of course, the fluid state should be approaching a state in S as t increases. However, if there is such a limit point, then that limit point itself typically will not be a stationary point, because if x did start at that limit point, then it will have to continue to move toward the final stationary point x * .
More generally, the failure of z 2,1 to actually reach 0 in finite time has practical implications for the FQR-T control in the original queueing system. It suggests that it should be beneficial to relax the one-way sharing rule, by introducing lower positive thresholds for z 1,2 and z 2,1 . For example, if z 2,1 (t) > 0 at some time t ≥ 0, and at the same time sharing should be done in the other direction (because of a new overload incident, with class 1 being more overloaded and needing to get help), then we will allow pool 2 to start helping class 1, provided that z 2,1 is smaller than some threshold s 2,1 > 0. In that case, if z 2,1 (t) > s 2,1 , then z 2,1 will cross the threshold s 2,1 in finite time, as can be seen from (B.7). It remains to examine the system performance in response to such more complex transient behavior.
For the cases covered by Theorem B.1, the system evolution over the entire halfline [0, ∞) is a continuous "soldering" of the different ODE's, but at the soldering points t i , the functions under consideration are typically not differentiable. Hence, there is no single ODE that captures the full dynamics of the system. To see why, consider the case in which s a 2 > 0 and κ > 0. Then, for t < t 1 , q 1 (t) = 0 andq 1 = 0, but for t 1 ≤ t < t 2 , q 1 (t) evolves according to (B.5), which typically has a strictly positive derivative at t 1 . Thus the left and right derivatives at t 1 are not equal. Similar arguments hold for all the other soldering points.

APPENDIX C: THE ALGORITHM AND MORE EXAMPLES
C.1. More on the Algorithm. Let {t m : m = 0, 1, 2, . . . , n} be the Euler steps, with t m+1 − t m = h. In our experiments we found h = 0.01 to be a good candidate for the step size since it is small enough to minimize numerical errors, while the number of iterations needed for the ODE to reach its stationary point, is just a few thousands. Hence the algorithm takes only a few seconds to terminate.
LetD(t) ≡ q 1 (t) − rq 2 (t), denote the weighted difference between the two fluid queues. The discretization of the ODE in the numerical algorithm means that if, at step k − 1,D(t k−1 ) / ∈ S b but is close to it, thenD(t k ) may miss the boundary, even though the (continuous) ODE is at the boundary at time t k . For that reason, if κ − h <D(t k ) < κ + h, then we force x(t k ) to be in S b , by takingD(t k ) = κ. Once we haveD(t k ) = κ we decide whether to keep staying on the boundary for the next Euler step, by checking whether (6.10) holds. According to the relation between the QBD drift rates at time t k , we decide whether we should apply the AP, in order to find π 1,2 (t k ), or rather set π 1,2 (t k ) to zero or one.
At any step in the algorithm, we must also decide which ODE to use.
That depends on the state of the system at each time, as described in §B.
If the fluid state is not in S, as in the initial period of the example in §11 and the example below, then we use the appropriate fluid model, as given in the proof of Theorem B.1.

C.2. An
Example with x * ∈ S + . We now consider the same example as in §11.3, except now we increase the arrival rate for class 1 substantially, so that x * ∈ S + . In particular, we let λ 1 = 3.0 instead of 1.3. Once again, the system is initialized empty. That means that the fluid solution in S is moving between the two regions S b and S + . In particular, the solution first hits S b , as was proved in Theorem B.1, but it stays there for a short amount of time, and then crosses to S + .
As before, we show the results multiplied by n = 1000 in the figures below. We see how z 2,2 starts increasing up to the time T in which z 1,2 (T )+z 2,2 (T ) = m 2 . At this time z 2,2 (T ) starts decreasing, and is replaced by class-1 fluid. Since no class-2 fluid is flowing to either of the service pool, all the class-2 fluid output is due to abandonment. We can also observe that z 2,2 eventually hits 0, even though z 2,2 satisfies the equation (B.7). This is due to the numerical errors, as described in §B.
In steady-state we have q * 2 = λ 2 /θ 2 = 900/0.3 = 3000 and q * 1 = (λ 1 − m 1 µ 1,1 − m 2 µ 1,2 )/θ 2 = 4000, as in Corollary 8.1 (ii). C.3. An Example With x * ∈ S − Moving Through S + And S b . The purpose of this example is to illustrate more complex dynamics. We make class 2 more overloaded than class 1, i.e., q a 1 < q a 2 , but we make the rates faster for class 1. Specifically, we considered the following model parameters: λ 1 = 13.0, λ 2 = 1.5, µ 1,1 = 10.0, µ 1,2 = 0.8, µ 2,2 = 1 θ 1 = 2, θ 2 = 0.2, r = 0.8 and κ = 0. Note that the arrival, service and abandonment rates are all substantially greater for class 1 than for class 2. Nevertheless, class 2 is more overloaded than class 1: q a 1 = 1.5 < 2.5 = q a 2 . For this example, x * = (1.5, 2.5, 0) ∈ S − . We applied our algorithm to this example, letting the system start empty, i.e., x(0) = 0. The results are shown in the remaining figures, where here the results are scale up by multiplying by n = 100. Since the class-1 arrival rate is so large, q 1 starts filling up rapidly, and becomes full first; see Figures 9 and 10. Since κ = 0, pool 2 starts helping class 1 as soon as pool 1 becomes full. At first, pool 2 has spare capacity. However, soon the spare capacity in pool 2 is exhausted. At that time, the solution hits S. Even at the time pool 2 becomes full, we have q 1 > rq 2 , so that the solution enters S via S + , Thus pool 2 continues to help class 1 even after it is fully occupied, causing a dip in z 2,2 ; see Figure 10. However, the ratio of the queue lengths q 1 /q 2 decreases from its peak of about 1.4 until it reaches the target ratio r = 0.8, producing the desired relation q 1 = rq 2 + κ; see Figure 12 At that time (about t = 1.15, the solution that was in S + hits the set A. At that time, π 1,2 (x) jumps from 1 down to a value about equal to 0.6; see Figure 11. For an interval of time, the solution remains in A with the queue ratio fixed at the target r = 0.8. However, the load imbalance cause the solution to move within A, causing π 1,2 (x) to decrease until it reaches 0 in the set A − , at about time t = 2.5. From A − , the solution moves immediately into S − , where it rapidly converges to its stationary point. Of course, the stationary point x * is not actually reached in finite time. Indeed, after S − is reached, z 1,2 decreases exponentially to 0, but z 1,2 (t) > 0 = z * for all t, consistent with Remark B.1.
From the figures, it is evident that the numerical solution is not identical to the real solution. Because of the discrete step sizes in the Euler steps, the numerical solution misses A initially. In fact, we have to design the algorithm such that it "discovers" when S b is missed, and then force it to hit S b . That is easy to do since, if x(t k ) ∈ S + and x(t k+1 ) ∈ S − , where t k is the time of the k th Euler step, k ≥ 1, then S b must have been missed. We can then compute q 1 (t k+1 ) and take rq 2 (t k+1 ) = q 1 (t k+1 ) − κ.
This discreteness of the numerical solution explains the erratic behavior of π 1,2 at the hitting time of A, shown in Figure 11. The thick vertical line just after time 1, exactly when r = 0.8 for the first time as can be seen from Figure 12, appears because π 1,2 jumps between 0 and 1 at each Euler step. These jumps are caused the solution missing S b at first.

APPENDIX D: MISSING PROOFS
Proof of Theorem 7.2. Since 0 ≤ z 1,2 ≤ m 2 and q i ≥ 0 in S, we only need to prove the upper bounds (7.5). For i = 1, 2, let u i (t) be the function describing the queue-length process (of queue i) in a modified system with no service processes (so that all the fluid output is due to abandonment). The queue-length process in the modified system evolves according to the ODEu  whose solution is It follows that u i (t) ≤ u i (0) ∨ λ i /θ i and, when u i (0) = q i (0), the the righthand side in (7.5) is an upper bound for u i (t). We now show that this is also a bound for q i (t). For that purpose, define the auxiliary function f i (t) ≡ q i (t) − u i (t), t ≥ 0, and observe that f i (0) = 0 andḟ i (0) < 0. Hence, f is decreasing at 0 with f (t) < f (0) for all t ∈ [0, δ) for some δ > 0. This implies that q i (t) < u i (t) for all t ∈ [0, δ).
We now want to show that q i (t) ≤ u i (t) for all t ≥ 0. For a proof by contradiction, assume that there exists some t 0 > 0 such that q i (t 0 ) > u i (t 0 ), and let t 1 ≡ sup{t < t 0 : q i (t) = u i (t)}, t 2 ≡ inf{t > t 0 : q i (t) = u i (t)}.
By the contradictory assumption and the continuity of q and u, we have 0 < t 1 < t 0 < t 2 . (t 2 may be infinite.) Then (D.1) q i (t) > u i (t) for all t 1 < t < t 2 .
Proof of Corollary 8.2. We prove (i) only. The proofs for (ii) and (iii) are similar. First assume that x * ∈ S b . Since z * 1,2 ≥ 0, It follows from the expression for z * 1,2 in (i) of Corollary 8.1 that if q a 2 ≥ 0 then q a 1 − κ ≥ rq a 2 .
x * , which was developed heuristically in [15] using flow-balance arguments, is indeed the unique stationary point for the ODE. Moreover, we provided mild conditions under which the solution x(t) converges to x * as t → ∞.
We also showed that the convergence to x * is exponentially fast, further justifying the steady-state analysis in [15]. We showed that the existence of a unique solution to the IVP (4.3) depends heavily on the characterization of the function Ψ in (4.1) and its topological properties. These properties, in turn, depend on the state space of Ψ, and the regions of the state space in which Ψ is continuous. These regions are further characterized by the probabilistic properties of the family of FTSP's {D t : t ≥ 0}.
To further relate to the model considered in our previous paper [15], in §B we considered the system at the time when the arrival rates first change. At that time, the system will typically be underloaded, so that the state space should not be S. After the change, we assume that the arrival rates are larger than the total service rate of the two pools. Specifically, we assumed Assumption A in §8. We then considered the first transient period [0, T ), where T is the time at which S b is hit. Using alternative fluid models (ODE's), we showed that T < ∞, under the conditions of Theorem B.1. The solutions to the fluid models during the first transient period are all exponential functions, so that this period also passes exponentially fast.
Finally, we developed an efficient algorithm to solve the IVP (4.3), based on the matrix geometric method. This algorithm solves the different fluid models described in §B, and combines these solutions with the solution to (4.2) once the set A is hit, where the AP takes place.
It remains to quantify or at least bound the number of times the fluid solution moves from one of the regions S + , S − or S b to one of the others. Of course, the complexity of a solution is constrained by the fact that the solution path cannot cross over itself. It also remains to consider more complicated dynamics than provided by a single change in the arrival rates. The numerical algorithm applies more generally, but it remains to establish mathematical results and examine the performance. For example, it remains to consider a second overload incident happening before the system has recovered from the first one. Finally, it remains to establish analogs of the results here for more complex models, e.g., with more than two classes and/or more than two service pools.