Joris Walraevens

J S H Van Leeuwaarden
email: j.s.h.v.leeuwaarden@tue.nl

Onno J Boxma
email: o.j.boxma@tue.nl.

Power series approximations for two-class generalized processor sharing systems

We develop power series approximations for a discrete-time queueing system with two parallel queues and one processor. If both queues are non-empty, a customer of queue 1 is served with probability β and a customer of queue 2 is served with probability 1β. If one of the queues is empty, a customer of the other queue is served with probability 1. We first describe the generating function U (z 1 , z 2 ) of the stationary queue lengths in terms of a functional equation, and show how to solve this using the theory of boundary value problems. Then, we propose to use the same functional equation to obtain a power series for U (z 1 , z 2 ) in β. The first coefficient of this power series corresponds to the priority case β = 0, which allows for an explicit solution. All higher coefficients are expressed in terms of the priority case. Accurate approximations for the mean stationary queue lengths are obtained from combining truncated power series and Padé approximation.

Introduction

Consider a discrete-time queueing model with two parallel queues that share a single processor. If both queues are non-empty at the beginning of a slot, a customer of queue 1 is served with probability (w.p.) β and a customer of queue 2 is served w.p. 1β. If one of the queues is empty, a customer of the other queue is served w.p. 1. This type of processor sharing occurs naturally in systems where different types of customers compete for resources. In telecommunication systems with integrated services, for instance, delay-sensitive streaming traffic shares resources with elastic traffic. So, if we consider the traffic arriving at queue 2 to be the delay-sensitive traffic in our model, a small β is necessary to limit the delay of this type of traffic. The exact choice of β should depend on the requirements (in terms of delay, loss, throughput, etc.) of both types of traffic.

The number of customers arriving at queue j (j = 1, 2) during slot k is denoted by a j,k . We assume that {a j,k , k > 0} forms a sequence of independent and identically distributed (i.i.d.) random variables. We denote the bivariate probability generating function (pgf) of a 1,k and a 2,k by A(z 1 , z 2 ) := E[z

a 1,k 1 z a 2,k 2 ]
. The mean number of arrivals in queue j is denoted by λ j . The customers from queue j need service for a geometrically distributed number of slots with mean 1/μ j . Since the service policy is work conserving, the stability condition is naturally given by ρ = λ 1 /μ 1 + λ 2 /μ 2 < 1.

The above queueing system gives rise to a random walk on the two-dimensional lattice in the quarter plane. The bivariate pgf of the stationary queue length distribution, denoted by U (z 1 , z 2 ), can be described in terms of a functional equation of the type K(z 1 , z 2 )U (z 1 , z 2 ) = K 00 (z 1 , z 2 )U (0, 0) + K 10 (z 1 , z 2 )U (z 1 , 0) + K 01 (z 1 , z 2 )U (0, z 2 ),

where K, K 00 , K 01 and K 10 depend on the input functions and input parameters; cf. ( 3)- [START_REF] Blanc | On a numerical method for calculating state probabilities for queueing systems with more than one waiting line[END_REF]. Using certain zero-tuples of the kernel K(z 1 , z 2 ), one can determine the functions U (z 1 , 0) and U (0, z 2 ) as solutions to a Riemann boundary value problem. This approach, developed in (10; 12) and surveyed in [START_REF] Cohen | Boundary value problems in queueing theory[END_REF], is for the model at hand outlined in Appendix A. The obtained formal solution, however, requires considerable numerical efforts, including the numerical determination of a conformal mapping. Alternative approaches exist. The most general method for obtaining the stationary distribution of nearestneighbor random walks in the quarter plane is developed by Fayolle, Iasnogorodski and Malyshev since the early seventies and summarized in the seminal book [START_REF] Fayolle | Random Walks in the Quarter Plane: Algebraic Methods, Boundary Value Problems and Applications[END_REF]. Their method also starts from [START_REF] Adan | Queueing models with multiple waiting lines[END_REF], but then builds on the analytic continuation of the functions U (z 1 , 0) and U (0, z 2 ). The rather ingenious continuation method, again combined with boundary value problems, then leads to solutions for U (z 1 , z 2 ) that are valid in the entire complex plane. Also, the conformal mapping that is required for our approach in Appendix A can in many cases be replaced by the use of Weierstrass elliptic functions, or by a Fredholm integral equation (as for Examples 1 and 2 in Subsection 5.1); for more details see [START_REF] Fayolle | Random Walks in the Quarter Plane: Algebraic Methods, Boundary Value Problems and Applications[END_REF]. We choose not to pursue the approach in [START_REF] Fayolle | Random Walks in the Quarter Plane: Algebraic Methods, Boundary Value Problems and Applications[END_REF], but rather to resort to a more numerically-oriented approach.

Other approaches for analyzing two-dimensional queueing models include the uniformization technique [START_REF] Kingman | Two similar queues in parallel[END_REF], the compensation method (3), and the power series approximation (PSA), see for instance (6; 7; 16). For a comparison of the approaches see [START_REF] Adan | Queueing models with multiple waiting lines[END_REF]. PSA is based on power series expansions of steady-state probabilities as functions of a certain parameter of the system, usually the load ρ, and was introduced in [START_REF] Hooghiemstra | Power series for stationary distributions of coupled processor models[END_REF]. By using the balance equations of the queueing system, the coefficients of the terms in the power series can be calculated iteratively. A disadvantage of this approach is the deterioration of the accuracy when ρ increases. We propose a novel version of PSA that differs from the conventional approach in two ways. Firstly, we construct a power series expansion for the bivariate pgf U (z 1 , z 2 ) directly from the functional equation [START_REF] Adan | Queueing models with multiple waiting lines[END_REF]. Secondly, we construct a power series in β rather than in ρ. This makes sense, since we are primarily interested in the results of our model for small values of β (and for all possible values of the load). Note further that the queueing system is symmetric in β in the sense that β = 0 means priority for queue 2, and β = 1 means priority for queue 1. This symmetry is not present for the parameter ρ. Therefore, our PSA approach leads to the most accurate approximations not only near β = 0 but also near β = 1 (by constructing the power series in 1β). This symmetry furthermore helps us in the construction of good approximations for all β.

Our PSA approach can be summarized as follows. For U (z 1 , z 2 ; β) := U (z 1 , z 2 ) we construct the power series

U (z 1 , z 2 ; β) = ∞ m=0 V m (z 1 , z 2 )β m , (2)
and we outline a procedure to determine the functions V m iteratively. The first term V 0 of this power series corresponds to the priority case β = 0, which is well studied and allows for an explicit solution, cf. [START_REF] Gakhov | Boundary Value Problems[END_REF]. The second term V 1 provides a first-order correction to the priority case for small β. All higher terms V m can be expressed in V 0 . A final remark concerns the chosen modeling of the service times. Although deterministic service times of exactly one slot come natural for discrete-time queueing systems, we have opted to extend it to geometrically distributed service times. This does not complicate the analysis significantly, while it allows us to derive results for the well-known continuous-time generalized processor sharing queueing system (see [START_REF] Fayolle | Two coupled processors: the reduction to a Riemann-Hilbert problem[END_REF]) directly from the discrete-time results. This is accomplished by letting the slot length go to zero and by scaling the arrival and service processes.

The paper is outlined as follows. In Section 2, we construct the functional equation for U (z 1 , z 2 ). An expression for U (z 1 , z 2 ) in terms of the solution of a boundary value problem is presented in Appendix A. In Section 3 we present, as our main contribution, the PSA approach for iteratively solving the functional equation. Approximations obtained from the PSA for the mean queue length are discussed in Section 4, along with some numerical validations in Section 5. In Section 6, we show how our discrete-time framework leads to results for the continuous-time counterpart. Some conclusions are presented in Section 7.

The functional equation

The length of queue j at the beginning of slot k is denoted by u j,k , j = 1, 2. We assume that the customer in service belongs to the queue it arrived in. We have the following system equations relating (u 1,k+1 , u 2,k+1 ) to (u 1,k , u 2,k ):

• If u 1,k = 0, u 2,k = 0: u j,k+1 = a j,k , j = 1, 2. • If u 1,k = 0, u 2,k > 0: u 1,k+1 = a 1,k and u 2,k+1 = u 2,k -1 + a 2,k , w.p. μ 2 , u 2,k + a 2,k w.p. 1 -μ 2 . • If u 1,k > 0, u 2,k = 0: u 2,k+1 = a 2,k and
u 1,k+1 = u 1,k -1 + a 1,k , w.p. μ 1 , u 1,k + a 1,k w.p. 1 -μ 1 . • If u 1,k > 0, u 2,k > 0: (u 1,k+1 , u 2,k+1 ) = ⎧ ⎪ ⎨ ⎪ ⎩ (u 1,k -1 + a 1,k , u 2,k + a 2,k ), w.p. βμ 1 , (u 1,k + a 1,k , u 2,k -1 + a 2,k ), w.p. (1 -β)μ 2 , (u 1,k + a 1,k , u 2,k + a 2,k ), w.p. β(1 -μ 1 ) + (1 -β)(1 -μ 2 ).
We define

U k (z 1 , z 2 ) as E z u 1,k 1 z u 2,k 2
, for all k, and we translate the system equations into pgfs:

U k+1 (z 1 , z 2 ) =A(z 1 , z 2 ) U k (0, 0) + 1 -μ 2 + μ 2 z 2 (U k (0, z 2 ) -U k (0, 0)) + 1 -μ 1 + μ 1 z 1 × (U k (z 1 , 0) -U k (0, 0)) + β 1 -μ 1 + μ 1 z 1 + (1 -β) 1 -μ 2 + μ 2 z 2 × (U k (z 1 , z 2 ) -U k (0, z 2 ) -U k (z 1 , 0) + U k (0, 0)) .
In steady-state, U k (z 1 , z 2 ) and U k+1 (z 1 , z 2 ) can be replaced by U (z 1 , z 2 ). By letting k → ∞ in the above equation, we find the functional equation ( 1) with kernel

K(z 1 , z 2 ) = z 1 z 2 -[(1 -βμ 1 -(1 -β)μ 2 )z 1 z 2 + (1 -β)μ 2 z 1 + βμ 1 z 2 ]A(z 1 , z 2 ) (3)
and

K 00 (z 1 , z 2 ) = [βμ 2 z 1 (z 2 -1) + (1 -β)μ 1 (z 1 -1)z 2 ]A(z 1 , z 2 ) (4) K 10 (z 1 , z 2 ) = (1 -β)[μ 2 z 1 (z 2 -1) -μ 1 (z 1 -1)z 2 ]A(z 1 , z 2 ) (5) K 01 (z 1 , z 2 ) = β[μ 1 (z 1 -1)z 2 -μ 2 z 1 (z 2 -1)]A(z 1 , z 2 ). (6)
The functional equation (1) relates U (z 1 , z 2 ) to U (z 1 , 0), U (0, z 2 ) and U (0, 0) and can be solved using the theory of boundary value problems. This is done in Appendix A for a generalization of (1) that includes the starting position and transient behavior.

The power series approximation

We introduce the notation U (z 1 , z 2 ; β) := U (z 1 , z 2 ) to express that this bivariate pgf is a function of β. First, we rearrange [START_REF] Adan | Queueing models with multiple waiting lines[END_REF] as

G(z 1 , z 2 )U (z 1 , z 2 ; β) -G 10 (z 1 , z 2 )U (z 1 , 0; β) -G 00 (z 1 , z 2 )U (0, 0; β) (7) =β • G 10 (z 1 , z 2 )[U (z 1 , z 2 ; β) -U (0, z 2 ; β) -U (z 1 , 0; β) + U (0, 0; β)],
where

G(z 1 , z 2 ) =z 2 -A(z 1 , z 2 )(μ 2 + (1 -μ 2 )z 2 ), (8)
G 10 (z 1 , z 2 ) =A(z 1 , z 2 )(μ 2 (z 2 -1) -μ 1 (1 -z -1 1 )z 2 ), G 00 (z 1 , z 2 ) =A(z 1 , z 2 )μ 1 (1 -z -1 1 )z 2 .
Note that one of the difficulties in solving the functional equation ( 7) is that it comprises both boundary functions U (0, z 2 ; β) and U (z 1 , 0; β). Our approach is based on the observation that only one of the two boundary functions appears at the left-hand side of [START_REF] Blanc | A numerical study of a coupled processor model[END_REF].

We assume, for the moment, that U (z 1 , z 2 ; β) is an analytic function of β in a neighborhood of 0. We will argue later on that this assumption is valid; see Subsection 3.1 below and Appendix B. We can then represent U (z 1 , z 2 ; β) by the power series expansion (2) for all z 1 and z 2 in the unit disk. Substitution of (2) into [START_REF] Blanc | A numerical study of a coupled processor model[END_REF] yields

G(z 1 , z 2 ) ∞ m=0 V m (z 1 , z 2 )β m -G 10 (z 1 , z 2 ) ∞ m=0 V m (z 1 , 0)β m -G 00 (z 1 , z 2 ) ∞ m=0 V m (0, 0)β m =G 10 (z 1 , z 2 ) ∞ m=0 [V m (z 1 , z 2 ) -V m (0, z 2 ) -V m (z 1 , 0) + V m (0, 0)]β m+1 .
Equating coefficients of corresponding powers of β at both sides results in the following functional equation for V m :

G(z 1 , z 2 )V m (z 1 , z 2 ) = G 10 (z 1 , z 2 )(V m (z 1 , 0) + P m-1 (z 1 , z 2 )) + G 00 (z 1 , z 2 )V m (0, 0), (9)
for all m ≥ 0, with

P m (z 1 , z 2 ) := V m (z 1 , z 2 ) -V m (0, z 2 ) -V m (z 1 , 0) + V m (0, 0),
for m ≥ 0 and P -1 (z 1 , z 2 ) := 0.

We shall now outline how to determine expressions for V m (z 1 , z 2 ). For a certain fixed m, we assume that P m-1 (z 1 , z 2 ) is known and we want to express V m in terms of P m-1 . One can prove by a generalization of Rouché's theorem (2) that G(z 1 , z 2 ) (equation ( 8)) has one zero in the unit disk of z 2 for an arbitrary z 1 in the unit disk. Denote this zero by Y (z 1 ). It is uniquely defined in the unit disk as G(z 1 , Y (z 1 )) = 0 and |Y (z 1 )| < 1. The implicit function theorem then says that Y (z 1 ) is an analytic function in the unit disk. In fact, Y (z 1 ) is the pgf of a random variable (see [START_REF] Walraevens | Delay characteristics in discrete-time GI-G-1 queues with non-preemptive priority queueing discipline[END_REF] for a similar example). Since U (z 1 , z 2 ) is analytic for all z 1 and z 2 in the unit disk, the V m (z 1 , z 2 ) are as well. Therefore, the right-hand side of (9) should equal zero for z 2 = Y (z 1 ). This gives

V m (z 1 , 0) = - G 00 (z 1 , Y (z 1 )) G 10 (z 1 , Y (z 1 )) V m (0, 0) -P m-1 (z 1 , Y (z 1 )). (10)
Upon substituting [START_REF] Cohen | Boundary Value Problems in Queueing System Analysis[END_REF] into [START_REF] Cohen | Boundary value problems in queueing theory[END_REF] we obtain

V m (z 1 , z 2 ) = 1 G(z 1 , z 2 ) G 00 (z 1 , z 2 )G 10 (z 1 , Y (z 1 )) -G 10 (z 1 , z 2 )G 00 (z 1 , Y (z 1 )) G 10 (z 1 , Y (z 1 ))
× V m (0, 0) + G 10 (z 1 , z 2 )Q m-1 (z 1 , z 2 ) (11)
with Q -1 (z 1 , z 2 ) := 0 and for m ≥ 0,

Q m (z 1 , z 2 ) :=P m (z 1 , z 2 ) -P m (z 1 , Y (z 1 )) =V m (z 1 , z 2 ) -V m (z 1 , Y (z 1 )) -V m (0, z 2 ) + V m (0, Y (z 1 )). (12)
The last step in finding an expression for V m in terms of P m-1 (or Q m-1 ) is the calculation of V m (0, 0). This constant is found from the normalization condition. Since U (1, 1; β) = 1 for all β, it follows that V 0 (1, 1) = 1 and V m (1, 1) = 0 for all m > 0. Setting z 1 = z 2 = 1 in (11) and using that Q m (1, 1) = 0 for all m ≥ 0, we find that V 0 (0, 0) = 1ρ and V m (0, 0) = 0 for m > 0. We finally arrive at the following relation between V m and Q m-1 for m > 0:

V m (z 1 , z 2 ) = A(z 1 , z 2 )[μ 2 z 1 (z 2 -1) -μ 1 (z 1 -1)z 2 ]Q m-1 (z 1 , z 2 ) z 1 [z 2 -A(z 1 , z 2 )(μ 2 + (1 -μ 2 )z 2 )] . ( 13
)
Here,

V 0 (z 1 , z 2 ) = μ 1 μ 2 (1 -ρ)A(z 1 , z 2 )(z 1 -1)(z 2 -Y (z 1 )) [z 2 -A(z 1 , z 2 )(μ 2 + (1 -μ 2 )z 2 )][μ 1 (z 1 -1)Y (z 1 ) -μ 2 z 1 (Y (z 1 ) -1)] . ( 14
)
Hence, starting from V 0 in [START_REF] Gakhov | Boundary Value Problems[END_REF], every function V m can be determined iteratively via [START_REF] Fayolle | Random Walks in the Quarter Plane: Algebraic Methods, Boundary Value Problems and Applications[END_REF] and [START_REF] Fayolle | Two coupled processors: the reduction to a Riemann-Hilbert problem[END_REF]. For β = 0, the second queue has strict priority over the first. The pgf of the queue length in this case equals U (z 1 , z 2 ; 0) = V 0 (z 1 , z 2 ) as given in [START_REF] Gakhov | Boundary Value Problems[END_REF]. This last expression is indeed the pgf in a discrete-time preemptive resume priority queueing system with geometric service times and the first queue having low priority (see e.g. ( 24)).

Analyticity of U

(z 1 , z 2 , β) in a neighborhood of β = 0
Proving the analyticity of U (z 1 , z 2 , β) in a neighborhood of β = 0 is not straightforward. There are however several (potential) approaches to tackle this problem.

A first approach is to combine the implicit function theorem with the balance equations of the Markov chain that describes the queue lengths. The implicit function theorem basically says that if a set of N equations in N + 1 variables (x 1 , , . . . , x N , y) satisfies certain conditions, and if (x

(0) 1 , . . . , x (0) N , y (0)
) is a known solution of the set of equations, then there are unique (analytic) functions u 1 (y), . . . , u N (y) such that u i (y (0) ) = x (0) i (i = 1, . . . , N) and (u 1 (y), . . . , u N (y), y) is a solution of the set of equations (see for instance [START_REF] Krantz | The Implicit Function Theorem: History, Theory and Applications[END_REF]). In our problem, x i would represent the mass functions of the queue lengths and y would be β. This approach is for instance taken in [START_REF] Hooghiemstra | Power series for stationary distributions of coupled processor models[END_REF] for a related problem. The main difficulty is that the implicit function theorem only applies to a finite set of equations. Therefore, in [START_REF] Hooghiemstra | Power series for stationary distributions of coupled processor models[END_REF], the proof consists of three steps: (i) construction of a finite Markov chain from the original infinite Markov chain, (ii) use of the implicit function theorem to prove that the stationary process of this Markov chain is analytic in the parameter and (iii) a proof that this analyticity is carried over to the (stationary process of the) infinite Markov chain. A fourth step in our case should be a proof that the analyticity of the probability mass function of the queue lengths of both queues is carried over to the bivariate pgf U (z 1 , z 2 ).

A second approach would be to reformulate (1) in terms of a Dirichlet problem for a circle, for which Schwarz's formula (see [START_REF] Brandt | On the sojourn times for many-queue head-of-the-line processor-sharing systems with permanent customers[END_REF]) yields the solution to the functional equation in terms of an elliptic integral. Such an integral would provide an explicit characterization, and for specific choices of A(z 1 , z 2 ) one could then investigate the analyticity of U (z 1 , z 2 , β). See also Remark 7 at the end of appendix B.

A third approach is the use of the implicit function theorem on the functional equation ( 1) directly. For this problem the 'traditional' implicit function theorem is of no use, since we do not have a set of equations. However, there are variants for the implicit function theorem in functional analysis (see for instance [START_REF] Dieudonné | Foundations of Modern Analysis[END_REF]).

For our problem, the proof could go as follows: the left hand side of the functional equation ( 1) is a mapping f (β, U (z 1 , z 2 )) from the product space of two Banach spaces to another Banach space. If this mapping satisfies certain conditions (most notably that the functional derivative D 2 (f )(β, U (z 1 , z 2 )) is a linear homeomorphism) and if there is a U 0 (z 1 , z 2 ) so that f (0, U 0 (z 1 , z 2 )) = 0 then there exists a unique (analytic) function W (z 1 , z 2 , β) so that f (β, W (z 1 , z 2 , β)) = 0. In our case, we have a solution U 0 (z 1 , z 2 ) = V 0 (z 1 , z 2 ), see expression [START_REF] Gakhov | Boundary Value Problems[END_REF].

Yet another approach is the one discussed in the PhD thesis of W.B. van den Hout (( 22), Section 2.5). Van den Hout first shows that the power series in, say, β as obtained in the PSA do not necessarily converge for all β inside the unit circle, and may not even be analytic at β = 0. He then provides assumptions which are sufficient for such analyticity in a neighbourhood of β = 0; see in particular Theorem 2.1 on p. 41 of [START_REF] Van Den Hout | A Numerical Approach to Markov Processes[END_REF]. Translated into our setting, we would have to pose conditions on

A(z 1 , z 2 ) such that (i) one can bound V m (z 1 , z 2 ) from above by some Vm (z 1 , z 2 ) for all |z 1 | ≤ 1, |z 2 | ≤ 1 and for m > M for some M > 0, and (ii) m β m Vm (z 1 , z 2 ) converges in a neighbourhood of β = 0. The latter holds in particular if the Vm (z 1 , z 2 ) decay geometrically fast in m. The recursive expression (13) of V m (z 1 , z 2 ) in Q m-1 (z 1 , z 2 ), and hence in V m-1 (z 1 , z 2 ), is the key to proving that bound for V m (z 1 , z 2 ).
We have chosen the third approach to prove the analyticity of U in a neighborhood of β = 0. A sketch of this proof is given in Appendix B.

Conservation laws

Let us first look at the special case μ = μ 1 = μ 2 . The pgf of the total number of customers in the queue can be found from (1) by setting z = z 1 = z 2 :

U (z, z) = A(z, z)μ(z -1)(1 -ρ) z -A(z, z)(μ + (1 -μ)z) . ( 15
)
Expression [START_REF] Guillemin | Analysis of generalized processor-sharing systems with two classes of customers and exponential services[END_REF] for U (z, z) is in fact the pgf of a discrete-time M X /Geo/1 queue with the pgf of the number of arrivals in a slot equal to A(z, z) and with geometrical service times with mean 1/μ. It is clear that the total system behaves as such a queueing system when μ 1 = μ 2 , since customers are served with rate μ when the system is busy, irrespective of β. As a result, V 0 (z, z) = U (z, z) and V m (z, z) = 0 for m ≥ 1, which also follows from the PSA approach.

When μ 1 = μ 2 , the total system content is no longer independent of β. However, the total amount of unfinished work expressed in number of slots does not depend on β. Since each customer of queue j present in the system at the beginning of a slot needs a geometrically distributed number of slots service time with mean 1/μ j , the pgf of the total unfinished work is given by:

U μ 1 z 1 -(1 -μ 1 )z , μ 2 z 1 -(1 -μ 2 )z = A μ 1 z 1 -(1 -μ 1 )z , μ 2 z 1 -(1 -μ 2 )z (z -1)U (0, 0) z -A μ 1 z 1 -(1 -μ 1 )z , μ 2 z 1 -(1 -μ 2 )z
, which is again found from (1). So we have that

V 0 μ 1 z 1 -(1 -μ 1 )z , μ 2 z 1 -(1 -μ 2 )z = U μ 1 z 1 -(1 -μ 1 )z , μ 2 z 1 -(1 -μ 2 )z and V m μ 1 z 1 -(1 -μ 1 )z , μ 2 z 1 -(1 -μ 2 )z = 0, m > 0.
This is also found from expressions ( 13) and ( 14). We will make use of this property when approximating the mean number of customers in the queues.

Approximations of performance measures

We have obtained an algorithm to determine the V m (z 1 , z 2 ) for each desired m, but the actual calculation is far from straightforward. The reason for this is that l'Hôpital's rule has to be used multiple times, which leads to expressions for V m (z 1 , z 2 ) that become more complex with m. For instance, determining

Q m (z 1 , z 2 ) requires V m-1 (z 1 , Y (z 1 )
) and V m-1 (0, z 2 ). The application of l'Hôpital's rule is necessary to obtain these latter functions from [START_REF] Fayolle | Random Walks in the Quarter Plane: Algebraic Methods, Boundary Value Problems and Applications[END_REF]. This problem is even more significant when calculating performance measures such as the moments; the generating functions then have to be differentiated in 1 which again requires multiple applications of l'Hôpital's rule.

We restrict the remaining discussion to the mean queue length in queue j given by

E[u j ] = ∞ m=0 β m ∂V m (z 1 , z 2 ) ∂z j z 1 =z 2 =1 . ( 16
) Since V m (μ 1 z/(1 -(1 -μ 1 )z), μ 2 z/(1 -(1 -μ 2 )z)) = 0 for m > 0, it is easily seen that ∂V m (z 1 , z 2 ) ∂z 1 z 1 =z 2 =1 = - μ 1 μ 2 ∂V m (z 1 , z 2 ) ∂z 2 z 1 =z 2 =1
, for m ≥ 1. Therefore, for m ≥ 1, we only need to calculate one of the two derivatives in the previous formula. The partial derivative ∂V 0 (z 1 , z 2 )/∂z j | z 1 =z 2 =1 , j = 1, 2 can be calculated easily from [START_REF] Gakhov | Boundary Value Problems[END_REF]. We get the results for the mean queue lengths in the low-priority and high-priority queue in the discrete-time M X /Geo/1 preemptive resume priority queue, as discussed before.

Let us now assume to have found the exact values v j,m for ∂V m (z 1 , z 2 )/∂z j | z 1 =z 2 =1 for j = 1, 2 and m = 0, . . . , M. Truncation of the power series [START_REF] Hooghiemstra | Power series for stationary distributions of coupled processor models[END_REF] leads to

E[u j ] = M m=0 v j,m β m + O(β M +1 ). ( 17
)
This truncation yields accurate approximations for small β.

The problem is symmetric in β in the sense that the PSA can also be constructed in β = 1β (instead of in β). If the v j,m are calculated for general A(z 1 , z 2 ), μ 1 and μ 2 , this second approximation can also be calculated directly by interchanging the roles of both queues. So, if one is interested in E[u j ] for β near 1, we can use that

E[u j ] = M m=0 ṽj,m (1 -β) m + O((1 -β) M +1 ), (18)
with ṽj,m = ∂ Ṽm (z 1 , z 2 )/∂z 3-j | z 1 =z 2 =1
, Ṽm given by V m (in ( 13)) with A(z 1 , z 2 ) replaced by A(z 2 , z 1 ), and μ 1 and μ 2 interchanged. Truncation yields approximations which are accurate near 0 or 1. In fact they provide the exact 0-to M -th order derivatives in 0 or 1. Padé approximants replace the power series ( 16) by a rational functional. We can hence approximate E[u j ] by

[L/N ] E[u j ] (β) = L l=0 u j,l β l N n=0 w j,n β n , ( 19
)
with L, N and the coefficients u j,l and w j,n chosen such that

[L/N ] E[u j ] (β) = M m=0 v j,m β m + O(β M +1 ), [L/N ] E[u j ] (β) = M m=0 ṽj,m (1 -β) m + O((1 -β) M +1 ),
i.e., such that the approximant has the correct derivatives up to order M in both 0 and 1. These 2(M + 1) derivatives to be matched by expression [START_REF] Krantz | The Implicit Function Theorem: History, Theory and Applications[END_REF] lead to a set of 2(M + 1) equations with L + N + 1 unknowns (coefficients u j,l (l = 0, . . . , L) and w j,n (n = 1, . . . , N); we use the normalization w j,0 = 1). Every choice of (L, N ) with L + N = 2M + 1 thus leads in general to a unique solution for the u j,l (l = 0, . . . , L) and w j,n (n = 1, . . . , N) in terms of v j,m and ṽj,m , m = 0, . . . , M. If one is interested in approximate formulas which are accurate for the whole range [0, 1] of β, the generalized Padé approximants are the best choice. Note that the [2M + 1/0] E[u j ] approximant is a polynomial, like the truncated version of the PSAs. We have observed that the [0/2M + 1] E[u j ] approximant is usually among the most accurate ones (see also Section 5). Note though that these Padé approximants should be used carefully, since the denominators can introduce poles in the approximation.

Remark 1. (First-order correction) The second term V 1 in the PSA equals the difference between the model described in this paper and the preemptive resume priority queue, for β going to 0. We find the following limits:

lim β→0 U (z 1 , z 2 ; β) -U (z 1 , z 2 ; 0) β = A(z 1 , z 2 )[μ 2 z 1 (z 2 -1) -μ 1 (z 1 -1)z 2 ] z 1 [z 2 -A(z 1 , z 2 )(μ 2 + (1 -μ 2 )z 2 )] × [U (z 1 , z 2 ; 0) -U (z 1 , Y (z 1 ); 0) -U (0, z 2 ; 0) + U (0, Y (z 1 ); 0)] lim β→0 E[u 2 ] -E[u 2 |β = 0] β = E[u 2 |β = 0] -E[u 2 1 u 1 =0 |β = 0] 1 -ρ 2 ( 20
)
lim β→0 E[u 1 ] -E[u 1 |β = 0] β = - μ 1 μ 2 E[u 2 |β = 0] -E[u 2 1 u 1 =0 |β = 0] 1 -ρ 2 ,
with 1 X the indicator function of the event X. These limits are important as they give first-order correction terms to the priority results for a near-priority queueing system.

Numerical examples

We now compare the PSA approximations to simulation results and investigate the influence of some parameters on the mean queue lengths. Throughout this section, we consider deterministic service times of one slot. We consider a generalized processor sharing discipline as analysed in the paper. Thus, when both queues are nonempty at the beginning of a slot, a customer of queue 1 (queue 2) is served w.p. β (1β) during that slot.

If one of the queues is empty, the other queue is served. Because of the work-conservation property, we can concentrate on the mean queue length of only one queue, say queue 2.

Validation of the approximations

Example 1. Assume the number of arrivals to both queues and the total number of arrivals in a slot to be binomially distributed. More precisely, assume the following bivariate pgf of the number of class-1 and class-2 arrivals during a slot:

A(z 1 , z 2 ) = 1 + λ 1 2 (z 1 -1) + λ 2 2 (z 2 -1) 2 . ( 21
)
Figure 1 depicts the approximations ( 17) and ( 18) as a function of β for increasing M . The arrival rates λ 1 and λ 2 are 0.7 and 0.1. We have simulated the system for β = 0, 0.1, . . . , 1 (crosses in the figure). The horizontal lines (M = 0) equal the values for the priority queues (β = 0 and β = 1 respectively). Figure 1 confirms that the PSA approximations are indeed accurate for β near 0 and near 1 and that more terms provide larger regions for β where the accuracy is good. However, as can also be seen from this figure, the approximations deteriorate (rapidly) for β away from 0 and 1. This can be dealt with using the Padé approximants introduced in Section 4. Figure 2 depicts the Padé approximants [START_REF] Krantz | The Implicit Function Theorem: History, Theory and Applications[END_REF] for the same example as in Figure 1 and for M = 3. It displays the approximants for N = 0, 2, 5 and 7 (and L = 2M + 1 -N ). Example 2. Assume the arrivals of both classes to be a sequence of two independent geometrically distributed random variables with means λ 1 and λ 2 respectively, i.e.

A(z 1 , z 2 ) = 1 -λ 1 1 -λ 1 z 1 1 -λ 2 1 -λ 2 z 2 .
We again concentrate on queue 2. Figure 3 depicts the mean class-2 content as a function of β for λ 1 = 0.7 and λ 2 = 0.1. The same conclusions as for Figure 1 can be drawn, albeit that it is more pronounced that the PSA for β works best for small β and the one for β = 1β for β near 1. We can again construct similar Padé approximants. 3.9567e-2 1.3237e-1 2.4944e-1 3.9111e-1 4.3844e-1 0.8 5.7460e-2 2.0558e-1 4.2880e-1 8.0254e-1 1.3614e+0 0.9 8.1576e-2 3.1367e-1 7.3150e-1 1.6875e+0 5.1812e+0

Table 1: First-order corrections for near-priority queues

Influence of input parameters on the mean system content for near-priority systems

Again consider the example of the binomially distributed arrival batch sizes (Expression ( 21)) and service times of one slot. The total arrival rate is given by λ = λ 1 + λ 2 and we define the fraction of class-2 arrivals α := λ 2 /λ. Table 1 displays some values of the first-order correction term (the first derivative of the mean class-2 queue content in β = 0, see Formula ( 20)) for particular values of λ and α. Some interesting observations can be made: Firstly, the correction term increases with increasing load for a constant value of α. This is expected since for high loads the number of customers in the class-1 buffer is usually non-zero, thus sharing the server with this buffer, even for a small percentage of the time can have a large influence. For small arrival rates, one of the buffers is almost always empty and the service discipline does not play an important role. A second observation is that for a given arrival rate the correction term first increases with α, reaches a maximum and then decreases. In α = 0 and α = 1, the correction term equals 0. Similar conclusions can be drawn from Table 2, where we show the relative first-order correction

∂V 1 (z 1 , z 2 ) ∂z 1 ∂V 0 (z 1 , z 2 ) ∂z 1 z 1 =z 2 =1
. This is a measure for the relative effect of giving a small share of the processor's time to the class-1 queue in the priority system. Here, the relative first-order correction does not go to 0 for α → 0 (in fact it tends to λ/(2λ) for this arrival process). The reason is that the mean class-2 queue length in the priority case equals 0 for α = 0. For low loads, the value of the relative first-order correction term is maximal for α = 0 and is a strictly decreasing function in α. For high loads, a maximum is reached for some α > 0.

QUES9188_source
We can conclude that an increase of the total arrival rate not only results in an absolute increase of the firstorder correction term, but also in an increase relative to the priority result. A GPS schedule has more impact when the arrival rate is high. Our results are therefore especially useful for these high arrival rates as they give a significant first-order (and even higher-order) correction term(s) to the priority result. The examples also show that these correction terms are sensitive to the parameters of the system, which in turn shows the necessity of the obtained formulas.

Continuous-time results

We now sketch how the discrete-time results lead to results for the continuous-time generalized processor sharing system. This continuous-time model is the most prominent subclass of the model in [START_REF] Fayolle | Two coupled processors: the reduction to a Riemann-Hilbert problem[END_REF], there called coupled processors; see also [START_REF] Konheim | Processor-sharing of two parallel lines[END_REF] for an analysis of the symmetric coupled processor model and see [START_REF] Osogami | Analysis of cycle stealing with switching times and thresholds[END_REF] for an approximate analysis of "cycle stealing" in coupled processors.

Assume that arrivals occur to both queues according to independent Poisson processes with arrival rates λ * 1 and λ * 2 respectively. The service times of customers in queue j are exponentially distributed with mean 1/μ * j . The processor is shared in the following way: when both queues are non-empty, queue 1 is served with rate β (0 ≤ β ≤ 1 without loss of generality) while queue 2 is served with rate 1β. When one of both queues is empty, the other queue is served with rate 1.

We outline in this section that we can find the distribution of the number of customers at a random point in time from the results of the discrete-time case. We therefore show that the pgfs of the numbers of arrivals and the service times as well as the functional equation ( 7) go to their continuous-time counterparts. Similar approaches were taken in (4; 21).

Let us divide the time axis into equal intervals of length Δ. We first define arrival and service processes in the discrete-time case that scale to the above continuous-time arrival and service processes. Define therefore

μ j =μ * j Δ, (22)
A(z 1 , z 2 ) =(1 -λ * 1 Δ + λ * 1 Δz 1 )(1 -λ * 2 Δ + λ * 2 Δz 2 ). (23)
The service times of class-j expressed in number of slots are geometrically distributed with parameter μ j , so their pgf is given by

S j (z) = μ * j Δ 1 -(1 -μ * j Δ)z . ( 24
)
The Laplace-Stieltjes transform of the service times in continuous time for the slot lengths going to zero equals

S * j (s) = lim Δ→0 S j (e -sΔ ) = μ * j s + μ * j . ( 25
)
The interarrival times between non-empty batches in the discrete-time model are geometrically distributed with parameter A(0, 0), i.e., their pgf is given by

I(z) = (1 -A(0, 0))z 1 -A(0, 0)z , with A(0, 0) = 1 -(λ * 1 + λ * 2 )Δ + O(Δ 2 ).
Similarly as for the service times, we get that the Laplace-Stieltjes transform of the interarrival times in continuous time for the slot length going to zero is exponential with mean 1/(λ * 1 + λ * 2 ). The pgf of the number of customers arriving in a batch is then given by

A(z 1 , z 2 ) -A(0, 0) 1 -A(0, 0) = (λ * 1 z 1 + λ * 2 z 2 )Δ + O(Δ 2 ) (λ * 1 + λ * 2 )Δ + O(Δ 2 )
.

Hence, for Δ going to zero, the interarrival times between batches are exponentially distributed with mean 1/(λ * 1 + λ * 2 ) and a batch consists of a customer arrival in the first queue with probability λ * 1 /(λ * 1 + λ * 2 ) or a customer arrival in the second queue with probability λ * 2 /(λ * 1 + λ * 2 ). We thus conclude that the arrival process into both queues converges to two independent Poisson processes with parameters λ * 1 and λ * 2 . With μ j and A(z 1 , z 2 ) as in ( 22)-( 23), the discrete-time arrival and service processes converge to Poisson arrivals and exponential service times for Δ → 0. Since β is the probability that a class-1 customer is served in a slot, this is the fraction of time class-1 customers are served when both queues are non-empty, in the limit for Δ going to zero. When we substitute μ j and A(z 1 , z 2 ) by ( 22)-( 23) in the functional equation [START_REF] Blanc | A numerical study of a coupled processor model[END_REF], we arrive at

F 1 (z 1 , z 2 )Δ + F 2 (z 1 , z 2 )Δ 2 + F 3 (z 1 , z 2 )Δ 3 = 0, (26)
with

F 1 (z 1 , z 2 ) = K * (z 1 , z 2 )U (z 1 , z 2 ) -K * 00 (z 1 , z 2 )U (0, 0) -K * 10 (z 1 , z 2 )U (z 1 , 0) -K * 01 (z 1 , z 2 )U (0, z 2 )
and

K * (z 1 , z 2 ) =(1 -β)μ 1 (z 1 -1)z 2 + βμ 2 z 1 (z 2 -1) -z 1 z 2 (λ 1 (z 1 -1) -λ 2 (z 2 -1)) K * 00 (z 1 , z 2 ) =βμ 1 (z 1 -1)z 2 + (1 -β)μ 2 z 1 (z 2 -1) K * 10 (z 1 , z 2 ) =β(μ 2 z 1 (z 2 -1) -μ 1 (z 1 -1)z 2 ) K * 01 (z 1 , z 2 ) =(1 -β)(μ 1 (z 1 -1)z 2 -μ 2 z 1 (z 2 -1)).
The precise expressions of F 2 (z 1 , z 2 ) and F 3 (z 1 , z 2 ) are not important for the further discussion. It suffices to know that -like F 1 (z 1 , z 2 ) -they are linear functions in U (z 1 , z 2 ), U (z 1 , 0), U (0, z 2 ) and U (0, 0), with coefficients analytic in the whole complex plane. Therefore, the three functions F i are analytic at least in the region |z 1 | < 1 and |z 2 | < 1. For Δ → 0 the first term in (26) is dominant. The functional equation thus converges to

F 1 (z 1 , z 2 ) = 0 (27)
for Δ → 0. This is indeed the functional equation of the continuous-time GPS system described in [START_REF] Fayolle | Two coupled processors: the reduction to a Riemann-Hilbert problem[END_REF]. Since there is only one normalized solution of this functional equation inside the unit disk, the solution described in this paper for the discrete-time functional equation evolves to the solution of ( 27) for Δ → 0 and we can thus directly find approximations for the continuous-time GPS queueing system from the results in sections 3-4.

Note that we can also find results for the continuous-time GPS system with batch arrivals, by choosing a less restrictive A(z 1 , z 2 ). Our discrete-time results can furthermore show the influence of discretizing time on the performance measures. When a continuous-time model is used to approximate a discrete-time model, our results can also quantify the error that is introduced, as a function of the slot length.

Conclusions

We studied a discrete-time two-queue Generalized Processor Sharing system in this paper. When customers are present in both queues, the queues are served with probability β and 1β respectively. We developed a novel technique based on Power Series Approximations of the joint probability generating function in the parameter β. The coefficients of the power terms are iteratively calculated starting from the constant term. This constant term is the joint pgf of a priority queue (β = 0).

By truncating the power series, we find good approximations for the means of the numbers of customers in both queues, for small β and, by symmetry, for β near 1. Interpolation techniques lead to more accurate approximations. A major advantage of the technique over the standard boundary value problem solution technique (described in the appendix) is that the formulas for say the mean queue lengths are explicit in the input parameters and require no additional numerical effort. Therefore they can be important in for instance control problems, where an optimal β is to be found given the delay requirements of both types of traffic.

The developed technique is promising to deal in general with the analysis of queueing systems with some sort of coupling. Examples are queues with a Packet-based Generalized Processor Sharing scheduling and a tandem queue where the two queues share a single processor. Another interesting topic is the Generalized Processor Sharing queue with three (or more) classes. The theory of boundary value problems has not been developed for problems with more than two dimensions.

A Solution to the boundary value problem

In this appendix, we present the solution to a more general form of the functional equation (1) in terms of a Riemann-Hilbert boundary value problem. Introducing U k;x,y (z 1 , z 2 ) := U k (z 1 , z 2 ) with initial condition u 1,0 = x, u 2,0 = y and

Φ x,y (r, z 1 , z 2 ) := ∞ k=0 r k U k;x,y (z 1 , z 2 ), (28)
we have for |r| < 1,

|z 1 | ≤ 1, |z 2 | ≤ 1: [z 1 z 2 -rψ(z 1 , z 2 )]Φ x,y (r, z 1 , z 2 ) =rA(z 1 , z 2 ){[βμ 2 z 1 (z 2 -1) + (1 -β)μ 1 (z 1 -1)z 2 ]Φ x,y (r, 0, 0) + (1 -β)[μ 2 z 1 (z 2 -1) -μ 1 (z 1 -1)z 2 ]Φ x,y (r, z 1 , 0) + β[μ 1 (z 1 -1)z 2 -μ 2 z 1 (z 2 -1)]Φ x,y (r, 0, z 2 )} + z x+1 1 z y+1 2 , ( 29
) with ψ(z 1 , z 2 ) := A(z 1 , z 2 )[(1 -βμ 1 -(1 -β)μ 2 )z 1 z 2 + (1 -β)μ 2 z 1 + βμ 1 z 2 ]. (30)
It should be noticed that ψ(z 1 , z 2 ) is the generating function of a pair of non-negative, integer-valued random variables. Formula (29) has the following global form:

K(r, z 1 , z 2 )Φ x,y (r, z 1 , z 2 ) (31) =rK 00 (z 1 , z 2 )Φ x,y (r, 0, 0) + rK 10 (z 1 , z 2 )Φ x,y (r, z 1 , 0) + rK 01 (z 1 , z 2 )Φ x,y (r, 0, z 2 ) + z x+1 1 z y+1 2 ,
with the kernel K(r, z 1 , z 2 ) being defined as

K(r, z 1 , z 2 ) := z 1 z 2 -rψ(z 1 , z 2 ). (32)
Random walks on the two-dimensional lattice in the first quadrant of the plane that give rise to functional equations and kernels of the type as in (31) and (32) are discussed by Cohen [START_REF] Cohen | Boundary value problems in queueing theory[END_REF]. In such random walks, the steps to the West, South-West and South are at most one. In [START_REF] Cohen | Boundary value problems in queueing theory[END_REF], it is sketched how this class of random walks, typically arising in queueing models, can be analysed via a transformation to a two-dimensional boundary value problem of mathematical physics, like a Riemann or Riemann-Hilbert problem. In Part II of (10), a much more detailed exposition of this approach is presented, for a slightly more restricted class of two-dimensional random walks. The kernel of the random walk that is studied there contains the kernel that features in (32). In this appendix, we sketch the way in which this boundary value approach can be used to determine the generating function Φ x,y (r, z 1 , z 2 ) of the two-dimensional queueing problem that was presented in Section 1. We distinguish between four steps.

A.1 The zerotuples of the kernel

Obviously, the generating function Φ x,y (r, z 1 , z 2 ) should be finite for all zerotuples

(ẑ 1 , ẑ2 ) of K(r, z 1 , z 2 ) with |ẑ 1 | ≤ 1, |ẑ 2 | ≤ 1. Let A(r) := {(z 1 , z 2 ) : K(r, z 1 , z 2 ) = 0, |z 1 | ≤ 1, |z 2 | ≤ 1}.
Then for all (ẑ 1 , ẑ2 ) ∈ A, one must require that

rK 00 (ẑ 1 , ẑ2 )Φ x,y (r, 0, 0) + rK 10 (ẑ 1 , ẑ2 )Φ x,y (r, ẑ1 , 0) + rK 01 (ẑ 1 , ẑ2 )Φ x,y (r, 0, ẑ2 ) + ẑx+1 1 ẑy+1 2 = 0. (33)
If one can construct functions Φ x,y (r, z, 0) and Φ x,y (r, 0, z) which are regular in |z| < 1, continuous in |z| ≤ 1, and satisfy (33) for every zerotuple (ẑ 1 , ẑ2 ) ∈ A, then Φ xy (r, z 1 , z 2 ) follows from (29) and the problem is solved. As outlined in Section 1 of ( 9), one may restrict oneself to consideration of a suitable subset S 1 × S 2 of A. Indeed, if there exist curves S 1 and S 2 , with S 1 ⊂ {z 1 :

|z 1 | ≤ 1} and S 2 ⊂ {z 2 : |z 2 | ≤ 1}
, and a one-to-one map z 1 = ω(z 2 ) from S 2 to S 1 such that (ω(ẑ 2 ), ẑ2 ) is a zerotuple of K(r, z 1 , z 2 ) for all ẑ2 ∈ S 2 , then the following holds. If functions Φ x,y (r, z, 0) and Φ x,y (r, 0, z) can be constructed that are regular for |z| < 1, continuous for |z| ≤ 1, and that satisfy (33) for all zerotuples (ẑ 1 , ẑ2 ) with ẑ1 = ω(ẑ 2 ), ẑ2 ∈ S 2 , then by analytic continuation these Φ x,y (r, z, 0) and Φ x,y (r, 0, z) satisfy (33) for all zerotuples (ẑ 1 , ẑ2 ) of A (9).

A.2 A suitable set of zerotuples

Let us now consider the construction of the curves S 1 and S 2 . While there are many possible choices here, it is important to make a choice that leads to tractable numerical techniques for obtaining the analytic continuations mentioned above. In [START_REF] Cohen | Boundary Value Problems in Queueing System Analysis[END_REF], the following choice is proposed. Let s be traversing the unit circle |s| = 1. Let z 1 = g(r, s)s and z 2 = g(r, s)s -1 , with g such that it makes the kernel zero: K(r, g(r, s)s, g(r, s)s -1 ) = 0. So in our case, with the kernel given by (32), g should satisfy g 2 = rψ(gs, gs -1 ).

(34)

Rouché's theorem implies that there are exactly two zeros of this equation satisfying |g| ≤ 1. In fact, in our case one of these two zeros is zero, due to ψ(0, 0) = 0. Now take

S 1 (r) := {z 1 : z 1 = g(r, s)s, |s| = 1}, S 2 (r) := {z 2 : z 2 = g(r, s)s -1 , |s| = 1}.

A.3 Preparing the ground for a Riemann boundary value problem

We next want to construct Φ x,y (r, z, 0) and Φ x,y (r, 0, z) that are regular in |z| < 1, continuous in |z| ≤ 1, and that satisfy (33) for all (ẑ 1 , ẑ2 ) with ẑ1 = ω(ẑ 2 ), ẑ2 ∈ S 2 . To accomplish this, one may solve a Riemann boundary value problem; see, e.g., Gakhov (14) for an extensive discussion of such boundary value problems, that aim to determine functions which are regular inside respectively outside a certain contour and satisfy a particular relation on that contour -the boundary. To be able to formulate such a boundary value problem, we need one more step. One can show (cf. [START_REF] Cohen | Boundary Value Problems in Queueing System Analysis[END_REF], Part II) that for the curves S 1 (r) and S 2 (r) there exists a unique simple closed contour L(r) in the p-plane, and functions z 1 (r, p) and z 2 (r, p) such that the following holds: Here C + and C -denote the interior and exterior of a closed contour C, and univalent means that z 1 (r, p 1 ) = z 1 (r, p 2 ) for p 1 = p 2 . Hence (i) and (ii) can be reformulated as:

(i) z 1 (r, p) : L + (r) → S + 1 (
• z 1 (r, p) is a conformal mapping of L + (r) into S + 1 (r), • z 2 (r, p) is a conformal mapping of L -(r) into S + 2 (r).
The curve L(r) can be determined by solving a particular integral equation, cf. Section II.3.6 of [START_REF] Cohen | Boundary Value Problems in Queueing System Analysis[END_REF].

A.4 The Riemann boundary value problem

Since z 1 (r, p) is regular and univalent for p ∈ L + (r) and z 2 (r, p) is regular and univalent for p ∈ L -(r), Φ x,y (r, z 1 (r, p), 0) also is a regular function for p ∈ L + (r) and continuous for p ∈ L + (r) L(r), and similarly Φ x,y (r, 0, z 2 (r, p)) is a regular function for p ∈ L -(r) and continuous for p ∈ L -(r) L(r). We are now ready to formulate a standard Riemann-type boundary value problem:

Determine two functions Ω 1 (r, p) := Φ x,y (r, z 1 (r, p), 0), p ∈ L + (r) L(r), and Ω 2 (r, p) := Φ x,y (r, 0, z 2 (r, p)), p ∈ L -(r) L(r), such that Ω 1 (r, p) is regular for p ∈ L + (r) and continuous for p ∈ L + (r) L(r), Ω 2 (r, p) is regular for p ∈ L -(r) and continuous for p ∈ L -(r) L(r), Ω 1 (r, 0) = Φ x,y (r, 0, 0), lim |p|→∞ Ω 2 (r, p) = Φ x,y (r, 0, 0), and on the boundary L(r), the functions Ω 1 (r, p) and Ω 2 (r, p) satisfy the relation (with appropriate functions H(r, p) and h(r, p), cf. (31)):

Ω 1 (r, p) = H(r, p)Ω 2 (r, p) + h(r, p), p ∈ L(r). ( 35
)
The solution of this boundary value problem may be found in [START_REF] Gakhov | Boundary Value Problems[END_REF], see also Section I.2 of [START_REF] Cohen | Boundary Value Problems in Queueing System Analysis[END_REF].

As observed before, now that Φ x,y (r, z 1 , 0) and Φ x,y (r, 0, z 2 ) are found for z 1 ∈ S + 1 (r) and z 2 ∈ S + 2 (r), one obtains Φ x,y (r, z 1 , z 2 ) first for z 1 ∈ S + 1 (r) and z 2 ∈ S + 2 (r), and finally for

|z 1 | ≤ 1, |z 2 | ≤ 1 via analytic continuation.
We end this appendix with several remarks.

Remark 2.

It should be noticed that the initial conditions u 1,0 = x, u 2,0 = y occur in the function h(r, p) in (35). Further observe that, if the initial state is a set of random variables (u 1,0 , u 2,0 ) = (X, Y ), with bivariate pgf

U 0 (z 1 , z 2 ) = ∞ x=0 ∞ y=0 P (X = x, Y = y)z x 1 z y 2
, then the last term in the right-hand side of (29) is replaced by z 1 z 2 U 0 (z 1 , z 2 ). This affects h(r, p), but not the solution approach. Remark 3. The kernel K(r, z 1 , z 2 ) has exactly the same form as the kernel that is considered throughout Part II of [START_REF] Cohen | Boundary Value Problems in Queueing System Analysis[END_REF], but the behaviour of the random walk in the interior of the first quadrant, that we consider, slightly deviates in the interior from the behaviour of the random walk that is being considered in that Part II. That implies that a slightly different Riemann boundary value problem results. Remark 4. It is easily seen that ψ(0, 0) = 0; indeed, at most one customer is served per time slot, and hence random walk transitions to the South-West are not possible. In Sections II.3.10-12 of (10), kernels with this special feature are treated in detail. It turns out that now S 1 (r) and S 2 (r) are traversed twice if s traverses the unit circle |s| = 1 once. Furthermore, one has to distinguish between the cases d dz ψ(z, 0) z=0 > (<, =) d dz ψ(0, z) z=0 , as these cases lead to different positions of the point z = 0 w.r.t. the contours S 1 (r) and S 2 (r), and consequently to slightly different analytic treatments. For the sake of exposition, let us restrict ourselves to the symmetric case that the random variables a 1,k and a 2,k are exchangeable, i.e., P (a 1,k = i, a 2,k = j) = P (a 1,k = j, a 2,k = i) for all i, j = 0, 1, . . . , and that μ 1 = μ 2 = 1 and β = 1 2 . Then ψ(z 1 , z 2 ) = A(z 1 , z 2 )(z 1 + z 2 )/2 and d dz ψ(z, 0) z=0 = d dz ψ(0, z) z=0 , leading to z 1 = 0 ∈ S 1 (r) and z 2 = 0 ∈ S 2 (r). The equation determining g(r, s) now is (cf. (34)):

g -rA(gs, gs -1 ) s + s -1 2 = 0.
The contour L(r) in this symmetric case turns out to be a circle, with center at 1 2 and radius 1 2 (see Section II.3.12 of (10)). If, moreover, only P (a 1,k = a 2,k = 0) = A 0,0 , P (a 1,k = 1, a 2,k = 0) = A 1,0 , P (a 1,k = 0, a 2,k = 1) = A 0,1 and P (a 1,k = 1, a 2,k = 1) = A 1,1 are possibly non-zero, then the equation determining g(r, s) reduces to a quadratic equation:

g -r[A 0,0 + A 1,0 g(s + s -1 ) + A 1,1 g 2 ] s + s -1 2 = 0,
from which the contours S 1 (r) and S 2 (r) are easily determined.

Remark 5. If the steady-state distribution of the joint queue-length distribution in our model exists, then its generating function Φ(z 1 , z 2 ) can be determined via Φ(z 1 , z 2 ) = lim r→1 (1r)Φ x,y (r, z 1 , z 2 ). In Section II.3.9 of (10), the details of this approach are presented; in Section II.2.16, for the case of a symmetric twodimensional random walk, it is outlined how one can directly handle the steady-state case. Here one considers the kernel K(1, z 1 , z 2 ), identifying a contour L(1) and zeros (z 1 (1, p), z 2 (1, p)) of the kernel, and formulating a Riemann boundary value problem for Φ(z 1 (1, p), 0) and Φ(0, z 2 (1, p)) with boundary L(1).

Remark 6. In order to determine performance measures like the mean steady-state queue lengths, one has to evaluate quantities like d dz Φ(z, 1) z=1 . Using our approach, this requires the numerical determination of L(1), and of (the analytic continuation of) the conformal mappings from S 1 (r) and S 2 (r), respectively, to L(r). Finally, a numerical analysis is required of the singular contour integrals that specify Φ(z, 0) and Φ(0, z). We refer to Part IV of [START_REF] Cohen | Boundary Value Problems in Queueing System Analysis[END_REF] for an extensive discussion of such numerical aspects. While it turns out to be possible to obtain numerical values of such performance measures, the analysis is quite involved. We have hence decided to concentrate in this paper on an alternative approach for the steady-state analysis of the queueing model under consideration -a method that has several novel aspects, and that allows us to obtain numerical values in a more straightforward manner. Finally, let us mention that these numerical problems could be circumvented by using the approach in (13) that strongly relies on analytic continuations, and that makes use of Weierstrass elliptic functions or a Fredholm integral equation instead of a conformal mapping.

B Proof of analyticity of U in a neighborhood of β = 0

In this appendix, we argue that U (z 1 , z 2 ; β) is an analytic function in a neighborhood of β = 0. We therefore use the implicit function theorem for Banach spaces (see theorem (10.2.3), page 272 in [START_REF] Dieudonné | Foundations of Modern Analysis[END_REF]), stated first. Lemma 1 (An implicit function theorem for Banach spaces). Let E, F , G be three Banach spaces, f a ptimes continuously differentiable mapping of an open subset A of E × F into G. Let (x 0 , y 0 ) be a point of A such that f (x 0 , y 0 ) = 0 and that the partial derivative D 2 f (x 0 , y 0 ) is a linear homeomorphism of F onto G. Then there is an open neighborhood Y 0 of x 0 in E such that, for every open connected neighborhood Y of x 0 , contained in Y 0 , there is a unique continuous mapping u of Y into F such that u(x 0 ) = y 0 , (x, u(x)) ∈ A and f (x, u(x)) = 0 for any x ∈ Y . Furthermore, u is p-times continuously differentiable in Y .

The proof of analyticity of U (z 1 , z 2 ; β) in a neighborhood of β = 0 evolves as follows. First, define B 2 as the Banach space comprising all bivariate analytic bounded functions in D 2 , with D the open complex unit disk. Similarly, define B 3 as the Banach space comprising all trivariate analytic bounded functions in D 3 , that have a limit of 0 for the first two arguments going to 1.

Define E, F and G in Lemma 1 as E = C, F = B 2 and G = B 3 × C.

We define the mapping f as:

f (β, U ) =[K(z 1 , z 2 , β)U (z 1 , z 2 ) -K 00 (z 1 , z 2 , β)U (0, 0) -K 10 (z 1 , z 2 , β)U (z 1 , 0) -K 01 (z 1 , z 2 , β)U (0, z 2 ), U(1, 1) -1]

with K, K 00 , K 10 and K 01 given in ( 3)-( 6). This mapping is defined from an open subset A of E × F to G, where A includes the point (0, V 0 ) (V 0 is the pgf of the system content of both classes for the HOL priority scheduling, i.e. for β = 0).

Since K, K 00 , K 10 and K 01 are bounded analytic functions in D 3 , and since f is affine in U and β, it is easily seen that f is p-times continuously differentiable for all p. We further observe that

f (0, V 0 ) = [0, 0].
Since f is affine in U , the (Banach space) derivative d 2 f (0, V 0 ) (see [START_REF] Dieudonné | Foundations of Modern Analysis[END_REF], section 8 for more explanation on differential calculus in Banach spaces) equals d 2 f (0, V 0 )(U ) = [K(z 1 , z 2 , 0)U (z 1 , z 2 ) -K 00 (z 1 , z 2 , 0)U (0, 0) -K 10 (z 1 , z 2 , 0)U (z 1 , 0), U(1, 1)].

We now prove that this mapping is a homeomorphism in four successive steps, namely that (i) it is continuous, that the mapping is (ii) injective and (iii) surjective (and thus bijective) and that (iv) the inverse (d 2 f (0, V 0 )) -1 is a continuous mapping.

• First of all, it is clear that d 2 f (0, V 0 ) is a continuous mapping for the same reasons as that the mapping f itself is continuous.

• Assume that d 2 f (0, V 0 )(U 1 ) equals d 2 f (0, V 0 )(U 2 ) for given U 1 and U 2 . Then K(z 1 , z 2 , 0)(U 1 (z 1 , z 2 ) -U 2 (z 1 , z 2 )) -K 00 (z 1 , z 2 , 0)(U 1 (0, 0) -U 2 (0, 0)) -K 10 (z 1 , z 2 , 0)(U 1 (z 1 , 0) -U 2 (z 1 , 0)) = 0, U 1 (1, 1) -U 2 (1, 1) = 0.

Or in other words, f (0, U 1 -U 2 ) = (0, -1).

However, this functional equation has the zero-solution as a unique solution (see Theorem 3.3 in ( 5)), and as a result U 1 = U 2 and d 2 f (0, V 0 ) is injective.

• To prove that the function d 2 f (0, V 0 ) is surjective, we solve the equation d 2 f (0, V 0 )(U ) = (g, c) with g a bivariate analytic bounded function in D 2 with limit 0 for its arguments going to 1, and c a complex number. We thus solve K(z 1 , z 2 , 0)U (z 1 , z 2 ) -K 00 (z 1 , z 2 , 0)U (0, 0) -K 10 (z 1 , z 2 , 0)U (z 1 , 0) = g(z 1 , z 2 ),

U (1, 1) = c.
This functional equation has the same form as functional equation [START_REF] Cohen | Boundary value problems in queueing theory[END_REF], and the same technique can be applied, leading to the solution

U (z 1 , z 2 ) = 1 - λ 1 μ 1 - λ 2 μ 2
c -g (1) (1, 1) μ 1 -g (2) (1, 1) μ 2 .

• The U obtained in the previous bullet is (d 2 f (0, V 0 )) -1 . This mapping is easily seen to be continuous.

We conclude that the mapping U → d 2 f (0, V 0 )(U ) is a linear homeomorphism. It thus follows from Lemma 1 that U (z 1 , z 2 ; β), as defined in the paper, is p-times differentiable at β = 0, and this for all p.

Remark 7. Although the particular two-dimensional process in this appendix, and hence our Riemann-Hilbert problem, is a special case of the more general class of problems treated in [START_REF] Cohen | Boundary Value Problems in Queueing System Analysis[END_REF], there is an alternate approach to solving the functional equation (31) by exploiting the special structure of our problem. From (31) we see that

βK 10 (z 1 , z 2 ) + (1 -β)K 01 (z 1 , z 2 ) = 0, (36)
so that the two functions K 10 and K 01 are equal up to a multiplicative constant. This feature makes that the Riemann-Hilbert problem can be reduced to a Dirichlet problem for a circle, for which Schwarz's formula yields the solution to the functional equation in terms of an elliptic integral. Such integrals provide efficient algorithms for computing performance characteristics by numerical integration. This approach was taken in (8; 12; 15). We choose not the pursue this approach here, as it requires a construction based on the analytic continuations of the functions K 10 and K 01 , and a detailed study of the involved conformal mappings, which is rather involved.

Ackowledgment

The first author wishes to thank Dieter Fiems for stimulating conversations and for pointing into the direction of Banach spaces for the analiticity proof in Appendix B. We thank the referee for his/her detailed reading of the manuscript and constructive comments.

Figure 1 :

 1 Figure 1: Truncation approximations of E[u 2 ] for binomially distributed arrival batch sizes with arrival rates λ 1 = 0.7 and λ 2 = 0.1.

Figure 2 :

 2 Figure 2: Padé approximations of E[u 2 ] for binomially distributed arrival batch sizes with arrival rates λ 1 = 0.7 and λ 2 = 0.1.

Figure 3 :

 3 Figure 3: Truncation approximations of E[u 2 ] for geometrically distributed arrival batch sizes with means λ 1 = 0.7 and λ 2 = 0.1.

  r) is regular and univalent for p ∈ L + (r),(ii) z 2 (r, p) : L -(r) → S + 2 (r) is regular and univalent for p ∈ L -(r), (iii) z 1 (r, p) = ω(z 2 (r, p)) for p ∈ L(r), ω(•) being a one-to-one map from S 2 (r) onto S 1 (r).

18 ×

 18 K 00 (z 1 , z 2 , 0)K 10 (z 1 , Y (z 1 ), 0) -K 10 (z 1 , z 2 , 0)K 00 (z 1 , Y (z 1 ), 0) K 10 (z 1 , Y (z 1 ), 0)K(z 1 , z 2 , 0) + g(z 1 , z 2 ) K(z 1 , z 2 , 0) -K 10 (z 1 , z 2 , 0)g(z 1 , Y (z 1 )) K 10 (z 1 , Y (z 1 ), 0)K(z 1 , z 2 , 0) , with g (j) (1, 1) = ∂g(z 1 , z 2 )/∂z j | z 1 =z 2 =1 and Y (z 1 ) = A(z 1 , Y (z 1 ))μ 2 (1 -A(z 1 , Y (z 1 ))(1μ 2 )

Table 2 :

 2 Relative first-order corrections for near-priority queues

	.tex; 21/07/2010; 15:23

QUES9188_source.tex; 21/07/2010; 15:23 p. 2

QUES9188_source.tex; 21/07/2010; 15:23 p. 3

QUES9188_source.tex; 21/07/2010; 15:23 p. 4

QUES9188_source.tex; 21/07/2010; 15:23 p. 6

QUES9188_source.tex; 21/07/2010; 15:23 p. 8

QUES9188_source.tex; 21/07/2010; 15:23 p. 9

QUES9188_source.tex; 21/07/2010; 15:23 p. 10

QUES9188_source.tex; 21/07/2010; 15:23 p. 15

QUES9188_source.tex; 21/07/2010; 15:23 p. 16

QUES9188_source.tex; 21/07/2010; 15:23 p. 17

Postdoctoral fellow with the Fund for The research was done during a stay of the author at the EURANDOM Research institute, and was supported by a travel grant of the FWO-Vlaanderen.


Alain Bonardi
email: alain.bonardi@univ-paris8.fr

Herman Akdag
email: herman.akdag@lip6.fr

Building Fuzzy Rules in an Emotion Detector

Keywords: Assistant, Fuzzy Rules, Intents, Emotions, Detection, Art Director, Performance

In this article, we propose a methodology to automatically construct fuzzy rules in a Fuzzy Rule-Based System that detects emotions from an actor's performance during a show.

The idea is to collect video data from a lot of performances of the same show from which it should be possible to visualize all the emotions and intents or more precisely "intent graphs". To do this, we analyze the collected data defining low-level descriptors that are converted into high-level characterizations thanks to aggregations. The following step is to partition correctly the universes in classes, depending on the retrieved data. Finally the fuzzy rules built from the classes are used to label the detected emotions.

Introduction

Entertainment productions are more and more numerous and directing is often a complex task, notably with the productions containing a part of improvisation. However the performers of these productions need to be directed since the improvisation often follows some rules.

In this context, it is necessary to conceive tools that could help the art director in his task of actors' performance supervision. The tool we propose here is a kind of assistant that gives a representation of a set of complex data for the exploration and observation of a show, and it is also a way of better understanding the creation and execution of the show. One important point is that, as every computer program, the assistant must be deterministic and systematic, i.e. it should always give the same results for a given entry and it must look over the whole data. Like high-level sportsmen that have tools to analyze, correct and improve their gesture, we want to propose a tool to assist the creative artists, to let them better understand the phases of a performance and to help them in their creation process.

We have chosen to work on a virtual opera called Alma Sola written by Bonardi and Zeppenfeld where a performer plays (sings and dances) different blocks from different universes (such as Prologue, Love, Pleasure, etc.). Alma Sola is an open form of opera [START_REF] Bonardi | New Approaches of Theatre and Opera Directly Inspired by Interactive Data-Mining[END_REF]. The performer embodies a feminine Faust and wanders through the various universes split into blocks. She therefore interprets an opera playlist that she selects during the show itself. For instance, a performance can be : Love-3, then Wealth-5, then Pleasure-3, etc. The computer offers continuations to the performer and suggests the next block to be performed (cf. figure 1). We retrieve the actor's performance on video files and look for emotions in the data. A distinction must be made between intent and emotion: indeed intent corresponds to the conscious part of the emotion. Thus the assistant enables the comparison between the performer's intent and the rendered emotions.

It is widely known that fuzzy logic offers good tools to deal with such subjective concepts [START_REF] Zadeh | Toward a generalized theory of uncertainty (GTU) -an outline[END_REF], emotions here, this is why we shall use a fuzzy rule-based system (FRBS) to detect the performer's emotions.

The article is organized as follows: first we give an overview of the existing research about assistants in art, then we describe our assistant and notably the way we retrieve emotion descriptors. In the third section we show how the descriptors are partitioned in classes. This step is necessary to build correctly the fuzzy rules that form the inference system. Finally, section 4 concludes this study.

Research about performer's assistants

In a way, research about performer's assistants has existed for centuries. One immediately thinks of the use of mirrors in dance.

At the time when mirrors become common place in Europe (Renaissance), the first treatise about dance is released by Thoinot Ar-beau in 1589. Dance has given the example for a long time among performing arts. Many notation systems (Feuillet, Laban, Benesh, etc.) have been developed along centuries. All these systems, as the musical notation, have in common to give through a score the causes of the phenomenon (and not its result) to try to access the author's intentions. For instance, scores encode descriptions either implicit (you have to play an F, but the score does not tell you how to do it) or explicit gestures to be achieved (cross hands when playing the piano, for instance), and from these structured indications, a dancer or a musician tries to infer some of the author's intentions [START_REF] T R I P A T H I&s | Automated Gesture Segmentation From Dance Sequences[END_REF].

Another approach consists in dealing with another set of intentions, the ones expressed by the performers such as we experience them, without direct reference to the author's. Hopefully these two systems of intention have a non-empty intersection (we sometimes feel this is not the case or on the contrary we sometimes feel a perfect adequacy between the work and its performers). In the field of music, for instance, one can "measure" a kind of difference between the prescription (typically the score) and its realization (from listening to musicians). In his approach of "Virtual scores", Manoury has considered [START_REF] Manoury | Les partitions virtuelles[END_REF] in his pieces for solo instrument and live electronics the computation of this difference (Jupiter for solo flute and live electronics, 1987, En Echo for voice and live electronics, 1991).

The first step in that approach consists in measuring various aspects of the performance, using captors in a broad meaning [START_REF] Zeppenfeld | L'acteur face aux technologies[END_REF]. Various captors are nowadays available: video camera, wireless microphone, ultrasound device, carpet detectors, digital compass, etc. Recently, Camurri and his team have developed the first robust platform for the analysis of gestures and consequently of performer's emotions. It is named EyesWeb [START_REF] C A M U R R I | EyesWeb -toward gesture and affect recognition in dance/music interactive systems[END_REF], [START_REF] Lindström | Affect, attitude and evaluation of multi-sensory performances[END_REF]. It is not based on captors implemented on the performer's body (with heavy batteries and radio transmission), but on video capture with a static shot. It is based on a graphical language that implements many descriptors of gestures: quantity of motion, stability, etc. EyesWeb has become a worldwide standard for performance analysis.

Roughly ten years before, the Ircam institute (and Cycling 74 company) had developed the real-time digital sound analysis and synthesis platform to complete Max software. It is named MSP (Methods for Sound Processing) and is now included in Max software (Max/MSP). This software is a worldwide standard in real-time sound analysis.

At the present time the state-of-the-art consists mainly in inferring a few emotional states from the raw data delivered by EyesWeb and/or Max/MSP. A project relatively similar to ours is developed by Friberg and his team. They have conceived a real-time algorithm to analyze emotional expression in musical performance and body movement [START_REF] Friberg | A fuzzy analyzer of emotional expression in music performance and body motion[END_REF]. In the framework of a game named "Ghost in the cave", the player has to express different emotions using his/her body or his/her voice, and these emotions are the input values of the software. They use EyesWeb to recover body movements and sound descriptors (sound level, instant tempo, articulation, attack rate, high-frequency content). But in this game, there is an immediate feedback, and the player has to move constantly and talk until the software reacts according to his/her wishes. Our aim is not the same, i.e.

the assistant adapts to the performer and not the contrary. There are also many works on facial expressions recognition with fuzzy techniques, such as the approach developed by Ralescu & Hartani, for example [START_REF] Ralescu | Fuzzy Modeling based Approach to Facial Expressions Understanding[END_REF].

Project Description

As explained above, we have chosen to work on the interactive opera Alma Sola written by A. Bonardi and C. Zeppenfeld. Two dissimilar scenes have been extracted to be used as "sample scenes": they are the chanted Prologue (which is improvised, unwritten) and the Love Universe (which is "strictly written"). The performer (a singer/dancer, here) is filmed by a camera in wide and static shot (cf. figure 2, left) and his/her voice is recorded in a separate file, in order to handle both sources of data separately.

The project is centered on two main phases: the acquired scene processing (with the performer) and the sound capture processing. In this article we focus on the first phase, where EyesWeb has been used (through a dedicated patch we have written) to analyze accurately, understand and exploit non-verbal expressive gestures. Several parameters can be extracted from the video file: quantity of motion, stability, motion duration, pause duration in a scene, contraction index and surface of the performer in the image, convex hull of the body silhouette, velocity, acceleration, etc.

This patch construction step is not trivial and represents hours of tests with sample videos and it implies a concertation with the art directors. First of all, the background of the image must be removed -using the difference between two frames -in order to keep only the performer's movement, since the camera is static. Next, the parameters can be easily extracted because we are sure they only concern the performer, and not the background.

Here is the description of some of the most interesting parameters for our problem (Table 1 shows some of them). These parameters can be seen as gesture primitives.

The quantity of motion is computed by the number of pixels changing position between two instants (the white pixels in figure 2); the convex hull of the body silhouette is the bounding rectangle of the white pixels (cf. figure 2); the stability is the ratio of the height of the silhouette's center of gravity on the length of the segment connecting the lower points of the silhouette; the contraction index is the ratio of the silhouette's surface over the surface of the convex hull. The stability is an important descriptor since it gives good insight on whether the performer is near the ground or not, i.e. whether the performer puts himself at risk or not. The contraction index is also very useful and reflects whether the performer is effusive or not.

Choosing these parameters judiciously (called video descriptors) allows us to compute various aggregations of each set of values for each descriptor (cf. a simple example in figure 3). The chosen aggregators are: partial and general means (computed on a part of the captured scene or on the whole scene), standard deviation, covariance, etc. denoted A 1 ,A 2 ...A 9 . Then, to a meta-level, we characterize and categorize the sequences of each scene thanks to an FRBS. The number of categories used in this opera is five (Sleepy, Angry, Happy, LoveBeliever (could be also called Effusive), LittleEffusive). Figure 4 sums the whole process up. (We capture 25 frames per second, i.e., 25 values per second for each gesture descriptor.)

As can be seen in Figure 4, once the descriptors are extracted from the video files, they are aggregated in order to give a description of each emotion we want to recognize. Figure 5 shows graphically the results for the Universe of Love. The y-axis represents the normalized (between 0 and 2) interval described in table 1, i.e. it is in ms, number of pixels or . . . The next step is to classify in partitions the aggregation results. For example, when the performer expresses an emotion such as happiness, he/she moves a lot. But this is not always easy to guess even if the emotions share some general patterns.

The fuzzy partitioning allowing the construction of the rules is now described.

Building the Fuzzy Rules

Fuzzy partitioning of the values taken by the aggregation results must cover the domain of discourse. However a simple uniform distribution of the classes on the axis is not appropriate since sometimes small, other times average or in other cases big variations lead to an emotion change, depending on the aggregated vectors.

Like Martínez & al. in [START_REF] Martínez | Fuzzy Tools to Deal with Heterogeneous Information in Engineering Evaluation Processes[END_REF], we propose a categorization depending on the data distribution. First of all, five classes are considered: Very Low, Low, Average, High and Very High values (denoted VL, L, A, H, VH ). For the moment, we consider the classes as simple intervals. The fuzzy sets construction will be explained below. For each aggregator and for all the records of each scene, the building of class Low is dynamically performed as follows: the interval's left bound is the minimum value of the whole data denoted x 0 .T h i s way, when analyzing a new record of the same scene, if there are smaller values than x 0 ,they will be categorized as Very Low values. In the same way, for each aggregator and for all the records of each scene, the building of class High is dynamically performed as follows: the interval's right bound is the maximum value of the whole data denoted x 2 . Then, for each aggregator and for all the records of each scene, the average value of the whole data is denoted x 1 .B o t hi n t e r v a l s[ x 0 ,x 1 ]and[x 1 ,x 2 ] are split into three equal sub-intervals, each.

Finally, the interval for class L is the concatenation of the first two sub-intervals, the inter-

x 0 x 2 x 1 (mean) (min) (max)
LA H

Figure 6: Splitting the intervals.

Emotions A1 A2 A3 A4 A5 A6 A7 A8 A9 LoveBel. L A L L A L L A L Sleepy H L H H L H H L H Angry A H L A A L A H L Happy L H L L H L L H L LittleEff. H L L H L L H L L
Table 2: Classification for the Universe of Love.

val for class A is the concatenation of the next two sub-intervals and the interval for class H is the concatenation of the last two subintervals. Figure 6 shows the intervals split. Thus it is possible to classify the aggregation values for the five emotions, as shown in table 2 where A i is the i th aggregator, in the case of the Love Universe. This will permit to establish the fuzzy rules.

The next step consists in building the fuzzy subsets representing the five classes for the FRBS. Figure 7 shows how this is performed.

The smallest class-interval, L in our example, is defined as an L-R fuzzy number and sets the construction unit for the other intervals.

If we denote by (a [START_REF] Bouchon-Meunier | La logique floue et ses applications[END_REF], the building of the five classes in the example of figure 7 is as follows:

F ,b F ,c F ,d F ) an L-R fuzzy subset F (with Supp(F )=[ a F ,d F ]a n d Ker(F )=[b F ,c F ])
VL : ⎧ ⎪ ⎨ ⎪ ⎩ a VL = b VL =inf x∈A i x c VL =2x 0 -x 0,1 d VL = x 0,1 L : ⎧ ⎪ ⎨ ⎪ ⎩ a L =2x 0 -x 0,1 = c VL b L = c L = x 0,1 = d VL d L = x 1 A : ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ a A = x 0,1 = b L b A = x 1 = d L c A = x 0,2 + x 1,1 -x 1 d A =2x 1,1 -c A H : ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ a H = c A b H = d A c H =2x 1,2 -b H d H =2x 2 -c H VH : ⎧ ⎪ ⎨ ⎪ ⎩ a VH = c H b VH = d H c VH = d VH =sup x∈A i x
This way, the smallest class-interval, either L or H, is an L-R fuzzy number. Generally speaking, at least, one of the two classes (L or H ) must be an L-R fuzzy number. Nevertheless A and H (or respectively A and L,d epending on which class contains the smallest sub-interval) may be either L-R fuzzy numbers or L-R fuzzy intervals. They are actually fuzzy numbers under a certain condition depending on a parameter that we call ε:

∀X ∈{L, A, H },w edenotebyX p the preceding class in the universe. After having computed a X ,b X ,c X ,d X such as in the example of figure 7 for class A or H,i fc X -b X <ε then c X and b X must be re-computed, i.e.

b ′ X = c ′ X =1/2(b X + c X )
The value taken by ε determines the membership value v to which the classes overlap.

In particular, if ε<c X -b X then v =0 .5a n dX is an L-R fuzzy interval, otherwise (i.e

. ε ≥ c X -b X ) v = 1 a X -b ′ X a X c X p -b ′ X d X p -a X + b ′ X -c X p + d X p + a X and
X is an L-R fuzzy number.

In the software we propose (see section 5), we also offer the possibility to use L-R fuzzy numbers exclusively (cf. figure 8) or L-R fuzzy intervals exclusively.

Finally, it is easy to establish the fuzzy rules according to Table 2, one rule per line. Here is an example for the Love Universe:

Rule 1:

If A 1 is L &A 2 is A &A 3 is L & A 4 is L &A 5 is A &A 6 is L & A 7 is A &A 8 is L &A 9 is L Then the emotion is LoveBeliever x 1 (mean) x 0 (min) x 2 (max) VL x 0 x 0,1 x 0,2 x 1 x 1,1 x 1,2 x 2 µ 1 A i LA H V H H LA

Application

The application we have developed implements the concepts explained above. We have worked on two universes: Love and Prologue. The performer has played several times each universe and the software has tried each time to detect the emotions perceived. The video files obtained for the performances are split into smaller files that are given to the software. After having chosen the fuzzy subsets that will be used for the partitioning (cf. section 4), the values of the aggregated descriptors are displayed and a graphical result is also proposed. Figure 9 shows that both Love Universe and Prologue have rather been per-formed with Sleepy emotion. results are better (i.e., the results for the emotions are more clearcut) using partitioning with both L-R fuzzy numbers and L-R fuzzy intervals because they adapt to the level of precision, depending on the values retrieved from the video file.

6C o n c l u s i o n

In this paper we have presented a system that is able to give clues to the art director in order to evaluate a performer's rendition. This is done thanks to a fuzzy rule-based system that detects the actor's emotions during a per-formance show. One originality is the way we construct the fuzzy classes when partitioning the universes before the rule construction. They are dynamically built according to the values characterizing the performance.

As a future work, we will evaluate our system with musical experts by asking them to fill in a perception grid for each scene. Moreover, we will try our assistant on other test sets, i.e. with more records from Alma Sola (with the same performer or not) but also with records from other shows. Moreover, it would be very interesting to include a back propagation in the software: when an unexpected emotion is detected, the assistant should suggest modifications of his behaviour to the performer in order to obtain best results during the next detection.

Concerning Alma Sola opera itself, we can also imagine in the future that the assistant could be used to classify the blocks performed according to the detected emotions and then contribute to the design of the open form.

Figure 1 :

 1 Figure 1: Alma Sola sequence diagram.

Figure 2 :

 2 Figure 2: Left, an image taken from the video file; right, the corresponding convex hull of the body silhouette (performer: Claire Maupetit). Descriptor Interval Intuitive description Movement duration [0, +∞[i nm s Pause duration [0, +∞[i nm s Contraction surface [0,nb of pixels in the video[ Quantity of pixels representing the performer in the image

Figure 3 :

 3 Figure 3: Some of the aggregations used.

Figure 4 :

 4 Figure 4: Our system components.

Figure 5 :

 5 Figure 5: Aggregated descriptors for each emotion in the Universe of Love.

Figure 7 :

 7 Figure 7: Building the fuzzy classes.

Figure 8 :

 8 Figure 8: The L-R fuzzy numbers used in the software.

Figure 9 :

 9 Figure 9: A screenshot of the software.

Table 1 :

 1 Descriptor definitions.

Acknowledgements

We thank Adrien Revault d'Allonnes and Murat Goksedef for their collaboration and help on this work. This project was funded by the Maison des Sciences de l'Homme Paris Nord.


Jean B Lasserre

RECOVERING AN HOMOGENEOUS POLYNOMIAL FROM MOMENTS OF ITS LEVEL SET

Keywords: homogeneous polynomials, sublevel sets, moments, inverse problem from moments

Let K := {x : g(x) ≤ 1} be the compact sub-level set of some homogeneous polynomial g. Assume that the only knowledge about K is the degree of g as well as the moments of the Lebesgue measure on K up to order 2d. Then the vector of coefficients of g is solution of a simple linear system whose associated matrix is nonsingular. In other words, the moments up to order 2d of the Lebesgue measure on K encode all information on the homogeneous polynomial g that defines K (in fact, only moments of order d and 2d are needed).

Introduction

The inverse problem of reconstructing a geometrical object K ⊂ R n from the only knowledge of moments of some measure µ whose support is K is a fundamental problem in both applied and pure mathematics with important applications in e.g. computer tomography, inverse potentials, signal processing, and statistics and probability, to cite a few. In computer tomography, for instance, the X-ray images of an object can be used to estimate the moments of the underlying mass distribution, from which one seeks to recover the shape of the object that appears on some given images. In gravimetry applications, the measurements of the gravitational field can be converted into information concerning the moments, from which one seeks to recover the shape of the source of the anomaly.

Of course, exact reconstruction of objects K ⊂ R n is in general impossible unless K has very specific properties. For instance, if K is a convex polytope then exact recovery of all its vertices has been shown to be possible via a variant of what is known as prony method. Only a rough bound on the number of vertices is required and relatively few moments suffice for exact recovery. For more details the interested reader is referred to the recent contribution of Gravin et al. [START_REF] Gravin | The Inverse Moment Problem for Convex Polytopes[END_REF] and the references therein. On the other hand, Cuyt et al. [START_REF] Cuyt | Multidimensional integral inversion, with applications in shape reconstruction[END_REF] have shown that approximate recovery of a general n-dimensional shape is possible by using an interesting property of multidimensional Padé approximants, analogous to the Fourier slice theorem for the Radon transform.

Contribution. From previous contributions and their references, it is transparent that exact recovery of an n-dimensional shape is a difficult problem that can be solved only in a few cases. And so identifying such cases is of theoretical and practical interest. The goal of this paper is to identify one such case as we show that exact recovery is possible when K ⊂ R n is the (compact) sublevel set {x ∈ R n : g(x) ≤ 1} associated with an homogeneous polynomial g. By exact recovery we mean recovery of all coefficients of the polynomial g. In fact, exact recovery is not only possible but rather straightforward as it suffices to solve a linear system with a nonsingular matrix! Moreover, only moments of order d and 2d of the Lebesgue measure on K are needed. As already mentioned, exact recovery is possible only if K has very specific properties and indeed, crucial in the proof is a property of levels sets associated with homogeneous polynomials (and in fact, also true for level sets of positively homogeneous nonnegative functions). A polynomial

Main result

f ∈ R[x] d is written x → f (x) = α∈N n f α x α , for some vector of coefficients f = (f α ) ∈ R s(d) . A real-valued polynomial g : R n → R is homogeneous of degree d (d ∈ N) if g(λx) = λ d g(x) for all λ and all x ∈ R. Given g ∈ R[x], denote by G ⊂ R n the sublevel set {x : g(x) ≤ 1}.
If g is homogeneous then G is compact only if g is nonnegative on R n (and so d is even). Indeed suppose that g(x 0 ) < 0 for some x 0 ∈ R n ; then by homogeneity, g(λx 0 ) < 0 for all λ > 0 and so G contains a half-line and cannot be compact.

Main result. The main result is based on the following result of independent interest valid for positively homogeneous functions (and not only homogeneous polynomials). A function

f : R n → R is positively homoge- neous of degree d ∈ R if f (λx) = λ d f (x)
for all λ > 0 and all x ∈ R n . Lemma 1. Let f : R n → R be a measurable, positively homogeneous and nonnegative function of degree 0 < d ∈ R, with bounded level set {x : f (x) ≤ 1}. Then for every k ∈ N and α ∈ N n :

(2.1)

{x : f (x)≤1} x α f (x) k dx = n + |α| n + kd + |α| {x : f (x)≤1}
x α dx.

Proof. To prove (2.1) we use an argument already used in Morosov and Shakirov [START_REF] Morosov | New and old results in resultant theory[END_REF][START_REF] Morosov | Introduction to integral discriminants[END_REF]. With α ∈ N n , let α := (α 2 , . . . , α n ) ∈ N n-1 and define z := (z 2 , . . . , z n ). Let φ : R + → R be measurable and consider the integral R n φ(g(x)) x α dx. Using the change of variable x 1 = t and x i = tz i for all i = 2, . . . , n, and invoking homogeneity, one obtains:

R n φ(f (x)) x α dx = R n φ(t d f (1, z 2 , . . . , z n )) t n+|α|-1 z α d(t, z) = d -1 ∞ 0 u (n+|α|)/d-1 φ(u) du × A α with A α = R n-1 z -αf (1, z) -(n+|α|)/d dz.
Hence the choices t → φ(t)

:= I [0,1] (t) and t → φ(t) := t k I [0,1] (t) yield d {x : f (x)≤1} x α dx = A α 1 0 u (n+|α|)/d-1 du = A α d n + |α| d {x : f (x)≤1} f (x) k x α dx = A α 1 0 u (n+kd+|α|)/d-1 du = A α d n + kd + |α| ,
respectively. And so (2.1) follows.

With g ∈ R[x] d being an homogeneous polynomial of degree d, consider now the matrix M d (λ) of moments of order 2d associated with the Lebesgue measure on G = {x : g(x) ≤ 1}, that is, M d (λ) is a real square matrix with rows and columns indexed by the monomials x α , α ∈ N n d , and with entries

(2.2) M d (λ)[α, β] = {x : g(x)≤1}
x α+β dx =:

λ α+β , ∀ α, β ∈ N n d .
So all entries of M d (λ) are moments of order 2d. Our main result is as follows:

Theorem 2. Let g ∈ R[x] d be homogeneous of degree d with unknown coefficient vector g ∈ R s(d) and with compact level set G = {x : g(x) ≤ 1}. Assume that one knows the moments λ = (λ α ) for the Lebesgue measure on G, for every α ∈ N n with |α| = 2d and |α| = d. Then:

(2.3) g = n + d n + 2d M d (λ) -1 λ (d)
where λ (d) = (λ α ), α ∈ N n d , is the vector of all moments of order d.

Proof. Use (2.1) with k = 1 and |α| = d to obtain

β∈N n : |β|=2d g β λ α+β = n + d n + 2d λ α , ∀ |α| = d, or in matrix form M d (λ) g = n + d n + 2d λ d ,
from which the desired result follows if M d (λ) is non singular. But this follows from the fact that G has nonempty interior.

There are alternative ways for obtaining g from the moments λ. It suffices to apply (2.1) for a family F of multi-indices α ∈ N n whose cardinal |F| matches the dimension n+d-1 d of the vector g. For instance, with n = 2 and d = 2 (g is a quadratic form with vector of coefficients g = (g 20 , g 11 , g 02 )), g can also be obtained by: Similarly, with n = 2 and d = 4 (g is a quartic form with vector of coefficients g = (g 40 , g 31 , g 22 , g 13 , g 04 )), g can also be obtained by

2. 1 .

 1 Notation and definitions. Let R[x] be the ring of polynomials in the variables x = (x 1 , . . . , x n ) and let R[x] d be the vector space of polynomials of degree at most d (whose dimension is s(d) := n+d n ). For every d ∈ N, let N n d := {α ∈ N n : |α| (= i α i ) = d}, and let v d (x) = (x α ), α ∈ N n , be the vector of monomials of the canonical basis (x α ) of R[x] d . Denote by S k the space of k × k real symmetric matrices with scalar product B, C = trace (BC); also, the notation B 0 (resp. B ≻ 0) stands for B is positive semidefinite (resp. positive definite).

λ

  20 λ 11 λ 02 λ 30 λ 21 λ 12 λ 21 λ 12 λ 30 above inverse matrix exists.

λ

  40 λ 31 λ 22 λ 13 λ 04 λ 50 λ 41 λ 32 λ 23 λ 14 λ 41 λ 32 λ 23 λ 14 λ 05 λ 60 λ 51 λ 42 λ 33 λ 24 λ 42 λ 33 λ 24 λ 15 λ 06 above inverse matrix exists.