# A multi-orthogonal polynomials' approach to bulk queueing theory

We consider a stationary Markov process that models certain queues with a bulk service of fixed number m of admitted customers. We find an integral expression of its transition probability function in terms of certain multi-orthogonal polynomials. We study the convergence of the appropriate scheme of simultaneous quadrature rules to design an algorithm for computing this integral expression.

## Authors

• 1 publication
11/21/2020

### Orthogonal polynomials on planar cubic curves

Orthogonal polynomials in two variables on cubic curves are considered, ...
01/25/2021

### There are EXACTLY 1493804444499093354916284290188948031229880469556 Ways to Derange a Standard Deck of Cards (ignoring suits) [and many other such useful facts]

In this memorial tribute to Joe Gillis, who taught us that Special Funct...
12/31/2021

### Subresultant of several univariate polynomials

Subresultant of two univariate polynomials is a fundamental object in co...
03/15/2021

### Spline quadrature and semi-classical orthogonal Jacobi polynomials

A theory of spline quadrature rules for arbitrary continuity class in a ...
11/01/2021

### On the Markov extremal problem in the L^2-norm with the classical weight functions

This paper is devoted to Markov's extremal problems of the form M_n,k=su...
04/08/2020

### PolyChaos.jl – A Julia Package for Polynomial Chaos in Systems and Control

Polynomial chaos expansion (PCE) is an increasingly popular technique fo...
01/05/2021

### Modified discrete Laguerre polynomials for efficient computation of exponentially bounded Matsubara sums

We develop a new type of orthogonal polynomial, the modified discrete La...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

In [6] N. Bailey considered a queueing model with a bulk service admitting batches with a maximum size . This model describes several situations where the customers are admitted in groups, for example transportation of people in vehicles with a maximum passenger capacity. In our case we assume that batches have a fixed size (see [24, Chapter 4]

). We have an example of these situations in some shuttle routes where their buses only start their trip if they complete the vehicle’s capacity. We can find more examples of this behavior in certain bureaucratic process, for instance this happens in some public museums where the guides only admit a fixed number of customers and if this number is not completed they don’t work. Let us assume that an individual arrives at each epoch of Poisson occurrence with rate

, and the service only starts if the queue size reaches or exceeds . The service time distribution of a batch is assumed to be exponential with parameter . In terms of Kendall’s notation [20] used in [24], we say that our system is .

This system corresponds to a stationary Markov process with path function whose transition probability function

 Pi,j(t)=Pr{X(t+s)=j|X(s)=i},

satisfies the following conditions

 (1.1) ⎧⎪⎨⎪⎩Pi,i(Δt)=1−λΔt+o(Δt),Pi,i+1(Δt)=λΔt+o(Δt),asΔt→0,ifi

When we are in the classical case of a Birth-and-Death process (see [19]). Birth and Death processes have been used to model any M/M/1 system (see [5, Chapter III]). As in [10, Chapter XVII]), considering (1.1) with the Chapman-Kolmogorov equations, we obtain the following system of infinite differential equation with :

 P′i,j(t)=−(λ+μ)Pi,j(t)+λPi,j−1(t)+μPi,j+m(t),(i,j)∈Z+×Z+=Z2+.

Let be a infinite matrix function , then is the solution of the following initial matrix value problem

 (1.2) ⎧⎪⎨⎪⎩P′(t)=P(t)A,t≥0,P(0)=I,

where is such that for each :

 ai,i+1=λ,ai,i=−(λ+μ),ai,i−m=μ,% andai,j=0,j∉{i+1,i,i−m},

schematically

 A⊤=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝−(λ+μ)000⋯0μ00⋯λ−(λ+μ)00⋯00μ0⋯0λ−(λ+μ)0⋯000μ⋯⋮⋱⋱⋱⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

denotes the transpose of matrix .

According to W. Sternberg (see [29]) the problem (1.2) has a unique solution. In [19], S. Karlin and J. McGregor arrived at an integral expression for the entries of in the case of , which is a Birth-and-Death process. They used some aspects of the orthogonal polynomial theory. Using Karlin and McGregor’s approach, T. S. Chihara and M. E. H. Ismail (see [9]) connected sequences of orthogonal polynomials with a queueing model. Later on in [17] M. E. H. Ismail revisited this connection. In this paper we also use Karlin and McGregor’s method to study our case for general , using now properties of multi-orthogonal polynomials.

Let be a sequence of polynomials generated by the recurrence relation

 (1.3) (λ+μ−x)Qn(x)=λQn+1(x)+μQn−m(x),n≥0

with initial conditions

 Q−m≡Q1−m≡⋯≡Q−1≡0andQ0=1,

which implies that

 Qn(x)=1λn(λ+μ−x)n,n=0,…,m.

We consider a family of vector polynomials

that satisfy the recurrence relation

 (1.4) q0=(1,0,0…,0),q1=(0,λ,0…,0),…,qm−1=(0,0,…,λm−1),(λ+μ−x)qr(x)=μqr+m(x)+λqr−1(x),r=1,2,….

Given two complex numbers and we denote by the section of straight line that connects with . containing the endpoints. The following Theorem is proved in Section 4:

###### Theorem 1.1.

There exists a system of weights , defined on the starlike set

 Σ0=m⋃k=0[λ+μ,λ+μ+(m+1)m+1√μλmmm/(m+1)exp(2πikm+1)],

such that:

 (1.5) −∫Qn(x)(q1,r(x)+q2,r(x)ρ2,2(x)+⋯+qm,r(x)ρ2,m(x))ρ1(x)dx=δn,r,

where

 δn,r=⎧⎪⎨⎪⎩0ifn≠r,1ifn=r,

denotes Kroneker’s delta.

In the proof of above Theorem 1.1 we follow ideas published in the papers [3, 4]. Let us state our main result:

###### Theorem 1.2.

The matrix function whose entries are

 (1.6) Pn,r(t)=−∫e−xtQn(x)(q1,r+q2,rρ2,2+⋯+qm,rρ2,m)(x)ρ1(x)dx

is the unique solution of the initial value problem (1.2).

In Section 3 we obtain expressions for the polynomials , and , and in Section 4 we give more details about the structure of the weights

. Such expressions are not easily manejable. We think that an extension of a Golub-Welsch algorithm would be very interesting to estimate the integral that appears in (

1.6) using the recurrence relations in (1.3) or (2.3), and so avoid the computation of the weights . Despite this work does not give a computation algorithm, in Section 6, we give an appropriate convergent scheme of quadrature rules.

## 2. Bi-orthogonality

In this Section we review some concepts and simplify results from [3, 4] to accommodate their statements and notation to our special case. They are necessary to prove Theorem 1.1 and Theorem 1.2 in Section 4. We consider the family of polynomials

 (2.1) Kn(z)=λnQn(λ+μ−z),n∈Z+.

Then the elements of the family satisfy the following recurrence relation

 (2.2) Kn+1(z)=zKn(z)−μλmKn−m(z),n≥0,

with initial conditions

 K−m≡K1−m≡⋯≡K−1≡0andK0=1.

We see

 Kn(z)=zn,n=0,1,…,m.

This kind of families of polynomials has been studied in several publications such as [3, 4, 21, 22, 23]. We also introduce the vectors of polynomials that satisfy the recurrence relation

 (2.3) p0=(1,0,…,0,0),p1=(0,1,…,0,0),…,pm−1=(0,0,…,0,1),zpr(z)=μλmpr+m(z)+pr−1(z),r=0,1,2,….

We have that

 (2.4) pr(z)=1λrqr(λ+μ−z).

The recurrence relations (2.2) and (2.3) are both associated to the infinite non-symmetric matrix

 ⎛⎜ ⎜ ⎜ ⎜ ⎜⎝000⋯0μλm00⋯100⋯00μλm0⋯010⋯000μλm\par⋮⋱⋱⋱⎞⎟ ⎟ ⎟ ⎟ ⎟⎠⊤,

which determines an operator in the space of sequences . Set the standard basis of . is defined as follows

 Le0=μλmemand% Lej=ej−1+μλmem+j,j=1,2,….

Let be a polynomial with an arbitrary degree . Then , denotes the operator:

 q(L)u=(bκLκ+⋯+b1L+b0)u=bκLκu+⋯+b1Lu+b0u,u∈ℓ2.
###### Lemma 2.1.

Consider the sequence of polynomials satisfying the recurrence relation (2.2). Then

 en=Kn(L⊤)e0,n∈N={1,2,…}.
###### Proof.

We recall that for each we have that , then . This means that

 K1(L⊤)e0=L⊤e0=e1,K2(L⊤)e0=(L⊤)2e0=L⊤e1=e2,⋮Km(L⊤)e0=(L⊤)me0=(L⊤)n−1e1=⋯=L⊤em−1=em.

We now consider . For each , let us assume that . By the recurrence relation (2.2), we have that

 Kn+1(L⊤)=L⊤Kn(L⊤)−μλmKn−m(L⊤),

hence

 Kn+1(L⊤)e0=L⊤Kn(L⊤)e0−μλmKn−m(L⊤)e0=L⊤en−μλmen−m
 μλmen−m+en+1−μλmen−m=en+1.

The proof is completed by induction. ∎

Given an arbitrary vector polynomial we denote for every system with , the ”dot” product

 q(L)⨀U=m∑j=1qj(L)uj.
###### Lemma 2.2.

Consider the sequence of vector polynomials satisfying the recurrence relation (2.3) and denote . Then

 er=pr(L⊤)⨀E,r∈Z+={0,1,2,…}.
###### Proof.

We observe that for each

 pr(L⊤)⨀E=(r−20,…,0,1,0,…,0)⨀(e0,…,em)=er.

Let us consider . For each , let us assume that . By the recurrence relation (2.3), we have that

 μλmpr+m(L)⨀E=Lpr(L)⨀E−pr−1(L)⨀E.

This implies that

 μλmpr+m(L)⨀E=Ler−er−1=λμmer+m+er−1−er−1=μλmer+m.

The proof is completed by induction. ∎

In we consider the regular inner product. This is that given two and elements of

 u⋅v=∞∑j=0uj¯¯¯vj=uv⋆,

where represents the complex conjugate transpose operator of . When is real .

###### Theorem 2.3.

Set . The polynomial in (2.2) and in (2.3) satisfy the orthogonality relations:

 (2.5) Kn(L)pr(L)⨀E⋅e0=δn,r.
###### Proof.

From Lemma 2.1 and Lemma 2.2 we have immediately

 Kn(L)pr(L)⨀E⋅e0=pr(L)⨀E⋅Kn(L⊤)e0=er⋅en=δr,n.

This proves (2.5). ∎

## 3. An algebraic equation: some properties

In [4] the following the algebraic equation is sudied:

 (3.1) ωm+1−zωm+μλm=0.

We enumerate some needed properties of its multivalued solution , which is an algebraic function of order and genus . Its inverse rational function

 (3.2) z(ω)=ω+μλmωm,

gives the composition of a conformal map of the Riemann sphere to a Riemann surface of the function and the projection of to the complex plane. In order to describe this surface we introduce certain concepts and notations, and also review some results stated in [4]. Given a complex numbers we denote by with , the ray starts at with slope . Set

 a=(m+1)m+1√μλmmm/(m+1).

Let us declare the following sets:

 (3.3) S0=m⋃k=0[0,aexp2πikm+1],Se∞=m⋃k=0[0,∞)×exp((2k+1)iπm+1),So∞=m⋃k=0[0,∞)×exp(2kiπm+1).

We also introduce

 Ω0=C∖S0,Ωj=C∖{So∞∪Se∞},1≤j

According to [4, Proposition 1] there exist global branches of in (3.1) such that

 ωj∈H(Ωj),j=0,1,…,m,

and

 (3.4) |ω0|≥|ω1|≥⋯≥|ωm|,on C.

In the inequalities are strict.

Let us describe then the Riemann surface associated to the multi-valued function such that the sheet structure corresponds to the choice of branches as above in (3.4). Let denote the projection of on , with (multivalued) representing its inverse branches. We have that

 R=¯¯¯¯¯¯¯¯¯¯¯¯¯¯m⋃j=0Rj,R0=π−10(Ω0∪{∞}),Rj=π−1j(Ωj),j=1,…,m.

We also denote a closed contour in which separates and . The contour is orientated such that lies on the left, this is positive side. We take

 π(∂R0,1)=∂Ω0=S0+∪S0−,π(∂Rj,j+1)=⎧⎪⎨⎪⎩Se∞+∪Se∞−,j=2l−1Se∞+∪So∞−,j=2lj=1,…,m−1.

Hence

 R={m⋃j=0Rj}⋃{m−1⋃j=0∂Rj,j+1}.

We now sketch the proof of [4, Proposition 1] where some properties of the multivalued function in (3.1), are stated and they are needed for our purposes. From Vieta’s relations in (3.1) we have that the solutions of satisfy that

 (3.5) m∑k=0ωk(z)=zandm∏k=0ωk(z)=μλm.

Combining (3.2) and (3.5) we have that (multivalued function) has a branch point at infinity of order . Then we suppose that ’s branches, denoted by , have the following behavior

 (3.6) ω0(z)=z+O(1)andωj(z)=O(1z1/m),j=1,…,m,asz→∞.

First we analyze for . Notice when . The inverse function (see (3.2)) is decreasing in the interval and increasing on . This means that there are two positive branches of in . We denote them by and such that

 ω0(x)>m+1√mμλm>ω1(x),x∈(a,+∞).

When is even there also is a negative branch of on , that we denote by . Then we can introduce the notation

 (3.7)

Let us define an analytic continuation of from to . Consider a solution of (3.1), then the following equality holds

 (e2πim+1ω(z))m+1−ze2πim+1(e2πim+1ω(z))m+μλm=0.

So for any fixed there exists a choice of branch , satisfying the symmetric relation:

 ωk(ze2πim+1)=e2πim+1ωj(z),z∈C.

We choose ,

 (3.8) ωk(ze2πim+1)=e2πim+1ωj(z),z∈C.

Using Monodromy Theorem (see [1, Chapter I 15]), we extend analytically the branch from a neighborhood of to the domain . We declare . Since (3.8) is satisfied near for , then

 ω0(ze2πim+1)=e2πim+1ω0(z),z∈Ω0.

Taking into account the behavior at infinity in (3.6), the branch has an analytic continuation from to the simple connected domain

 ~A={z:|argz|<πm+1}∖[0,a].

Starting from above set and using (3.8) with we extend analytically as follows

 ω1(ze2πim+1)=e2πim+1ω1(z),z∈Ω1={Ω0m=1Ω0∖Se∞m>1.

Since , we have that

 ω1(z)=¯¯¯¯¯¯¯¯¯¯¯¯ω1(¯¯¯z),z∈C.

So combining the above identity with (3.8) we have that satisfies the jump condition

 ¯¯¯ω1+(z)=e2πim+1ω1−(z),z∈Se∞,

where

 2πiω1±(x)=limz→x,z∈Se∞±ω1(z),x∈∘Se∞,

and is the closure of . Hence .

We proceed analogously with the remaining branches of . Consider the fundamental sector of

 A={z:|argz|<πm+1}.

We extend the branches to using the symmetry relations in (3.8):

 (3.9) ωj(ze2πim+1)=e2πim+1ωj(z),z∈Ωj,j=2,…,m.

A branch with even index, , is defined in as a direct analytic continuation of the branch to the domain from the upper part of the boundary of , and to the domain from the lower part of the boundary of . While a branch with odd index, , is defined on taking the analytic continuation of the branch to from the lower part of the boundary of , and to from the upper part of the boundary of .

Consider the domain

 B={z:0

and define the function as follows

 ω(z)=⎧⎪ ⎪⎨⎪ ⎪⎩ω1(z),πm+1

Let us recall (3.7) and (3.8). If the maximum principle to the holomorphic function we can extend the inequality in to and therefore to the whole complex plane . Analogously we can extend all the relations stated in (3.7) to and obtain (3.4).

Once we have sketched the proof of the relations (3.4), we define

 (3.10) 2πiω0±(x)=limz→x,z∈S0±ω0(z),x∈∘S0=m⋃k=0(0,aexp2πikm+1).

Using (3.4) we observe that the following functions defined on , are real valued and never change sign.

 (3.11) ˜ρj(x)=12πi⎡⎣1ωj0+(x)−1ωj0−(x)⎤⎦=1ωj0(x)−1ωj1(x),j=1,…,m.

So given ,

 ˜ρj(x)=˜ρ1(x)r2,j(x),j=1,…,m,

where

 r2,1≡1andr2,j(x)=j−1∑k=01ωj−1−k0(x)ωk1(x),x∈S0,j=2,…,m.

Since never vanishes in and taking into account (3.6) we have that

 (3.12) zj−1ωj0(z)∈H(Ω0)andzj−1ωj0(z)=O(1z)asz→∞,j=1,…,m.

Now using (3.11) and Sokhotski–Plemelj (see [14, Chapter I, equality (4.8)]) we obtain that:

 (3.13) 1ωj0(z)=1zj−1∫xj−1r2,j(x)˜ρ1(x)dxx−z,j=1,…,m.

Let us express the polynomials satisfying the recurrence relation (2.2) and the vector polynomials in (2.3), in terms of the solutions of the equation (3.1).

###### Theorem 3.1.

Fix . The polynomial , can be expressed

 (3.14) Kn(z)=m∑j=0ωn+m+1jωm+1j−mμλm.

Given , the vector polynomial has the form

 pr=(p0,r,p1,r,…,pm,r),

with

 (3.15) pj,r=μλmωm+1j−mμλm(1ωr0+ωjωr1+⋯+ωmjωrm),j=0,1,…,m.
###### Proof.

In (3.14) we follow the steps of [3]. We observe that the following function satisfies the recurrence relation (2.2)

 Kn=A0ωn0+A1ωn1+⋯+Amωnm.

The coefficients are defined from the initial conditions:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩A0+A1+A2+⋯+Am=1A0ω0+A1ω1+A2ω2+⋯+Amωm=zA0ω20+A1ω21+A2ω22+⋯+Amω2m=z2⋯⋯⋯⋯⋯⋯A0ωm0+A1ωm1+A2ωm2+⋯+Amωmm=zm

From Cramer’s rule we find

 (3.16) Ai=m∏k=0,k≠j(z−ωk)m∏k=0,k≠j(ωj−ωk).

Using now Vieta’s formulas in (3.5), we arrive at

 m∏k=0(z−ωk)=zm+1−m∑j=0ωjzm+⋯+(−1)m+1m∏j=0ωj=μλm.

We also have that

 ωm+1j−zωmj+μλm=0⟹μλmz−ωj=ωm.

Then

 m∏k=0,k≠j(z−ωk)=μλmz−ωj=ωmj.

Consider the functions

 Pz(ω)=m∏k=0(ω−ωk)=ωm+1−zωm+μλn.

Then the denominator of in (3.16) can be written as:

 m∏k=0,k≠j(ωj−ωk)=P′z(ωj)=(m+1)ω