A Constraint Programming Approach for Solving a Queueing Control Problem
aa r X i v : . [ c s . A I] O c t Journal of Artificial Intelligence Research 32 (2008) 123–167 Submitted 09/07; published 05/08
A Constraint Programming Approach for Solving aQueueing Control Problem
Daria Terekhov [email protected]
J. Christopher Beck [email protected]
Department of Mechanical & Industrial EngineeringUniversity of Toronto, Canada
Abstract
In a facility with front room and back room operations, it is useful to switch workersbetween the rooms in order to cope with changing customer demand. Assuming stochasticcustomer arrival and service times, we seek a policy for switching workers such that theexpected customer waiting time is minimized while the expected back room staffing issufficient to perform all work. Three novel constraint programming models and severalshaving procedures for these models are presented. Experimental results show that a modelbased on closed-form expressions together with a combination of shaving procedures is themost efficient. This model is able to find and prove optimal solutions for many probleminstances within a reasonable run-time. Previously, the only available approach was aheuristic algorithm. Furthermore, a hybrid method combining the heuristic and the bestconstraint programming method is shown to perform as well as the heuristic in termsof solution quality over time, while achieving the same performance in terms of provingoptimality as the pure constraint programming model. This is the first work of which weare aware that solves such queueing-based problems with constraint programming.
1. Introduction
The original motivation for the study of scheduling and resource allocation problems withinartificial intelligence (AI) and constraint programming (CP) was that, in contrast to Op-erations Research (OR), the full richness of the problem domain could be represented andreasoned about using techniques from knowledge representation (Fox, 1983). While much ofthe success of constraint-based scheduling has been due to algorithmic advances (Baptiste,Le Pape, & Nuijten, 2001), recently, there has been interest in more complex problemssuch as those involving uncertainty (Policella, Smith, Cesta, & Oddi, 2004; Sutton, Howe,& Whitley, 2007; Beck & Wilson, 2007). In the broader constraint programming commu-nity there has been significant work over the past five years on reasoning under uncertainty(Brown & Miguel, 2006). Nonetheless, it is recognized that “constraint solving under changeand uncertainty is in its infancy” (Brown & Miguel, 2006, p. 754).Queueing theory has intensively studied the design and control of systems for resourceallocation under uncertainty (Gross & Harris, 1998). Although much of the study has beendescriptive in the sense of developing mathematical models of queues, there is prescriptivework that attempts to develop queue designs and control policies so as to optimize quantitiesof interest (e.g., a customer’s expected waiting time) (Tadj & Choudhury, 2005). One ofthe challenges to queueing theory, however, is that analytical models do not yet extend to c (cid:13) erekhov & Beck richer characteristics encountered in real-world problems such as rostering for call centres(Cezik & L’Ecuyer, 2008).Our long-term goal is to integrate constraint programming and queueing theory with twoends in mind: the extension of constraint programming to reason better about uncertaintyand the expansion of the richness of the problems for which queueing theory can be broughtto bear. We do not achieve this goal here. Rather, this paper represents a first step wherewe solve a queueing control problem with constraint programming techniques. Specifically,we develop a constraint programming approach for a queueing control problem which arisesin retail facilities, such as stores or banks, which have back room and front room operations.In the front room, workers have to serve arriving customers, and customers form a queueand wait to be served when all workers are busy. In the back room, work does not directlydepend on customer arrivals and may include such tasks as sorting or processing paperwork.All workers in the facility are cross-trained and are assumed to be able to perform backroom tasks equally well and serve customers with the same service rate. Therefore, it makessense for the managers of the facility to switch workers between the front room and theback room depending both on the number of customers in the front room and the amountof work that has to be performed in the back room. These managers are thus interested infinding a switching policy that minimizes the expected customer waiting time in the frontroom, subject to the constraint that the expected number of workers in the back room issufficient to complete all required work. This queueing control problem has been studied indetail by Berman, Wang and Sapna (2005), who propose a heuristic for solving it.Our contributions are twofold. Firstly, constraint programming is, for the first time,used to solve a stochastic queueing control problem. Secondly, a complete approach for aproblem for which only a heuristic algorithm existed previously is presented.The paper is organized as follows. Section 2 presents a description of the problemand the work done by Berman et al. (2005). In the next section, three CP models for thisproblem are proposed. Sections 4 and 5 present methods for improving the efficiency of thesemodels, focusing on dominance rules and shaving procedures, respectively. Section 6 showsexperimental results comparing the proposed CP models and combinations of inferencemethods. The performance of the CP techniques is contrasted with that of the heuristicmethod of Berman et al. Based on these results, a hybrid method is proposed and evaluatedin Section 7. In Section 8, a discussion of the results is presented. Section 9 describes relatedproblems and states some directions for future work. Section 10 concludes the paper. Anappendix containing the derivations of some expressions used in the paper is also included.
2. Problem Description
Let N denote the number of workers in the facility, and let S be the maximum number ofcustomers allowed in the front room at any one time. When there are S customers present,arriving customers will be not be allowed to join the front room queue and will leave withoutservice. Customers arrive according to a Poisson process with rate λ . Service times in thefront room follow an exponential distribution with rate µ . The minimum expected numberof workers that is required to be present in the back room in order to complete all of thenecessary work is assumed to be known, and is denoted by B l , where l stands for ‘lower
1. The notation used by Berman et al. (2005) is adopted throughout this paper.
CP Approach for a Queueing Control Problem bound’. Only one worker is allowed to be switched at a time, and both switching time andswitching cost are assumed to be negligible. The goal of the problem is to find an optimalapproach to switching workers between the front room and the back room so as to minimizethe expected customer waiting time, denoted W q , while at the same time ensuring that theexpected number of workers in the back room is at least B l . Thus, a policy needs to beconstructed that specifies how many workers should be in the front room and back room ata particular time and when switches should occur. The term policy is used in the queueing control literature to describe a rule which prescribes,given a particular queue state, the actions that should be taken in order to control the queue.Most of the research on the optimal control of queues has focused on determining when aparticular type of policy is optimal, rather than on finding the actual optimal values ofthe parameters of the policy (Gross & Harris, 1998). The term optimal policy is used inthe literature to mean both the optimal type of policy and the optimal parameter values for a given policy type. The distinction between the two is important since showing thata particular policy type is optimal is a theoretical question, whereas finding the optimalvalues for a specific policy type is a computational one.The policy type adopted here is the one proposed by Berman et al. (2005). A policy isdefined in terms of quantities k i , for i = 0 , . . . , N and states that there should be i workersin the front room whenever there are between k i − + 1 and k i customers (inclusive) inthe front room, for i = 1 , , . . . , N . As a consequence of this interpretation, the followingconstraints have to hold: k i − k i − ≥ k ≥ k N = S . For example, consider a facilitywith S = 6 and N = 3, and suppose the policy ( k , k , k , k ) = (0 , , ,
6) is employed.This policy states that when there are k + 1 = 1 or k = 2 customers in the front room,there is one worker in the front room; when there are 3 customers, there are 2 workers; andwhen there are 4, 5, or 6 customers, all 3 workers are employed in the front. Alternatively, k i can be interpreted as an upper bound on the number of customers that will be servedby i workers under the given policy. Yet another interpretation of this type of switchingpolicy comes from noticing that as soon as the number of customers in the front room isincreased by 1 from some particular switching point k i , the number of workers in the frontroom changes to i + 1. This definition of a policy forms the basis of the model proposed byBerman et al., with the switching points k i , i = 0 , . . . , N −
1, being the decision variablesof the problem, and k N being fixed to S , the capacity of the system.This policy formulation does not allow a worker to be permanently assigned to the backroom: the definition of the k i ’s is such that every worker will, in some system state, beserving customers. Due to this definition, there exist problem instances which are infea-sible with this policy type yet are feasible in reality. Consequently, the proposed policyformulation is sub-optimal (Terekhov, 2007). However, the goal of the current work is todemonstrate the applicability of constraint programming to computing the optimal valuesfor the given policy type and not to address theoretical optimality questions. Therefore,the term optimal policy is used throughout the paper to refer to the optimal numericalparameters for the policy type proposed by Berman et al. (2005). erekhov & Beck In order to determine the expected waiting time and the expected number of workers in theback room given a policy defined by particular values of k i , Berman et al. first define a setof probabilities, P ( j ), for j = k , k + 1 , . . . , S . Each P ( j ) denotes the steady-state (long-run) probability of the queue being in state j , that is, of there being exactly j customersin the facility. Based on the Markovian properties of this queueing system (exponentialinter-arrival and service times), Berman et al. define a set of detailed balance equations forthe determination of these probabilities: P ( j ) λ = P ( j + 1)1 µ j = k , k + 1 , . . . , k − P ( j ) λ = P ( j + 1)2 µ j = k , k + 1 , . . . , k − P ( j ) λ = P ( j + 1) iµ j = k i − , k i − + 1 , . . . , k i − P ( j ) λ = P ( j + 1) N µ j = k N − , k N − + 1 , . . . , k N − . The probabilities P ( j ) also have to satisfy the equation P Sj = k P ( j ) = 1. Intuitively, insteady-state, the average flow from state j to state j + 1 has to be equal to the averageflow from state j + 1 to state j (Gross & Harris, 1998). Since P ( j ) can be viewed as thelong-run proportion of time when the system is in state j , the mean flow from state j to j + 1 is P ( j ) λ and the mean flow from state j + 1 to j is P ( j + 1) iµ for some i between 1and N depending on the values of the switching points.From these equations, Berman et al. derive the following expressions for each P ( j ): P ( j ) = β j P ( k ) , (2)where β j = j = k (cid:16) λµ (cid:17) j − k (cid:0) i (cid:1) j − k i − X i if k i − + 1 ≤ j ≤ k i i = 1 , . . . , N , (3) X i = i − Y g =1 (cid:18) g (cid:19) k g − k g − (4)( X ≡ , i = 1 , . . . , N.P ( k ) can be calculated using the following equation, which is derived by summing bothsides of Equation (2) over all values of j : P ( k ) S X j =0 β j = 1 . (5) CP Approach for a Queueing Control Problem
All quantities of interest can be expressed in terms of the probabilities P ( j ). Expectednumber of workers in the front room is F = N X i =1 k i X j = k i − +1 iP ( j ) , (6)while the expected number of workers in the back room is B = N − F . (7)The expected number of customers in the front room is L = S X j = k jP ( j ) . (8)Expected waiting time in the queue can be expressed as W q = Lλ (1 − P ( k N )) − µ . (9)This expression is derived using Little’s Laws for a system of capacity k N = S .Given a family of switching policies K = { K ; K = { k , k , ..., k N − , S } , k i integers, k i − k i − ≥ k ≥ , k N − < S } , the problem can formally be stated as: minimize K ∈ K W q (10) s.t. B ≥ B l P Sj = k P ( j ) = 1 equations (1) , (6) , (7) , (8) , (9) . Berman et al. (2005) refer to this problem as problem P . It is important to note that B , F and L are expected values and can be real-valued. Consequently, the constraint B ≥ B l states that the expected number of workers in the back room resulting from the realizationof any policy should be greater than or equal to the minimum expected number of backroom workers needed to complete all work in the back room. At any particular time point,there may, in fact, be fewer than B l workers in the back room.As far as we are aware, the computational complexity of problem P has not beendetermined. Berman et al. (p. 354) state that solving problem P “exactly is extremelydifficult since the constraints set (the detailed balance equations) changes when the policychanges.” Berman et al. (2005) propose a heuristic method for the solution of this problem. Thismethod is based on two corollaries and a theorem, which are stated and proved by theauthors. These results are key to the understanding of the problem, and are, therefore,repeated below. erekhov & Beck
Theorem 2.1 (Berman et al.’s Theorem 1)
Consider two policies K and K ′ which areequal in all but one k i . In particular, suppose that the value of k ′ J equals k J − , for some J from the set { , ..., N − } such that k J − k J − ≥ , while k ′ i = k i for all i = J . Then (a) W q ( K ) ≥ W q ( K ′ ) , (b) F ( K ) ≤ F ( K ′ ) , (c) B ( K ) ≥ B ( K ′ ) . As a result of Theorem 2.1, it can be seen that two policies exist which have specialproperties. Firstly, consider the policyˆ K = { k = 0 , k = 1 , k = 2 , ..., k N − = N − , k N = S } . This policy results in the largest possible F , and the smallest possible B and W q . Becausethis policy yields the smallest possible expected waiting time, it is optimal if it is feasible.On the other hand, the smallest possible F and the largest possible W q and B are obtainedby applying the policyˆˆ K = { k = S − N, k = S − N + 1 , ..., k N − = S − , k N = S } . Therefore, if this policy is infeasible, the problem (10) is infeasible also.Berman et al. propose the notions of eligible type 1 and type 2 components. An eligibletype 1 component is a switching point k i satisfying the condition that k i − k i − > < i < N or k i > i = 0. A switching point k i is an eligible type 2 component if k i +1 − k i > ≤ i < N . More simply, an eligible type 1 component is a k i variablewhich, if decreased by 1, will still be greater than k i − , while an eligible type 2 componentis a k i variable which, if increased by 1, will remain smaller than k i +1 . Eligible type 1components and eligible type 2 components will further be referred to simply as type 1 andtype 2 components, respectively.Based on the definitions of policies ˆ K and ˆˆ K , the notions of type 1 and 2 components,and Theorem 2.1, Berman et al. propose a heuristic, which has the same name as theproblem it is used for, P :1. Start with K = ˆˆ K .2. If B ( K ) < B l , the problem is infeasible. Otherwise, let imb W q = W q ( K ) and imb K = K . Set J = N .3. Find the smallest j ∗ s.t. 0 ≤ j ∗ < J and k j ∗ is a type 1 component. If no such j ∗ exists, go to 5. Otherwise, set k j ∗ = k j ∗ −
1. If B ( K ) < B l , set J = j ∗ and go to 5.If B ( K ) ≥ B l , go to 4.4. If W q ( K ) < imb W q , let imb W q = W q ( K ) and imb K = K . Go to 3.5. Find the smallest j ∗ s.t. 0 ≤ j ∗ < J and k j ∗ is a type 2 component. If no such j ∗ exists, go to 6. Otherwise, set k j ∗ = k j ∗ + 1. If B ( K ) < B l , repeat 5. If B ( K ) ≥ B l ,go to 4.6. Stop and return imb K as the best solution. CP Approach for a Queueing Control Problem
Parameter Meaning S front room capacity N number of workers in the facility λ arrival rate µ service rate B l expected number of workers required in the back roomTable 1: Summary of problem parameters.Notation Meaning Definition F expected number of workers in the front room Equation (6) B expected number of workers in the back room Equation (7) L expected number of customers in the front room Equation (8) W q expected customer waiting time Equation (9)Table 2: Summary of quantities of interest.Limiting the choice of j ∗ to being between 0 and J , and resetting J every time aninfeasible policy is found, prevents the heuristic from entering an infinite cycle. The heuristicguarantees optimality only when the policy it returns is ˆ K or ˆˆ K .Empirical results regarding the performance of heuristic P are not presented in thepaper by Berman et al. (2005). In particular, it is not clear how close policies provided by P are to the optimal policies. Berman et al.’s model of problem P requires five input parameters (Table 1) and hasexpressions for calculating four main quantities of interest (Table 2), most of which arenon-linear.
3. Constraint Programming Models
Some work has been done on extending CP to stochastic problems (Tarim, Manandhar,& Walsh, 2006; Tarim & Miguel, 2005; Walsh, 2002). Our problem is different from theproblems addressed in these papers because all of the stochastic information can be explicitlyencoded as constraints and expected values, and there is no need for either stochasticvariables or scenarios.It is known that creating an effective constraint programming model usually requiresone to experiment with various problem representations (Smith, 2006). Consequently, herewe present, and experiment with, three alternative models for problem P (Equation (10)).Although the first model is based directly on the formulation of Berman et al., further erekhov & Beck models were motivated by standard CP modelling techniques, such as the use of dualvariables (Smith, 2006). We investigate the following three CP models: • The
If-Then model is a CP version of the formal definition of Berman et al. • The
PSums model uses a slightly different set of variables, and most of the constraintsare based on closed-form expressions derived from the constraints that are used in the
If-Then model. • The
Dual model includes a set of dual decision variables in addition to the variablesused in the
If-Then and
PSums models. Most of the constraints of this model areexpressed in terms of these dual variables.
There are a number of modelling components that are common to each of the constraintmodels. Before presenting the models, we therefore present their common aspects.
All three of the proposed models have a set of decision variables k i , i = 0 , , . . . , N , repre-senting the switching policy. Each k i from this set has the domain [ i, i + 1 , . . . , S − N + i ]and has to satisfy the constraint k i < k i +1 (since the number of workers in the front room, i ,increases only when the number of customers, k i , increases). Due to Berman et al.’s policydefinition, k N must equal S . All models include variables and constraints for the representation of the balance equations(Equation (1)), and expressions for F , the expected number of workers in the front room,and L , the expected number of customers in the front room. However, these representationsdiffer slightly depending on the model, as noted below in Sections 3.2, 3.3 and 3.4.A set of auxiliary variables, βSum ( k i ), defined as P k i j = k i − +1 β j , for all i from 1 to N −
1, is included in each of the models (see Equations (2)–(4) for the definition of β j ).These are necessary for representing Equation (11), which relates these variables to P ( k ), afloating point variable with domain [0 ..
1] representing the probability of having k customersin the facility. These auxiliary variables and constraint ensure that an assignment of alldecision variables leads to a unique solution of the balance equations. We discuss the formaldefinition of these auxiliary variables in Section 3.2.4. P ( k ) N X i =0 βSum ( k i ) = 1 (11)The back room constraint, B ≥ B l , is stated in all models as N − F ≥ B l . The equationfor W q is stated in all models as Equation (9). If-Then
Model
The initial model includes the variables P ( j ) for j = k , k + 1 , . . . , k , k + 1 , . . . , k N − , k N , each representing the steady-state probability of there being j customers in the front CP Approach for a Queueing Control Problem minimize W q subject tok i < k i +1 ∀ i ∈ { , , . . . , N − } ; k N = S ;( k i ≤ j ≤ k i +1 − → P ( j ) λ = P ( j + 1)( i + 1) µ, ∀ i ∈ { , , . . . , N − } , ∀ j ∈ { , , . . . , S − } ;( j < k ) ≤ ( P ( j ) = 0) , ∀ j ∈ { , , . . . , S − N − } ; S X j =0 P ( j ) = 1;( k = j ) → P ( j ) N X i =0 βSum ( k i ) = 1 , ∀ j ∈ { , , . . . , S } ; L = S X j =0 jP ( j );( k i − + 1 ≤ j ≤ k i ) → r ( i, j ) = iP ( j ) , ∀ i ∈ { , , . . . , N } , ∀ j ∈ { , , . . . , S } ;( k i − + 1 > j ∨ j > k i ) → r ( i, j ) = 0 , ∀ i ∈ { , , . . . , N } , ∀ j ∈ { , , . . . , S } ; F = N X i =1 S X j =0 r ( i, j ); W q = Lλ (1 − P ( k N )) − µ ; N − F ≥ B l ; auxiliary constraints. Figure 1: Complete
If-Then
Modelroom. These floating point variables with domain [0 ..
1] have to satisfy a system of balanceequations (Equation (1)) and are used to express L and F .The complete If-Then model is presented in Figure 1.
The balance equations are represented by a set of if-then constraints. For example, the firstbalance equation, P ( j ) λ = P ( j + 1) µ for j = k , k + 1 , ..., k −
1, is represented by theconstraint ( k ≤ j ≤ k − → P ( j ) λ = P ( j + 1) µ . Thus, somewhat inelegantly, an if-then erekhov & Beck constraint of this kind has to be added for each j between 0 and S − k i , k i +1 for i from 0 to N − N × S if-then constraints.The probabilities P ( j ) also have to satisfy the constraint P Sj = k P ( j ) = 1. The difficultywith this constraint is the fact that the sum starts at j = k , where k is a decision variable.In order to deal with this complication, we add the meta-constraint (( j < k ) ≤ ( P ( j ) = 0))for each j from the set { , , . . . , S − N − } . This implies that all values of P ( j ) with j less than k will be 0 and allows us to express the sum-of-probabilities constraint as P Sj =0 P ( j ) = 1. A set of if-then constraints also has to be included in order to represent Equation (6) asa constraint in our model. This is due to the dependence of this constraint on sums ofvariables between two switching points, which are decision variables. More specifically,we add a set of variables r ( i, j ) for representing the product of i and P ( j ) when j isbetween k i − + 1 and k i , and the constraints ( k i − + 1 ≤ j ≤ k i ) → r ( i, j ) = iP ( j ) and( k i − + 1 > j ∨ j > k i ) → r ( i, j ) = 0 for all i from 1 to N and for all j from 0 to S . Thetotal number of these if-then constraints is 2 N ( S + 1). F can then be simply stated as asum over the indices i and j of variables r ( i, j ). L is defined according to Equation (8). Since the meta-constraint (( j < k ) ≤ ( P ( j ) = 0))has been added to the model in order to ensure that P ( j ) = 0 for all j < k , the constraintfor L can be simply stated as the sum of the products of j and P ( j ) over all j from 0 to S : L = S X j =0 jP ( j ) . (12) The model includes a set of N − X i for all i from 3 to N ( X and X are always equal to 1). Instead of including S β j variables (refer to Equation (2) forthe definition of β j ), we use N + 1 continuous auxiliary variables βSum ( k i ) with domain[0 . . . S (cid:16) λµ (cid:17) S ], which represent sums of β j variables between j = k i − + 1 and j = k i (inclusive). βSum ( k ) is constrained to equal 1, while the rest of these variables are defined
2. We do not need to add this constraint for all j from 0 to S because the upper bound of the domain of k is S − N . CP Approach for a Queueing Control Problem according to Equation (13). The validity of this equation is proved in the appendix. βSum ( k i ) = k i X j = k i − +1 β j = X i (cid:18) λµ (cid:19) k i − − k +1 (cid:18) i (cid:19) − λiµ ! ki − ki − − λiµ ! if λiµ = 1 X i (cid:18) λµ (cid:19) k i − − k +1 (cid:18) i (cid:19) ( k i − k i − ) otherwise. (13) ∀ i ∈ { , . . . , N } ;The sum P k N j = k β j can then be expressed as P Ni =0 βSum ( k i ). The requirement that P ( k ) × P k N j = k β j has to equal 1 is stated as a set of if-then constraints ( k = j ) → P ( j ) × P Ni =0 βSum ( k i ) = 1 for all j ∈ { , , . . . , S } .In summary, the auxiliary constraints that are present in the model are: βSum ( k ) =1, Equation (13) and X i = Q i − g =1 (cid:16) g (cid:17) k g − k g − ∀ i ∈ { , . . . , N } .The If-Then model includes a total of 3
N S + 2 N + S + 1 if-then constraints, which areoften ineffective in search because propagation only occurs either when the left-hand side issatisfied or when the right-hand side becomes false. Consequently, the next model attemptsto avoid, as much as possible, the use of such constraints. PSums
Model
The second CP model is based on closed-form expressions derived from the balance equations(the details of the derivations are provided in the appendix). The set of P ( j ) variablesfrom the formulation of Berman et al. is replaced by a set of P Sums ( k i ) variables for i = 0 , . . . , N −
1, together with a set of probabilities P ( j ) for j = k , k , k , . . . , k N . Notethat P ( j ) is defined for each switching point only, not for all values from { , , . . . , S } . The P Sums ( k i ) variable represents the sum of all probabilities between k i and k i +1 − PSums model is presented in Figure 2. The remainder of this sectionprovides details of this model.
Balance equations are not explicitly stated in this model. However, expressions for P ( k i )and P Sums ( k i ) are derived in such a way that the balance equations are satisfied. The P Sums ( k i ) variables are defined in Equation (14). Equation (15) is a recursive formula forcomputing P ( k i +1 ). P Sums ( k i ) = k i +1 − X j = k i P ( j ) erekhov & Beck = P ( k i ) 1 − (cid:20) λ ( i + 1) µ (cid:21) k i +1 − k i − λ ( i + 1) µ if λ ( i +1) µ = 1 P ( k i )( k i +1 − k i ) otherwise. (14) P ( k i +1 ) = (cid:18) λ ( i + 1) µ (cid:19) k i +1 − k i P ( k i ) , ∀ i ∈ { , , . . . , N − } . (15)Additionally, the probability variables have to satisfy the constraint P N − i =0 P Sums ( k i ) + P ( k N ) = 1. F , the expected number of workers in the front room, can be expressed in terms of P ( k i )and P Sums ( k i ) as shown in Equation (16): F = N X i =1 i [ P Sums ( k i − ) − P ( k i − ) + P ( k i )] . (16) The equation for L is L = N − X i =0 L ( k i ) + k N P ( k N ) (17)where L ( k i ) = k i P Sums ( k i ) + P ( k i ) λ ( i + 1) µ × (cid:16) λ ( i +1) µ (cid:17) k i +1 − k i − ( k i − k i +1 ) + (cid:16) λ ( i +1) µ (cid:17) k i +1 − k i ( k i +1 − k i −
1) + 1 (cid:16) − λ ( i +1) µ (cid:17) . The auxiliary constraints are exactly the same as in the
If-Then model. However, therequirement that P ( k ) × P k N j = k β j has to equal 1 is stated as P ( k ) P Ni =0 βSum ( k i ) = 1,rather than as a set of if-then constraints, because this model has an explicit closed-formexpression for the variable P ( k ). Dual
Model
The problem can be alternatively formulated using variables w j , which represent the numberof workers in the front room when there are j customers present. The w j variables can be CP Approach for a Queueing Control Problem minimize W q subject tok i < k i +1 ∀ i ∈ { , , . . . , N − } ; k N = S ; P ( k i +1 ) = (cid:18) λ ( i + 1) µ (cid:19) k i +1 − k i P ( k i ) , ∀ i ∈ { , , . . . , N − } ; P Sums ( k i ) = P ( k i ) 1 − (cid:20) λ ( i + 1) µ (cid:21) k i +1 − k i − λ ( i + 1) µ if λ ( i +1) µ = 1 P ( k i )( k i +1 − k i ) otherwise , ∀ i ∈ { , , . . . , N − } ; N − X i =0 P Sums ( k i ) + P ( k N ) = 1; P ( k ) × N X i =0 βSum ( k i ) = 1; F = N X i =1 i [ P Sums ( k i − ) − P ( k i − ) + P ( k i )] ; L = N − X i =0 L ( k i ) + k N P ( k N ) ,L ( k i ) = k i P Sums ( k i ) + P ( k i ) λ ( i + 1) µ × (cid:16) λ ( i +1) µ (cid:17) k i +1 − k i − ( k i − k i +1 ) + (cid:16) λ ( i +1) µ (cid:17) k i +1 − k i ( k i +1 − k i −
1) + 1 (cid:16) − λ ( i +1) µ (cid:17) , ∀ i ∈ { , , . . . , N − } ; W q = Lλ (1 − P ( k N )) − µ ; N − F ≥ B l ; auxiliary constraints. Figure 2: Complete
PSums
Model erekhov & Beck referred to as the dual variables because, compared to the k i ’s, the roles of variables andvalues are reversed (Hnich, Smith, & Walsh, 2004; Smith, 2006). As stated by Smith, theuse of dual variables in a constraint programming model is beneficial if some constraints ofthe problem are easier to express using these new variables. This is the case for our problem– in fact, the use of dual variables allows us to significantly reduce the number of if-thenconstraints necessary for stating relations between the probabilities.In our Dual model, there are S + 1 w j variables, each with domain [0 , , . . . , N ]. Thesevariables have to satisfy the following equations: w = 0, w S = N and w j ≤ w j +1 for all j from 0 to S −
1. Additionally, the complete set of k i variables is included in this model,since some constraints are easier to express using the k i ’s rather than the w j ’s.The complete Dual model is presented in Figure 3, while the details are discussed below.
Given the dual variables, the balance equations can be restated as P ( j ) λ = P ( j + 1) w j +1 µ, ∀ j ∈ { , , . . . , S − } . (18)This formulation of the balance equations avoids the inefficient if-then constraints. The restof the restrictions on the probability variables are stated in terms of the k i variables, as in the If-Then model. In particular, the constraints P Sj =0 P ( j ) = 1 and (( k > j ) ≤ ( P ( j ) = 0)) ∀ j ∈ { , . . . , S − N − } are present in this model. In order to use redundant variables, a set of channelling constraints has to be added to themodel to ensure that an assignment of values to one set of variables will lead to a uniqueassignment of variables in the other set. The following channelling constraints are included: w j < w j +1 ↔ k w j = j ∀ j ∈ { , , . . . , S − } , (19) w j = w j +1 ↔ k w j = j ∀ j ∈ { , , . . . , S − } , (20) w j = i ↔ k i − + 1 ≤ j ≤ k i ∀ j ∈ { , , . . . , S } , ∀ i ∈ { , . . . , N } . (21)Constraints (19) and (20) are redundant given the constraint w j ≤ w j +1 . However, suchredundancy can often lead to increased propagation (Hnich et al., 2004). One direction forfuture work is examining the effect that removing one of these constraints may have on theperformance of the program. The expression for the expected number of workers in the front room is F = P Sj =0 w j P ( j ). The constraint used to express L is identical to the one used in the If-Then model: L = P Sj =0 jP ( j ). This equation is valid because all P ( j ) for j < k are constrained to be 0. CP Approach for a Queueing Control Problem minimize W q subject tow = 0 , w S = Nw j ≤ w j +1 ∀ j ∈ { , , . . . , S − } ; k i < k i +1 ∀ i ∈ { , , . . . , N − } ; k N = S ; w j < w j +1 ↔ k w j = j ∀ j ∈ { , , . . . , S − } ; w j = w j +1 ↔ k w j = j ∀ j ∈ { , , . . . , S − } ; w j = i ↔ k i − + 1 ≤ j ≤ k i ∀ i ∈ { , , . . . , N } , ∀ j ∈ { , , . . . , S } ; P ( j ) λ = P ( j + 1) w j +1 µ ∀ j ∈ { , , . . . , S − } ;( k = j ) → P ( j ) N X i =0 βSum ( k i ) = 1 , ∀ j ∈ { , , . . . , S } ;( j < k ) ≤ ( P ( j ) = 0) , ∀ j ∈ { , , . . . , S − N − } ; S X j =0 P ( j ) = 1; F = S X j =0 w j P ( j ); L = S X j =0 jP ( j ); W q = Lλ (1 − P ( k N )) − µ ; N − F ≥ B l ; auxiliary constraints. Figure 3: Complete
Dual
Model erekhov & Beck
Statistic/Model
If-Then PSums Dual N + 1 N + 1 N + S + 2 S + 1 2 N + 1 S + 1 N ( S −
1) + 2 S + 2 2 N + 1 3 S − N + 2 F N ( S + 1) + 1 1 1 L N + 1 1 N S + 2 N + S + 1 0 S + 1total N S + 4 N + 2 S + 6 6 N + 5 N S + 3 N + 6 S + 8Table 3: Summary of the main characteristics of the three proposed CP models. The auxiliary constraints present in this model are exactly the same as in the
If-Then model. The requirement that P ( k ) × P k N j = k β j = 1 is also stated as a set of if-thenconstraints ( k = j ) → P ( j ) P Ni =0 βSum ( k i ) = 1 for all j ∈ { , , . . . , S } . Table 3 presents a summary of the number of variables and constraints in each of the threeproposed models. It can be seen that the
PSums model has a smaller number of probabilityvariables and constraints but a slightly larger number of constraints for representing L thanthe other two models, and no if-then constraints. The Dual model has a larger number ofdecision variables than the
If-Then and
PSums models. This does not imply that the searchspace is bigger in this model because the two sets of variables are linked by channellingconstraints and so are assigned values via propagation. The
Dual allows the simplestrepresentations of F and L , each requiring only one constraint. The number of probabilityconstraints in the Dual is smaller than or equal to the number of such constraints in the
If-Then model and greater than in the
PSums model. However, the actual representationof these constraints is the most straightforward in the
Dual since it neither requires if-thenconstraints nor closed-form expressions.It is hard to determine, simply by looking at Table 3, which of the models will be moreefficient since, in CP, a larger number of constraints and/or variables may actually lead tomore propagation and to a more effective model (Smith, 2006). However, it is known thatif-then constraints do not propagate well, and, since the difference in the number of theseconstraints between the
PSums model and both the
If-Then model and the
Dual model isquite significant, one may expect the
PSums model to have some advantage over the othertwo models.In fact, preliminary experiments with all three models showed poor performance (seeTable 4 in Section 6). Due to the complexity of constraints relating the decision variables tovariables representing probabilities, there was little constraint propagation, and, essentially,search was required to explore the entire branch-and-bound tree. As a consequence, in thefollowing two sections we examine dominance rules (Beck & Prestwich, 2004; Smith, 2005)and shaving (Caseau & Laburthe, 1996; Martin & Shmoys, 1996), two stronger inference
CP Approach for a Queueing Control Problem forms used in CP. In Section 8, we investigate why the models without dominance rulesand shaving need to search the whole tree in order to prove optimality, and also discussdifferences in the performance of the models based on our experimental results.
4. Dominance Rules
A dominance rule is a constraint that forbids assignments of values to variables whichare known to be sub-optimal (Beck & Prestwich, 2004; Smith, 2005). For problem P , adominance rule states that, given a feasible solution, K , all further solutions have to haveat least one switching point assigned a lower value than the value assigned to it in K . Inother words, given two solutions K and K ′ , if the W q value resulting from policy K ′ issmaller than or equal to the W q value resulting from K , then there have to exist switchingpoints k ′ i and k i in K ′ and K , respectively, satisfying the condition k ′ i < k i . The followingtheorem states the dominance rule more formally. Theorem 4.1 (Dominance Rule)
Let K = ( k , k , . . . , k N ) and K ′ = ( k ′ , k ′ , . . . , k ′ N ) betwo policies such that k = k ′ = 0 , k = k ′ = 1 , . . . , k J − = k ′ J − = J − and k J = k ′ J (i.e. at least one of k J , k ′ J is strictly greater than J ) for some J ∈ { , , . . . , N − } . Let W q ( K ) and W q ( K ′ ) denote the expected waiting times resulting from the two policies K and K ′ , respectively. If W q ( K ′ ) < W q ( K ) , then there exists i ∈ { J, J + 1 , J + 2 , . . . , N − } forwhich k ′ i < k i . Proof: [By Contraposition] We prove the contrapositive of the statement in Theorem4.1: we assume that there does not exist i ∈ { J, J + 1 , . . . , N − } such that k ′ i < k i and show that, given this assumption, W q ( K ′ ) has to be greater than or equal to W q ( K ).Assume no i ∈ { J, J + 1 , . . . , N − } exists for which k ′ i < k i . Then one of thefollowing is true:(a) k n = k ′ n for all n ∈ { J, J + 1 , . . . , N − } , or(b) there exists at least one j ∈ { J, J + 1 , . . . , N − } such that k ′ j > k j , and thevalues of the rest of the switching points are the same in the two policies.Case (a) implies that K ′ and K are the same policy, and so W q ( K ) = W q ( K ′ ).To prove (b), suppose there exists exactly one j ∈ { J, J + 1 , J + 2 , . . . , N − } suchthat k ′ j > k j , and k n = k ′ n for all n ∈ { J, . . . , N − } \ { j } . Then K and K ′ aredifferent in the value of exactly one switching point. Consequently, by Theorem 2.1, W q ( K ′ ) ≥ W q ( K ). Similarly, by applying Theorem 2.1 several times, the resultgeneralizes to cases when there exists more than one j such that k ′ j > k j .Therefore, if no i ∈ { J, J + 1 , . . . , N − } exists for which k ′ i < k i , it follows that W q ( K ′ ) ≥ W q ( K ). In other words, if W q ( K ′ ) < W q ( K ), then there exists i ∈ { J, J + 1 , . . . , N − } for which k ′ i < k i . ✷ Theorem 4.1 implies that, given a feasible policy (0 , , . . . , k ∗ i , k ∗ i +1 , . . . , k ∗ N − , S ) whereall switching points with index i or greater are assigned values that are strictly greater thantheir lower bounds, we know that a solution with smaller W q has to satisfy the constraint(( k i < k ∗ i ) ∨ ( k i +1 < k ∗ i +1 ) ∨ . . . ∨ ( k N − < k ∗ N − )). Therefore, in order to implement erekhov & Beck the dominance rule, we add such a constraint during search every time a feasible policy isfound, which should lead to a reduction in the size of the search space. Section 6 presentsexperimental results regarding the usefulness of this technique.
5. Shaving
Shaving is an inference method that temporarily adds constraints to the model, performspropagation, and then soundly prunes variable domains based on the resulting state ofthe problem (Demassey, Artigues, & Michelon, 2005; van Dongen, 2006). For example, asimple shaving procedure may be based on the assignment of a value a to some variable x .If propagation following the assignment results in a domain wipe-out for some variable, theassignment is inconsistent, and the value a can be removed from the domain of x (Demasseyet al., 2005; van Dongen, 2006). In a more general case, both the temporary constraint andthe inferences made based on it can be more complex. Shaving has been particularly usefulin the job-shop scheduling domain, where it is used to reduce the domains of start and endtimes of operations (Caseau & Laburthe, 1996; Martin & Shmoys, 1996). For such problems,shaving is used either as a domain reduction technique before search, or is incorporated intobranch-and-bound search so that variable domains are shaved after each decision (Caseau& Laburthe, 1996).In a shaving procedure for our problem, we temporarily assign a particular value to aswitching point variable, while the rest of the variables are assigned either their maximumor their minimum possible values. Depending on whether the resulting policies are feasibleor infeasible, new bounds for the switching point variables may be derived.The instance N = 3, S = 6, λ = 15, µ = 3, B l = 0 .
32 is used below for illustrationpurposes. Policy ˆ K , which always yields the smallest possible W q , for this instance is( k , k , k , k ) = (0 , , ,
6) and policy ˆˆ K , which always yields the greatest possible W q , is(3 , , , .. .. ..
5] and [6]for k , k , k and k , respectively. At any step, shaving may be able to reduce the domainsof one or more of these variables. B l -based Shaving Procedure The initial shaving procedure consists of two cases in which either the upper or the lowerbounds of variables may be modified. In the first case, the constraint k i = min( k i ), wheremin( k i ) is the smallest value from the domain of k i , is temporarily added to the problemfor some particular value of i between 0 and N . All other switching points are assignedtheir maximum possible values using the function gMax . Given an array of variables, thefunction gMax assigns the maximum possible values to all of the variables that do not yethave a value, returning true if the resulting assignment is feasible, and false otherwise. Themaximum possible values are not necessarily the upper bound values in the domains of thecorresponding variables, rather they are the highest values in these domains that respectthe condition that k n < k n +1 , ∀ n ∈ { , ..., N − } . In the example, if k is assigned thevalue 1, while the rest of the variables are unbound, gMax would result in policy (0 , , , B value of 0 . CP Approach for a Queueing Control Problem
Recall that an assignment is infeasible when it yields a B value which is smaller than B l . When the policy resulting from the addition of k i = min( k i ) and the use of gMax isinfeasible, and min( k i ) + 1 ≤ max( k i ), the constraint k i > min( k i ) can be added to theproblem: if all variables except k i are set to their maximum values, and the problem isinfeasible, then in any feasible policy k i must be greater than min( k i ). Such reasoning isvalid since Theorem 2.1 states that increasing the value of a switching point will increase B . Note that if the solution is feasible, it should be recorded as the best-so-far solution ifits W q value is smaller than the W q value of the previous best policy. For easier reference,this part of the shaving procedure will be referred to as the gMax case.In the second case, the constraint k i = max( k i ) is added to the problem for some i between 0 and N , where max( k i ) is the maximum value from the domain of k i . The restof the variables are assigned the minimum values from their domains using the function gMin . These assignments are made in a way that respects the constraints k n < k n +1 , ∀ n ∈ { , ..., N − } . If the resulting policy is feasible, then the constraint k i < max( k i ) canbe permanently added to the problem, assuming max( k i ) − ≥ min( k i ). Since all variablesexcept k i are at their minimum values already, and k i is at its maximum, it must be true,again by Theorem 2.1, that in any better solution the value of k i has to be smaller thanmax( k i ). This case will be further referred to as the gMin case.In both cases, if the inferred constraint violates the current upper or lower bound ofa k i , then the best policy found up to that point is optimal. Whenever the domain of aswitching point is modified as a result of inferences made during the gMax or gMin case, allof the switching points need to be re-considered. If the domain of one variable is reducedduring a particular shaving iteration, some of the temporary constraints added in the nextround of shaving will be different from the ones used previously, and, consequently, newinferences may be possible. Thus, the shaving procedure terminates when optimality isproved or when no more inferences can be made.Consider the example mentioned above. Suppose the constraint k = 0 is added to theproblem, and all the rest of the variables are assigned their maximum possible values (using gMax ). The resulting policy is ( k , k , k , k ) = (0 , , , B valueof 0.63171, which implies that this policy is feasible because 0 . > .
32 = B l , and sono domain reductions can be inferred. The constraint k = 0 is then removed. Since thedomain of k has not been modified, the procedure considers the next variable. Thus, theconstraint k = 1 is added, and all the other variables are again set to their maximum values.The resulting policy is (0 , , , B value is 0 . k = 1 is then removed, and k = 2 is added. When all the other variables areset to their maximum values, the resulting policy is (0 , , , B valueof 0.1116577, which is smaller than B l . Thus, this policy is infeasible, and the constraint k > k to [3 .. k = 3.Now, consider the gMin case and assume that all variables have their full initial domains.Suppose the constraint k = 3 is added to the problem. All the rest of the variables areassigned their smallest possible values consistent with k = 3. Thus, the policy (3 , , ,
6) isconsidered. This policy has a B value of 0 . K , so if it wereinfeasible, the problem would have been infeasible). The value of k in any better solution erekhov & Beck has to be smaller than 3, and so the domains of the variables become [0 .. .. .. k = 3 is removed, and, since the domain of k has been modified, theconstraint k = 2 is added next. The policy that is considered now is (2 , , , .. .. .. k = 2 is removed, and the next one added is k = 1. The corresponding policy assigned by gMin is infeasible, and no domain reductions are made. Since the addition of k = 1 did notresult in any domain reductions, there is no need to reconsider the variable k until after allother switching points have been looked at. Consequently, the next temporary constraintto be added is k = 4.In the complete B l -based shaving procedure, we can start either with the gMin or the gMax case. Since policies considered in the gMin case will generally have smaller waitingtime than ones considered in the gMax case, it may be beneficial to start with the gMin case. This is the approach we take.Our complete B l -based shaving algorithm is presented in Figure 4. It is assumed in all ofthe algorithms presented that the functions “add( constraint )” and “remove( constraint )”add and remove constraint to and from the model, respectively.Upon the completion of this shaving procedure, the constraint W q ≤ bestW q , where bestW q is the value of the best solution found up to that point, is added ( W q ≤ bestW q rather than W q < bestW q is added because of numerical issues with testing equality offloating point numbers). However, although such a constraint rules out policies with higher W q as infeasible, it results in almost no propagation of the domains of the decision variablesand does little to reduce the size of the search tree. In order to remedy this problem, anothershaving procedure, this time based on the constraint W q ≤ bestW q is proposed in the nextsub-section. The issue of the lack of propagation of the domains of k i from the addition ofthis constraint will be discussed in more detail in Section 8.1. W q -based Shaving Procedure The W q -based shaving procedure makes inferences based strictly on the constraint W q ≤ bestW q : the constraint B ≥ B l is removed prior to running this procedure in order toeliminate the possibility of incorrect inferences. As in B l -based shaving, a constraint of theform k i = max( k i ), where max( k i ) is the maximum value in the domain of k i , is addedtemporarily, and the function gMin is used to assign values to the rest of the variables.Because the B l constraint has been removed, the only reason for the infeasibility of thepolicy is that it has a W q value greater than the best W q that has been encountered so far.Since all switching points except k i are assigned their smallest possible values, infeasibilityimplies that in any solution with a smaller expected waiting time, the value of k i has to bestrictly smaller than max( k i ). This shaving procedure is stated in Figure 5. W q -based and B l -based shaving will result in different domain reductions since they arebased on two different constraints. Moreover, using the two together may cause moredomain modifications than when either is used by itself. Therefore, it makes sense to run the B l -based and W q -based shaving procedures alternately (with W q and B l constraints added CP Approach for a Queueing Control Problem
Algorithm 1: B l -based Shaving Input: S , N , µ , λ , B l (problem instance parameters); bestSolution (best solution foundso far) Output: bestSolution with (possibly) modified domains of the variables k i , or bestSolution with proof of its optimality while (there are domain changes) for all i from 0 to N − while (shaving successful for Domain ( k i ))add( k i = max( Domain ( k i )) ) if ( gMin ) if (new best solution found) bestSolution = currentSolution ; if ( max( Domain ( k i )) − ≥ min( Domain ( k i )) )add( k i < max( Domain ( k i )) ) else return bestSolution ; stop, optimality has been provedremove( k i = max( Domain ( k i )) ) while (shaving successful for Domain ( k i ))add( k i = min( Domain ( k i )) ) if ( gMax ) if (new best solution found) bestSolution = currentSolution ; elseif ( min( Domain ( k i )) + 1 ≤ max( Domain ( k i )) )add( k i > min( Domain ( k i )) ) else return bestSolution ; stop, optimality has been provedremove( k i = min( Domain ( k i )) )Figure 4: B l -based shaving algorithm erekhov & Beck Algorithm 2: W q -based Shaving Input: S , N , µ , λ , B l (problem instance parameters); bestSolution (best solution foundso far) Output: bestSolution with (possibly) modified domains of the variables k i , or bestSolution with proof of its optimality while (there are domain changes) for all i from 0 to N − while (shaving successful for Domain ( k i ))add( k i = max( Domain ( k i )) ) if (! gMin ) if ( max( Domain ( k i )) − ≥ min( Domain ( k i )) )add( k i < max( Domain ( k i )) ) else return bestSolution ; stop, optimality has been provedremove( k i = max( Domain ( k i )) )Figure 5: W q -based shaving algorithmand removed appropriately) until no more domain pruning is possible. Such a combinationof the two shaving procedures will be referred to as AlternatingShaving .The
AlternatingShaving procedure can be effectively combined with search in the follow-ing manner.
AlternatingShaving can be run initially, until no further domain modificationsare possible. Search can then be performed until a better solution is found, at which point
AlternatingShaving can be applied again. Subsequently, search and shaving can alternateuntil one of them proves optimality of the best solution found. Such an approach maybe successful because if search finds a new best solution, a new constraint on W q will beadded, and so W q -based shaving may be able to reduce the upper bounds of the switchingpoint variables. This way of combining search and shaving will be further referred to as AlternatingSearchAndShaving .Other variations of shaving are also possible. In particular, both B l -based and W q -basedshaving procedures can be extended to make inferences about values of two switching points.For example, one can assign maximum values to a pair of switching point variables, whileassigning minimum values to the rest. If the resulting policy is feasible, then a constraintstating that at least one variable from this pair has to be assigned a smaller value can beadded to the problem. Preliminary experiments indicated that shaving procedures basedon two switching points do not, in general, result in more effective models. Such proceduresdo not explicitly reduce the domains of the switching point variables but rather add a setof constraints to the model which do not appear, in practice, to significantly reduce the sizeof the search space. One possible direction for future work may be to further investigatethese variations of shaving. CP Approach for a Queueing Control Problem
6. Experimental Results
Several sets of experiments were performed in order to evaluate the efficiency of the pro-posed models and the effectiveness of dominance rules and shaving procedures, as well as tocompare the performance of the best CP model with the performance of heuristic P . Allconstraint programming models were implemented in ILOG Solver 6.2, while the heuristicof Berman et al. was implemented using C++.We note that numerical results obtained in the experiments are sensitive to the level ofprecision that is set. In all constraint programming models, we set the default ILOG Solverprecision to 0.000001. Doing so implies that all floating point variables in our model areconsidered bound when the maximum ( max ) and minimum ( min ) values of their intervalsare such that (( max − min ) / ( max { , | min |} ) ≤ . The information gained from policies ˆ K and ˆˆ K is explicitly used in the implementation ofall three models. If ˆˆ K is infeasible, then the program stops as there is no feasible solutionfor that instance. Otherwise, if ˆ K is feasible, then it is optimal. The two cases when ˆ K isoptimal or ˆˆ K is infeasible are therefore trivial and are solved easily both by the CP modelsand by Berman et al.’s heuristic. Although instances in which ˆˆ K is optimal are very hard tosolve without shaving, using the elementary B l -based shaving procedure will always resultin a (usually fast) proof of the optimality of this policy. This case is also trivial for Bermanet al.’s heuristic. Consequently, the experimental results presented here are based only onthe instances for which the optimal solution is between ˆ K and ˆˆ K .Preliminary experiments indicated that the value of S has a significant impact on theefficiency of the programs since higher values of S result in larger domains for the k i variablesfor all models and a higher number of w j variables for the Dual model. As indicated inTable 3 in Section 3.5, S also has a big impact on the number of constraints in the If-Then and
Dual models. Therefore, we consider instances for each value of S from theset { , , . . . , } in order to gain an accurate understanding of the performance of themodel and the heuristic. We note that for most instances with S greater than 100, neitherour method nor Berman et al.’s heuristic P may be used due to numerical instability.Thirty instances were generated for each S in such a way as to ensure that the instanceis feasible and that the optimal policy is neither ˆˆ K nor ˆ K . We created random combinationsof parameter values for each S and chose the instances for which policy ˆˆ K was found tobe feasible, but not optimal, and ˆ K was determined to be infeasible. In order to checkthat ˆˆ K is not optimal, it is sufficient to find a feasible solution that has one switchingpoint assigned a value lower than its upper bound. When generating combinations ofparameters, the values of N were chosen with uniform probability from the set {
2, . . . , 38 } ,
3. Numerical values in some of the results are slightly different from the ones presented in our previouswork (Terekhov & Beck, 2007) due to some minor errors discovered after the publication of that paper.The main conclusions and analysis of the previous work remain valid, however.4. Jean-Fran¸cois Puget - personal communication. erekhov & Beck the values of λ from { } , the values of µ from {
1, . . . , 49 } and the values of B l from {
1, . . . , 4 } . There appears to be no easy way of determining whether a given instancewill have ˆˆ K or ˆ K as the optimal solution based only on the parameter values. Moreover,preliminary experiments indicated that problem difficulty depends on a combination ofproblem parameters (especially S , N and B l ) rather than on one parameter only.A 10-minute time limit on the overall run-time of the program was enforced in theexperiments. All experiments were performed on a Dual Core AMD 270 CPU with 1 MBcache, 4 GB of main memory, running Red Hat Enterprise Linux 4. In order to perform comparisons among the CP models, and between the CP models andBerman et al.’s heuristic, we look at mean run-times, the number of instances in whichthe optimal solution was found, the number of instances in which optimality was proved,the number of instances in which the best-known solution was found and the mean relativeerror (MRE). MRE is a measure of solution quality that allows one to observe how quicklya particular algorithm is able to find a good solution. MRE is defined as
M RE ( a, M ) = 1 | M | X m ∈ M c ( a, m ) − c ∗ ( m ) c ∗ ( m ) (22)where a is the algorithm that is used to solve the problem, M is the set of problem instanceson which the algorithm is being tested, c ( a, m ) is the cost of a solution found for instance m by algorithm a , and c ∗ ( m ) is the best-known solution for instance m . As we generatedthe instances, c ∗ ( m ) is the best solution we found during our these experiments. Each CP model was tested with and without shaving and dominance rules. A total of 30CP-based methods were therefore evaluated. A model with B l -based shaving is a modelwhich runs the B l -based shaving procedure until no more domain changes are possible, addsa constraint on the value of W q based on the best solution found during the shaving proce-dure and runs search for the rest of the time. Similarly, models with W q -based shaving and AlternatingShaving are models which run the W q -based shaving procedure and the Alter-natingShaving procedure, respectively, until it is no longer possible to reduce the domainsof the switching point variables, add a constraint requiring W q to be less than the expectedwaiting time of the best solution found by the shaving procedure and use search for therest of the time. As described previously, AlternatingSearchAndShaving alternates betweensearch and the
AlternatingShaving procedure. In all models, search assigns switching pointsin increasing index order. The smallest value in the domain of each variable is tried first.
Table 4 presents the number of instances, out of 300, for which the optimal solution wasfound and proved by each of the 30 proposed CP-based methods. This table indicatesthat the
PSums model is the most effective of the three, proving optimality in the largestnumber of instances regardless of the use of dominance rules and shaving. With
Alter-
CP Approach for a Queueing Control Problem
No Shaving B l -based W q -based Alternating AlternatingSearch
Shaving Shaving
Shaving AndShaving
D ND D ND D ND D ND D ND
If-Then
105 105 192 191 105 105 219 218 234 234
PSums
126 126 202 201 126 126 225 225 238 238
Dual
105 105 191 191 105 105 218 218 232 232Table 4: Number of instances for which the optimal solution is found and its optimality isproved within 10 CPU-minutes out of a total of 300 problem instances (D - withdominance rules, ND - without dominance rules). natingSearchAndShaving , PSums proves optimality in the largest number of instances: in79 .
3% of all instances, 238 out of 239 instances for which optimality has been proved byany model.Figure 6 shows how the MRE changes over the first 50 seconds of run-time for
If-Then , PSums and
Dual models with
AlternatingSearchAndShaving , and for Berman’s heuristic(we comment on the performance of the heuristic in Section 6.4).
PSums is, on average,able to find better solutions than the other two models given the same amount of run-time.In Table 5, additional statistics regarding the performance of the three models with
AlternatingSearchAndShaving and without dominance rules are presented (we comment onthe same statistics for P in Section 6.4). In particular, for each model, the number ofinstances in which it finds the best solution (out of 300), the number of instances in whichit finds the optimal solution (out of 239 cases for which optimality has been proved) andthe number of times it proves optimality (out of 300) are presented. It can be seen thatall models find the optimal solution in the 239 instances for which it is known. However,the PSums model proves optimality in 4 more instances than the
If-Then model and in 6more instances than the
Dual . PSums also finds the best-known solution of any algorithmin 97 .
6% of all the instances considered. A more detailed discussion of the differences inthe performance of the CP models is presented in Section 8.2.
From Table 4, it can be observed that the CP models without shaving and with the W q -based shaving procedure prove optimality in the fewest number of cases. The similarityin performance of models without shaving and with W q -based shaving is not surprisingbecause the W q -based procedure is able to start pruning domains only when the value ofthe best policy found prior to this procedure is quite good. When the W q -based procedureis used alone, it only has one solution to base its inferences on, namely ˆˆ K . Since all policieswill result in a smaller expected waiting time than ˆˆ K , this procedure by itself is useless.Employing the B l -based shaving procedure substantially improves the performance ofall models: without dominance rules, the If-Then , the
PSums and the
Dual models proveoptimality in 86, 75 and 86 more instances, respectively, than the corresponding modelswithout shaving or with only W q -based shaving; with dominance rules, the situation is erekhov & Beck . . . . . . Run Time (seconds) M ean R e l a t i v e E rr o r . . . . . . . . . . . . If−Then ModelPSums ModelDual ModelBerman et al.’s Heuristic P1
Figure 6: Comparison of MRE of three CP models with
AlternatingSearchAndShaving withBerman’s heuristic P .equivalent. These results imply that inferences made based on the B l constraint are effectivein reducing the domains of the decision variables.Models with AlternatingShaving and
AlternatingSearchAndShaving perform even bet-ter than models employing only the B l -based shaving procedure. The real power of W q -based shaving only becomes apparent when it is combined with B l -based shaving because B l -based shaving often finds a good solution and a good value of bestW q , allowing the W q -based procedure to infer more domain reductions. This observation explains why Al-ternatingSearchAndShaving performs better than
AlternatingShaving . In particular, in
Al-ternatingSearchAndShaving , the W q -based procedure is used both after a new best solution CP Approach for a Queueing Control Problem
PSums
293 239 238
If-Then
282 239 234
Dual
279 239 232 P
282 238 0Table 5: Comparison of three CP models (with
AlternatingSearchAndShaving and withoutdominance rules) with Berman’s heuristic P .is found during shaving and after one is found during search. Therefore, if a higher qualitysolution is found during search, it will be used by the W q -based procedure to further prunethe domains of the switching point variables.In Figure 7, the average run-times for each value of S from 10 to 100 are presented for thefour shaving procedures with the PSums model. Since a run-time limit of 600 seconds wasused throughout the experiments, we assumed a run-time of 600 seconds for all instancesfor which optimality has not been proved within this limit. Therefore, the mean run-timesreported throughout this paper are underestimates of the true means. Figure 7 shows that,for each value of S , the AlternatingSearchAndShaving procedure gives the best performance.It can also be seen that, as S increases, it becomes increasingly difficult to prove optimalityand average run-times increase. As stated previously, this is due to larger domains of theswitching point variables. The AlternatingSearchAndShaving procedure, however, is ableto significantly reduce the domains of the k i variables and therefore provides an effectivemethod for instances with higher values of S as well. Table 4 indicates that there is rarely any difference in the number of instances solved tooptimality between models with and without dominance rules. No difference at all is visiblefor any model without shaving, with W q -based shaving and AlternatingSearchAndShaving .Recall that dominance rules are implemented by the addition of a constraint on thevalues of the switching point variables after a solution is found. Such a constraint will bemore effective when more of the switching point variables are assigned their minimum valuesin the current solution. Usually, such policies are also the ones which result in a smallerexpected waiting time. Similarly, W q -based shaving is useful only when a solution withsmall expected waiting time is found. This leads to the conjecture that dominance rulesmay be effective only for the same instances for which the W q -based shaving procedure iseffective. This conjecture is supported by the results of Table 4. In particular, when the W q -based shaving procedure is used by itself, it makes inferences based only on the policyˆˆ K , a solution that is generally of poorest quality for an instance. The method with a singlerun of W q -based shaving therefore heavily relies on search. Since search takes a long time tofind a feasible solution of good quality, the effectiveness of dominance rule-based constraintsis also not visible within the given time limit. erekhov & Beck
20 40 60 80 100
Values of S A v e r age R un − T i m e s ( s e c ond s )
20 40 60 80 100
20 40 60 80 100
20 40 60 80 100
Wq−based ShavingBl−based ShavingAlternatingShavingAlternatingSearchAndShaving
Figure 7:
PSums model with various shaving techniques: average run-times for each valueof S . Average run-times for PSums without shaving are not shown in this graphsince the resulting curve would be indistinguishable from the one for W q -basedshaving.On the other hand, in AlternatingSearchAndShaving , the W q -based procedure plays akey role because it makes domain reductions based on high quality solutions produced by B l -based shaving, and later, search. Dominance rules do not play a role in this proceduresince shaving is used after every new solution is found. However, even if dominance ruleconstraints were explicitly incorporated in the procedure (i.e. if they were added beforeeach new run of search), they would be redundant since they serve essentially the samepurpose as the W q -based shaving procedure.When shaving is not used, the results are equivalent to those achieved when W q -basedshaving is employed. The explanation for the absence of a difference between models with CP Approach for a Queueing Control Problem
20 40 60 80 100 e − e − e − e − e − Values of S M ean R e l a t i v e E rr o r
20 40 60 80 100 e − e − e − e − e −
20 40 60 80 100 e − e − e − e − e −
20 40 60 80 100 e − e − e − e − e − PSums Model At 10 secondsPSums Model At 150 secondsPSums Model At 500 secondsBerman et al.’s Heuristic P1
Figure 8: MRE for each value of S for P and the P Sums model.and without dominance rules is therefore also the same. In particular, it takes a long timefor a solution to be found whose quality is such that it allows the dominance rule constraintto effectively reduce the size of the search tree.When B l -based shaving and AlternatingShaving are used, dominance rules are some-times helpful. In both cases, this is because after these two shaving procedures, subsequentsearch usually finds a good solution quickly, and, since W q -based shaving is not used againat that point, the dominance rule constraint that is added can be effective in reducing thesize of the search tree.Overall, it can be observed that using AlternatingSearchAndShaving without dominancerules is more effective than using B l -based shaving or AlternatingShaving with dominance erekhov & Beck rules. Therefore, in further comparisons, the focus is only on models with
AlternatingSearchAndShaving without dominance rules. P vs. Best Constraint Programming Approach Empirical results regarding the performance of heuristic P are not presented by Bermanet al. (2005), and so the ability of P to find good switching policies has not been explicitlyevaluated in previous work. We wanted to find out how well the heuristic actually performsby comparing it to our CP methods.In Table 5, we present several measures of performance for the three proposed modelswith AlternatingSearchAndShaving and for the heuristic P . The heuristic performs well,finding the best-known solution in only eleven fewer instances than the PSums model, inthree more instances than the
Dual model and in the same number of instances as the
If-Then model. Moreover, the heuristic finds, but, of course, cannot prove, the optimalsolution in 238 out of 239 instances for which the optimal is known. The three CP modelsfind the optimal solution in all 239 of these. The run-time of the heuristic is negligible,whereas the mean run-time of the
PSums model is approximately 130 seconds (the meanrun-times of the other two models are slightly higher: 141 seconds for the
If-Then modeland 149 seconds for the
Dual model).Table 5 also shows that the
PSums model is able to find the best-known solution in 11more instances than the heuristic. Closer examination reveals that there are 275 instancesin which both the
PSums model and P find the best-known solution, 18 instances in whichonly PSums is able to do so and 7 in which only the heuristic finds the best-known.From Figure 6, it can be observed that the heuristic achieves a very small MRE in anegligible amount of time. After about 50 seconds of run-time, the MRE over 300 instancesresulting from
PSums with
AlternatingSearchAndShaving becomes comparable to that ofthe heuristic MRE. In Figure 8, the MRE over 30 instances for each value of S is presentedfor the heuristic and for PSums with
AlternatingSearchAndShaving at 10, 150 and 500seconds of run-time. After 10 seconds, the performance of
PSums is comparable to thatof the heuristic for values of S smaller than or equal to 40, but the heuristic appears tobe quite a bit better for higher values of S . At 150 seconds, the performance of PSums iscomparable to that of the heuristic except at S values of 50 and 80. After 500 seconds, PSums has a smaller MRE over the 300 instances and also a lower (or equal) MRE for eachvalue of S except 50 and 100.Overall, these results indicate that the heuristic performs well–its run-time is negligible,it finds the optimal solution in all but one of the cases for which it is known, and it findsthe best solution in 94% of all instances. Moreover, it results in very low MRE. Although PSums with
AlternatingSearchAndShaving is able to achieve slightly higher numbers inmost of the performance measures, it is clear that these improvements are small given thatthe
PSums run-time is so much higher than the run-time of the heuristic. PSums - P Hybrid
Naturally, it is desirable to create a method that would be able to find a solution of highquality in a very short amount of time, as does Berman’s heuristic, and that would alsohave the same high rate of being able to prove optimality within a reasonable run-time as
CP Approach for a Queueing Control Problem
PSums
293 239 238 P
282 238 0
PSums - P
300 239 238HybridTable 6: Comparison of
PSums model with
AlternatingSearchAndShaving and Berman etal.’s heuristic P with the Hybrid model.does PSums with
AlternatingSearchAndShaving . It is therefore worthwhile to experimentwith a
PSums - P Hybrid, which starts off by running P and then, assuming the instance isfeasible, uses the PSums model with
AlternatingSearchAndShaving to find a better solutionor prove the optimality of the solution found by P (infeasibility of an instance is proved ifthe heuristic determines that policy ˆˆ K is infeasible).Since it was shown that heuristic P is very fast, running it first incurs almost nooverhead. Throughout the analysis of experimental results, it was also noted that theperformance of the W q -based shaving procedure depends on the quality of the best solutionfound before it is used. We have shown that the heuristic provides solutions of very highquality. Therefore, the first iteration of the W q -based procedure may be able to significantlyprune the domains of switching point variables because of the good-quality solution foundby the heuristic. Continuing by alternating the two shaving techniques and search, whichhas also been shown to be an effective approach, should lead to good results.The proposed Hybrid algorithm was tested on the same set of 300 instances that wasused above. Results illustrating the performance of the Hybrid as well as the performance of P and PSums with
AlternatingSearchAndShaving are presented in Table 6. The Hybrid isable to find the best solution in all 300 cases: in the 275 instances in which both the heuristicand
PSums find the best-known solution, in 18 in which only
PSums finds the best-knownand in 7 in which only the heuristic does so. The Hybrid finds the optimal solution (for thoseinstances for which it is known) and proves optimality in as many instances as the
PSums model. The mean run-time for the Hybrid is essentially identical to the mean run-time of
PSums with
AlternatingSearchAndShaving , equalling approximately 130 seconds.Thus, the Hybrid is the best choice for solving this problem: it finds as good a solutionas the heuristic in as little time (close to 0 seconds), it is able to prove optimality in asmany instances as the best constraint programming method, and it finds the best-knownsolution in all instances considered. Moreover, all these improvements are achieved withoutan increase in the average run-time over the
PSums model.
8. Discussion
In this section, we examine some of the reasons for the poor performance of the CP modelswithout shaving, suggest reasons for the observed differences among the CP models, discussthe performance of the Hybrid and present some perspectives on our work. erekhov & Beck
In our experiments, we have some instances for which even the
PSums - P Hybrid with
AlternatingSearchAndShaving is unable to find and prove the optimal solution within the10-minute time limit. In fact, in many of these instances, the amount of time spent duringsearch is higher than the time spent on shaving, and the run-time limit is usually reachedduring the search, rather than the shaving, phase. Further analysis of the algorithms’behaviour suggests that this poor performance of search can be explained by the lack of back-propagation. Back-propagation refers to the pruning of the domains of the decision variablesdue to the addition of a constraint on the objective function: the objective constraintpropagates “back” to the decision variables, removing domain values and so reducing search.In the CP models presented above, there is very little back-propagation.Consider a model without shaving. Throughout search, if a new best solution is found,the constraint W q ≤ bestW q , where bestW q is the new objective value, is added to themodel. However, the domains of the switching point variables are usually not reduced inany way after the addition of such a constraint. This can be illustrated by observing theamount of propagation that occurs in the PSums model when W q is constrained.For example, consider an instance of the problem with S = 6, N = 3, λ = 15, µ = 3,and B l = 0 .
32 (this instance is used in Section 5 to illustrate the shaving procedures). Theinitial domains of the switching point variables are [0 .. .. ..
5] and [6]. The initialdomains of the probability variables P ( k i ) for each i , after the addition of W q boundsprovided by ˆ K and ˆˆ K , are listed in Table 7. The initial domain of W q , also determined bythe objective function values of ˆ K and ˆˆ K , is [0 . .. . L and F , are [2 . e − ..
6] and [0 .. . W q ≤ . . W q is reduced to [0 . .. . L becomes [1 . ..
6] and thedomain of F remains [0 .. . P ( k i ) after this addition are listed inTable 7. The domains of both types of probability variables are reduced by the additionof the new W q constraint. However, the domains of the switching point variables remainunchanged. Therefore, even though all policies with value of W q less than 0 . W q to be less than or equal to this value does not result in anyreduction of the search space. It is still necessary to search through all policies in order toshow that no better feasible solution exists.One of the reasons for the lack of pruning of the domains of the k i variables due to the W q constraint is likely the complexity of the expression W q = Lλ (1 − P ( k N )) − µ . In the exampleabove, when W q is constrained to be less than or equal to 0.306323, we get the constraint0 . ≥ L − P ( k N )) − , which implies that 9 . − P ( k N )) ≥ L . This explains whythe domains of both L and P ( k N ) change upon this addition to the model. The domains ofthe rest of the P ( k i ) variables change because of the relationships between P ( k i )’s (Equation(15)) and because of the constraint that the sum of all probability variables has to be 1.Similarly, the domains of P Sums ( k i )’s change because these variables are expressed interms of P ( k i ) (Equation (14)). However, because the actual k i variables mostly occur asexponents in expressions for P Sums ( k i ), P ( k i ), and L ( k i ), the minor changes in the domainsof P Sums ( k i ), P ( k i ), or L ( k i ) that happen due to the constraint on W q have no effect onthe domains of the k i . This analysis suggests that it may be interesting to investigate a CP CP Approach for a Queueing Control Problem
Before addition of W q ≤ . W q ≤ . j P ( j ) P Sums ( j ) P ( j ) P Sums ( j ) k [4 . e − .. . ..
1] [4 . e − .. . .. . k [1 . e − ..
1] [0 ..
1] [0 . ..
1] [0 .. . k [2 . e − .. .
6] [2 . e − ..
1] [0 . .. . . .. . k [4 . e − ..
1] N/A [0 . .. . P ( j ) and P Sums ( j ) variables for j = k , k , k , k , before and afterthe addition of the constraint W q ≤ . If-Then and
Dual models, the domains of the decision variables are notreduced when a bound on the objective function value is added, although the domains of allprobabilities, L and F are modified. In both of these models, the constraints relating F , L and the probability variables to the variables k i are the balance equations, which are quitecomplex. The domains of the probability variables do not seem to be reduced significantlyenough due to the new W q bound so as to result in the pruning of k i domains because ofthese constraints.These observations served as the motivation for the proposed shaving techniques. Inparticular, the W q -based shaving procedure reduces the domains of switching point variableswhen it can be shown that these values will necessarily result in a higher W q value than thebest one found up to that point. This makes up for some of the lack of back-propagation.However, even when this procedure is used after each new best solution is found, as in AlternatingSearchAndShaving , it is not always able to prune enough values from the domainsof the k i ’s so as to be able to prove optimality within 10 minutes of run-time. It cantherefore be seen that inferences based on the value of W q are very limited in their powerand, therefore, if the domains of switching point variables are large after shaving, then itwill not be possible to prove optimality in a short period of time. Experimental results demonstrate that the
PSums model is the best out of the three modelsboth with and without shaving. In this section, we examine the models in more detail inan attempt to understand the reasons for such differences.
PSums and the other two models
In order to analyze the performance of the models without shaving, we look at the meannumber of choice points statistics, which give an indication of the size of the search spacethat is explored. To compare our three models, we look at the mean number of choicepoints considered before the first feasible solution is found and the mean total number ofchoice points explored within 600 seconds of run-time. erekhov & Beck
21 Instances Solved by
PSums
Only 105 Instances Solved by All ModelsFirst Solution Total First Solution Total
If-Then
PSums
Dual
PSums , the latter statisticcorresponds to the total number of choice points needed to prove optimality.In Table 4, it is shown that, without shaving, the
PSums model proves optimality in21 more instances than the other two models and that there are 105 instances in which allthree models prove optimality. In Table 8, we present the mean number of choice pointsstatistics for all three models for both of these sets of instances. It can be seen that themean number of choice points that need to be explored by the
PSums model in order tofind an initial solution is smaller than for the other two models, both for the 105 instancesthat are eventually solved to optimality by all models and for the 21 instances that are onlysolved to optimality by
PSums . Because the same variable and value ordering heuristicsare used in all models, this observation implies that more propagation occurs during searchin the
PSums model than in the other two models. This claim is further supported by thefact that the mean total number of choice points for the 105 instances that are solved byall models is smaller for
PSums than for the other two models.Table 8 also shows that, for the 21 instances that are only solved to optimality by
PSums , the mean total number of choice points is the highest for the
PSums model. Since
PSums is the only one out of the three models to solve these instances, this implies thatpropagation is happening faster in this model. This observation is confirmed by the resultsfrom the 105 instances that are solved by all three models: for these instances, the
Dual explores an average of 713 choice points per second, the
If-Then model explores an averageof 895 choice points per second and the
PSums model explores an average of 1989 choicepoints per second. In other words, it appears that propagation in the
PSums model is morethan twice as fast as in the other two models.A more detailed examination of the results showed that in 82 out of the 105 instancesthat were solved by all models, the number of choice points explored, for a given instance, isthe same for all models. Moreover, for instances which have the same value of S and N , thenumber of choice points explored is equal. In Figure 9, the run-times of the three modelsas the number of choice points increases is illustrated. In order to create this graph, weaveraged the run-times of all instances for which the number of choice points examined isthe same. Some points in the figure are labeled as ( S, N ) in order to show the relationshipbetween the number of choice points, the values of S and N , and the run-times. We notethat there is one instance, when S = 10 and N = 6, for which the number of choice pointsis the same as for the instances when S = 10 and N = 4. However, for all other instances CP Approach for a Queueing Control Problem
Choice Points R un − T i m e ( s e c ond s ) If−Then ModelPSums ModelDual Model(50,3)(60,3)(20,6)(70,2) (70,3) (20,7)(40,4) (20,8)(30,5) (20,9)
Figure 9: Run-times averaged over all instances with equal number of choice points ex-plored, for 82 instances for which the number of choice points is the same for allmodels. The labels in the points indicate the ( S , N ) values for the correspondinginstances.(out of 82), there is a one-to-one correspondence between ( S, N ) and the number of choicepoints.Several observations can be made from Figure 9. Firstly, this graph demonstrates thatpropagation in the
PSums model is faster than in the other models. Secondly, the behaviourof the
PSums model appears to be quite different from that of both the
If-Then and the
Dual models. The run-times of the
PSums model seem to be significantly influenced by thevalue of N . For example, when S = 20, the run-times of this model increase as N increases erekhov & Beck from 6 to 9. Moreover, given two instances of which one has a high S and a low N , andthe other has a low S and a high N , the PSums model generally needs a longer time tosolve the instance with a low S and a high N (e.g., compare the points (20, 7) and (40, 4),or (20, 7) and (70, 3)). For the If-Then and the
Dual models, there are several cases whenthe run-times for instances with a high S and a low N are higher than for instances with alow S and a high N (e.g., compare the run-time at (70, 3) with that of (20, 7)), althoughthe opposite happens as well (e.g., compare (50, 3) and (40, 4)). Thus, it appears that N is the parameter influencing the run-times of PSums the most, while for other two models,both S and N are influential, with S having a greater effect. Although these characteristicsrequire additional investigation, one possible reason for such differences in model behaviourcould be the relationship between the number of constraints in the models and the problemparameters. From Table 3, it is known that the number of constraints is mostly determinedby the value of S in the If-Then and
Dual models (since S is typically larger than N ), andby the value of N in the PSums model. Combining observations from Table 3 and Figure9, it appears that the effect of N and S on the run-times is due to their influence on thenumber of constraints in the models.Overall, this examination indicates that the superiority of the PSums model withoutshaving is caused by its stronger propagation (Table 8) and by the fact that propagation isfaster (Figure 9).When shaving is employed, the
PSums model also performs better than the
Dual and
If-Then models, proving optimality in a greater number of instances (Table 4) and findinggood-quality solutions faster (Figure 6). In all models, the shaving procedures make thesame number of domain reductions because shaving is based on the W q and B l constraints,which are present in all models. However, the time that each shaving iteration takes isdifferent in different models. Our empirical results show that each iteration of shavingtakes a smaller amount of time with the PSums model than with the
If-Then and
Dual models. Thus, with shaving, the
PSums model performs better than the other two, bothbecause shaving is faster and because subsequent search, if it is necessary, is faster.
If-Then and the
Dual models
A comparison of the
If-Then model with
AlternatingSearchAndShaving with the
Dual with
AlternatingSearchAndShaving using Figure 6 shows that the
If-Then model is usually ableto find good solutions in a smaller amount of time. Moreover, as shown in Table 5, the
If-Then model with
AlternatingSearchAndShaving finds the best solution in three moreinstances, and proves optimality in two more instances, than the
Dual model with the sameshaving procedure. With other shaving procedures and without shaving, the same statisticsshow almost no difference in the performance of the two models. It was expected thatthe
Dual would outperform the
If-Then model because it uses a simpler representation ofthe balance equations and expressions for F and L , and has a smaller number of if-thenconstraints. (there are Table 8 shows that the Dual has to explore a smaller number ofchoice points to find the initial solution. In the 105 instances that both of these modelssolve, the total number of choice points explored by the
Dual is also smaller. However, the
If-Then model is faster, exploring, on average, 895 choice points per second compared tothe average of 713 choice points per second explored by the
Dual . One possible explanation
CP Approach for a Queueing Control Problem for the
Dual being slower is the fact that it has to assign more variables (via propagation)than the other models. In particular, in order to represent a switching policy, the
Dual hasto assign S + 1 w j variables in addition to N + 1 k i variables, with S usually being muchlarger than N . PSums - P Hybrid
Experimental results demonstrate that the
PSums - P Hybrid finds good solutions quicklyand is able to prove optimality in a large number of instances. It should be noted however,that no synergy results from the combination: the number of instances for which optimalityis proved does not increase and the run-times do not decrease. Moreover, the hybrid modelfinds the best-known solution in all test instances simply because there are cases in whichonly the
PSums model or only the heuristic is able to do so. No new best solutions areobtained by using the
PSums - P Hybrid to solve the problem. It appears that startingthe
PSums model from the solution found by the heuristic does not lead to a significantincrease in the amount of propagation. Also, the fact that the heuristic finds a good-qualitysolution does not improve the overall performance since, if search has to be used, placing aconstraint on W q that requires all further solutions to be of better quality has little effecton the domains of the decision variables. These observations imply that in order to createa more effective model for the problem, one would need to improve back-propagation byadding new constraints or reformulating the existing ones. If back-propagation is improved,the good-quality heuristic solution may result in better performance of the hybrid approach. The CP methods that we have developed are, in some ways, non-standard. More commonapproaches when faced with the poor results of our three basic CP models (without shaving)would have been to create better models, to develop a global constraint that could representand efficiently reason about some relevant sub-structure in the problem, and/or to inventmore sophisticated variable- and value-ordering heuristics. Shaving is more a proceduraltechnique that must be customized to exploit a particular problem structure. In contrast, abetter model or the creation of a global constraint are more in-line with the declarative goalsof CP. Our decision to investigate shaving arose from the recognition of the need to moretightly link the optimization function with the decision variables and the clear structure ofthe problem that both appeared and proved to be ideal for shaving.We believe that there is scope for better models and novel global constraints. Modellingproblem P using CP was not straightforward because the formulation proposed by Bermanet al. contains expressions such as Equation (6), where the upper and lower limits of a sumof auxiliary variables are decision variables. Such constraints are not typical for problemsusually modelled and solved by CP and there appear to be no existing global constraintsthat could have been used to facilitate our approach. In spite of these issues, our modelsdemonstrate that CP is flexible enough to support queueing constraints. However, webelieve it is likely that the generalized application of CP to solve a larger class of queueingcontrol problems will require global constraints specific to expressions commonly occurringin queueing theory. Given the back-propagation analysis and the fact that the problem is erekhov & Beck to find and prove optimality, we are doubtful that, for P , sophisticated search heuristicswill perform significantly better than our simple heuristics.While this is both the first time that CP has been used to solve any queueing controlproblem and the first time that instances of P have been provably solved to optimality,the work in this paper can be viewed as somewhat narrow: it is a demonstration that aparticular queueing control problem can be solved by constraint programming techniques.The work does not immediately deliver solutions to more general problems, however, webelieve that it does open up a number of directions of inquiry into such problems.1. It appears that there is no standard method within queueing theory to address queue-ing control optimization problems. This first application opens the issue of whetherCP can become an approach of choice for such problems.2. As noted in Section 1, there is increasing interest in incorporating reasoning aboutuncertainty into CP-based problem solving. Queueing theory can provide formula-tions that allow direct calculation of stochastic quantities based on expectation. Thechallenge for CP is to identify the common sub-structures in such formulations anddevelop modelling, inference, and search techniques that can exploit them.3. Challenging scheduling problems, such as staff rostering at call centres (Cezik &L’Ecuyer, 2008), consist of queues as well as rich resource and temporal constraints(e.g., multiple resource requirements, alternative resources of different speeds, taskdeadlines, precedence relations between tasks). We believe that a further integrationof CP and queueing theory could prove to be a promising approach to such problems.4. The ability to reason about resource allocation under uncertainty is an importantcomponent of definitions of intelligent behaviour such as bounded rationality (Simon,1997). While we cannot claim to have made significant contribution in this direction,perhaps ideas from queueing theory can serve as the inspiration for such contributionsin the future.
9. Related Work and Possible Extensions
Several papers exist that deal with similar types of problems as the one considered here.For example, Berman and Larson (2004) study the problem of switching workers betweentwo rooms in a retail facility where the customers in the front room are divided into twocategories, those “shopping” in the store and those at the checkout. Palmer and Mitrani(2004) consider the problem of switching computational servers between different types ofjobs where the randomness of user demand may lead to unequal utilization of resources.Batta, Berman and Wang (2007) study the problem of assigning cross-trained customerservice representatives to different types of calls in a call centre, depending on estimateddemand patterns for each type of call. These three papers provide examples of problems forwhich CP could prove to be a useful approach. Investigating CP solutions to these problemsis therefore one possible direction of future work.Further work may also include looking at extensions of the problem discussed in thispaper. For example, we may consider a more realistic problem in which there are resourceconstraints for one or both of the rooms, or in which workers have varying productivity.
CP Approach for a Queueing Control Problem
Another direction for further work is improvement of the proposed models. In particu-lar, in all models, but especially in
PSums , there are constraints with variable exponents.One idea for improving the performance of these constraints is to explicitly represent thedifferences between switching points (i.e., k i +1 − k i ) as variables. Another idea is to in-vestigate a model based on the logarithms of probabilities rather than the probabilitiesthemselves. Additionally, ways of increasing the amount of back-propagation need to beexamined.The goal of this paper was to demonstrate the applicability of constraint programmingto solving a particular queueing control problem. The main direction for future work is,therefore, to explore the possibility of further integrating CP and queueing theory in anattempt to address stochastic scheduling and resource allocation problems. Such problemsare likely to involve complex constraints, both for encoding the necessary stochastic infor-mation and for stating typical scheduling requirements such as task precedences or resourcecapacities. Combining queueing theory with CP may help in solving such problems.
10. Conclusions
In this paper, a constraint programming approach is proposed for the problem of findingthe optimal states to switch workers between the front room and the back room of a re-tail facility under stochastic customer arrival and service times. This is the first work ofwhich we are aware that examines solving such stochastic queueing control problems usingconstraint programming. The best pure CP method proposed is able to prove optimalityin a large proportion of instances within a 10-minute time limit. Previously, there existedno non-heuristic solution to this problem aside from naive enumeration. As a result of ourexperiments, we hybridized the best pure CP model with the heuristic proposed for thisproblem in the literature. This hybrid technique is able to achieve performance that isequivalent to, or better than, that of each of the individual approaches alone: it is able tofind very good solutions in a negligible amount of time due to the use of the heuristic, and isable to prove optimality in a large proportion of problem instances within 10 CPU-minutesdue to the CP model.This paper demonstrates that constraint programming can be a good approach for solv-ing a queueing control optimization problem. For queueing problems for which optimality isimportant or heuristics do not perform well, CP may prove to be an effective methodology.
Appendix A. Constraint Derivations
In this section, the derivations of constraints in the
PSums model and expressions for theauxiliary variables and constraints are presented.
A.1 Closed-form Expressions in the
PSums model
The
PSums model has two sets of probability variables, P ( k i ), i = 0 , , . . . , N , the proba-bility of there being k i customers in the front room, and P Sums ( k i ), i = 0 , , . . . , N −
5. Thanks to an anonymous reviewer for this suggestion. erekhov & Beck explicitly stated in this model. However, expressions for P ( k i ) and P Sums ( k i ) are de-rived in such a way that the balance equations are satisfied. The technique used for thesederivations is similar to that used by Berman et al. (2005) to simplify the calculation ofprobabilities.Consider the balance equation P ( j ) λ = P ( j + 1) iµ , which is true for j = k i − , k i − +1 , . . . , k i − i ∈ { , . . . , N } . In particular, this subset of the balance equations is P ( k i − ) λ = P ( k i − + 1) iµP ( k i − + 1) λ = P ( k i − + 2) iµ ... P ( k i − λ = P ( k i ) iµ. These equations imply the following expressions: P ( k i − ) λiµ = P ( k i − + 1) (23) P ( k i − + 1) λiµ = P ( k i − + 2)... P ( k i − λiµ = P ( k i ) . Combining these together, we get P ( k i ) = (cid:16) λiµ (cid:17) k i − k i − P ( k i − ) for all i from 1 to N , or P ( k i +1 ) = (cid:18) λ ( i + 1) µ (cid:19) k i +1 − k i P ( k i ) , ∀ i ∈ { , , . . . , N − } . (24)This equation is included in the PSums model and has previously been stated as Equation(15).Similarly, from Equation (23), we see that P ( k i − + 1) = λiµ P ( k i − ) for all i from 1 to N , or P ( k i + 1) = λ ( i + 1) µ P ( k i ) , ∀ i ∈ { , , . . . , N } . (25)Using Equation (25), an expression for P Sums ( k i ), the sum of probabilities P ( j ) for j between k i and k i +1 −
1, can be derived as follows:
P Sums ( k i ) = k i +1 − X j = k i P ( j )= P ( k i ) + P ( k i + 1) + P ( k i + 2) + . . . + P ( k i +1 − P ( k i ) + P ( k i ) λ ( i + 1) µ + P ( k i ) (cid:18) λ ( i + 1) µ (cid:19) CP Approach for a Queueing Control Problem + . . . + P ( k i ) (cid:18) λ ( i + 1) µ (cid:19) k i +1 − k i − = P ( k i ) " λ ( i + 1) µ + (cid:18) λ ( i + 1) µ (cid:19) + . . . + (cid:18) λ ( i + 1) µ (cid:19) k i +1 − k i − = P ( k i ) 1 − (cid:20) λ ( i + 1) µ (cid:21) k i +1 − k i − λ ( i + 1) µ if λ ( i +1) µ = 1 P ( k i )( k i +1 − k i ) otherwise. (26)The last step in the derivation is based on the observation that the expression [1 + λ ( i +1) µ + (cid:16) λ ( i +1) µ (cid:17) + . . . + (cid:16) λ ( i +1) µ (cid:17) k i +1 − k i − ] is a geometric series with the common ratio λ ( i +1) µ .When λ ( i +1) µ is 1, the expression is simply a sum of k i +1 − k i ones. The expression for P Sums ( k i ) has been previously stated as Equation (14). A.1.1 Expected Number of Workers Constraint F can be expressed in terms of P ( k i ) and P Sums ( k i ) using the following sequence of steps: F = N X i =1 k i X j = k i − +1 iP ( j )= N X i =1 i [ P ( k i − + 1) + P ( k i − + 2) + . . . + P ( k i −
1) + P ( k i )]= N X i =1 i [ P Sums ( k i − ) − P ( k i − ) + P ( k i )] . (27) A.1.2 Expected Number of Customers Constraints
The equation for L can be derived in a similar manner: L = k N X j = k jP ( j )= k − X j = k jP ( j ) + k − X j = k jP ( j ) + . . . + k N − X j = k N − jP ( j ) + k N P ( k N )= L ( k ) + L ( k ) + . . . + L ( k N − ) + k N P ( k N )= N − X i =0 L ( k i ) + k N P ( k N ) (28) erekhov & Beck where L ( k i ) = k i P ( k i ) + ( k i + 1) P ( k i + 1) + ( k i + 2) P ( k i + 2) + . . . + ( k i +1 − P ( k i +1 − k i P ( k i ) + k i P ( k i + 1) + k i P ( k i + 2) + . . . + k i P ( k i +1 −
1) + P ( k i + 1)+ 2 P ( k i + 2) + . . . + ( k i +1 − k i − P ( k i +1 − k i [ P ( k i ) + P ( k i + 1) + P ( k i + 2) + . . . + P ( k i +1 − P ( k i + 1)+ 2 P ( k i + 2) + . . . + ( k i +1 − k i − P ( k i +1 − k i P Sums ( k i ) + P ( k i ) λ ( i + 1) µ + 2 P ( k i ) (cid:18) λ ( i + 1) µ (cid:19) + . . . + ( k i +1 − k i + 1) P ( k i ) (cid:18) λ ( i + 1) µ (cid:19) k i +1 − k i − = k i P Sums ( k i ) + P ( k i ) λ ( i + 1) µ × " λ ( i + 1) µ + 3 (cid:18) λ ( i + 1) µ (cid:19) + . . . + ( k i +1 − k i − P ( k i ) (cid:18) λ ( i + 1) µ (cid:19) k i +1 − k i − = k i P Sums ( k i ) + P ( k i ) λ ( i + 1) µ k i +1 − k i − X n =0 n (cid:18) λ ( i + 1) µ (cid:19) n − (29)= k i P Sums ( k i ) + P ( k i ) λ ( i + 1) µ × (cid:16) λ ( i +1) µ (cid:17) k i +1 − k i − ( k i − k i +1 ) + (cid:16) λ ( i +1) µ (cid:17) k i +1 − k i ( k i +1 − k i −
1) + 1 (cid:16) − λ ( i +1) µ (cid:17) . A.2 Auxiliary Variables
All constraint programming models contain Equation (13) (restated below as Equation(30)) for defining the βSum ( k i ) variables, which are necessary for expressing an auxiliaryconstraint that ensures that the balance equations have a unique solution. The validity ofthis equation is proved by the following derivation, which uses the formula for the sum of afinite geometric series in the last step: βSum ( k i ) = k i X j = k i − +1 β j = (cid:18) λµ (cid:19) k i − +1 − k (cid:18) i (cid:19) k i − +1 − k i − X i + (cid:18) λµ (cid:19) k i − +2 − k (cid:18) i (cid:19) k i − +2 − k i − X i + . . . + (cid:18) λµ (cid:19) k i − k (cid:18) i (cid:19) k i − k i − X i CP Approach for a Queueing Control Problem = X i (cid:18) λµ (cid:19) k i − − k +1 (cid:18) i (cid:19) × " λµ i + (cid:18) λµ (cid:19) (cid:18) i (cid:19) + . . . + (cid:18) λµ (cid:19) k − k − ( k i − − k +1) (cid:18) i (cid:19) k i − k i − − = X i (cid:18) λµ (cid:19) k i − − k +1 (cid:18) i (cid:19) k i − k i − − X n =0 (cid:18) λiµ (cid:19) n = X i (cid:18) λµ (cid:19) k i − − k +1 (cid:18) i (cid:19) − λiµ ! ki − ki − − λiµ ! if λiµ = 1 X i (cid:18) λµ (cid:19) k i − − k +1 (cid:18) i (cid:19) ( k i − k i − ) otherwise. (30) Acknowledgments
This research was supported in part by the Natural Sciences andEngineering Research Council and ILOG, S.A. Thanks to Nic Wilson and Ken Brown fordiscussions and comments on this work, and to Tom Goguen for careful proofreading ofthe final copy. A preliminary version of parts of this work has been previously published(Terekhov & Beck, 2007).
References
Baptiste, P., Le Pape, C., & Nuijten, W. (2001).
Constraint-based Scheduling . KluwerAcademic Publishers.Batta, R., Berman, O., & Wang, Q. (2007). Balancing staffing and switching costs in aservice center with flexible servers.
European Journal of Operational Research , ,924–938.Beck, J. C., & Prestwich, S. (2004). Exploiting dominance in three symmetric problems. In Fourth International Workshop on Symmetry and Constraint Satisfaction Problems .Beck, J. C., & Wilson, N. (2007). Proactive algorithms for job shop schedulng with proba-bilistic durations.
Journal of Artificial Intelligence Research , , 183–232.Berman, O., & Larson, R. (2004). A queueing control model for retail services having backroom operations and cross-trained workers. Computers and Operations Research , (2), 201–222.Berman, O., Wang, J., & Sapna, K. P. (2005). Optimal management of cross-trained workersin services with negligible switching costs. European Journal of Operational Research , (2), 349–369. erekhov & Beck Brown, K. N., & Miguel, I. (2006). Uncertainty and change. In Rossi, F., van Beek, P.,& Walsh, T. (Eds.),
Handbook of Constraint Programming , chap. 21, pp. 731–760.Elsevier.Caseau, Y., & Laburthe, F. (1996). Cumulative scheduling with task intervals. In
Proceed-ings of the Joint International Conference and Symposium on Logic Programming ,pp. 363–377. MIT Press.Cezik, M. T., & L’Ecuyer, P. (2008). Staffing multiskill call centers via linear programmingand simulation.
Management Science , (2), 310–323.Demassey, S., Artigues, C., & Michelon, P. (2005). Constraint-propagation-based cuttingplanes: An application to the resource-constrained project scheduling problem. IN-FORMS Journal on Computing , (1), 52–65.Fox, M. S. (1983). Constraint-Directed Search: A Case Study of Job-Shop Scheduling . Ph.D.thesis, Carnegie Mellon University, Intelligent Systems Laboratory, The Robotics In-stitute, Pittsburgh, PA. CMU-RI-TR-85-7.Gross, D., & Harris, C. (1998).
Fundamentals of Queueing Theory . John Wiley & Sons,Inc.Hnich, B., Smith, B. M., & Walsh, T. (2004). Dual modelling of permutation and injectionproblems.
Journal of Artificial Intelligence Research , , 357–391.Martin, P., & Shmoys, D. B. (1996). A new approach to computing optimal schedules forthe job shop scheduling problem. In Proceedings of the Fifth Conference on IntegerProgramming and Combinatorial Optimization , pp. 389–403.Palmer, J., & Mitrani, I. (2004). Optimal server allocation in reconfigurable clusters withmultiple job types. In
Proceedings of the International Conference on ComputationalScience and Its Applications (ICCSA’04) , pp. 76–86.Policella, N., Smith, S. F., Cesta, A., & Oddi, A. (2004). Generating robust schedulesthrough temporal flexibility. In
Proceedings of the Fourteenth International Conferenceon Automated Planning and Scheduling (ICAPS’04) , pp. 209–218.Simon, H. A. (1997).
Models of Bounded Rationality
Handbookof Constraint Programming , chap. 11, pp. 377–406. Elsevier.Solver (2006).
ILOG Scheduler 6.2 User’s Manual and Reference Manual . ILOG, S.A.Sutton, A. M., Howe, A. E., & Whitley, L. D. (2007). Using adaptive priority weighting todirect search in probabilistic scheduling. In
Proceedings of the Seventeenth Interna-tional Conference on Automated Planning and Scheduling , pp. 320–327.Tadj, L., & Choudhury, G. (2005). Optimal design and control of queues.
TOP , (2),359–412. CP Approach for a Queueing Control Problem
Tarim, S. A., Manandhar, S., & Walsh, T. (2006). Stochastic constraint programming: Ascenario-based approach.
Constraints , (1), 53–80.Tarim, S. A., & Miguel, I. (2005). A hybrid Benders’ decomposition method for solvingstochastic constraint programs with linear recourse.. In Joint ERCIM/CoLogNETInternational Workshop on Constraint Solving and Constraint Logic Programming ,pp. 133–148.Terekhov, D. (2007). Solving queueing design and control problems with constraint pro-gramming. Master’s thesis, Department of Mechanical and Industrial Engineering,University of Toronto.Terekhov, D., & Beck, J. C. (2007). Solving a stochastic queueing control problem withconstraint programming. In
Proceedings of the Fourth International Conference onIntegration of AI and OR Techniques in Constraint Programming for CombinatorialOptimization Problems (CPAIOR’07) , pp. 303–317. Springer-Verlag.van Dongen, M. R. C. (2006). Beyond singleton arc consistency. In
Proceedings of theSeventeenth European Conference on Artificial Intelligence (ECAI’06) , pp. 163–167.Walsh, T. (2002). Stochastic constraint programming. In
Proceedings of the FifteenthEuropean Conference on Artificial Intelligence , pp. 111–115., pp. 111–115.