{'fraction_non_alphanumeric': 0.04860963325530255, 'fraction_numerical': 0.010906575867205788, 'mean_word_length': 4.180431786862655, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reward Maximization in General Dynamic Matching Systems

We consider a matching system with random arrivals of items of different types. The items wait in queues -- one per each item type -- until they are"matched."Each matching requires certain quantities of items of different types; after a matching is activated, the associated items leave the system. There exists a finite set of possible matchings, each producing a certain amount of"reward". This model has a broad range of important applications, including assemble-to-order systems, Internet advertising, matching web portals, etc. We propose an optimal matching scheme in the sense that it asymptotically maximizes the long-term average matching reward, while keeping the queues stable. The scheme makes matching decisions in a specially constructed virtual system, which in turn control decisions in the physical system. The key feature of the virtual system is that, unlike the physical one, it allows the queues to become negative. The matchings in the virtual system are controlled by an extended version of the greedy primal-dual (GPD) algorithm, which we prove to be asymptotically optimal -- this in turn implies the asymptotic optimality of the entire scheme. The scheme is real-time, at any time it uses simple rules based on the current state of virtual and physical queues. It is very robust in that it does not require any knowledge of the item arrival rates, and automatically adapts to changing rates. The extended GPD algorithm and its asymptotic optimality apply to a quite general queueing network framework, not limited to matching problems, and therefore is of independent interest.


Introduction
We consider a dynamic matching system with random arrivals. Items of different types arrive in the system according to a stochastic process and wait in their dedicated queues to be "matched." Each matching requires certain quantities of items of different types; after a matching is activated, the associated items leave the system. There exists a finite number of possible matchings, each producing a certain amount of "reward". The objective is to maximize long-term average rewards, subject to the constraint that the queues of currently unmatched items remain stochastically stable. In this paper we propose a dynamic matching scheme and prove its asymptotic optimality. (In fact, the policy works for a more general objective, being a concave function of the long-term rates at which different matchings are used.) Figure 1 shows an example of a matching system with 4 item types. The items arrive as a random process, as individual items or in batches. The average arrival rate of type i items is α i . There exist 3 possible matchings; e.g. 1, 2 is a matching which matches one item of type 1 with one item of type 2. 2, 3, 4 is another matching which matches one item of types 2, 3 and 4. (In general, unlike in this example, a matching may require more than one item of any given type.) A matching can only be applied if all contributing items are present in the system; and if it is applied, the contributing items instantaneously leave the system. The analysis of static matching has a large literature (see, e.g., [9]). The dynamic model, which we focus on, has attracted a lot of attention recently, due to large variety of new (or relatively new) important applications. One example is assemble-to-order systems (see e.g. [12] and references therein), where randomly arriving Figure 1 An example of the matching model product orders are "matched" with sets of parts required for the product assembly. Another application is to Internet advertising [11], where the problem is to find appropriate matchings between the ad slots and the advertisers. Web portals as places for business and personal interactions is an important application; the problem in these portals (such as dating websites, employment portals, online games) is to match people with similar interests [3]. Matching problems also arise in systems with random arrival of customers and servers; for example, in taxi allocation, where matched "items" are passengers and taxis [8]. Further applications also can be found in [5,6].
Different control objectives may be of interest for matching systems. Gurvich and Ward [7] study the problem of minimizing finite-horizon cumulative holding costs for a model very close to ours. Plambeck and Ward [12], in the context of assemble-to-order systems, consider a model where item arrival rates can be controlled via a pricing mechanism; the objective includes queueing holding costs in addition to rewards/costs associated with order fulfillments, parts salvaging and/or expediting. Paper [12], in particular, proposes and studies a discrete-review policy; it involves solving an optimization problem at each review point.
A special case of the matching system, which received considerable attention, is where customers and servers are randomly arriving in the system and each server can be matched with one customer from a certain subset. This model, also known as the (stochastic) bipartite matching system, was initially studied by [6]. Majority of the previous research for this model was focused on finding the stationary distribution [1,2] and stability issues [3,4,10]. Bušić et al. [4] established the necessary and sufficient conditions for stabilizability of such systems, and have shown that the well known MaxWeight algorithm achieves maximum stability region. The problem of minimizing the long-term average holding cost for the bipartite matching system is studied by [5]. They have shown that with known arrival rates (and some other conditions on the problem structure), a threshold-type policy is asymptotically optimal in the (appropriately defined) heavy traffic regime.
In this paper, we show that the reward-maximizing optimal control of the matching model can be obtained by putting it into a typical queueing network framework. Our scheme uses a specially constructed virtual system, whose state, along with the state of the physical system, determines control decisions via a simple rule. In the virtual system any matching can be applied at any time and the queues are allowed to be negative. The matchings in the virtual system are controlled by (an extended version of) the Greedy Primal-Dual (GPD) algorithm [13], which maximizes a queueing network utility subject to stability of the queues. Negative queues in the virtual system can be interpreted as the shortages of physical items of the corresponding types. The GPD algorithm in [13] does not allow negative queues, so it is insufficient for the control of our virtual system. The main theoretical contribution of this paper is that we introduce and study an extended version of GPD, labeled EGPD, which does allow negative queues, and prove its asymptotic optimality under non-restrictive conditions that we specify. The approach of using a virtual system to control the original one has been used before, e.g. in [15], but the virtual system employed in this paper is substantially different, primarily because it allows negative queues.
Our proposed scheme is very robust in that it does not require a priori knowledge of item arrival rates, and automatically adjusts if/when the arrival rates change. It also covers a wide range of applications and control objectives. For example, in the context of assemble-to-order systems, the objective can include rewards/costs associated with order fulfillments, parts salvaging and/or expediting.
Although our scheme is designed (and proved asymptotically optimal) for the reward maximization objective, which does not include holding costs, we will discuss heuristic approaches to how the scheme can be used to achieve good performance in terms of a more general objective (including holding costs).
The paper is organized as follows. Section 2 contains notation used throughout the paper. In Section 3 we formally introduce the matching model and the reward maximization problem; here we also formally define the corresponding virtual system and the overall control scheme, in which the matching algorithm for the virtual system is a key part. In Section 4, we introduce the Extended Greedy Primal-Dual (EGPD) algorithm for a general network model, with queues that may be negative, and prove asymptotic optimality of EGPD; here we also show that the virtual system algorithm (in Section 3) is a special case of EGPD and thus is asymptotically optimal. (A reader interested mostly in applications of our proposed scheme may skip Section 4, at least at first reading.) We evaluate the performance of our scheme via simulations in Section 5. Finally, in Section 6, we discuss heuristics on how a more general objective, including holding costs, can be addressed by tuning EGPD parameters. Some conclusions are given in Section 7.

Basic Notation
We denote by R, R + and R − the set of real, real non-negative and real non-positive numbers, respectively. R N , R N + and R N − are the corresponding N -dimensional vector spaces. A vector x ∈ R N is often written as For a vector function f : For differentiable functions f : R → R and g : R N → R, we use f (t) (or (d/dt)f (t)) to denote the derivative with respect to t and ∇g(x For a set V and a real-valued function g(v), v ∈ V , denotes the subset of vectors v ∈ V which maximizes g(v).

Optimal Control of the Matching System
The outline of this section is as follows. First, we formally define the physical matching system in Section 3.1 and discuss the flexibility of this model to include a large variety of practical systems in Section 3.2. In Section 3.3 we introduce a virtual system, corresponding to the physical one. In Section 3.4 we define a control scheme, such that a certain algorithm runs on the virtual system, and control decisions for the physical system depend on those in the virtual one. We propose a specific algorithm for the virtual system in Section 3.5; this algorithm is asymptotically optimal in the sense that, under certain non-restrictive conditions, when the algorithm parameter (β) goes to zero, our entire physical/virtual control scheme maximizes average matching reward in the physical system. (The asymptotic optimality will be proved later, in Section 4.) We discuss features of the virtual system algorithm, and the conditions for its asymptotic optimality in Section 3.6.

Definition of the Physical Matching System
Consider a matching system with I item types forming set I = {1, · · · , I}. The customers arrive in batches, consisting of items of same or different types. To simplify exposition, assume that batches arrive as Poisson process, with each batch type chosen upon arrival, independently, according to some fixed distribution. There is a finite number of possible batch types. The average rate at which type i customers arrive into the system is α i > 0.
There is a finite set J = {1, · · · , J} of possible matchings. Let µ(j) = (µ i (j), i ∈ I), where µ i (j) ≥ 0 is the required number of type i items to form matching j ∈ J . Without loss of generality, we can and do assume that the "empty" matching, with all µ i = 0, is an element of J ; the empty matching is denoted ∅ . If a matching requires either zero or one item of each type, it is denoted by the subset of the required item types; say, 1, 2 denotes the matching requiring one item of type 1 and one item of type 2.
Without loss of generality, we can and do assume that the matching decisions are made only at the times of batch arrivals into the system. Essentially without loss of generality, we also assume that at those times at most m ≥ 1 matchings can be done. To simplify exposition, we further assume that m = 1 -it will be clear from our analysis that all results and (with very minor adjustments) proofs hold for arbitrary fixed m. Therefore, from now on we consider the system as operating in discrete (slotted) time t = 0, 1, 2, . . ., with i.i.d. batches arriving at those times, and exactly one (possibly empty) matching activated at each t.
Further, without loss of generality, we adopt the convention that the items arrived at time t are only available for matching at time t+1. (If items arriving at time t are immediately available for matching, the convention still holds if we simply pretend that they arrived at time t − 1, after the matching decision at time t − 1 was made.) Type i ∈ I items waiting to be matched form a first-come-first-served (FCFS) queue; its length is denoted Q i . At any time t, any one matching j ∈ J can be activated subject to the constraint that all the required items must be available in the system. With activation of matching j ∈ J , (i) Certain (real-valued) reward w j is generated; (ii) Number µ i (j) of items is removed from the queues of the corresponding types i.
Let X j be the long-term average reward generated by matching j, under a given control policy. We are interested in finding a dynamic matching policy, which maximizes a continuously differentiable concave utility function G(X 1 , · · · , X J ) subject to the constraint that all queue lengthsQ i (·) remain stochastically stable. Informally speaking, stochastic stability means that as time goes to infinity the queues do not "run away" to infinity, i.e. remain O(1). Formally, by stochastic stability we will understand positive recurrence of the underlying Markov process, describing the system evolution. (For example, if the process is a countablestate-space irreducible Markov chain, positive recurrence is equivalent to the existence of unique stationary probability distribution and to ergodicity.) Therefore, stochastic stability ensures that all arriving items are matched, without the backlogs and waiting times of unmatched items building up to infinity over time. Remark 1. Stability and long-term averages. We will give a specific definition of long-term average rewards X j later. When the process is Markov, positive recurrent, then X j can be thought of as the steadystate average reward u j due to type j matchings -we will elaborate on the relation between X j and u j later. Remark 2. More general µ i (j). Our model and the results hold -as is -in the case when the values of µ i (j) can be real numbers of any sign. A negative µ i (j) means that matching j adds |µ i (j)| items toQ i , and by convention any negative number of items of any type is always available for matching completion. We assume in this paper that µ i (j) are non-negative integers to keep the exposition intuitive.

Model Flexibility
The matching model defined in Section 3.1 is very flexible to include a variety of systems and their features. Let us consider assemble-to-order systems as an example. In such systems, orders for multiple products arrive as a random process. Each product requires a certain number of components of each type to be assembled. Components also arrive into the system as a random process. A product can only be assembled when all necessary parts are available; in which case it brings a certain reward (profit). This is a matching system where the components and product-orders of different types are "items", a completed product is a matching comprising one corresponding product-order and the required number of parts. Salvaging and/or disposing of the components is easily accommodated; namely, salvaging/disposing of one component, labeled as a type i item, can be treated as a matching i , with a reward that might be negative (as well as non-negative). Similarly with orders: discarding an order for a product, which is labeled as item type , is a matching with the corresponding (most likely, negative) reward. Expediting component delivery can be included as well. Suppose matching 1, 2, 3, 9 corresponds to product 9 assembled from (one unit of) parts 1, 2, 3, with the reward 20. However, the system has an option of expediting component 2, and receive it immediately, at the cost of 15. Then, assembling product 9 from already available components 1 and 3, and expedited component 2, can be modeled as a matching 1, 3, 9 with reward 20 − 15 = 5. (Another, more natural, way to model expediting of item 2 is to treat it as a "matching," requiring −1 type-2 items, with the reward −15. See Remark 2 above.) This discussion illustrates the flexibility of our model as long as the objective is to maximize average rewards associated with actions, such as matching, salvaging, expediting, etc. The model does not explicitly include holding costs. In Section 6 we propose and discuss heuristic extensions of our scheme which do implicitly take holding costs into account.

Virtual Matching System
We will propose a matching control scheme in Section 3.4, which in parallel to the physical system "runs" a virtual system, which determines the matching decisions for the physical one. The virtual matching system is defined as follows.
The virtual system has the same item types, set of matchings and arrival flows as the physical system. It is only different in that any matching can be activated at any time and the queues of the virtual system can be negative, as well as positive. Matchings in the virtual system are activated based on its own state, regardless of the state of physical system. The activated matchings in the virtual system become actual matchings in the physical system either immediately, or later in time, depending on availability of physical items. The virtual matchings, until they become actual ones, are called incomplete matchings. Incomplete matchings wait in a queue, which lists the incomplete matchings (their identities j) in the order of arrival; we denote the length of this queue byQ 0 . An incomplete matching becomes an actual one and leaves this queue when it is "completed" by all required physical items. (Incomplete matchings' queue, as we will see shortly, serves as the "interface" between the virtual and physical systems. In our figures and plots it is shown as part of the physical system.)

Control of the Physical Matching System via Virtual System
Denote by Q(t) = (Q i (t), i ∈ I) andQ(t) = (Q i (t), i ∈ I) the vectors of queue lengths in the virtual and physical systems, respectively, at time t. In this paper we always assume that the system is initialized in a state such that all physical and virtual queues are zero, Q i (0) =Q i (0) = 0, ∀i ∈ I, and there are no incomplete matchings,Q 0 (0) = 0. This means that the only feasible system states are those reachable from this "zero-state." At time t the following occurs sequentially: (i) A new matching is chosen in the virtual system based on Q(t). (We will give a specific rule in Section 3.5.) If it is a non-empty matching j, then the virtual queues are updated as Q := Q − µ(j), and a new type j incomplete matching is created and placed at the end of the (incomplete matchings') queue; so thatQ 0 :=Q 0 + 1.
(ii) The incomplete matchings' queue is scanned in FCFS order, to find the first incomplete matching j , which can be completed, i.e. such thatQ(t) ≥ µ(j ). If such matching j is found, it is completed, i.e. it is removed from the incomplete matchings' queue (so thatQ 0 :=Q 0 − 1), a physical matching j is created, and the corresponding number of physical items leaves the system,Q :=Q − µ(j ).
According to steps (i)-(iii) above, if matching j ∈ J is chosen in the virtual system at time t, the virtual queues change as follows: The evolution of the physical queues, if matching j ∈ J is completed is: Recall that we only consider feasible states of the queues -those reachable from the state where all virtual and physical queues are zero. Then we can make the following observations for the control scheme described above. For illustration, we will use Figure 2 showing a physical matching system with two item types and one possible matching and its corresponding virtual system. The system state shown on Figure 2 is such that: (a) in the physical system there are two type 1 items and no type 2 items; (b) the queue lengths in the virtual system are Q 1 (t) = 1, Q 2 (t) = −1; (c) there is one incomplete matching 1, 2 , which is incomplete because, while there is a type 1 item in the physical system (to complete it), there is no available (physical) type 2 item. (Note that at this point we did not specified yet the matching rule(s) for the virtual systemthis will be done in Section 3.5. So, the state on Figure 2 only illustrates the relation between virtual and physical systems, not a specific matching rule.) is the current surplus of type i items, beyond what is needed for completing all incomplete matchings. (On Figure 2, Q 1 = 1 indicates that there is one type 1 item in addition to one type 1 item which can be used for completion of the incomplete matching 1, 2 .) In addition to notations Q(t) andQ(t), let us denote byQ 0 (t) the state (list) of all incomplete matchings at time t.
The following simple Proposition 3 gives a total queue length bound (2) for the physical system in terms of the virtual one. This bound does not require any additional assumptions. Statements (ii) and (iii) of the proposition involve the notion of stochastic stability, which means positive recurrence of a Markov process. To keep the exposition simple, assume that the process (Q(t), t ≥ 0), describing the evolution of the virtual system, and the process [(Q(t), t ≥ 0), (Q(t), t ≥ 0), (Q 0 (t), t ≥ 0)] describing the evolution of the entire system, are countable-state-space Markov chains. (This is the case, for example, under the virtual system matching algorithm that we propose below in Section 3.5, and under linear utility function G.) Proposition 3. (i) At any t ≥ 0, the following relation between physical and virtual queues holds: where µ * .
is stochastically stable, then the steady-state average rates at which different matchings are activated are the same in the physical and virtual systems.
Proof. (i) Clearly, for all i ∈ I at all times, Q i (t) ≤Q i (t). Note that the total shortage of items of all types for the completion of all incomplete matchings is i Q − i (t); this means, in particular, that the total number of incomplete matchings is upper bounded asQ The total number of physical items in the system, iQ i (t), can be partitioned into those that are ready to be used for completion of incomplete matchings and the "surplus" items; the number of the former is upper bounded by µ * Q 0 (t); the number of the latter is equal to i Q + i (t). This implies the second part of (2). (ii) Follows from (i).
Remark 4. If m ≥ 1 matchings can be done after each arrival, the sequence of steps (i)-(iii) above is repeated m times.

Asymptotically Optimal Matching Algorithm for The Virtual System
We now specify the algorithm to be used for the control of the virtual system. This algorithms will be proved to be asymptotically optimal for the virtual system, and then (by Proposition 3) for the physical system as well -see Remark 6 below.

Algorithm 1 Matching Algorithm for the Virtual System
Let a (small) parameter β > 0 be fixed. At each time t = 1, 2, · · · , activate matching where running average values X j (t) (of the rewards obtained by activation of different matchings j) are updated as follows: and Q i (t) is updated according to rule (1) for all i ∈ I.
Note that if the function G is linear, say G(X) = j X j , then the partial derivatives in (3) are constant, and rule (3) becomes simply Moreover, in this case the algorithm does not need to keep track of the averages X j (t). As a result, both Consider the following Assumption 5 on the model structure. It is stated informally -its precise meaning will be given (in a more general context) later in Assumption 8 (Section 4). Also, in Section 3.6.2 we explain why this assumption is non-restrictive.
Assumption 5. For any subsetĪ ⊆ I, there exists a matching activation strategy, under which the longterm average drift of queues i ∈Ī is strictly positive and the long-term average drift of queues i ∈Ī is strictly negative.
When parameter β is small, then the running average X j (t) is (one notion of ) a long-term average rate at which rewards due to matching j are generated. (See Section 4.4.) We will prove in Section 4 (as a corollary of Theorem 10) that, under Assumption 5, Algorithm 1 is asymptotically optimal in the following sense. (It is described here informally -the formal result is Theorem 10, for the more general model in Section 4.) Let V be the set of those long-term rate vectors X that are achievable (by some control strategy) subject to the stability of the queues, and let V * be its optimal subset, V * = arg max X∈V G(X). Then, when β is small, Suppose now that the system process is Markov under Algorithm 1 (as is the case when function G is linear). Then Assumption 5 ensures the process stability (for example, by the argument described in Section 4.9 in [13]). In this case the steady-state average rewards (due to different matchings) u = (u 1 , . . . , u J ) are well defined. If the process is in stationary regime, then obviously EX(t) = u. Furthermore, the asymptotic optimality of Algorithm 1 in the sense described above, can be used to show that, as β → 0, the vector u converges to the optimal set V * (see Section 4.9 in [13]).
Remark 6. If Algorithm 1 is asymptotically optimal for the virtual system, then under our scheme it is also asymptotically optimal for the physical system. Indeed, the physical and virtual systems have the same set V of achievable long-term rate vectors X (subject to the stability of the queues). This is because any X achievable in the virtual system is achievable in the physical system as well (by our scheme, for which we have Proposition 3), and vice versa because obviously any control of the physical system can be applied to the virtual system. Therefore, under our scheme, if the virtual system produces (in the asymptotic limit) the optimal long-term rates X ∈ V * , the same optimal rates are produced (by Proposition 3) in the physical system.
3.6 Discussion of Algorithm 1

Basic intuition
The key feature of the virtual system is that it has an option of creating matchings "in advance," before all required physical items have arrived. These "advance" matchings are the ones we called incomplete. Virtual queues keep track of the items' availability: recall that if Q i < 0, |Q i | is the shortage of type i items, and if Q i ≥ 0, it is the surplus of type i items.
The intuition behind Algorithm 1 is the same as for the GPD algorithm in [13] (and other related works -see, e.g., [14] and references therein), but our model is more general in that the queues may have any sign. For simplicity of discussion, suppose the objective function is linear, G(X) = j X j , in which case Algorithm 1 specializes to (6). The rule "tries" to choose a matching j which brings large reward w j , but at the same time it "encourages" the drift of the queues towards 0. Indeed, recall that activation of any matching can only decrease the virtual queues. This means that the rule "encourages" the use of matchings that decrease positive Q i 's as much as possible and decrease negative Q i 's as little as possible; in other words, the rule encourages matchings requiring items of which there is a large surplus, and discourages matchings requiring items of which there is already a large shortage -this guarantees stability of the queues. When parameter β is small, the virtual queues "stabilize around correct levels" -positive or negative -which allows rule (6) to make "correct" decisions maximizing the average rewards.

Assumption 5 is non-restrictive
We now describe two common cases, in which Assumption 5 holds. These two cases cover a very large number of applications.
Case 1. Assumption 5 holds automatically in the special case when for each item type i there exists at least one matching requiring only type i items (namely, with µ i ≥ 1 and µ = 0 for = i). In this case it suffices to pick any parameter m (the number of matchings per each batch arrival) which is greater than µ * . = max j i µ i (j). This special case is very common for the following reason, which we illustrate using the simple model in Figure 2. If matching 1, 2 is the only possible (besides the empty matching), the system is unbalanced when the arrival rates are unequal, α 1 = α 2 , and cannot be stable. (If items arrive one-by-one, this particular system obviously cannot be stable even if α 1 = α 2 . More generally, any system with one-byone arrivals cannot be stable if its "matching graph" is bi-partite, see [10].) This shows that many practical systems typically need the option of using "single" matchings i anyway (salvaging or discarding individual items), to ensure stability, and then Assumption 5 holds.
Case 2. This case is more subtle. Suppose a system can potentially be made stable without requiring single-type matchings. For example, consider the system in Figure 2 in which the arrivals occur only in pairs (1,2). Suppose also that up to two matchings can be done upon each arrival (m = 2). On the face of it, Assumption 5 does not hold for this system. Indeed, the linear relation Q 1 (t) = Q 2 (t) holds at all times and, therefore, it is impossible for Q 1 and Q 2 to have different average drifts, which is required under Assumption 5. However, consider the orthogonal change of coordinates,Q 1 = Q 1 + Q 2 ,Q 2 = Q 1 − Q 2 , with λ(·) and µ(·) transformed accordingly. Then,Q 2 (t) ≡ 0, and the system can be considered as having only one queueQ 1 . For the latter system Assumption 5 does hold. Note that the algorithm itself does not need to perform any change of coordinates -it remains as is. This situation is generic: if there is an inherent linear dependence between the queues, Assumption 5 often holds for the system after an appropriate orthogonal change of coordinates. This is, in fact, the case for many bi-partite matching systems (with items arriving in pairs), including the one we consider later in Section 6.2.
To summarize the discussion in this subsection, Assumption 5 is essentially the assumption that the system can be made stable, plus a very common condition that the queues "can be moved in any direction" within the subspace of feasible queue states.

A General Network Model and EGPD Algorithm
In this section we introduce the Extended Greedy Primal-Dual (EGPD) algorithm for a general network model, which includes the matching system as a special case. This algorithm is a generalization of the GPD algorithm of [13] in the sense that queues at some network nodes, we call them free nodes, are allowed to have any sign; as they evolve, these queues are "free" to change from positive to negative and vice versa. The model in [13] is such that queues at all nodes are constrained to be non-negative -in our model we call such nodes constrained. First, we will formally define the model and the underlying optimization problem in Sections 4.1-4.3. The optimization problem determines the best possible (under any control algorithm) long-term drifts of the queues, which maximize the network "utility" subject to the condition that queuedrifts are zero at free nodes and are non-positive at constrained nodes; the optimal solutions to this problem give the maximum possible network utility that can be achieved by any network control strategy subject to stability of the queues. We define the EGPD algorithm in Section 4.4. In Section 4.5, we show that, as the algorithm parameter β goes to 0, the "fluid scaled" version of the process converges to a random process with sample paths being what we define as EGPD-trajectories. In Section 4.6 we prove asymptotic optimality of the EGPD-algorithm, in the sense that EGPD-trajectories converge to the optimal set of the underlying optimization problem while keeping all queues uniformly bounded; in other words, EGPDalgorithm maximizes the system utility subject to stability. Finally, in Section 4.7 we show that Algorithm 1 (Section 3.5) for the virtual system of Section 3.3 is a special case of EGPD.
A reader interested only in the application of EGPD algorithm to the dynamic matching model of Section 3 may wish to skip at first reading the proofs in Sections 4.5-4.6.

The Model
Consider a network consisting of a finite set of nodes N = {1, 2, · · · , N }, N ≥ 1. The nodes are of two different types: N 1 constrained nodes form the set N c = {1, 2, · · · , N 1 } and N 2 = N − N 1 free nodes form N f = {N 1 + 1, N 1 + 2, · · · , N }. Either N c or N f is allowed to be an empty set. There is a queue associated with each node, where we denote by Q n (t) the queue length of node n ∈ N at time t and we will denote Q(t) = (Q n (t), n ∈ N ). The queue length of node n ∈ N c is always non-negative, but node n ∈ N f can have queue length of any sign.
The system operates in discrete time t = 1, 2, · · · . (By convention, we identify an integer time t with unit time interval [t,t+1), which is usually referred to as time slot t.) A finite number of controls is available, where we denote by K the set of controls. With activation of control k ∈ K at time t, the following occurs sequentially: (i) A certain (non-random) real amount ("number") µ n (k) ≥ 0 of items is removed from queue n and leaves the network. Queues in constrained nodes cannot go below zero; so if Q n (t) ≤ µ n (k), the entire content of queue n is removed.
According to steps (i) and (ii), the queue update rules for constrained and free nodes, given control k is chosen at time t, are as follows:

System Rate Region
For each k ∈ K and time t, consider the random vector b(k, t) = (b n (k, t), n ∈ N ) equal in distribution to λ(k) − µ(k). Clearly, b(k, t) is equal to the random vector of queue increments Q(t + 1) − Q(t) provided that control k is chosen at time t and assuming Q n (t) ≥ µ n (k) for all n ∈ N c . We call components of b(k, t) the nominal increments of queues upon control k at time t. Let k(t) denote the control chosen at time t by a given control policy.
Informally speaking, the finite-dimensional convex compact rate region V ⊂ R N is defined as the set of all possible long-term average values of b(k(t), t), which can be induced by different control policies. Formal definition of the rate region is as follows.
For each k ∈ K, denote by b(k) = Eb(k, t) the drift of queue lengths upon control k (at any time t when control k is activated). For a fixed probability distribution φ = (φ k , k ∈ K) (with φ k ≥ 0 and k∈K φ k = 1) consider the vector v(φ) = k∈K φ k b(k).
If we interpret φ k as the long-term average fraction of time slots when control k is chosen from the set of controls K, then v(φ) corresponds to the vector of long-term average drifts of Q(t), assuming that the queues in the constrained nodes never hit zero. Then the system rate region V is defined as the set of all possible vectors v(φ) corresponding to all possible φ.

Underlying Optimization Problem
Consider an open convex setṼ ⊆ R N such thatṼ ⊇ V . Consider a concave continuously differentiable utility function H :Ṽ → R and the following optimization problem: {v ∈ V : v n ∈ R − , ∀n ∈ N c and v n = 0, ∀n ∈ N f } = ∅.
If Assumption 7 holds, we denote by V * ⊆ V the set of optimal solutions of (10). The dual to optimization problem (10) is and we denote by Q * the closed convex set of optimal solutions q * ∈ R N1 + × R N2 of problem (12). For any v * ∈ V * and any q * ∈ Q * , the compementary slackness condition holds: In Section 4.4, we will introduce an algorithm, which is asymptotically optimal under the following assumption, which is stronger than Assumption 7.

Assumption 8.
For any subsetN f ⊆ N f , there exists v ∈ V such that v n > 0 for n ∈N f and v n < 0 for n ∈N f . Assumption 8 means that there always exists a control policy which provides, simultaneously, a strictly negative average drift to all the constrained node queues and non-zero average drifts toward zero for all free node queues.
Note that under Assumption 8, the set Q * is compact. Indeed, the optimal value of the problem (10) is equal to for any v * ∈ V * and any q * ∈ Q * . Set Q * must be bounded, because otherwise, from Assumption 8, there would exist v ∈ V such that v n < 0 for all nodes with q n ≥ 0, and v n > 0 for all nodes with q n < 0. Then we can arbitrarily increase the RHS of (14) by choosing q * ∈ Q * with large |q * n |. The problem that we are going to address is as follows. Let X denote a long-term average value of b(k(t), t) under a given dynamic control policy, that is, a policy of choosing k(t) depending on the system state. We are interested in finding a dynamic control policy such that when optimization problem (10) is feasible, and moreover, the stronger Assumption 8 holds, the corresponding X is close to V * , while the system queues remain stochastically stable.

Extended Greedy Primal-Dual Algorithm
Consider the following control policy: Algorithm 2 EGPD algorithm for the general network model At time t = 1, 2, · · · , choose a control where β > 0 is a small parameter. Here X(t) is the running average of b(k(t), t), updated as follows: and Q(t) is updated according to (7) and (8).
The initial condition is X(0) ∈Ṽ . Note that such initial condition and update rule (16) imply that X(t) ∈Ṽ for all t ≥ 0. Hence the system evolution is well-defined for all t ≥ 0, since the gradient and argmax in (15) are well-defined.
Also note that, if 0 < β < 1, then for t ≥ 1 Therefore, when t is large, X(t) is essentially the geometric average of values of b(k(τ ), τ ) up to time t − 1. When t is large and β > 0 is small, X(t) is (one notion of) the long-term average of values of b(k(τ ), τ ) up to time t − 1.

Asymptotic Regime and Fluid Limit
We define EGPD-trajectory as a pair of absolutely continuous functions (x, q) = ((x(t), t ≥ 0), (q(t), t ≥ 0)), each taking values in R N and satisfying the following conditions: and for almost all t ≥ 0, where Functions x(t) and q(t) are dynamically changing primal and dual variables, respectively, for problems (10) and (12), which arise as asymptotic limits of the fluid scaled version of the process as described next.
Consider a sequence of processes (X β , Q β ), indexed by a parameter β, where β ↓ 0 along a sequence B = {β j } ∞ j=1 with β j > 0 for all j. The initial state (X β (0), Q β (0)) ∈Ṽ is fixed for each β ∈ B. (The processes and variables associated with a fixed parameter β will be supplied by superscript β.) We need to augment the definition of the process. Let us assume X β (t) and Q β (t) are functions defined on t ∈ R + and constant within each time slot [l, l + 1), l = 0, 1, 2, · · · . Thus for each β, consider the (continuous-time) process Z β = (X β , Q β ), where For each β, is the fluid scaled version of process Z β , obtained by The following theorem is straightforward modification of Theorem 3 in [13], which we present without proof.
Theorem 9. Consider a sequence of process {z β } with β ↓ 0 along set B. Each process is considered as a random element in the Skorohod space of RCLL ("right continuous with left limits") functions. Assume that z β (0) → z(0), where z(0) is a fixed vector in R 2N such that X(0) ∈Ṽ . Then, the sequence {z β } is relatively compact and any weak limit of this sequence (i.e a process obtained as the weak limit of a subsequence of {z β }) is a process with sample paths z being EGPD-trajectories (with initial state z(0)) with probability 1.

Global Attraction Result
The following theorem is the main result of this section which shows the convergence of EGPD-trajectories to the saddle set V * × Q * .
The proof of Theorem 10 is a generalization of that of Theorem 2 in [13] -all steps of the latter are extended to our more general setting. For this reason we will not give a complete proof of Theorem 10 in this paper, because it is lengthy. Instead, we demonstrate the key points involved in the generalization, by proving in this section the convergence (27) for the special case when x(0) ∈ V and H(·) is strictly concave.
A time point t ≥ 0 is called "regular " if conditions (17)-(19) are satisfied and proper derivatives x (t), q (t) and f (t) exist. Almost all t are regular.
Let us introduce the following function: n∈N y 2 n , v ∈Ṽ , y n ∈ R + for n ∈ N c , y n ∈ R for n ∈ N f . Lemma 11. The trajectory (q(t), t ≥ 0) is such that Proof. Within this proof, we say that a vector-function (or scalar function) α(t), t ≥ 0, is uniformly bounded if sup t≥0 α(t) < ∞. By Assumption 8, the following holds for some fixed number δ > 0. For any t ≥ 0, there exists ξ = (ξ n , n ∈ N ) ∈ V such that for any n, |ξ n | ≥ δ, ξ n > 0 if q n < 0, and ξ n < 0 if q n ≥ 0. Then for any regular t ≥ 0 (and a corresponding ξ) we have: Since ∇H(x(t)) and x(t) are uniformly bounded, so it the second term in (32). When q(t) is large, the first term in (32) is large positive. We conclude that (d/dt)F (x(t), q(t)) ≥ 1 > 0 as long as q(t) ≥ C 1 > 0, for some fixed constants 1 and C 1 . Since H(x(t)) is uniformly bounded, we can pick C 2 > 0 sufficiently large so that F (x(t), q(t)) ≤ −C 2 implies q(t) ≥ C 1 and then (d/dt)F (x(t), q(t)) ≥ 1 > 0. We then conclude that lim inf t→∞ F (x(t), q(t)) ≥ −C 2 and, therefore, inf t≥0 F (x(t), q(t)) > −∞. The latter (along with the uniform boundedness of H(x(t))) implies that q(t) is uniformly bounded.
Lemma 12. For any EGPD-trajectory, at any regular time t ≥ 0, and Furthermore, if Assumption 7 holds, Proof. Noting q n (t) = v n (t) and v * n = 0, for any n ∈ N f , every step of the proof is analogous to that of Lemma 3 in [13].
Select an arbitrary point q * ∈ Q * and associate it with the following function is the Lagrangian of problem (10) with the dual variable equal to q * ∈ Q * . Having strictly concave H(·) implies that H * (·) is also a strictly concave function and is the unique optimal solution. Note that ∇H * (v) = ∇H(v) − q * .
Lemma 13. Consider F * (·, ·) associated with an arbitrary q * ∈ Q * . Then for all (regular) t ≥ 0, and Proof. The proof is analogous to that of Lemma 5 in [13]. The only difference is the existence of free nodes, where we can easily validate this Lemma by using q n (t) = v n (t) and v * n = 0 for any n ∈ N f .

Mapping of the Virtual Matching System of Section 3.3 into EGPD Framework
Now we are in position to show that Algorithm 1 for the control of the virtual system in the original matching model in Section 3 is a special case of EGPD Algorithm. The mapping of the virtual system of Section 3.3 into the more general model of Section 4.1 is as follows. Consider the following system, which we refer to as a modification of the virtual system. Suppose the item types I are modelled as free nodes and let the set of matchings J be the set of controls K. Let us add one constrained node per each matching j ∈ J .
(These additional nodes are the utility nodes in the terminology of GPD algorithm [13].) From this point on, for convenience of the notations, we replace the set of indices of item types I with {J + 1, · · · J + I} and denote by I c = {1, · · · , J} the set of all constrained nodes. For the constrained nodes we adopt the convention that they never receive any inputs, i.e. λ j (t) ≡ 0, j ∈ I c . We also fix a sufficiently large c > 0, so that w j − c < 0 for all constrained nodes, and for each constrained node (or, matching) j ∈ I c we set by convention µ j (j) = c − w j > 0 and µ i (j) = c > 0, i ∈ I c \ {j}. These conventions about the constrained nodes guarantee that under any control strategy, their queues are automatically stable. In fact, for any i ∈ I c and any initial value Q i (0), the queue length Q i (t) will decrease until it hits 0 within a finite time and then it will remain at 0. This allows to assume, without loss of generality, that Q i (t) ≡ 0 for all constrained nodes . Finally, we define the utility function H(v) as follows: Given these conventions, it is easy to see that the problem of maximizing G(X 1 , . . . , X J ) (subject to the stability of the queues) in the original matching system is equivalent to the problem of maximizing H(X 1 , . . . , X J , X J+1 , . . . , X J+I ) (subject to the stability of the queues) in the modified system defined in this subsection. The latter system is a special case of the general system of Section 4.1. If we specialize Algorithm 2 to the modified system, and then rewrite it in terms of the original virtual system, we obtain Algorithm 1. Assumption 8, specialized to the modified system and expressed in terms of the original virtual system, gives the formal meaning of Assumption 5 (which is stated informally).
The mapping described in this section and the asymptotic optimality of the EGPD algorithm under Assumption 8 imply the asymptotic optimality of Algorithm 1 under Assumption 5.

Simulations
In this section, we evaluate the performance of EGPD algorithm via simulations. Consider the system described in Section 1. We extend the set of possible matchings by including "single" matchings (see Section 3.6.2): The reward vector is w = (0, −1, −1, 1, 2, 5, 4, 7) where its j-th component corresponds to the jth element of the matchings' set. We consider a linear utility function, namely the sum of average rewards due to different matchings. The vector of arrivals rates is α = (1.2, 1.5, 2, 0.8).
For our linear utility function, the EGPD algorithm for the virtual system is given by rule (6).
A. Average reward maximization. We use parameter β = 0.01. Figure 3 shows the queue trajectories of the virtual and physical systems under the EGPD algorithm. All queues are initially empty. We observe that all queues are quickly "converging". Nearly all type 2 and 4 items are matched right after they enter the system, while there exist a queue of around 100 items of types 1 and 3. The rates at which matchings are activated under EGPD algorithm are provided in table 1, which shows that these rates are close to the optimal ones, obtained by solving the underlying optimization problem (which is a linear program in this case). Therefore, as expected, the algorithm yields near optimal performance for small β. Note that solving the optimization problem requires the knowledge of arrival rates (as well as other system parameters), while our algorithm need not know arrival rates.   Figure 4 Average matching reward under the EGPD algorithm.
graph, the running average reward under EGPD algorithm is getting very close to optimal objective value and this convergence is sufficiently fast.
B. Effect of parameter β. In order for βQ in the virtual system to "stay close" to some q * ∈ Q * , parameter β should be small. Therefore, as long as parameter β is sufficiently small, the algorithm is nearly optimal and the virtual queue lengths are roughly of the order 1/β. As β is increasing, the accuracy of the algorithm in terms of average reward maximization decreases, while the queues become smaller.
The dependence of the average reward on β for the considered scenario is shown on Figure 5. First, we note that the average reward remains nearly optimal for values of β almost as large as 1 (i.e. not even very small in absolute terms). Then, as β changes from 1 to about 10, the average reward decreases and reaches the lower "plateau," and then remains constant for β ≥ 10. Thus, as expected, the algorithm is effective in terms of reward maximization when β is sufficiently small (less than 1 in our scenario); when β is sufficiently large (greater than 10 in our scenario), the average reward is also roughly independent of β, but is at a lower, suboptimal level. We note that larger values of β have the benefit of reducing the queues and, as a result, reducing (as we will see next) the algorithm response (or, adaptation) time to changes of the items' arrival rates. (Shorter queue also mean lower holding costs, if such are a part of the model. This will be discussed in Section 6.) Therefore, the value of parameter β should be chosen, very informally speaking, "as large as possible, but not larger".
C. Automatic adaptation to changes in arrival process. An important robustness issue is how quickly the EGPD algorithm responds to the changes in the arrival process. In the following experiment, the arrival rates are changed to α = (1.8, 0.8, 1.4, 1) at time 2000. This change leads to different optimal matching rates and thus different optimal value. If quick response to arrival rate changes is important, a larger β is preferable. Here we use β = 0.1. Figure 6 shows the queue trajectories of the virtual and physical systems.
We observe that EGPD automatically adapts to the new arrival rates and reaches the new "right" queue lengths, without using any a priori information on this change. 6 Heuristics for the Objective Including both Matching Rewards and Holding Costs

General discussion
The scheme we proposed in Section 3.4 for the matching model is asymptotically optimal for the reward maximization problem. (We will refer to this entire scheme as EGPD, because EGPD is its key part, applied to the virtual system and determining matching choices.) In practical systems, the objective may be more general, namely maximizing the average "profit" defined as average reward minus average queue holding cost. We now informally discuss how EGPD can be used to achieve better profit in the system (even though it is not specifically designed for that).
For the purposes of the discussion below, we assume linear holding costs with rate vector c = (c i , i ∈ I); that is the average holding cost over interval [0, T ] is Suppose the arrival rates are scaled up by a factor r > 0. This simply speeds up the process r times, so that the average reward increases r times, while the holding cost remains same. Thus for systems with "high" arrival rates, the rewards dominate the profit objective and we expect the average profit obtained by the EGPD algorithm to be "close" to the optimal one. In other cases, holding costs may dominate, for example when the system is in (appropriately defined) heavy traffic (see [7,5]) -this makes the queues necessarily large. When the optimal average rewards and optimal holding cost are on the same scale, the EGPD parameter settings can be used to control the tradeoff between these two performance measures, thus potentially improving the average profit. We now briefly discuss different heuristic approaches for profit improvement within the framework of our scheme.
Choice of parameter β. As discussed in Section 5, as long as parameter β is sufficiently small, the virtual queue lengths under EGPD are large, roughly of the order 1/β. To see how this affects the holding cost, consider two cases: (i) If Q i (t) ≥ 0, thenQ i (t) will also be large (of the order of at least 1/β) since the inequalityQ i (t) ≥ Q i (t) holds for all i ∈ I at all t.
(ii) If Q i (t) < 0, this has an indirect impact on the holding cost. In particular, large |Q i (t)| in this case would imply more incomplete matchings. This subsequently results in a higher holding cost.
Therefore, parameter β should be chosen as large as possible, but not to exceed the level beyond which the average rewards start to be significantly (negatively) affected.
Additional queue scaling. Consider arbitrary positive weights γ i , i ∈ I. All the results for the EGPD algorithm hold if we use more general rule instead of (3). In this case, it is the weighted vector (γ i βQ i (t), i ∈ I) (not βQ(t)) that will be close to an optimal dual solution q * . This property may be used to reduce the holding cost by giving higher weights to more "expensive" queues (with large c i ), thus making them relatively smaller.
Matching completion order. There is a flexibility in choosing which incomplete matching to complete first. For the average matching reward maximization this does not matter (so, earlier we specified FCFS rule for concreteness). However, if the holding costs are a consideration, one may pick incomplete matchings with higher associated holding cost to be completed first.

Simulation: Average Profit in a Bipartite Matching System
Consider a bipartite matching system, where items arrive in pairs, and the matchings are pairs as well. It is depicted in Figure 7. There are 8 item types {1, 2, 3, 4, 1 , 2 , 3 , 4 }. The arrival graph is on the left, where each edge shows a possible arrival pair, and the plot in the right hand side is the matching graph with edges representing the possible matchings. Up to two matchings can be done per each arrival (m = 2).

Figure 7
Illustration of the matching system.
We consider the process in discrete time t = 1, 2, · · · . The arrival process is i.i.d. across time. Specifically, at each time t, a pair of items enters the system. The probabilities (rates) of different arrival pairs are specified in Table 2. It is easy to check that this system satisfies necessary and sufficient condition [4] for bi-partite matching systems to be stabilizable. The condition, called NCond in [4], is as follows. Suppose the matching graph is connected. Consider a subset T of items from the top part of the bi-partite graph, and denote by α T the total arrival rate of all items in T . Denote by B(T ) the subset of items from the bottom part of the graph that can be matched with at least one item in T , and by α B(T ) the total arrival rate of all items in B(T ). Then, the system is stabilizable if and only if α T < α B(T ) holds for any strict subset T of "top" items.
Now, let us see if our Assumption 5 holds. Since this is a bipartite matching system, with items arriving and departing in pairs, virtual queues satisfy the following linear relation Q 1 (t) + Q 2 (t) + Q 3 (t) + Q 4 (t) − (Q 1 (t) + Q 2 (t) + Q 3 (t) + Q 4 (t)) ≡ 0. However, given that NCond condition holds, it is easy to see that our Assumption 5 (formally given by Assumption 8) holds for this system in the sense described in Section 3.6.2, namely after an orthogonal change of coordinates. (We emphasize again that the algorithm itself remains as is, it does not need to do any change of coordinates.) Therefore, the EGPD algorithm is asymptotically optimal for this system for the average reward maximization objective.
Assume linear holding costs, c ·Q(t), with the cost rate vector c = (0.1, 0.2, 0.3, 0.4, 0.4, 0.3, 0.2, 0.1). The matching rewards for different matchings are given in Table 3. We simulated this system under EGDP scheme. Figure 8 shows the dependence of EGPD average performance metrics on the parameter β. The range of β is shown within which the average reward declines from its optimal (largest) value to the "plateau" it reaches when β is large. Parts (a), (b) and (c) show average holding cost, reward and profit, respectively; the average profit is the average reward minus the average holding cost.  We see that the average profit is maximized within a certain range of values of β, where, roughly speaking, the average reward is "still" close to optimal and the average holding cost is "already" close to the best achievable by EGPD. We conjecture that the average profit with such choice of β is reasonably close to the optimal profit under any control algorithm. Verifying and quantifying this informal conjecture is an interesting subject for future research.

Conclusions
In this paper we have proposed an approach for optimal dynamic control of general matching systems. The central idea is using a virtual matching system allowing negative (as well as positive) queues, as part of the overall control scheme. The virtual system fits into a queueing network framework, except the queues may be negative, and it is controlled by an extended version of the GPD algorithm, called EGPD. We prove EGPD asymptotic optimality. The approach is very generic, not restricted to special cases, such as bipartite matching. The proposed scheme is also very robust in the sense that it does not require the knowledge of input rates, and automatically adapts to changing input rates. Simulations demonstrate good performance of the algorithm.
Although the scheme that we develop has the average reward maximization as its objective, the parameter setting can be used to achieve good performance in terms of the more general objective, which includes holding costs. Addressing this and other more general objectives within a dynamic control framework, not requiring a priori knowledge of the item arrival rates, is an important future subject.