An Optimal Control Approach to Optimization of Multiclass Queueing Networks

Florin Avram ^{∗}  Dimitris Bertsimas ^{†}  Michael Ricard ^{‡}

June 1994

Abstract

We propose an optimal control approach to the optimization of fluid models of open multiclass queueing networks. Using Pontryagin’s maximum principle we show that the optimal policies of specific examples of networks are of the threshold type, the parameters of which we calculate explicitly. These examples offer insights to the structure of the optimal policy. We also propose an optimal numerical approach to the problem that dicretizes the problem and solves it as a linear programming problem. We also propose a heuristic algorithm for the problem that has the attractive feature that it learns by solving smaller instances. The heuristic leads to a very fast and completely automatic process for solving the problem.

∗Florin Avram, Department of Mathematics, Northeastern University.

†Dimitris Bertsimas, Sloan School of Management and Operations Research Center, MIT, Cambridge, MA 02139. Research partially supported by a Presidential Young Investigator Award DDM-9158118 with matching funds from Draper Laboratory.

‡Operations Research Center, MIT.

1 Introduction

A multiclass queueing network is one that services multiple types of customers which may differ in their arrival processes, service requirements, routes through the network as well as costs per unit of waiting time. The fundamental optimization problem that arises in open networks is to determine an optimal policy for sequencing and routing customers in the network that minimizes a linear combination of the expected sojourn times of each customer class. There are both sequencing and routing decisions involved in these optimization problems. A sequencing policy determines which type of customer to serve at each station of the network, while a routing policy determines the route of each customer.

There are several important applications of such problems: packet-switching communication networks with different types of packets and priorities, job shop manufacturing systems, scheduling of multi-processors and multi-programmed computer systems, to name a few.

The control of multiclass queueing networks is a mathematically challenging problem. In order to achieve optimality, stations have to decide how to sequence competing customer types at each point in time, based on information about the load conditions of various other stations. These interactions between various stations create serious dependencies among them and prevent not only optimization but even performance analysis of a given policy. To indicate the difficulty of the problem, it is worth mentioning that even with Poisson arrivals and class dependent exponential service times, and for the simplest possible policy, FCFS, product form or analytical solutions are not available. Naturally, optimizing a multiclass queueing network is an even harder problem. Thus, not surprisingly, simulation is the most common practice among researchers and practitioners as a tool of evaluating heuristic policies.

Optimizing a multiclass queueing network has three characteristics that add to the complexity of the problem: a) The problem is dynamic, i.e., sequencing decisions need to be taken dynamically in time, b) the problem is combinatorial, i.e., there are many possible discrete solutions and c) the problem is stochastic, i.e., there are uncertainties regarding the arrival and service processes.

This paper takes the philosophical view that the most important components of the problem are the first two, i.e., the dynamic and combinatorial character of the problem. In other words, while the stochastic component of the problem is important, the structure of the optimal policies can be revealed by a deterministic and dynamic model. In particular we study in this paper a deterministic fluid model for multiclass networks that has also been studied by Chen and Yao [6], who propose a suboptimal myopic solution using a sequence of linear programing problems (see also Atkins and Chen [3]). Our contibutions in the present paper are as follows:

1. Since mathematically, the problem is a linear control problem over a polyhedral region, we use Pontryagin’s maximum principle to solve specific examples optimally. These examples offer significant qualitative insights as they suggest that the optimal policy is of the threshold type. We calculate exactly the parameters of the optimal policies. In contrast with the work of Chen and Yao [6], our approach offers explicit threshold policies. Several of the examples we

solve have not been solved previously in the literature.

2. We also view the problem as an infinite linear programming problem (see Anderson and Nash [2]). We solve then the problem exactly based on a discretization scheme. Although this approach has dimensionality problems and only offers numerical insights, it provides us with a yardstick for evaluating the performance of heuristic policies.

3. Based on the exact solution of simpler cases we propose a learning heuristic that solves the problem both symbolically (thus offering qualitative insight on near optimal policies) and numerically extremely fast. Comparing the exact and approximate solutions we found that the current implementation of the heuristic provides very accurate policies that are close (and often identical) to the optimal ones even for large scale problems. An attractive feature of the heuristic is that in principle it has a learning capability, in the sense that the more problems one can solve optimally, the better the heuristic can be (we have not, however, implemented this aspect of the heuristic).

Literature Review Optimal control approaches to fluid models in the context of routing in communication systems have been studied by several researchers. For an informative review of this work see Filipiak [7]. Local optimal solutions can be found by gradient search techniques or global optimal solutions for specific problems may sometimes be found by exploiting the structure of the problems at hand. Solution techniques of both types can be found in [12], [10], [17], [11] or [5].

Multiclass queueing networks operating in heavy traffic conditions have been successfully modeled via Brownian motions. Introduced by Harrison [?] and further explored by Wein, this approach proposes heuristic policies which typically outperform more traditional ones. This approach has been more successful in closed networks [?] and networks with controllable input [?, ?], but has not been as successful in scheduling open networks. In particular, Harrison and Wein show in [?] that a threshold policy is consistent with the optimality conditions for a Brownian two-station, three-class network which we also consider in this paper (we explicitly characterize the optimal policy, which is of hte threshold type in some scenaria). Wein [?, ?] proposes priority rules and admission control policies in open networks where admission control is allowed. For a nice survey of the heavy-traffic approach for optimization of multiclass networks the reader, is referred to Kelly and Laws [?].

In studies that concern lower bounds for general networks, Ou and Wein [?] derive pathwise lower bounds for general open queueing networks with deterministic routing. Bertsimas, Paschalidis and Tsitsiklis [4] (see also Kumar and Kumar [?]) propose linear and nonlinear programming problems that provide lower bounds for the optimal solution values. Structure of the paper The rest of the paper is organized as follows: In Section 2, we formally define the sequencing problem for multiclass open networks, describe the mathematical model, briefly review the maximum principle, and describe the results we obtained

[figure: λ_{1}; λ_{3}; μ_{1}; μ_{3}; μ_{2}]

Figure 1: An example queueing network

in solving exactly some simple examples of queueing networks, with the goal of obtaining insights on the structure of the optimal policy. In Section 3 we propose an exact optimal algorithm that solves the problem numerically. In Section 4 we propose a heuristic algorithm for the problem that solves the problem very fast, both analytically and numerically. In Section 5 we report extensive computational results with the goal of obtaining insights into the structure of the optimal policy and the performance of the heuristic. The final section contains some concluding remarks.

2 Problem Formulation

Before presenting the general problem formulation we examine an example of a queueing network composed of three classes and two stations; classes one and two are served at stations one and two respectively while class three competes with class one for service at station one (see Figure 1). Jobs for classes one and three arrive at rates λ_{1} and λ_{3} per time period and jobs of class i are processed at rate μ_{i} per time period.

Assuming that there are jobs in the system for all three classes, the decision of the problem is if station one should service class one or three (it is clear that neither station will be idle). This decision will be represented by a control variable that is associated with each class. The control variable, u_{i} denotes the fraction of time that the system is servicing jobs in class i. In this example, it might be optimal to initially service class three at station one if either the initial queue for the second station is large or it is very expensive to keep jobs in the second queue.

This problem can be formulated as follows: For i = 1, 2, 3, let x_{i}(t) denote the number of jobs of class i at time t. The cost per unit time for holding a job of class i is c_{i}. Finally, the control variables, u_{i}(t), will denote the fraction of effort being spent processing jobs of class i at time t (0 ≤ u_{i}(t) ≤ 1). Let x(t), u(t), c denote the corresponding vectors. Let T a large enough time so the system will empty btime T. Let x̃ be the given vector of the number of customers at time 0. The control problem is:

min  ∫_{0}^{T} c′ x(t)dt
ẋ_{1}(t) = λ_{1} − u_{1}(t)μ_{1}
ẋ_{2}(t) = u_{1}(t)μ_{1} − u_{2}(t)μ_{2}
ẋ_{3}(t) = λ_{3} − u_{3}(t)μ_{3}
x(0) = x̃
u_{1}(t) + u_{3}(t) ≤ 1
u_{2}(t) ≤ 1
x(t) ≥ 0
u(t) ≥ 0.

Our primary goal is to provide insight to the question of what the optimal policy u(0) is as a function of the initial state vector x̃. In other words we would like to find u(0) for all x̃.

A general queueing network will have m stations and n job classes (n ≥ m). A job class is specific to a station; as a particular job flows through the network, it defines a different class in each station. A single station may service multiple job classes, each with its own service rate. The arrivals for each job class either come from another station or from outside the system. Each class has a control variable associated with it specifying the fraction of time that the station spends processing the jobs in this class. As the jobs are processed, they either become members of a new class or leave the system. The data of the problem is as follows: A vector b of external arrivals, a vector μ of service rates and the routing is given.

The relationship between classes and stations can be represented in an n × n matrix, A. Each row of the matrix represents a class of jobs and each column represents a control in the system. Since there is exactly one control for each class, assume that control i is used for class i. The i^{th} row of matrix A contains information for class i. Positive elements represent arrivals for class i and there is one negative element specifying the service rate for the class. Specifically, the elements of the matrix A are nonzero in the following cases: A_{ij} = μ_{j}, if jobs of class j become jobs of class i, A_{ii} = −μ_{i}.

Using the same notation as the example above [e.g. x_{i}(t) is the number of jobs of class i in the system at time t, and u_{i}(t) is the fraction of time spent processing class i jobs at time t], the differential equations for x(t) can be written:

ẋ(t) = Au + b.

In the example above,

b = (λ_{1}, 0, λ_{3})^{t}

A = [ \begin{matrix}−μ_{1} & 0 & 0 \\ μ_{1} & −μ_{2} & 0 \\ 0 & & −μ_{3} \end{matrix} ]

The control variables u_{i} represent the fraction of effort devoted to the i’th class by the station where it is processed.Thus, the sum of the control variables for all the classes processed at the same station must be less than one. This can be represented with an m × n binary matrix, D, where

D_{ij} = \{ \begin{matrix}1 & if station i services class j \\ 0 & otherwise. \end{matrix}

This constraint can then be written as

Du(t) ≤ e.

Then, the general problem is

(CONT ROL) min  ∫_{0}^{T} c′ x(t)dt  (1)
ẋ(t) = Au(t) + b  (2)
Du(t) ≤ e  (3)
x(0) = x̃  (4)
x(t) ≥ 0  (5)
u(t) ≥ 0.  (6)

Problem (CONT ROL) can be viewed as a linear program in infinite dimensions. The number of rows is finite, but the number of variables is infinite, since we have an x(t) vector for each infinitesimal time period [t, t + Δ]. Such problems could be solved by using a discrete approximation and solving a linear program in a finite (albeit large) number of variables; this approach will be discussed in section 4.

Problem (CONT ROL) can also be found in the field of calculus of variations that offers techniques such as the Pontryagin maximum principle to yield necessary conditions for an optimal solution. These techniques provide insights on the structure of the optimal policies, and lead sometimes also to explicit solutions. We will solve below certain examples by the maximum principle. But first we will give a ”myopic” solution to the problem (CONT ROL), similar to that of Chen and Yao [?] for comparison. Since the instantaneous costs incurred at a fixed time t are c′ x(t), a myopic policy will choose u(t) to minimize c′ ẋ(t). Using the special structure of the matrix A yields the problem

max ∑(f_{i}u_{i})

Du ≤ e

u ≥ 0

where f_{i} = μ_{i} (c_{i} − c_{n}),where n = n(i) is the class following class i (and c_{n} = 0 if i is the last class on the route). The solution of the above problem is simply to serve at each station among the classes with positive f_{i} that one for which f_{i} is maximal. The myopic solution turns out to be optimal in the case of one server networks. However, in general it can miss important features of the optimal solution.

2.1 Review of the Pontryagin Maximum Principle

In this subsection, the Pontryagin Maximum Principle [15] is reviewed. The maximum principle provides necessary conditions for an optimal trajectory and its associated control. Define the problem as

min  ∫_{0}^{T} c′ x(t)dt
ẋ(t) = Au(t) + b
Du(t) ≤ e
x(t) ∈ X
u(t) ≥ 0
x(0) = x̃

where x̃ is the initial condition and X = {x : x ≥ 0}. We define the Hamiltonian function as

H(x, u, y, t) = c′ x(t) + y(t)′ [Au(t) + b].

, where the multipliers y(t) satisfy the ”adjoint” system

\frac{dy_{i}}{dt} = −\frac{dH}{dx_{i}} = −c_{i}

, with certain end conditions called transversality conditions. Then, it turns out that the segment of the optimal trajectory that lies in the interior of X may be found by minimizing for each fixed t the Hamiltonian with respect to u.This may again be carried out explicitly, like in the myopic solution. The optimal policy is to serve at each station the class with maximum index P_{i}, where P_{i}(t) = μ_{i}(y_{i} − y_{n}) (with n denoting the class following i) is this time however a dynamic index (changing with time). Actually, the constraints x_{i} ≥ 0 may in general complicate considerably the solution of the unconstrained problem. A further simplification may be obtained here noting that under the additional assumption b = 0 (no arrivals from the outside) , once a quantity i is depleted at t_{i},it cannot be optimal to allow it to become positive again. This extra constraint determines u_{i}(t) for t > t_{i};also, the cofactor y_{i}(t) becomes unnecessary for t > t_{i}, and this reduces the constrained problem to a sequence of unconstrained problems of decreasing dimension. The connection between the problems is given by end conditions for the cofactors, which are obtained from the fact that along the optimal trajectory x̄ and ȳ, ū the Hamiltonian must satisfy

max_{u} H(x̄, ȳ, u, t) = H(x̄, ȳ, ū, t) = 0

This ”continuity equation” provides us alsowith the general form of the equation of a ”switch” curve, by which we mean the locus of points where a certain fraction d of servicing effort is transferred from a class i to a class c = c(i) (c for competitor). Since of d(P_{i}−P_{c}) in the Hamiltonian accompanies the transfer of effort between the classes i and c, but the Hamiltonian is continuous, it follows that a SWITCH at t may occur iff P_{i}(t) = P_{c}(t). This equation must hold both at forced switches, called from now on DEPLETIONS (when a class is depleted), and at switches occurring in the interior of the state space, called from now on INTERRUPTS , and provides us both with the end conditions necessary to solve the y multipliers, and with the complete solution of the control problem, provided that we can guess the time order of the INTERRUPTS and DEPLETIONS. WE illustrate this in a couple of examples below. Meanwhile, we note some general facts implied by the maximum principle.

Proposition:a) THe cofactors y_{i}(t) are decreasing linear functions in their domain t ≤ t_{i}. b) The dynamic indices P_{i}(t) are piecewise linear functions with at most N_{i} + 1 changes of slope, where N_{i} denotes the number of stations following i on its route (thus , the next of i, n_{i}, may go through at most N_{i} + 1 changes).

We outline now the results , obtained in the appendix by the above method, for two particular networks:

1)Two queues in tandem. The Myopic solution would be to serve at the first station iff c_{1} ≥ c_{2}, and to never serve (with disastruous long term consequences) if c_{1} < c_{2}) . However, the true solution in this last case turns out to be a HOLD policy: hold x_{1} until\frac{x̃_{1}}{x̃_{2}} < \frac{c_{2}}{c_{1}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}} The equation for the hold point is obtained by determining the two points where the dynamic index P_{1}(t) becomes positive (the INTERRUPT), and ceases to be positive (the DEPLETION. This geometry requires that P_{1} has a ”tooth” shape (i.e. its first slope , c_{2} − c_{1} has to be positive, and also t_{2} < t_{1}, since t_{2}, the depletion of x_{2} and the ”peak of the tooth”, has to precede t_{1}, the ”end of the tooth”.

2)The fork intersection with K routes. In this case, all the cofactors of the classes served at the first station may be ”tooth” shaped. Thus, the optimal solution, obtained by figuring out which of the K ”teeth ” is currently higher, may involve up to 2K − 2 ”INTERRUPTS”, which occur at the points where one of the ”teeth” ”rises” above the others. The occurrence of the ”teeth” is caused, like in the previous example, when the time of depletion of n(i) precedes that of i.

Chen and Yao have pointed out already that the myopic policy may be changed adaptively, as soon as certain classes are depleted.In the above examples we see moreover that the correct optimal solution ”foresees” future depletions on the route, and changes before the depletions actually occur.

Finally, we provide the complete results obtained in the appendix for a particular type of fork intersection, the example from the beginning of this section.

Theorem 1 The optimal control for the problem at t = 0 is given by the following five cases.

Conditions  u(0)
Case 1  c_{1}μ_{1} < c_{3}μ_{3}  (0,1,1)
Case 2  c_{1}μ_{1} > c_{3}μ_{3}  (1,1,0)
c_{1}μ_{1} − c_{3}μ_{3} > c_{2}μ_{1}
Case 3 c_{1}μ_{1} > c_{3}μ_{3} (0,1,1)
c_{1}μ_{1} − c_{3}μ_{3} < c_{2}μ_{1}
μ_{1} > μ_{2}
Case 4 c_{1}μ_{1} > c_{3}μ_{3} (1,1,0)
c_{1}μ_{1} − c_{3}μ_{3} < c_{2}μ_{1}
μ_{1} < μ_{2}
\frac{x̃_{1}}{x̃_{2}} > \frac{c_{2}μ_{1}}{c_{1}μ_{1}−c_{3}μ_{3}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}}
Case 5 c_{1}μ_{1} > c_{3}μ_{3} (0,1,1)
c_{1}μ_{1} − c_{3}μ_{3} < c_{2}μ_{1}
μ_{1} < μ_{2}
\frac{x̃_{1}}{x̃_{2}} < \frac{c_{2}μ_{1}}{c_{1}μ_{1}−c_{3}μ_{3}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}}

3 Insights From Examples

In this section we find the optimal policy for specific examples using the maximum principle. Our goal is to offer insights into the structure of the optimal policy. All of the problems considered assume that the system starts with jobs of every class in the system (x̃ > 0).

3.1 Competing Queues

Consider the problem in Figure 2 with a single server and n classes. Jobs for class i arrive at a rate of λ_{i} per unit time and can be serviced at a rate μ_{i} per unit time. This problem is solved by the well known “μc” rule, but it will be illustrative to derive the “μc” rule using the maximum principle.

[figure: λ_{1}; μ_{1}; λ_{2}; μ_{3}; ⋮; λ_{n}; μ_{n}]

Figure 2: Multiple classes being served by one station

The problem is written as:

min  ∫_{0}^{T} c′ x(t)dt
ẋ_{1}(t) = λ_{1} − u_{1}(t)μ_{1}
⋮

ẋ_{n}(t) = λ_{n} − u_{n}(t)μ_{n}
x(0) = x̃
∑_{i}^{n} u_{i}(t) ≤ 1
x(t) ≥ 0
u(t) ≥ 0

The Hamiltonian function is

H = y_{0}c′ x(t) + ∑_{i}^{n} y_{i}(t)(λ_{i} − μ_{i}u_{i})

= y_{0}c′ x(t) − ∑_{i}^{n} μ_{i}y_{i}u_{i} + ∑_{i}^{n} y_{i}λ_{i}

At time t = 0, class i will be serviced if μ_{i}y_{i} = min_{1≤j≤n} μ_{j}y_{j}. Class i will be served until x_{i}(t) = 0 or μ_{i}y_{i}(t) = μ_{j}y_{j}(t) for some j.

Theorem 2 The optimal policy is u_{K}(0) = 1 and u_{i}(0) = 0 if i ≠ K where K = arg max_{1≤j≤n} c_{j}μ_{j}.

Proof It is clear from the maximum principle that exactly one u_{i}(0) will equal 1 and the rest 0. Define P_{i}(t) = −μ_{i}y_{i}(t) which gives Ṗ_{i}(t) = μ_{i}c_{i}y_{0}. All of the derivatives Ṗ_{i}(t) are negative and Ṗ_{K}(t) is the most negative. Since P_{i}(T) = 0, P_{K}(0) = max_{1≤j≤n} P_{i}(0) (if there are no jumps). The sufficient conditions of the maximum principle are satisfied so this is also optimal. □ This example supports the philosophical view of the paper that it is the dynamic and combinatorial character of the problem and not the stochastic character that is the most important. In this example the same policy is optimal for the stochastic model as well.

3.2 Tandem Queues

Consider the queueing network in Figure 3 with two stations in tandem. Jobs enter the system with rate λ and are serviced at station one at rate μ_{1}. All jobs then pass to station two where they are served at rate μ_{2} and then exit the system.

[figure: λ; μ_{1}; μ_{2}]

Figure 3: Two stations in tandem

The problem is written as:

min  ∫_{0}^{T} c′ x(t)dt

ẋ_{1}(t) = λ − u_{1}(t)μ_{1}
ẋ_{2}(t) = u_{1}(t)μ_{1} − u_{2}(t)μ_{2}
x(0) = x̃
u_{i}(t) ≤ 1
x(t) ≥ 0
u_{i}(t) ≥ 0

It is clear that λ < μ_{1} and λ < μ_{2} are required for stability and that u_{2}(t) = 1 as long as x_{2}(t) > 0. The problem is concerned with finding the optimal u_{1}(0) as a function of c, λ, μ and x̃. The maximum principle offers insights in some cases and leaves only one case to be solved analytically.

Theorem 3 The optimal control for the tandem problem at t = 0 is given by the following four cases.

Conditions u(0)
Case 1  c_{1} > c_{2}  (1,1)
Case 2  c_{1} < c_{2}  (0,1)
μ_{1} > μ_{2}
Case 3 c_{1} < c_{2} (0,1)
μ_{1} < μ_{2}
\frac{x̃_{1}}{x̃_{2}} < \frac{c_{2}}{c_{1}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}}
Case 4 c_{1} < c_{2} (1,1)
μ_{1} < μ_{2}
\frac{x̃_{1}}{x̃_{2}} > \frac{c_{2}}{c_{1}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}}

Proof The hamiltonian function for this problem is

H = y_{0}[c_{1}x̃_{1} + c_{2}x̃_{2}] + y_{1}(t)[λ − u_{1}μ_{1}] + y_{2}(t)[u_{1}μ_{1} − u_{2}μ_{2}]
= y_{0}[c_{1}x̃_{1} + c_{2}x̃_{2}] + u_{1}[y_{2}(t)μ_{1} − y_{1}(t)μ_{1}] + y_{1}(t)λ − u_{2}[y_{2}(t)μ_{2}]

Case 1

While x(t) > 0, u_{1}(t) = 0 if and only if y_{2}(t) − y_{1}(t) + ty_{0}(c_{1} − c_{2}) ≤ 0, since u must maximize H. If u_{1}(0) = 0 and c_{1} > c_{2}, then u_{1}(t) = 0 for at least as long as x_{2}(t) > 0. During this time the expensive, class one jobs are stored at station one. It would clearly cost less if as many of these jobs as possible were stored at station two, where the holding cost is cheaper.

Case 2

Assume that u(0) = (1, 1). Since u must maximize H, this implies that y_{1}(0) ≤ y_{2}(0) ≤ 0. Since c_{1} < c_{2}, y_{1}(t) ≤ y_{2}(t) so u_{1}(t) = 1 at least as long as x(t) > 0. Using the fluid model, it the jobs could be stored at station 1 until x_{2}(t) = 0 and then immediately have station 2 working at full capacity since μ_{1} > μ_{2}. It would clearly be cheaper to use this second alternative.

[figure: x_{2}; x_{1}; ẋ = (λ − μ_{1}, μ_{1} − μ_{2}); ẋ = (λ, −μ_{2})]

Figure 4: Switching curve for case 3

Case 3

This case is the most difficult. The maximum principle shows that while x(t) > 0, there can be at most one switch in u_{1}. Geometrically, the state space is partitioned by a single line (Figure 4). This makes it possible to find the optimal switching time directly. If u_{1}(0) = 0, then it switches to 1 at some time, t, when y_{1}(t) = y_{2}(t) and then remains 1 until x(t) is no longer positive. Assume u(0) = (0, 1) (y_{2}(0) ≤ y_{1}(0) and y_{2}(0) ≤ 0) and c_{1} < c_{2} and μ_{1} < μ_{2}. It is possible to find the relationship between x_{1} and x_{2} when u_{1} is 0 or 1 by using ideas from [11]. When u_{1} = 1, ẋ_{1} = λ − μ_{1} and ẋ_{2} = μ_{1} − μ_{2}

x_{1}(t) = (λ − μ_{1})t + x_{1}(0)

x_{2}(t) = (μ_{1} − μ_{2})t + x_{2}(0)

x_{2}(t) = \frac{μ_{1} − μ_{2}}{λ − μ_{1}} x_{1}(t) + K

where K = x_{2}(0) − x_{1}(0)\frac{μ_{1}−μ_{2}}{λ−μ_{1}} . Similarly, when u_{1} = 0, ẋ_{1} = λ and ẋ_{2} = −μ_{2} or

x_{2}(t) = \frac{−μ_{2}}{λ} x_{1}(t) + K̄

where K̄ = x_{2}(0) + x_{1}(0)\frac{μ_{2}}{λ}.

It is now easy to compute the total cost of starting with u = (1, 1). This cost can then be compared with the total cost of starting with u = (0, 1) and switching to u_{1} = 1 after a some time. If the second cost is lower, it is better to hold station one.

The cost incurred without holding station 1 will be computed in two stages. In the first stage, both x_{1}(t) and x_{2}(t) are positive. In the second stage, x_{2}(t) = 0 and x_{1}(t) > 0. In the first stage, u(0) = (1, 1).

In the integral below, t̄ is the time when x_{2}(t) first becomes 0:

t̄= \frac{x̃_{2}}{μ_{2} − μ_{1}}

∫ t̄ ; 0[c_{1}x_{1}(t) + c_{2}x_{2}(t)] dt

= [\frac{c_{1}}{2} (λ − μ_{1}) + \frac{c_{2}}{2} (μ_{1} − μ_{2})] t̄^{2} + [c_{1}x̃_{1} + c_{2}x̃_{2}] t̄

The cost for stage two will be considered next. After time t̄, the control is u = (1, μ_{1}/μ_{2}). The system will be able to maintain x_{2}(t) = 0 since μ_{1} < μ_{2} so it will only accrue charges for this stage until x_{1}(t) = 0 at time t = x̃_{1}/(μ_{1} − λ). Define

T_{1} ≡ \frac{x̃_{1}}{μ_{1} − λ} − \frac{x̃_{2}}{μ_{2} − μ_{1}}

The cost for stage two is

∫ T_{1}+t̄ ; t̄c_{1}x_{1}(t)dt

= ∫ T_{1}+t̄ ; t̄[c_{1}(λ − μ_{1})t + c_{1}x̃_{1}] dt
= \frac{c_{1}x̃_{2}^{2}(μ_{1} − λ)}{2(μ_{2} − μ_{1})^{2}} − \frac{c_{1}x̃_{1}x̃_{2}}{μ_{2} − μ_{1}} + \frac{c_{1}x̃_{1}^{2}}{2(μ_{1} − λ)}

Therefore, if the system does not hold station one, the total cost is:

[\frac{c_{1}}{2} (λ − μ_{1}) − \frac{c_{2}}{2} (μ_{2} − μ_{1})] \frac{x̃_{2}^{2}}{(μ_{2} − μ_{1})^{2}} + [c_{1}x̃_{1} + c_{2}x̃_{2}] \frac{x̃_{2}}{μ_{2} − μ_{1}} +

\frac{c_{1}x̃_{2}^{2}(μ_{1} − λ)}{2(μ_{2} − μ_{1})^{2}} − \frac{c_{1}x̃_{1}x̃_{2}}{μ_{2} − μ_{1}} + \frac{c_{1}x̃_{1}^{2}}{2(μ_{1} − λ)}

= \frac{c_{2}x̃_{2}^{2}}{2(μ_{2} − μ_{1})} + \frac{c_{1}x̃_{1}^{2}}{2(μ_{1} − λ)}
≡ F_{N}(x̃_{1}, x̃_{2})

The cost of holding station one must now be computed to compare with F_{N}. The cost will be the cost of holding station one for some small time, Δ and then switching to u = (1, 1) (which carries a cost of F_{N}(x(Δ)). First, compute the cost for 0 ≤ t ≤ Δ (u = (0, 1)),

∫_{0}^{Δ}(c_{1}x_{1}(t) + c_{2}x_{2}(t))dt
= \frac{Δ^{2}}{2} (c_{1}λ − c_{2}μ_{2}) + [c_{1}x̃_{1} + c_{2}x̃_{2}]Δ

And, at t = Δ,

x_{1}(Δ) = λΔ + x̃_{1}

and

x_{2}(Δ) = x̃_{2} − μ_{2}Δ

If Δ is small enough so x(Δ) > 0, then the cost for time t > Δ is F_{N}(x_{1}(Δ), x_{2}(Δ)). It is now possible to compare the cost of not holding at all with the cost of holding for some small time period. It will be beneficial to hold if

\frac{Δ^{2}}{2} (c_{1}λ − c_{2}μ_{2}) + (c_{1}x̃_{1} + c_{2}x̃_{2})Δ + F_{N}(x_{1}(Δ), x_{2}(Δ)) − F_{N}(x̃_{1}, x̃_{2})

is negative.

\frac{Δ^{2}}{2} (c_{1}λ − c_{2}μ_{2}) + (c_{1}x̃_{1} + c_{2}x̃_{2})Δ + F_{N}(x_{1}(Δ), x_{2}(Δ)) − F_{N}(x̃_{1}, x̃_{2})

= \frac{Δ^{2}}{2} [c_{1}λ (1 + \frac{λ}{μ_{1} + λ}) + c_{2}μ_{2} (\frac{μ_{2}}{μ_{2} − μ_{1}} − 1)] + Δ [c_{1}x̃_{1} + c_{2}x̃_{2} − \frac{c_{2}μ_{2}x̃_{2}}{μ_{2} − μ_{1}} + \frac{c_{1}λx̃_{1}}{μ_{1} − λ}]
= Δ^{2}A + ΔB

Since A > 0, then Δ^{2}A + ΔB < 0 if and only if 0 < Δ < −B/A. Since Δ > 0, it will be beneficial to hold only if B < 0 The expression for B can be rewritten as

[c_{1}x̃_{1} + c_{2}x̃_{2} − \frac{c_{2}μ_{2}x̃_{2}}{μ_{2} − μ_{1}} + \frac{c_{1}λx̃_{1}}{μ_{1} − λ}] = x̃_{1} [c_{1} + \frac{c_{1}λ}{μ_{1} − λ}] − x̃_{2} [−c_{2} + \frac{c_{2}μ_{2}}{μ_{2} − μ_{1}} ]

Therefore, it is optimal to hold when

x̃_{1} [c_{1} + \frac{c_{1}λ}{μ_{1} − λ}] < x̃_{2} [−c_{2} + \frac{c_{2}μ_{2}}{μ_{2} − μ_{1}} ]

or

\frac{x̃_{1}}{x̃_{2}} < \frac{c_{2}}{c_{1}} \frac{μ_{1} − λ}{μ_{2} − μ_{1}}

This completes the proof of this case.

Case 4

This case is the complement of case 3 so the proof is very similar to the derivation of that case. □ Note that cases 3 and 4 exhibit an optimal policyof the threshold type, in which we have explicitly computed the threshold curve.

3.3 Tandem Queues with Competition in the First Node

Consider the queueing network in Figure 5 with two stations in tandem and two queues competing for service in the first node. Class one jobs enter the system with rate λ_{1} and are serviced at station one at rate μ_{1}. All class one jobs then pass to station two where they are served at rate μ_{2} and then exit the system. Class three

[figure: λ_{1}; λ_{3}; μ_{1}; μ_{3}; μ_{2}]

Figure 5: Tandem system with competition

jobs enter the system with rate λ_{3} at station one and are serviced at rate μ_{3}. The problem is written as:

min  ∫_{0}^{T} c′ x(t)dt
ẋ_{1}(t) = λ_{1} − u_{1}(t)μ_{1}
ẋ_{2}(t) = u_{1}(t)μ_{1} − u_{2}(t)μ_{2}
ẋ_{3}(t) = λ_{3} − u_{3}(t)μ_{3}
x(0) = x̃
u_{1}(t) + u_{3}(t) ≤ 1
u_{2}(t) ≤ 1
x(t) ≥ 0
u_{i}(t) ≥ 0

It is clear that λ_{1}/μ_{1}+λ_{3}/μ_{3} < 1 and λ_{1} < μ_{2} are required for stability, u_{2}(t) = 1 as long as x_{2}(t) > 0 and as long as x_{3}(t) > 0, u_{3}(t) = 1 − u_{1}(t). The problem is to choose between u(0) = (1, 1, 0) and u(0) = (0, 1, 1) This problem is conceptually a combination of the previous two networks studied since x_{1} and x_{3} compete at station one and x_{1} and x_{2} are related in a tandem system.

Theorem 4 The optimal control for the problem at t = 0 is given by the following five cases.

Conditions  u(0)
Case 1  c_{1}μ_{1} < c_{3}μ_{3}  (0,1,1)
Case 2  c_{1}μ_{1} > c_{3}μ_{3}  (1,1,0)
c_{1}μ_{1} − c_{3}μ_{3} > c_{2}μ_{1}
Case 3 c_{1}μ_{1} > c_{3}μ_{3} (0,1,1)
c_{1}μ_{1} − c_{3}μ_{3} < c_{2}μ_{1}
μ_{1} > μ_{2}
Case 4 c_{1}μ_{1} > c_{3}μ_{3} (1,1,0)
c_{1}μ_{1} − c_{3}μ_{3} < c_{2}μ_{1}
μ_{1} < μ_{2}
\frac{x̃_{1}}{x̃_{2}} > \frac{c_{2}μ_{1}}{c_{1}μ_{1}−c_{3}μ_{3}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}}
Case 5 c_{1}μ_{1} > c_{3}μ_{3} (0,1,1)
c_{1}μ_{1} − c_{3}μ_{3} < c_{2}μ_{1}
μ_{1} < μ_{2}
\frac{x̃_{1}}{x̃_{2}} < \frac{c_{2}μ_{1}}{c_{1}μ_{1}−c_{3}μ_{3}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}}

Proof The hamiltonian function for this problem is

H = y_{0}[c_{1}x̃_{1} + c_{2}x̃_{2} + c_{3}x̃_{3}] + y_{1}(t)[λ_{1} − u_{1}μ_{1}] + y_{2}(t)[u_{1}μ_{1} − u_{2}μ_{2}] + y_{3}(t)[λ_{3} − u_{3}μ_{3}]
= y_{0}[c_{1}x̃_{1} + c_{2}x̃_{2} + c_{3}x̃_{3}] + u_{1}[y_{2}(t)μ_{1} − y_{1}(t)μ_{1}] − u_{2}[y_{2}(t)μ_{2}] − u_{3}[y_{3}(t)μ_{3}]
+y_{1}(t)λ_{1} + y_{3}(t)λ_{3}

For the first case, note that if the system only consisted of classes one and three, priority would clearly go to class three and class one would hold. The current system only modifies the path of class one, so there is no reason to change the policy at station one.

For the remaining cases, note that class one will be processed if and only if μ_{1}(y_{2} − y_{1}) > −μ_{3}y_{3} or y_{2}μ_{1} > y_{1}μ_{1} − y_{3}μ_{3}. Compare this with the two class tandem where class one was processed if and only if y_{2} > y_{1}. This, along with the fact that class three is processed whenever class one is not, means that this system is a “modified tandem” where the cost used for the first station is c_{1}μ_{1} − c_{3}μ_{3} while the cost for the second station is c_{2}μ_{1}. The second through fifth cases now follow. □

3.4 Tandem Queues with Competition at the Second Node

Consider the queueing network in Figure 6 with two stations in tandem and two queues competing for service in the second node. Class one jobs enter the system with rate λ_{1} and are serviced at station one at rate μ_{1}. All class one jobs then pass to station two where they are served at rate μ_{2} and then exit the system. Class three jobs enter the system with rate λ_{3} at station two and are serviced at rate μ_{3}.

[figure: λ_{1}; λ_{3}; μ_{2}; μ_{3}; μ_{1}]

Figure 6: Tandem system with competition

The problem is written as:

min  ∫_{0}^{T} c′ x(t)dt
ẋ_{1}(t) = λ_{1} − u_{1}(t)μ_{1}
ẋ_{2}(t) = u_{1}(t)μ_{1} − u_{2}(t)μ_{2}
ẋ_{3}(t) = λ_{3} − u_{3}(t)μ_{3}
x(0) = x̃
u_{1}(t) ≤ 1
u_{2}(t) + u_{3}(t) ≤ 1
x(t) ≥ 0
u(t) ≥ 0

It is clear that λ_{1} < μ_{1} and λ_{1}/μ_{2} + λ_{3}/μ_{3} < 1are requi red for stability and as long as x_{3}(t) > 0, u_{3}(t) = 1 − u_{2}(t). The problem is to decide if class one should be held or not and if station two should process class two or three. Specifically, the four possible choices for u(0) are (1,1,0), (1,0,1), (0,1,0) and (0,0,1). This problem is also a combination of the competition and tandem problems since x_{1} and x_{2} are in tandem and x_{2} competes with x_{3}. The proposed optimal policy is presented in theorem 5 without proof.

Theorem 5 The optimal control for the problem at t = 0 is given by the following eight cases.

Conditions  u(0)
Case 1 c_{1} > c_{2} (1,1,0)
c_{2}μ_{2} > c_{3}μ_{3}
Case 2 c_{1} > c_{2} (1,0,1)
c_{3}μ_{3} > c_{2}μ_{2}
Case 3 c_{1} < c_{2} (0,0,1)
μ_{1} > μ_{2}
c_{3}μ_{3} > c_{2}μ_{2}
Case 4 c_{1} < c_{2} (0,1,0)
μ_{1} > μ_{2}
c_{2}μ_{2} > c_{3}μ_{3}
Case 5 c_{1} < c_{2} (1,1,0)
μ_{1} < μ_{2}
c_{2}μ_{2} > c_{3}μ_{3}
\frac{x̃_{1}}{x̃_{2}} > \frac{c_{2}}{c_{1}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}}
Case 6 c_{1} < c_{2} (0,1,0)
μ_{1} < μ_{2}
c_{2}μ_{2} > c_{3}μ_{3}
\frac{x̃_{1}}{x̃_{2}} < \frac{c_{2}}{c_{1}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}}
Case 7 c_{1} < c_{2} (1,0,1)
μ_{1} < μ_{2}
c_{3}μ_{3} > c_{2}μ_{2}
\frac{x̃_{1}}{x̃_{2}} > \frac{c_{2}}{c_{1}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}}
Case 8 c_{1} < c_{2} (0,0,1)
μ_{1} < μ_{2}
c_{3}μ_{3} > c_{2}μ_{2}
\frac{x̃_{1}}{x̃_{2}} < \frac{c_{2}}{c_{1}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}}

4 A Numerical Approach to the General Problem

As mentioned in Section 2, the control problem can be viewed as an infinite dimensional linear program. Such problems have been studied by many different researchers, yet there does not exist an effcient algorithm for the general infinite dimensional LP [2]. Some results are known for specific problem forms, one of which is the problem at hand. The scheduling problem being considered is known in the literature as a separated continuous linear program (SCLP). Anderson discusses basic solutions for this problem in [1] and duality theorems for this problem are presented in [2],[?] and [16]. Pullan presents an iterative algorithm to solve a SCLP in [16]. The algorithm uses discrete approximations to find an upper bound and then generates an improvement step. This section presents a numerical method to generate exact solutions for the problem that also uses discrete approximation. As opposed to using an improvement step to improve the approximation, the discretization is continually refined until the solution no longer improves.

Recall that the problem at hand is written

min  ∫_{0}^{T} c′ x(t)dt
ẋ(t) = Au + b
Du ≤ e
x(0) = x̃
x(t) ≥ 0
u(t) ≥ 0

where A is an n × n matrix and D is an m × n matrix. The problem has a finite number of rows and an infinite number of variables. An approximation to the problem can be obtain by discretizing the time interval and solving the resulting finite linear program. Divide the time interval [0, T] into T/h intervals of length h and only allow the controls to change value at the ends of the intervals. The derivative of the state variables is now written ẋ(j) = Au(t) + b for t = 0 … T/h − 1 and t ≤ j ≤ t + 1. This leads to a linear program with 2nT/h variables:

min  ∑_{t=0}^{T/h} c′ x(t)
x(t) = x(t − 1) + [Au(t − 1) + b]h for t = 1 … T/h
Du(t) ≤ e  for t = 0 … T/h − 1
x(t) ≥ 0  for t = 1 … T/h
u(t) ≥ 0  for t = 0 … T/h − 1
x(0) = x̃

Pullan has shown that the optimal solution to this problem can be used to derive a feasible solution to the original problem [16] and a solution to this problem would clearly give an upper bound for the original problem.

It is intuitively clear that the smaller the stepsize, h, the better the approximation to the original problem. Of course, decreasing h has a sharp tradeoff in problem size. Moreover, it is entirely possible that a moderate stepsize may produce an acceptable approximation for some problems so it would not be worthwhile to use a small stepsize. The algorithm used begins with a moderate stepsize and gets an initial solution. It then keeps dividing the stepsize in half and solving the new problem. The process continues until the optimal objective value decreases less than 5% that of the previous iteration and the optimal u(0) is binary.

While this may appear to be a great deal of computational effort, there is a reason to divide the stepsize in half at each iteration. Given the optimal solution with stepsize h, it is a straightforward process to derive a feasible solution to a problem with stepsize h/2:

x(t + h/2) = x(t) + [Au(t) + b]h/2
u(t + h/2) = u(t)

After the stepsize is divided in half, a new feasible solution is computed so the new problem does not need to be solved from the beginning; the simplex method can begin from the feasible solution in hand.

c = (3, 5), λ = 1, μ = (2, 3)

1 2 3 4 5 6 7 8 9 10  1  2  3  4  5  6  7  8  9  10 X2  X1 •  •  • •  • • •  • • •  • • • •  • • • •  • • • • •  • • • • • • ∘ ∘ ∘ ∘ ∘ ∘ ∘ ∘ ∘ ∘  ∘ ∘ ∘ ∘ ∘ ∘ ∘ ∘ ∘  ∘ ∘ ∘ ∘ ∘ ∘ ∘ ∘ ∘  ∘ ∘ ∘ ∘ ∘ ∘ ∘ ∘  ∘ ∘ ∘ ∘ ∘ ∘ ∘  ∘ ∘ ∘ ∘ ∘ ∘ ∘  ∘ ∘ ∘ ∘ ∘ ∘  ∘ ∘ ∘ ∘ ∘ ∘  ∘ ∘ ∘ ∘ ∘  ∘ ∘ ∘ ∘

∘ = hold, • = serve, — = Actual threshold

Figure 7: Use of discrete LP with tandem problem

As expected, the numerical solutions generated with this method agree with the analytic solutions derived in the previous section. As an example, we solve numerically the two queue tandem problem for multiple values of x(0). The solution given by the LP is compared with the threshold found in theorem 3 in Figure 7. The numerical appraoch correctly reproduces the threshold character of the optimal policy that we have derived analytically in the previous section.

5 A Learning Heuristic that Captures Interactions Among Classes

Given the complicated character of the optimal policy in the examples we have examined and given the dimensionality difficulties that the exact approach of the previous section has, it would be desirable to design an efficient approximate algorithm for the problem. This section develops such an algorithm.

Assuming A is nonsingular, we rewrite equation (2) as

u(t) = A^{−1}(ẋ(t) − b).  (7)

The general problem is then

min  ∫_{0}^{T} c′ x(t)dt
DA^{−1}(ẋ(t) − b) ≤ e
A^{−1}(ẋ(t) − b) ≥ 0
x(0) = x̃
x(t) ≥ 0

We define the (m + n) × n matrix Ā as

Ā = [ \begin{matrix}DA^{−1} \\ −A^{−1} \end{matrix} ]

and the m + n vector b̄ as

b̄ = [ \begin{matrix}e + DA^{−1}b \\ −A^{−1}b \end{matrix} ]

The problem can now be written as

min  ∫_{0}^{T} c′ x(t)dt
Āẋ(t) ≤ b̄
x(0) = x̃
x(t) ≥ 0

The maximum principle says that while x > 0, the optimal control will be an extreme point of {u : Du(t) ≤ e, u ≥ 0}. But these points are also the extreme points of the polyhedron Āẋ(t) ≤ b̄. It is a well known fact from linear programming that any extreme point of Āẋ ≤ b̄ is also the unique optimal solution of

min  w′ ẋ  (8)
Āẋ ≤ b̄

for some appropriately chosen objective function, w which is a function of x̃ (see, for example, [13]). For the remainder of this section, the time parameter will not be specified and t = 0 will be assumed.

The heuristic operates by choosing w and then solving the above linear program to find the optimal ẋ. Given ẋ, it is a simple matter to recover u using (7). We next propose a heuristic waof choosing w that takes into account pairwise interactions between classes. Pairwise interactions In a general network, two classes i and j can only be related in one of three ways: either compete in the same station or are in tandem in consecutive stations or they can not related at all. Having solved the problems of competing classes and two classes in tandem exactly in Sections 3.1 and 3.2 respectively, we can calculate the corresponding w, so that the reformulation (8) gives the exact optimal policy. We then choose the ratio w_{i}/w_{j} in the network to match the corresponding ratio in which two classes compete or they are in tandem respectively.

In this way we form a potentially overdetermined system of equations of the form

Mw = 0
w′ e = 1.

We claculate the least squares solution of the system. Note that software packages like Maple can calculate the solution either numerically or symbolically which might be useful in order to find a policy described in algebraic terms.

Philosophically the heuristic attacks the general problem by learning from special cases. The previous method only takes into account pairwise interactions between classes. Offcourse it is possible to include higher order interactions. Or more generally, when we solve exactly a specific network, we can compute the corresponding w and then look for classes in the network that are related in the same structure as in hte special network, specifying one extra equation. In this respect the heuristic can be seen as a learning heuristic. Given that we have implemented the heuristic using only pairwise interactions we will call the heuristic as: pairwise interaction heuristic.

We could extend the power of the method byconsidering three-way interactions (a three-wayinteraction heuristic) and examine the possible ways three classes can interact: Three classes can compete in the same station, they can be involved in two stations with two classes competing in the first (Section 3.3) or in the second station (Section 3.4), or finally thecan be involved in three stations in tandem. We have not implemented this three wayinteraction heuristic given that the pairwise interaction heuristic behaves very well.

In the next two subsections we find explicitly the w corresponding to the two building blocks of the pairwise interaction heuristic (two classes competing in a single station and two classes in tandem), describe an example of how the heuristic is used and describe the implementation of the heuristic.

5.1 Competing Classes

When two classes compete for service at the same station,

Ā = [ \begin{matrix}−1/μ_{1} & −1/μ_{2} \\ 1/μ_{1} & 0 \\ 0 & 1/μ_{2} \end{matrix} ]

and

b̄ = [ \begin{matrix}1 − (\frac{λ_{1}}{μ_{1}} + \frac{λ_{2}}{μ_{2}} ) \\ λ_{1}/μ_{1} \\ λ_{2}/μ_{2} \end{matrix} ]

The polyhedron Āẋ ≤ b̄ is shown in Figure 8. The three extreme points labeled A,

[figure: ẋ_{2}; ẋ_{1}; ẋ = (λ_{1} − μ_{1}, λ_{2}); u = (1, 0); u = (1, 1); ẋ = (λ_{1}, λ_{2} − μ_{2}); A; B; C; ◇; ◇; ◇]

Figure 8: The feasible region of competing classes

B and C correspond to u = (1, 0), (0, 1) and (0,0) respectively, but only the first two are of interest. From (7) u = (1, 0) corresponds to ẋ = (λ_{1} − μ_{1}, λ_{2}) and u = (1, 1) corresponds to ẋ = (λ_{1}, λ_{2} − μ_{2}).

For the heuristic to function correctly, w needs to be chosen so that point A is the optimal solution of (8) when c_{1}μ_{1} > c_{2}μ_{2} and point B is the optimal if c_{1}μ_{1} < c_{2}μ_{2} (if c_{1}μ_{1} = c_{2}μ_{2}, the problem is indifferent).

Theorem 6 If w = c, then (8) correctly determines the optimal u(0) for the problem with two competing classes.

Proof The optimal solution for (8) will be attained at the point A, B or C. Then the objective function value z_{i} corresponding to point i (i = A, B, C) for the linear program (8) is:

z_{A} = w_{1}λ_{1} − w_{1}μ_{1} + w_{2}λ_{2}
z_{B} = w_{1}λ_{1} + w_{2}λ_{2} − w_{2}μ_{2}
z_{C} = w_{1}λ_{1} + w_{2}λ_{2}

By choosing w = c, it is clear that z_{C} will never be the minimum and z_{A} < z_{B} if and only if c_{1}μ_{1} > c_{2}μ_{2}, which is exactly the optimal policy. . □

[figure: ẋ_{2}; ẋ_{1}; ẋ = (λ, −μ_{2}); u = (0, 1); ẋ = (λ − μ_{1}, μ_{1} − μ_{2}); u = (1, 1); A; D; B; C; ◇; ◇; ◇; ◇]

Figure 9: Two classes in tandem

5.2 Two Classes in Tandem

When two classes are in tandem,

Ā = [ \begin{matrix}−1/μ_{1} & 0 \\ −1/μ_{2} & −1/μ_{2} \\ 1/μ_{1} & 0 \\ 1/μ_{2} & 1/μ_{2} \end{matrix} ]

and

b̄ = [ \begin{matrix}1 − λ/μ_{1} \\ 1 − λ/μ_{2} \\ λ/μ_{1} \\ λ/μ_{2} \end{matrix} ]

The polyhedron Āẋ ≤ b̄ is shown in Figure 9. Extreme points C and D correspond to u = (0, 0) and u = (1, 0) respectively, which will never be optimal solutions to the scheduling problem. We will choose a w that correctly chooses between points A and B under the same conditions as theorem 3. The next theorem describes how to choose such a w.

Theorem 7 If w is defined as

w = \{ \begin{matrix}c & if μ_{1} > μ_{2} \\ (c_{1}, min (\frac{x_{2}}{x_{1}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}} c_{2}, c_{2})) & otherwise. \end{matrix}

then (8) correctly determines the optimal u(0) for the problem of two classes in tendem.

Proof Assume μ_{1} > μ_{2}, then theorem 3 states that point A should be optimal if c_{1} > c_{2} and point B otherwise. Let z_{i} be the objective value corresponding to point i (i = A, B, C, D).

z_{A} = c_{1}λ − c_{1}μ_{1} + c_{2}μ_{1} − c_{2}μ_{2}
z_{B} = c_{1}λ − c_{2}μ_{2}
z_{C} = c_{1}λ
z_{D} = c_{1}λ − c_{1}μ_{1} + c_{2}μ_{1}

Point C will clearly never be optimal. Since, c and μ are positive, z_{A} < z_{D} so point D will not be optimal either. Point A will be chosen over point B when z_{A} < z_{B} or c_{1} > c_{2}.

Now assume μ_{1} < μ_{2}. Theorem 3 states that point A should be optimal if either c_{1} > c_{2} or both c_{1} < c_{2} and x_{1}/x_{2} > c_{2}(μ_{1} − λ)/c_{1}(μ_{2} − μ_{1}) or, to be concise, point A should be optimal when

c_{1} > min (\frac{x_{2}}{x_{1}} \frac{μ_{1} − λ}{μ_{2} − μ_{1}} c_{2}, c_{2})

If \frac{x_{2}}{x_{1}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}} < 1, then

z_{A} = c_{1}(λ − μ_{1}) − \frac{x_{2}}{x_{1}} (μ_{1} − λ)c_{2}
z_{B} = c_{1}λ − \frac{x_{2}}{x_{1}} \frac{μ_{1} − λ}{μ_{2} − μ_{1}} c_{2}μ_{2}
z_{C} = c_{1}λ
z_{D} = c_{1}(λ − μ_{1}) + \frac{x_{2}}{x_{1}} \frac{μ_{1} − λ}{μ_{2} − μ_{1}} c_{2}μ_{1}

Once again, points C and D will never be optimal. Point will be optimal when

z_{A} < z_{B}
−c_{1}μ_{1} − \frac{x_{2}}{x_{1}} (μ_{1} − λ)c_{2} < −\frac{x_{2}}{x_{1}} \frac{μ_{1} − λ}{μ_{2} − μ_{1}} μ_{2}
\frac{x_{2}}{x_{1}} c_{2}(μ_{1} − λ) [\frac{μ_{2}}{μ_{2} − μ_{1}} − 1] < c_{1}μ_{1}
\frac{x_{2}}{x_{1}} c_{2} \frac{μ_{1} − λ}{μ_{2} − μ_{1}} < c_{1}
\frac{μ_{1} − λ}{μ_{2} − μ_{1}} \frac{c_{2}}{c_{1}} < \frac{x_{1}}{x_{2}}

Alternatively, point B will be optimal if

c_{1} < \frac{x_{2}}{x_{1}} \frac{μ_{1} − λ}{μ_{2} − μ_{1}} c_{2} < c_{2}.

Finally, assume \frac{x_{2}}{x_{1}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}} > 1 and μ_{1} < μ_{2}. In this case w is the same as the μ_{1} > μ_{2} case and point A will be chosen if c_{1} > c_{2}. Point B will be optimal if c_{1} < c_{2} < \frac{x_{2}}{x_{1}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}} c_{2}. □

5.3 Example

Having specified the building blocks of the pairwise interaction heuristic by analyzing an example studied in section 3.3. We will see that the heuristic correctly determines the optimal policy.

The problem is defined by

min  ∫_{0}^{T} c′ x(t)dt
ẋ_{1}(t) = λ_{1} − u_{1}(t)μ_{1}
ẋ_{2}(t) = u_{1}(t)μ_{1} − u_{2}(t)μ_{2}
ẋ_{3}(t) = λ_{3} − u_{3}(t)μ_{1}
x(0) = x̃
u_{1}(t) + u_{3}(t) ≤ 1
u_{2}(t) ≤ 1
x(t) ≥ 0
u(t) ≥ 0

Assume that μ_{1} < μ_{2}.

Ā = [ \begin{matrix}−1/μ_{1} & 0 & −1/μ_{3} \\ −1/μ_{2} & −1/μ_{2} & 0 \\ 1/μ_{1} & 0 & 0 \\ 1/μ_{2} & 1/μ_{2} & 0 \\ 0 & 0 & 1/μ_{3} \end{matrix} ]

and

b̄ = [ \begin{matrix}1 − (λ_{1}/μ_{1} + λ_{3}/μ_{3}) \\ 1 − λ/μ_{2} \\ λ_{1}/μ_{1} \\ λ_{1}/μ_{2} \\ λ_{3}/μ_{3} \end{matrix} ]

Classes one and three compete at station one so the first row of the matrix M is

\frac{w_{1}}{w_{3}} = \frac{c_{1}}{c_{3}} .

Classes one and two are in tandem so the second row is

\frac{w_{1}}{w_{3}} = \frac{c_{1}}{min (\frac{x_{2}}{x_{1}} \frac{μ_{1}−λ}{μ_{2}−μ_{1}} c_{2}, c_{2})}.

Classes two and three are not related, so that pair does not contribute to the system of equations.

If x_{2}(μ_{1} − λ)/x_{1}(μ_{2} − μ_{1}) < 1, then w would be

w = (c_{1}x_{1}(μ_{2} − μ_{1}), c_{2}x_{2}(μ_{1} − λ), c_{3}x_{1}(μ_{2} − μ_{1})).

Note that w can be scaled so w′ e = 1 does not need to be enforced – it is only used to insure a nontrivial solution for w. As mentioned in section 3.3, u_{2}(0) = 1 and u_{3}(0) = 1 − u_{1}(0) so the only extreme points to be considered correspond to u = (1, 1, 0) and u = (0, 1, 1). The extreme points are ẋ = (λ_{1} − μ_{1}, μ_{1} − μ_{2}, λ_{3}) and ẋ = (λ_{1}, −μ_{2}, λ_{3} − μ_{3}). Direct calculation reveals that the control (1,1,0) will be chosen if

c_{2}x_{2}(μ_{1} − λ_{1})μ_{1} < (c_{1}μ_{1} − c_{3}μ_{3})(μ_{2} − μ_{1})μ_{1}x_{1}

or

\frac{μ_{1} − λ}{μ_{2} − μ_{1}} \frac{c_{2}μ_{1}}{c_{1}μ_{1} − c_{3}μ_{3}} < \frac{x_{1}}{x_{2}}

if c_{1}μ_{1} > c_{3}μ_{3}. If c_{1}μ_{1} < c_{3}μ_{3}, the above condition will never be satisfied. This, along with the assumption at the start of this paragraph, agree with cases 2 and 4 of theorem 4.

If x_{2}(μ_{1} − λ)/x_{1}(μ_{2} − μ_{1}) > 1, then w = c. By computing the objective function values, it is easy to show that control (0,1,1) will be used when c_{1}μ_{1} − c_{3}μ_{3} < c_{2}μ_{1} (case 5 of theorem 4). This condition will always be satisfied when c_{3}μ_{3} > c_{1}μ_{1} so the control will always be (0,1,1) in this case (case 1 of theorem 4).

When μ_{1} > μ_{2}, then w = c once again and the result agrees with (case 2 of theorem 4).

5.4 Implementation of the Heuristic

We have implemented the pairwise interaction heuristic using the symbolic manipulation program Maple. The user inputs the topology of a network along with parameters λ and μ which the user can specify either numerically or symbolically if he is interested in the qualitative behavior of the problem. Depending on how the data is specified (numerically or symbolically) the implementation works in two different modes: numerical and symbolic. We next describe the implementation of each mode.

Under the symbolic mode, the program computes the objective function, w, by deriving a matrix M where w is the least squares solution to Mw = 0, w′ e = 1. For a particular network, there may be multiple instances of M, each one being used if the parameters satisfy certain conditions. In this stage, all of the candidate extreme points are also determined. Note that during this stage, both M and the extreme points are determined symbolically; this is based only on the structure of the network, not the network’s parameters. This information is stored in an intermediate file and then a Maple program is generated that will compute, for each set of conditions, the algebraic solution for w and the algebraic objective value for each extreme point. This Maple program creates a LaTeX file that prints, for each set of conditions, the conditions in effect, the objective value for each extreme point and the control value associated with that extreme point. This information can then be used to derive the complete policy; after determining which conditions are being used, the correct policy is the one with the lowest objective value. As an example, if the problem only had two classes, class 0 and class 1, competing at a single station, the output is included in Figure 10. This problem does not have any conditions Conditions

• none

Objective Values (use minimum)

• Control (1 0)

−c_{0}μ_{0}

• Control (0 1)

−c_{1}μ_{1}

Figure 10: Heuristic Output for Competing Queues

(there is only one version of the M matrix). Given the two objective functions, the policy is to use the minimum. After simplifying the expressions, it is easy to see that the output is correct. A current shortcoming of the software is that the expressions for the objective function often require “manual” simplification. It is important to note, however, that the policy generation is completely automated and that the user is left with a Maple program that can be used as a base for further simplification that is more easily undertaken by a human.

Under the numerical mode, the intermediate file is read along with instance data that specifies the parameters for the network. This information is used to determine which conditions are in effect and then the exact M matrix is calculated. The value of w is then found numerically and the linear program is solved to produce the policy to be used for this instance of the network. Finally the heuristic policy is simulated to produce an estimate of the cost. In this case the software package outputs a policy and its cost automatically without user intervention.

We have already seen that the pairwise interaction heuristic gives the correct optimal policy in one example. Naturally we would to compare the performance of the heuristic with the optimal policy (found using the discretization technique of Section 4). We perform this comparison in the next section.

6 Performance of the Heuristic

To test the performance of the heuristic many problems were processed both symbolically and numerically and the numerical results were compared with the optimal solution from the discrete approximation. The results of some sample networks are presented in this section.

6.1 A 4 class-2 station network

We first consider the network in Figure 11 that has four classes and two stations operating as two tandem queues. Class zero jobs arrive at the first station with rate λ_{0} and then become class one jobs at the second station before leaving the system. Similarly, class two jobs arrive at the second station with rate λ_{2} and become class three jobs at the first station before leaving the system. The parameters for λ, μ and c are shown in the figure. The values of x_{0} and x_{3} were fixed at 4 and 2 respectively while x_{1} and x_{2} varied from 1 to 10. In order to obtain insight to the difference of the optimal policy and heuristic policy we have calculated both policies for all combinations of initial conditions for x_{1} and x_{2}. We observe that the heuristic policy correctly identifies the existence of a threshold curve. The calulation of the slope of the curve is not exact but is very close to the optimal one.

The same network was also processed when x_{3} was fixed at 7. In this case the optimal and heuristic policies were identical. As another example, the network was examined with the parameters shown in Figure 13. Once again the two polic were quite close.

6.2 A 4 class-3 station network

The next network tested is shown in Figure 14. Class zero jobs arrive at station one with rate λ_{0} and then become class two jobs at the second station before leaving. Class one jobs arrive at the first station with rate λ_{1} and are then processed at the third station before leaving. The values of x_{2} and x_{3} were fixed at 2 and 6 respectively while x_{0} and x_{1} took on the values 2, 4, …, 10. The control for classes two and three will always be 1 so the only decision is which class to process at the first station. The control for time zero is shown in the graphs. As before the heuristic identifies correctly the threshold phenomenon and misses the exact slope slightly.

6.3 Reentrant lines

We next consider the reentrant network shown in Figure 15. The values of x_{0}, x_{1}, x_{3}andx_{4} were fixed at 2,4,3 and 5 respectively. The values of x_{2} and x_{5} were in the set 2, 4, …, 10. In this example, the control chosen by the heuristic and the optimal control was always U(0) = (0,0,0,1,1,1).

As a final example, consider the network in Figure 16 that was studied by Perkins and Kumar [14]. Because of the difficulty to solve the discrete LP to optimality only two instances of this network were processed. In the first instance, the parameters were

λ = 1, μ = (2, 3, 2, 3, 2, 5, 4, 5, 4)
c = (2, 3, 3, 2, 4, 3, 4, 5, 3), x = (5, 6, 4, 7, 6, 5, 8, 5, 2)

and the control U(0) was

U(0) = (0, 0, 0, 0, 0, 0, 0, 1, 1)

for both the LP and the heuristic. In the second instance, the parameters were

λ = 1, μ = (2, 3, 2, 3, 2, 5, 4, 5, 4)
c = (2, 3, 6, 2, 10, 3, 4, 5, 3), x = (5, 6, 6, 7, 6, 5, 8, 5, 2)

[figure: λ_{0}; λ_{2}; μ_{0}; μ_{1}; μ_{2}; μ_{3}]

λ_{0} = 1, λ_{2} = 2, μ = (4, 8, 3, 4), c = (4, 3, 2, 3)

2 4 6 8 10  2  4  6  8 10 x2  x1 Discrete LP  ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ • • • • • • • • • • • • • • • • • • • • • • • • •  2 4 6 8 10  2  4  6  8 10 x2  x1 Heuristic  ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ • • • • • • • • • • • • • • • • • • • • • • • • • • • ∘ ∘ ∘ ∘

∘ = (1,0,1,0) • = (1,1,0,0) ⋄ = (0,1,0,1) Figure 11: Example with two competing tandem queues

Once again U(0) was the same for both the heuristic and the LP; the value was U(0) = (0, 0, 0, 0, 1, 0, 0, 0, 1).

2 4 6 8 10  2  4  6  8 10 x2  x1 Discrete LP  ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ • • • • • • • • • • • • • • • • • • • •  2 4 6 8 10  2  4  6  8 10 x2  x1 Heuristic  ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ • • • • • • • • • • • • • • • • • • • •
• = (1,1,0,0) ⋄ = (0,1,0,1)

Figure 12: The same network with x_{3} = 7

λ_{0} = 1, λ_{2} = 2, μ = (2, 3, 4, 8), c = (6, 2, 1, 1) Figure 13: Results for the third instance

[figure: λ_{1}; λ_{0}; μ_{1}; μ_{0}; μ_{3}; μ_{2}]

λ_{0} = 1, λ_{1} = 1, μ = (4, 2, 6, 3), c = (3, 2, 4, 3)

2 4 6 8 10  2  4  6  8 10 x1  x0 Discrete LP ⋄ ⋄ ⋄ ⋄ ⋄  • • • • • • • • • • • • • • • • • • • •  2 4 6 8 10  2  4  6  8 10 x1  x0 Heuristic ⋄ ⋄  • • • • • • • • • • • • • • • ∘ ∘ ∘ ∘ ∘ ∘ ∘ ∘

∘ = (0,0,1,1) • = (1,0,1,1) ⋄ = (0,1,1,1) Figure 14: Results for the fork network

[figure: λ; μ_{0}; μ_{5}; μ_{1}; μ_{3}; μ_{2}; μ_{4}]

λ = 1, μ = (5, 3, 4, 3, 4, 3), c = (4, 5, 6, 7, 6, 4) Figure 15: Results for a reentrant Network

[figure: λ; μ_{6}; μ_{3}; μ_{0}; μ_{7}; μ_{4}; μ_{1}; μ_{5}; μ_{8}; μ_{2}]

Figure 16: Network from Kumar and Perkins

7 Concluding Remarks and Open Problems

We presented exact and heuristic approaches for a fluid model of a multiclass queueing network that offers, we believe, qualitative insight to the optimal policy. While the exact discretization linear programming algorithm solves the problem numerically, it has dimensionallity problems so that we can solve problems with up to ten stations. Obviously deriving an optimal algorithm for the problem is a challenging open problem.

We proposed a hierarchy of heuristics that capture interactions between various classes and is able to obtain solutions for networks with hundreds of classes of the type encountered in applications very efficiently. We have seen that the pairwise interaction heuristic produces near optimal policies. It would be desirable to determine what would be the benefits of a higher interaction heuristic.

Finally, perhaps the most challenging problem is to determine whether the addition of stochastic arrivals and services in the problem definition produces new insights in the problem. Under a fluid model the threshold curves are hyperplanes. We conjecture that the corresponding threshold curves are higher order curves, i.e., we suspect that qualitatively the optimal policies are unchanged.

References

[1] E.J. Anderson, Basic Solutions and a ‘Simplex’ Method for a Class of Continuous Linear Programs, working paper.

[2] Edward J. Anderson and Peter Nash, Linear Programming in Infinite-Dimensional Spaces, John Wiley & Sons, New York, 1987.

[3] Derek Atkins and Hong Chen, Dynamic Scheduling Control for a Network of Queues, The Proceedings of the 12th World Congress International Federation of Automatic Control (Sidney).

[4] Dimitris Bertsimas, Ioannis Ch. Paschalidis and John N. Tsitsiklis, Optimization of Multiclass Queueing Networks: Polyhedral and Nonlinear Characterizations of Achievable Performance, Annals of Applied Probability, 4, (1994), no. 1, 43-75.

[5] Arthur E. Bryson, Jr. and Yu-Chi Ho, Applied Optimal Control, Blaisdell, Waltham, MA, 1969.

[6] Hong Chen and David D. Yao, Optimal Scheduling of a Multi-Class Fluid Network, 1989.

[7] Filipiak, Modelling and Control of Dynamic Flows in Communication Networks, Springer Verlag, Berlin, 1988.

[8] J.M. Harrison, Brownian Models of Queueing Networks with Heterogeneous Customers, Presented at the IMA Workshop on Stochastic Differential Systems, 1986.

[9] J.M. Harrison and L.M. Wein, Scheduling Networks of Queues; Heavy Traffic Analysis of a Two-Station Closed Network, Operations Research, 38, (1990), no. 2, 1052-1064.

[10] R. Kalaba and R. Sridhar, Initial-Value Methods in Optimal Control Theory, The RAND Corporation, Santa Monica, CA, 1968.

[11] Morton I. Kamien and Nancy L. Schwartz, Dynamic Optimization: The Calculus of Variations and Optimal Control in Economics and Management, North-Holland, Amsterdam, 1991.

[12] I.H. Mufti, Computational Methods in Optimal Control Problems, Lecture Notes in Operations Research and Mathematical Systems, Springer Verlag, Berlin, 1970.

[13] Katta G. Murty, Linear Programming, John Wiley & Sons, New York, 1983.

[14] J.R. Perkins and P.R. Kumar, Optimal Control of Pull Manufacturing Systems, University of Illinois at Urbana, working paper, 1994.

[15] L.S. Pontryagin, V.G. Boltyanskii, R.V. Gamkrelidze, and E.F. Mishchenko, The Mathematical Theory of Optimal Processes, Interscience Publishers, New York, 1962.

[16] Malcolm C. Pullan, An Algorithm for a Class of Continuous Linear Programs, SIAM J. Control and Optimization, 31, (1993), no. 6, 1558-1577.

[17] Atle Seierstad and Knut Sydsæter, Sufficient Conditions in Optimal Control Theory, International Economic Review, 18, (1977), no. 2, 367-391.

[18] L.M. Wein, Optimal Control of a Two-Station Brownian Network, Mathematics of Operations Research, 15, (1990), 215-242.

[19] L.M. Wein, Scheduling Networks of Queues; Heavy Traffic Analysis of a Two Station Network with Controllable Inputs, 1990.