UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Reparameterization of special contours for boundary value problems in queueing systems Epema, Cody 2015

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2015_september_Epema_Cody.pdf [ 683.79kB ]
JSON: 24-1.0166496.json
JSON-LD: 24-1.0166496-ld.json
RDF/XML (Pretty): 24-1.0166496-rdf.xml
RDF/JSON: 24-1.0166496-rdf.json
Turtle: 24-1.0166496-turtle.txt
N-Triples: 24-1.0166496-rdf-ntriples.txt
Original Record: 24-1.0166496-source.json
Full Text

Full Text

Reparameterization of SpecialContours for Boundary ValueProblems in Queueing SystemsbyCody EpemaB.Sc., The University of British Columbia, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE COLLEGE OF GRADUATE STUDIES(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)July 2015c© Cody Epema, 2015AbstractWe present several tests to detect whether a simple reparameterizationof the contour is possible for the Riemann-Hilbert Boundary Value Problemassociated with the stationary distribution of a two dimensional queueingsystem. We cover several special cases, some of which have been previouslycovered using different methods, and present an expansion on the existingtheory to the case of elliptical contours. Two fast heuristic tests are alsopresented. Additionally, we construct a formal BVP solution for queueingsystems with homogeneous fundamental equations.iiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . ixChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 11.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Previous Work in the Field . . . . . . . . . . . . . . . . . . . 1Chapter 2: Background . . . . . . . . . . . . . . . . . . . . . . . 32.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2 Queueing Systems . . . . . . . . . . . . . . . . . . . . . . . . 42.2.1 Two Dimensional Queuing Systems . . . . . . . . . . . 42.2.2 Jump Probability Representation for a Random Walk 52.2.3 Stationary Distribution . . . . . . . . . . . . . . . . . 62.3 Generating Function Approach for Queueing Systems . . . . 102.3.1 Generating Functions . . . . . . . . . . . . . . . . . . 102.3.2 Indexes of Functions . . . . . . . . . . . . . . . . . . . 142.3.3 Riemann Hilbert BVPs . . . . . . . . . . . . . . . . . 152.4 The Beginning of the BVP Approach . . . . . . . . . . . . . . 172.4.1 The Kernel Method and Contours . . . . . . . . . . . 17iiiTABLE OF CONTENTS2.4.2 Counterparts . . . . . . . . . . . . . . . . . . . . . . . 19Chapter 3: Special Contours . . . . . . . . . . . . . . . . . . . . 213.1 Circles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.1.1 The Case for a Circle with a Center on the Origin . . 213.1.2 The Case for a Circle with a Center on the Real axis . 223.2 Ellipses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.2.1 Determining Parameters . . . . . . . . . . . . . . . . . 243.2.2 The Case for an Ellipse with a Major Axis on the RealAxis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.2.3 The Case for an Ellipse with a Minor Axis on the RealAxis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.3 A Fast Informal Test for Circular Contours . . . . . . . . . . 29Chapter 4: Homogeneous Queueing BVPs . . . . . . . . . . . . 334.1 General Theory for Homogeneous Fundamental Equations . . 344.1.1 Construction of New Branches . . . . . . . . . . . . . 344.1.2 Bounds for X0(y) and Y0(x) . . . . . . . . . . . . . . . 364.1.3 Solution via BVP . . . . . . . . . . . . . . . . . . . . . 39Chapter 5: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 43Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45ivList of TablesTable 2.1 Legend and Table of Counterparts . . . . . . . . . . . 20vList of FiguresFigure 2.1 The regions Sr of the state space. . . . . . . . . . . . 7Figure 2.2 A common representation scheme for the jump prob-abilities for a random walk in the quarter plane. . . . 8Figure 3.1 The two cases for an elliptical contour. . . . . . . . . 26Figure 3.2 The test for a circle with a center at the origin. . . . 30Figure 3.3 The test for a general circular contour. . . . . . . . . 31viAcknowledgmentsI would like to thank my supervisor Dr. Javad Tavakoli, first for his sup-port over the whole of my stay in this program, and second for introducingme to Dr. Yiqiang Zhao of Carleton University, who I would also like tothank. Both have provided useful feedback and invaluable help on the pre-sentation and scope of this thesis, and both have been wonderful to workwith.I would also like to thank the other members of my supervisory commit-tee, Dr. Wayne Broughton, and Dr. Rebecca Tyson, for being patient wherethey had more reasons not to be and donating their valuable time to helpme revise and polish this work.viiChapter 1Introduction1.1 OverviewIn this thesis we develop several original tests for finding special casesfor the contours employed in the Boundary Value Problem (BVP) approachto Queueing Systems, and extract a formal BVP solution for the stationaryprobability distribution on the state space of a queueing system that presentsthe case of a homogeneous equation.Some of the special contours in this text have appeared in queueing BVPsbefore, including Fayolle and Iasnogorodski’s 1979 paper [3], Cohen andBoxma’s 1983 book [2], and Avrachenkov, Nain, and Yechali’s 2014 paper[1]. Our goal in creating tests for other elegantly parametrizable contours isto allow further work to be done in this line, create easy tests to see whatcan be connected to the methods applied in similar problems, and hopefullyto allow for analytic solutions of certain BVPs rather than the formal oneswhich are most commonly seen. We will address these new results for specialcontours in chapter 3. Chapter 4 details a formal solution for a broad class ofqueueing BVPs, those with homogeneous fundamental equations. The restof chapter 1 provides a timeline for historically important texts in this fieldand a few recent queueing papers. Chapter 2 provides the mathematicalbackground for the work contained in this thesis, and the BVP methodsused in our approach.1.2 Previous Work in the FieldThe genesis of the BVP approach to queueing systems can be foundin the 1979 paper authored by G. Fayolle and R. Iasnogorodski [3], who11.2. Previous Work in the Fieldpioneered the method, and dealt with one case. Of particular interest hereis that this paper involved a circular contour.A text was published in 1988 by J. Cohen and O. Boxma [2], which dealtwith several cases of queueing systems and explained the required back-ground information for the Riemann-Hilbert BVPs in detail. Several caseswere also treated in detail here, including the very important ‘Join ShortestQueue’ model, which covers one of the most common arrival disciplines.In 1999, G. Fayolle, R. Iasnogorodski, and V. Malyshev published intheir own book [4] a broad class of random walks and the associated queue-ing systems were dealt with via BVP methods, by transforming the initialproblem into a Riemann-Hilbert-Carelman BVP with a shift, the ‘shift’ re-ferring to an automorphism of the independent variable appearing in someof the terms of the boundary condition. These authors also developed an-other approach, using covering spaces for the algebraic curve defined by thekernel equation. Said kernel equation (in this thesis, equation (2.43)) is de-termined by the fundamental equation for the random walk, so this algebraicapproach starts in the same way as the BVP approach. However, this alge-braic approach has a disadvantage compared to BVP methods, because itonly works for certain fundamental equations, depending on the propertiesof what they declare to be the group of the random walk, which is definedusing Galois automorphisms related to the algebraic curve associated withthe kernel equation.More recent papers in BVP queueing include [12], by J. Resing andL. O¨rmeci, published in 2003, [8], by F. Guillemin, C. Knessl, and S. vanLeeuwaarden, published in 2013, and [1], by K. Avrachenkov, P. Nain, andU. Yechiali, published in 2014. Avrachenkov et al. study a system wherethe state of the server is unknown to the queues and vice versa, so an emptyserver could take a job from either one of the orbiting queues, which retryjobs that were sent to wait, or from a new incoming job. Guillemin etal. model a chain of wireless transceivers whose transmissions can interferewith one another. The paper [8] posits a slightly different approach toqueueing BVPs, but the paper [1], while having a very interesting setup forthe fundamental equation, follows the BVP setup found in [3] quite closely.2Chapter 2Background2.1 NotationIn the interest of being unambiguous we will define our notations forseveral common operations:− For vectors v and w ∈ Cn we define the exponentiation of vectors:vw =n∏i=1vwii .− For a boolean statement V , we denote the indicator function 1 as1V ={1 if V0 if not V,− For sets A and B, where A is a subset of an additive group G, and foran element b ∈ G, we denote cosets and set difference operation:A+ b = {a+ b : a ∈ A},A− b = {a− b : a ∈ A},A\B = {e : e ∈ A, e 6∈ B}.− We define a neighborhood of a point z in C as a bounded simplyconnected open subset of C containing that point.− We also denote for a set A ⊆ C the closure of A as A, the interior of32.2. Queueing SystemsA as int A, and the boundary of A as bdy A. Specifically:A = A ∪ {l : ∃(zn)∞n=1 ⊆ A, limn→∞zn = l}.int A = {z ∈ A : z is in a neighborhood that is a subset of A}.bdy A = A \ int A.− For a closed simple contour B ⊂ C, we will denote the interior of thecontour under a positive (counter-clockwise) traversal as GB.− For a complex number z = a + ib, a, b ∈ R, we denote the complexconjugate z = a− ib.2.2 Queueing Systems2.2.1 Two Dimensional Queuing SystemsThe application at hand is a queueing system wherein jobs arrive intovarious queues where they wait for service by a processor. A processor maysend a job into a different queue, changing the job type, or remove it fromthe system. Each processor may only handle one type of job. The processorsdon’t have all the information about the system; they can only tell if thereis a job waiting in queue or not.We assume that the arrival rate of jobs and the service rate of the pro-cessor are independent of the number of jobs waiting in the system, of oneanother, and of the current time. Specifically, we can treat the lengths ofthe queues as the state space of a time-independent Markov chain, with adiscrete time variable, so that if we write down the length of queue i at timen as the random variable Xin, the following equations in the expected valuesof these variables represent the stated independence, given that n > m+ t,and that m, t ≥ 0, and that v, vi are some nonnegative integers:E[Xin|Xim+t = vm+t]= E[Xin|Xim = vm, . . . , Xim+t = vm+t], (2.1)E[Xin|Xim+t = v]= E[Xin+s|Xim+t+s = v]. (2.2)42.2. Queueing Systemswhere E[A|B] is the expected value of a function A given that B.Equation (2.1) states that the information of how one got to the currentstate is unimportant for predicting the future behavior of the system, andequation (2.2) states that the amount of time that the system took to arriveat its current state is similarly unimportant. These conditions are necessaryfor the approaches to analyzing this system that appear in the next threesubsections, as they rely heavily on the theory of Markov chains.Since a negative queue length is not meaningful in most applications, werestrict Xin to be non-negative. The focus of our approach will be on systemsthat either start as or reduce to two coupled queues, and so in all the workthat follows we will have Z+ × Z+, all pairs of nonnegative integers, as ourstate space, and denote the Markov chain for the state of the random walkat time n as the vector of random variables Xn =〈X1n, X2n〉.2.2.2 Jump Probability Representation for a Random WalkUp to this point, our time variable could have taken any real value, butwe want to work with a discrete time Markov chain, so we note that evenin a continuous time Markov chain, two events happening at exactly thesame time is analogous to randomly picking two real numbers and havingthem come out identical, a probabilistic impossibility which can be safelyneglected. We subdivide the real line that is the domain of our time variableinto intervals of equal size small enough that it is safe to say only one eventwill happen in each. In this way we will claim the random walk has boundedjumps, so the position in the state space for queue length will only everincrease or decrease by one at a time. This formulation of the problemallows us to work with a discrete space, discrete time, time independentMarkov chain that tracks both queues, both of which are bounded to havenon-negative queue lengths, which corresponds to a random walk in thequarter plane.All the processors in our queueing system can only tell if their queueis empty or nonempty, which implies a condition referred to as spatial ho-mogeneity. The state space can be divided up into regions where the jump52.2. Queueing Systemsprobabilities at any given point in a given region is the same as any otherpoint in the same region. These regions for the random walk in the quarterplane are depicted in figure 2.1, and are defined as follows:S = {1, 2, 3, . . . }2, (2.3)S′ = {1, 2, 3, . . . } × {0}, (2.4)S′′ = {0} × {1, 2, 3, . . . }, (2.5)S0 = {(0, 0)}. (2.6)It follows that we can define prij as the probability of transitioning fromstate (a, b) to state (a+ i, b+ j) first, for (a, b) in a given region r:pri,j = P(X1t+1 −X1t = i ∧X2t+1 −X2t = j|(X1t , X2t ) ∈ r). (2.7)Rather than write pSi,j or pS′′i,j , we will sometimes write pi,j or p′′i,j , and so on.Spatial homogeneity is equivalent to saying that we could instead condi-tion equation (2.7) on (X1t , X2t ) being any specific element of a given regionr, and it wouldn’t change the result, which means that knowing the proba-bility of being in any specific state is not needed to calculate pri,j .It is common in the literature to create a diagram that represents thetransition probabilities from one state to another for the associated Markovchain for the queueing system. Since we have spatial homogeneity here, it iscommon to represent these transition probabilities with diagrams similar tofigure 2.2, omitting those arrows that are associated with a jump probabilityof zero. Probabilities of zero commonly occur in this application, as queueingtheorists tend to study special cases.2.2.3 Stationary DistributionNow that the rudiments are available, we can state the problem aroundwhich this thesis is focused. We want to find the long-run or stationarydistribution over the state space, that is, the spatial distribution of theprobability of the system being in a specific state after the queue has run62.2. Queueing SystemsFigure 2.1: The regions Sr of the state space.for an arbitrarily long period of time. We define the stationary probabilityfor a point (i, j) in the state space aspiij = limn→∞P (Xn = (i, j)). (2.8)The stationary distribution encoded in equation (2.8) is the most importantdefinition in this section. It does not depend on the time variable n, andsatisfies the property that for an arbitrary time t0, if P (Xt0 = (i, j)) =piij , ∀(i, j) ∈ Z2+ then P (Xt0+1 = (i, j)) = piij ,∀(i, j) ∈ Z2+.From the theory of discrete space discrete time Markov chains [13] it isknown that for stable stationary distribution to exist, a Markov chain musthave all of the following properties:− Irreducible: There is a nonzero probability of eventually getting fromany one state to any other state.Irreducibility dictates that the Markov chain cannot get ‘stuck’ in aproper subset of the state space.72.2. Queueing SystemsFigure 2.2: A common representation scheme for the jump probabilities fora random walk in the quarter plane.− Aperiodic: The greatest common divisor of the set of the number ofsteps taken to go from any one specific state to itself with nonzeroprobability is one.Aperiodicity dictates that the state space cannot be divided into propersubsets such that the Markov chain always arrives in a specific subsetgiven time arguments in arithmetic progression.− Ergodic: The mean number of steps taken to return to any given statefrom itself is finite. 1We will primarily assume irreducibilty and aperiodicity for our pur-poses, although in the application of two-dimensional random walks, these1In [4], Fayolle et al. characterize ergodicity with stochastic matrices: If P is thestochastic transition matrix for an irreducible aperiodic Markov chain, then the saidMarkov chain is ergodic if and only if there exists a row vector pi with non-negativeentries that sum to 1 such that piP = pi. This is in fact how stationary distributions arecharacterized later in this chapter. We feel that this alternate definition offers no positivecharacterization of what ergodicity actually means.82.2. Queueing Systemsare easily characterized. For example, for there to be no chance of get-ting from any specific state to another, the queueing system would haveto have at least one nonfunctioning queue, (a queue with no arrivals or noservice,) or both queues would behave like one queue, with linked arrivalsand service. An example case of a periodic Markov chain would be whenp01,0 = 1, p′−1,1 = 1, p′′0,−1 = 1, resulting in an endless cycle, where jobs areinstantly served after arrival. In short, if our type of analysis is required,then the Markov chain is almost certainly irreducible and aperiodic.Ergodicity is dealt with in [4, page 3] by looking at the mean jump sizesin the x and y directions: if we writeMx =∑ipij , My =∑jpij ,M ′x =∑ip′ij , M′y =∑jp′ij ,M ′′x =∑ip′′ij , M′′y =∑jp′′ij ,(2.9)then, given that Mx and My are both not identically zero, the random walkis ergodic if and only if one of the following three sets of conditions holds:Condition 1:Mx < 0,My < 0,MxM′y −MyM′x < 0,M ′′xMy −M′′yMx < 0.(2.10)Condition 2:Mx < 0,My ≥ 0,M ′′xMy −M′′yMx < 0.(2.11)92.3. Generating Function Approach for Queueing SystemsCondition 3:Mx ≥ 0,My < 0,MxM′y −MyM′x < 0.(2.12)These conditions essentially embody the idea that if a random walk tendsto move away from the origin, it will grow past any bound, and will likelynever return to any given finite state.One needs to verify that the queueing system that one wants to modelis representable with a irreducible, aperiodic, ergodic Markov chain beforeusing the work that appears in the rest of the thesis. These conditions implythe existence of a unique stationary probability distribution for the Markovchain, and the existence of this distribution underpins most of the work thatfollows in this thesis.Calculation of the stationary distribution is desirable as it allows thecalculation of valuable statistics for the queueing system, and on its owngives how one should expect the queueing system to behave in the casewhere you have no information about the system’s exact behavior at orbefore the current time, relying instead and solely on the jump probabilitiesto inform your expectations.2.3 Generating Function Approach for QueueingSystems2.3.1 Generating FunctionsA generating function is a means of encoding a data set into a function,specifically, by making its Taylor coefficients store the data set. If we have fin the polynomial ring R(x, y) with x, y ∈ C, and f(x, y) =∑∞i,j=0,0 ai,jxiyj ,then f ’s derivatives contain the information to recover ai,j , as∂m+nf∂mx∂ny(0, 0) = m!n!am,n. (2.13)102.3. Generating Function Approach for Queueing SystemsThus, by solving for the function f(x, y), either analytically or by numer-ically computing it and its derivatives, one can recover the data set thefunction encodes.We can use our knowledge of how the stationary probability distributionmust behave in different regions of the state space to define a functionalequation in the generating functions for these distributions. The study ofthis functional equation will govern the rest of this thesis. Our approachhere comes from the introduction of [4].We denote the complex variables that our generating functions willtake as arguments as x, y ∈ C, defining for brevity u = 〈x, y〉, and R ={S, S′, S′′, S0}.First, we define the jump probability generating function over region rasPr(u) = E[u(Xn+1−Xn)1Xn∈r]. (2.14)We may read equation (2.14) with the information that we have aboutour application, givingPr(u) = ExChange in X1 at n︷ ︸︸ ︷(X1n+1 −X1n)yChange in X2 at n︷ ︸︸ ︷(X2n+1 −X2n)∣∣∣∣∣∣∣∣Xn ∈ rwhich, by the definition of expected value, and the definition in equation(2.7), implies thatPr(u) =∑i,jprijxiyj . (2.15)The generating function for the stationary probability distribution in theregion r is defined as follows:pir(u) = limn→∞E[uXn1{Xn∈r}]=∑(i,j)∈rpiijxiyj (2.16)We now start fresh and begin to construct our desired result, a rela-tionship tying all of the relevant generating functions for the stationary112.3. Generating Function Approach for Queueing Systemsdistributions together through generating functions for the jump probabil-ities. Note that since {Xn}∞n=1 is a Markov chain, the variables Xn+1 andXn+1 −Xn are independent, allowing us to factor E[uXn+1]as follows:E[uXn+1]= E[uXn+1−XnuXn]= E[uXn+1−Xn]E[uXn], (2.17)Since the expected value is just a summation over the state space, and (Xn)is a Markov chain with each step forward only dependent on the currentstate, we can split up equation (2.17) over the different regions and intoindependent parts:E[uXn+1]=∑r∈RE[u(Xn+1−Xn)1{Xn∈r}]E[uXn1{Xn∈r}], (2.18)and by equation (2.14), we can simply identify Pr(u), givingE[uXn+1]=∑r∈RE[uXn1{Xn∈r}]Pr(u) (2.19)Now, all we need to do is note thatE[uXn+1]=∑r∈RE[uXn+11{Xn+1∈r}], (2.20)and take the limit as n → ∞ of both sides of equation (2.18) to get ourdesired relationship in the generating functions, applying equations (2.20)and (2.16) as needed:limn→∞[∑r∈RE[uXn+11{Xn+1∈r}]]= limn→∞[∑r∈RE[uXn1{Xn+1∈r}]]Pr(u)=∑r∈Rpir(u) =∑r∈Rpir(u)Pr(u),(2.21)122.3. Generating Function Approach for Queueing Systemswhich simplifies to∑r∈R[1− Pr(u)]pir(u) = 0. (2.22)This equation (2.22) characterizes how the stationary distribution be-haves for random walks in the quarter plane explicitly, and independently ofthe division the quarter plane into regions with distinct jump probabilities.However, equation (2.22) is rather abstract, and we will want to use theinformation that we have for our queueing system to make it more concrete,using the regions and jump probabilities as defined above. After specifyingthese, equation (2.22) is almost always factored as in the following equation(2.23), referred to as the fundamental equation (F.E.) for the random walk.Q(x, y)pi(x, y) = q1(x, y)pi1(x) + q2(x, y)pi2(y) + q0(x, y)pi0. (2.23)We define the all of the functions present in equation 2.23 (with referencesto equation (2.22)) as follows:Q(x, y) = xy−1 +∑1≥i≥−11≥j≥−1pijxiyj = −(1− PS(u)), (2.24)q1(x, y) = x1−∑1≥i≥−11≥j≥0p′ijxiyj = 1− PS′(u), (2.25)q2(x, y) = y1−∑1≥i≥01≥j≥−1p′′ijxiyj = 1− PS′′(u), (2.26)q0(x, y) = 1−∑1≥i≥01≥j≥0p0ijxiyj = 1− PS0(u), (2.27)pi(x, y) =∞∑i,j=1piijxi−1yj−1 = piS(u), (2.28)132.3. Generating Function Approach for Queueing Systemspi1(x) =∞∑i=1pii0xi−1 = piS′(u), (2.29)pi2(y) =∞∑j=1pi0jyj−1 = piS′′(u), (2.30)pi0 = pi0,0 = piS0(u). (2.31)This fundamental equation is so named because solving for all of theunknown functions (pi(x, y), pi1(x), pi2(y), pi0,) in it allows recovery of thestationary distribution encoded in the generating functions. Chapter 4 isprimarily focused on a formal integral solution to this equation, which isuseful in numerical work.In order to solve for the unknown functions in this equation, we willbe using the theory of Riemann Hilbert boundary value problems, whichrequire the following:2.3.2 Indexes of FunctionsIf we have a complex function G defined over a simple closed piecewisesmooth positively oriented contour L, the index χ of the function over thecontour is a quantity we will define asχ =12piargLG =12pii∫Ld(Log(G(t)) =12pii∫LG′(t)G(t)dt, (2.32)assuming without loss of generality that 0 is in the interior of L. [4]Interpreting the second formula above, this index counts the total num-ber of times the function G goes over the branch cut of the logarithm underarguments coming from a positively oriented traversal of L, with a negativeindex corresponding to a net clockwise traversal in the image, and a positiveindex corresponding to a net counterclockwise traversal in the image.Cauchy’s argument principle [11] also grants another interpretation ofthe index, as the third formula in equation (2.32) gives the difference betweenthe number of roots and poles of G inside L, with roots increasing theindex and poles reducing it, with the multiplicity of both contributors being142.3. Generating Function Approach for Queueing Systemsaccounted for.2.3.3 Riemann Hilbert BVPsThe Riemann Hilbert boundary value problem was used by Cohen andBoxma to analyze queueing systems in their book [2]. This particular BVPwas also employed in the papers [1, 3, 8] to find formal solutions. Theformulation given here was provided in [2].The Riemann Hilbert boundary value problem is as follows: We wish tofind a function Φ defined over the complex plane holomorphic on the interiorand exterior of a given closed simple contour L excepting perhaps a pole offinite order at infinity. Additionally, Φ has limiting values coming from theinside and outside of L, labeled Φ+ and Φ− respectively.Defining the notation GA as the interior of a given closed simple contourA without including the contour itself, we can define these explicitly asΦ+ : L→ C, Φ+(x) = limt→xt∈GLΦ(t), (2.33)Φ− : L→ C, Φ−(x) = limt→xt6∈GL∪LΦ(t). (2.34)The values of Φ are unknown at the start, but if we have a relationbetween the limiting values Φ+ and Φ− given byG(x)Φ−(x) + g(x) = Φ+(x), x ∈ L, (2.35)the Riemann Hilbert boundary value problem is solvable given the fol-lowing conditions:1. G is continuous on L,2. G is nonzero on L,3. g satisfies the Ho¨lder condition |g(t1) − g(t2)| ≤ C|t1 − t2|µ on L, forreal constants C and 0 ≤ µ < 1.152.3. Generating Function Approach for Queueing SystemsIf we defineχ = the index of G over L,Γ(x) =12pii∫LLog(t−χG(t))dtt− x,(2.36)S+(x) = eΓ(x), x ∈ GL,S−(x) = x−χeΓ(x), x 6∈ GL ∪ L,(2.37)the solution to the homogeneous case of equation 2.35, where in the equationg ≡ 0, is given byΦ(x) ={S+(x) if x ∈ GL,S−(x) if x 6∈ GL ∪ L.(2.38)For the nonhomogeneous case, we have two types of solutions, both usingthe following definition:γ(x) =12pii∫Lg(s)dsS+(s)(s− x). (2.39)Which type of solution we use depends on the index:− When χ ≥ 0:Φ(x) ={S+(x)(γ(x) + Pχ(x)) if x ∈ GL,S−(x)(γ(x) + Pχ(x)) if x 6∈ GL ∪ L.(2.40)where Pχ is an arbitrary polynomial of degree χ over C.− When χ < 0:Φ(x) ={S+(x)γ(x) if x ∈ GL,S−(x)γ(x) if x 6∈ GL ∪ L.(2.41)but note that when χ ≤ −2, there are −χ − 1 additional conditions162.4. The Beginning of the BVP Approachthat must be satisfied for the BVP to have a solution, given by∫Lg(t)tk−1dtS+(t)= 0, for k ∈ {1, 2, . . . ,−χ− 1}. (2.42)The Riemann Hilbert BVP turns out to be just the right tool for the laststep in solving for one of the unknown functions in fundamental equation(2.23) in our application. We will use Φ to recover pi1 and pi2 individuallyand independently, by finding a G which satisfies the properties required bythe BVP setup for each case, the construction of which will leave g = 0.For more details on the construction of the solution to this boundary valueproblem, see [2] or [4].2.4 The Beginning of the BVP Approach2.4.1 The Kernel Method and ContoursNow equipped with information on the tool we will later use, we canbegin work on the problem proper. The typical starting point that a greatnumber of papers (for example, [1, 3, 5, 8]) use to tackle related problemsis here referred to as the kernel method, wherein we restrict the values of xand y to the algebraic curveQ(x, y) = 0. (2.43)This restriction does two jobs: First, the complexity of the fundamentalequation (2.23) is reduced by eliminating one unknown function, leavingus with only two unknown functions, which is consistent with the BVPapproach, as long as those functions satisfy the required properties. Second,it allows the construction of a contour for the BVP. If we factor Q asQ(x, y) = a(x)y2 + b(x)y + c(x), (2.44)172.4. The Beginning of the BVP Approachwe can apply the quadratic formula on the coefficients a(x), b(x), and c(x).We label the resulting two-branched function asY (x) =−b(x)±√b2(x)− 4a(x)c(x)2a(x), (2.45)and label the individual branches Y + and Y − corresponding to the choiceof sign for the ± sign in the right hand side of equation (2.45) asY +(x) =−b(x) +√b2(x)− 4a(x)c(x)2a(x), (2.46)Y −(x) =−b(x)−√b2(x)− 4a(x)c(x)2a(x). (2.47)Often newcomers to this discipline wonder whether restricting the vari-ables x and y by the relationship Q(x, y) = 0 implies that we have a re-stricted solution, and that perhaps what we have generated by the BVP isonly valid on this algebraic curve. The proper way of thinking of this re-striction is to recognize that the BVP only requires its boundary condition(equation (2.35)) to be met on the boundary, and that setting Q(x, y) = 0gives us this boundary. Using Q(x, y) = 0 to meet the condition on theboundary is fine, since the BVP solution allows us to use the informationon the boundary to move past the boundary by design. As long as theconditions are met the BVP gives us an unrestricted analytic solution.In order to get a contour in the complex plane we investigate the branchpoints of the function Y , that is to say, the points where both branches areequal, or, equivalently where the discriminant D(x) = b2(x)−4a(x)c(x) = 0.The degree of D(x) is either three or four 2 , and so there are at most fourof these branch points, call them xi, i = 1..4. If the degree of D(x) is onlythree, we write x4 = ∞. We know from the literature [4] that these xi2The degree of D is derived directly from the degrees of a(x), b(x), and c(x), which inturn are defined by the boundedness of the jumps in the queueing system’s state space toone away in any given direction, which in turn is derived from our queueing system onlyhaving one job entering the system, switching queue, or leaving the system at any giventime.182.4. The Beginning of the BVP Approachhave to be real, with the exception of x4 which is possibly∞, and that theyadmit the ordering 0 ≤ x1 < x2 ≤ 1 ≤ x3 < x4, with x2 6= x3. Restricted tothe real line, D(x) changes sign at each of its branch points, being negativebetween x1 and x2, and also between x3 and x4.From here, there are two equivalent ways of proceeding. For a definitionthat allows easy numerical computation and parameterization we can definethe contour L as the image of the line segment x1x2 under Y +(x), joinedwith the image of the line segment x2x1 under Y −(x):L = {Y +(x) : x ∈ [x1, x2]} ∪ {Y−(x) : x ∈ [x1, x2]} (2.48)Alternatively, we can define L as the limiting contour of the image undera branch of Y of a family of contours in the x-plane cut between x1 and x2,i.e. C\x1x2, that goes around the cut. This gives the same result as theabove procedure, with the advantage that we can treat that branch of Yas analytic in that cut plane, and state L in concrete terms with respect tothat branch of Y . This is mainly of use in proofs and analysis rather thannumerical work. Looking ahead, and using a specific branch of Y given inthe future equation (4.1.1), we can give the alternate characterization:L = lim↓0{Y0(x) : x ∈ {[x1, x2] + } ∪ {[x1, x2]− }} (2.49)Additionally, we define Lext identically, except using x3 and x4 insteadof x1 and x2.While we remain on L, we know that Q(x, y) is 0, allowing us to usethe modified F.E. as a relation on that boundary. From this point, BVPapproaches diverge in how they attack the problem, but most follow thiscontour and fundamental equation setup very closely.2.4.2 CounterpartsWe note that there is more than one way to factor Q(x, y): the alternatefactorization Q(x, y) = a˜(y)x2 + b˜(y)x+ c˜(y) allows us to proceed in exactlythe same way as before, and define a new contour, a new discriminant, and192.4. The Beginning of the BVP Approachso on. Rather than double the length of this thesis, we will introduce a tableof counterparts (table 2.4.2) giving names for the constructions implied bythe work we do following our admittedly arbitrary choice of factorization.Table 2.1: Legend and Table of Counterpartsa(x) path Counterpart c.f.a(x) a˜(y) (2.44)b(x) b˜(y) (2.44)c(x) c˜(y) (2.44)Y (x) X(y) (2.45)Y +(x) X+(y) (2.46)Y −(x) X−(y) (2.47)D(x) D˜(y) §2.4.1xi yi §2.4.1L M (2.45)y0 x0 (3.4)mr m˜r (3.11)ymax xmax (3.11)mi m˜i (3.11)f f˜ (3.11)h h˜ (3.11)t, u t˜, u˜ (3.16)Y0(x) X0(y) §4.1.1Y1(x) X1(y) §4.1.1Object c.f.piri,j , pii,j , pi′0,j , etc. (2.15)pi(x, y) (2.31)Q(x, y) (2.31)pi1(x) (2.31)q1(x, y) (2.31)pi2(y) (2.31)q2(x, y) (2.31)pi0 (2.31)q0(x, y) (2.31)From here, the path forward forks. Both chapters 3 and 4 use the setupthat appears here, but they proceed independently of one another. A smalldiscussion of the overlap (or the current triviality thereof) will appear in theconclusion, after the work has been discussed.20Chapter 3Special ContoursIn this section we will examine special cases of contours as used in theBVP analysis of queueing systems, wherein a simple reparameterization ofthe contour may be created. We analyze these special cases in the hopesthat they may be leveraged to find explicit solutions for the BVP integral.3.1 Circles3.1.1 The Case for a Circle with a Center on the OriginThis case has shown up in at least two papers [1, 6], where the contourfor the BVP is a circle. Also of interest is the characterization appearing in([4], 111), wherein it is shown that the contour is a circle when∣∣∣∣∣∣∣p1,−1 p1,0 p1,1p0,−1 p0,0 p0,1p−1,−1 p−1,0 p−1,1∣∣∣∣∣∣∣= 0. (3.1)In this chapter, we will focus on techniques that are extensible into othercontours.One characterization of a circle in the complex plane is {z : |z| = r},where r is a positive constant. We will try to find queueing BVP contoursthat are circles by setting|Y (x)| = r (3.2)for either branch.2Parts of this chapter were presented as a part of the 2014 CORS Annual GeneralConference, under the talk title ‘Special contours in queueing BVPs’.213.1. CirclesDeveloping this, we recall that the discriminant has negative value overthe line segment in question, so the branches of Y (x) are complex conju-gates of one another, seeing as a(x), b(x), and c(x) are polynomials with realcoefficients. This result allows the following computation for either branch,(with the notation Y ± being either branch of Y (x), it doesn’t matter which,)|Y ±(x)| =√Y +(x)Y −(x)=√−b(x) +√b2(x)− 4a(x)c(x)2a(x)·−b(x)−√b2(x)− 4a(x)c(x)2a(x)=√b2(x)− b2(x) + 4a(x)c(x)4a2(x)=√c(x)a(x).(3.3)Since equation (3.3) is just an identity coming from the form of Y , itholds regardless of whether or not the contour is a circle. So c(x)/a(x) is aconstant iff L is a circle with a center at the origin.3.1.2 The Case for a Circle with a Center on the Real axisOur first generalization is to account for circles that do not have a centerat the origin. In our application, the only restriction on where the contourmust lie comes from the fact that the branches of Y (x) are conjugates of oneanother- as the top part and the bottom part of the contour are the samedistance away from the real axis, if the contour is a circle, the center mustlie on the real axis.Moreover in this case, this center can easily be determined by takingY (x1) and Y (x2) to find the ends of the contour and settingy0 =Y (x2) + Y (x1)2. (3.4)To characterize a contour that is a circle with center y0 in the complexplane, we can set|Y ±(x)− y0| = r. (3.5)223.1. CirclesLemma. Independent of the shape of the contour,∣∣Y ±(x)− y0∣∣2 =c(x) + y0b(x)a(x)+ y20for real y0.We may expand |Y ±(x)− y0|:∣∣Y ±(x)− y0∣∣2 = (Y +(x)− y0)(Y−(x)− y0)= Y +(x)Y −(x)− y0Y−(x)− y0Y+(x) + y20,(3.6)since y0 = y0. Now, we can rewrite the right hand side of equation (3.6) asfollows:Y +(x)Y −(x)︸ ︷︷ ︸See ((3.3))−y0(Y−(x) + Y +(x)︸ ︷︷ ︸Real part of 2Y) + y20, (3.7)so equation (3.7) is equal to∣∣Y ±(x)− y0∣∣2 =c(x)a(x)− 2y0−b(x)2a(x)+ y20, (3.8)which simplifies to the statement of the lemma.In summary, if the contour is a circle with center y0 and radius r, theny0b(x) + c(x)a(x)= r2 − y02. (3.9)If the contour is a circle, equation (3.9) will hold as designed. In the interestof establishing equivalence, say we start with equation (3.9) assigning arbi-trary real y0 and r. That would imply that the circle characterization wouldhold for that y0 and |r|, the squaring step being reversible because we knowwhat sign a modulus has to take. So equation (3.9) holds for any y0 and riff the contour is a circle with those parameters.233.2. Ellipses3.2 EllipsesThe characterization of a circle in the previous section (equation (3.5))generalizes to a very similar characterization with ellipses, so we will extendour work to include this case. The locus characterization of an ellipse is thatfor any point y on the ellipse, the sum of the lengths of the lines from y to thefoci f1 and f2 is a constant, h. In the complex plane, this characterizationis written as|y − f1|+ |y − f2| = h.As before, the symmetry along the real axis strongly restricts which ellipsesmight appear in our application. There are two cases:− The major axis of the ellipse is on the real line, and we write the centerof the ellipse as y0, and the foci as y0 + f and y0− f , for a real y0 anda real f .− The minor axis of the ellipse is on the real line, and we write the centerof the ellipse as y0, and the foci as y0 + f and y0− f , for a real y0 anda strictly imaginary f .In both of these cases we can write the equation corresponding to the ellipseas follows:|y − y0 − f |+ |y − y0 + f | = h. (3.10)3.2.1 Determining ParametersIn order for the tests we develop to be a usable tool, we need to beable to determine what y0, f , and h should be before we start applyingthe equations in the next sections 3.2.2 and 3.2.3. Fortunately, this is easyenough, after computing the length of the axes of the ellipse, noting herethat mr is the length of the axis of the ellipse parallel to the real axis, and243.2. Ellipsesmi is the length of the axis of the ellipse parallel to the imaginary axis:y0 =Y +(x2) + Y +(x1)2,mr = Y+(x2)− Y+(x1),ymax = Y+(z) where[ddx√4a(x)c(x)− b2(x)2a(x)]x=z∈R= 0,mi = 2|Im(ymax)|,f =√m2r −m2i2,h = max(mr,mi).(3.11)Note that the case where mr = mi is in fact the degenerate case where the‘ellipse’ is a circle. The constant h defined here is twice what r would bein section 3.1, but otherwise what follows is identical to the stated result,allowing the results of the following sections to be applied freely withouttesting for the cases from section 3.1.ParameterizationsIf we have a contour L that is an ellipse of one of the forms admissibleto our application, with y0, f and h as defined as in equation (3.11), it canbe endowed with the parameterizationL ={mr2cos(t) +imi2sin(t) + y0 | t ∈ [0, 2pi]}. (3.12)The choice of notation for mr and mi allows equation (3.12) to work forboth cases. Sines and cosines behave extremely well in symbolic integrationand will hopefully allow non-formal, non-numerical solutions to be derivedfor specific queueing BVPs in future work.253.2. EllipsesFigure 3.1: The two cases for an elliptical contour.On the left, the major axis is on the real axis. On the right the minor axisis on the real axis.3.2.2 The Case for an Ellipse with a Major Axis on theReal AxisIn the first case, the case for an ellipse with a major axis on the realaxis, we can expand the focus equations as follows:|Y (x)− y0 − f |︸ ︷︷ ︸t+ |Y (x)− y0 + f |︸ ︷︷ ︸u= h (3.13)which expands to4t2u2 =(h2 − t2 − u2)2. (3.14)Using the substitution coming from the lemma above, we replace t2 and u2witht2 = |Y (x)− (y0 + f)|2 =(y0 + f)b(x) + c(x)a(x)+ (y0 + f)2 (3.15)263.2. Ellipsesandu2 = |Y (x)− (y0 − f)|2 =(y0 − f)b(x) + c(x)a(x)+ (y0 − f)2. (3.16)At this point, we can simply substitute equations (3.15) and (3.16) intoequation (3.14), and have a usable test for computationally determiningwhether a given contour is an ellipse. In the interest of having a final resultof a similar form to (3.9), we can simplify the result of said substitution to−h2a(x)b(x)y0 + 4 f2a(x)b(x)y0 − h2a(x)c(x) + f2b2(x)a2(x)= −1/4(h2 − 4 y02) (h2 − 4 f2).(3.17)This final formula allows for easy checking whether this case of the ellipseholds, as all the contents of the equation can be quickly determined usingthe definitions in equation (3.11).3.2.3 The Case for an Ellipse with a Minor Axis on theReal AxisFor this case, we proceed mostly as in the other, but without the luxuryof having real foci. Instead of being able to reuse our previous work directly,we have, for a general case of a complex z,|Y ±(x)− z|2 = (Y ±(x)− z)(Y ±(x)− z)= Y ±(x)Y ±(x)− Y ±(x)z − Y ±(x)z + zz.(3.18)What equation (3.18) means for this problem, is that there are two casesthat need to be treated, depending on the choice of branch we use, as thechoice of branch determines whether we are on the top half of the ellipseor the bottom. This division into cases makes sense, as the foci are not onthe real line, so for each individual focus, the modulus may vary based onchoice of branch, but as we have two conjugate foci in our equation and twoconjugate branches, this division into cases should cancel out in the end.The necessary arithmetic in these cases is analogous to the derivation of273.2. Ellipsesequation (3.8). We obtain, for z = w + iv, (w, v ∈ R, )|Y +(x)− (w + iv)|2 = −v√4 a(x)c(x)− b2(x)a(x)+ w2 +wb(x)a(x)+ v2 +c(x)a(x),|Y −(x)− (w + iv)|2 =v√4 a(x)c(x)− b2(x)a(x)+ w2 +wb(x)a(x)+ v2 +c(x)a(x).(3.19)Finally, we can substitute w = y0, iv = f into the above equations (3.19)and substitute each one in turn into equation (3.14), noting that the branch-dependent terms containingv√4 a(x)c(x)− b2(x)a(x)cancel as expected. Fi-nally, we isolate constants as in section 3.2.2 to get our analogous result toequations (3.9) and (3.17):h2y0a(x)b(x) + h2a(x)c(x) + 4f2a(x)c(x)− f2b2(x)a2(x)=h24(h2 + 4 f2 − 4y02) .(3.20)Combined with our work in the previous section, the two formulae (3.17)and (3.20) exhaust all possible cases in our application for an elliptical con-tour.Important NoteAll of the work in this chapter so far has been operating from the pathstarting from the factorization Q(x, y) = a(x)y2 + b(x)y + c(x), and all ofthe tests give information on the contour L. If we want information on thecontour M associated with the alternate factorization, we simply replace allof the terms that appear here with their counterparts that appear in table2. A Fast Informal Test for Circular Contours3.3 A Fast Informal Test for Circular ContoursThe degenerate cases corresponding to our work in section 3.1, where thecontour associated with the BVP approach to a random walk is a circle, havean interesting characteristic for the jump probabilities for that random walk.The arrow representation for the jump probabilities in a random walk in thequarter plane (similar to figure 2.2, but with arrows corresponding to jumpprobabilities of zero not drawn,) is often used in discussion and informalwork on the random walk. In this section we will present two weak, fast,informal tests that can be used on these diagrams.Observing closely the definitions for a(x), b(x), and c(x),Q(x, y) = xy∑i,j(pijxiyj)− 1= p−1,−1 + xp0,−1 + yp−1,0 + xy(p0,0 − 1) + x2p1,−1+ y2p−1,1 + x2yp1,0 + xy2p0,1 + x2y2p1,1,which is factorable to be= (x2p1,1 + xp0,1 + p−1,1)y2 + (x2p1,0 + x(p0,0 − 1) + p−1,0)y+ (x2p1,−1 + xp0,−1 + p−1,−1),wherein we name the coefficients= a(x)y2 + b(x)y + c(x),(3.21)we may note that the condition given in equation (3.3) is equivalent, bymatching like terms, to the vector equation, for k = r2:kp1,1p0,1p−1,1 =p1,−1p0,−1p−1,−1 . (3.22)293.3. A Fast Informal Test for Circular ContoursEquation (3.22) means that we can automatically disqualify a subset ofrandom walks from having a circular contour centered at the origin. Thevalue of k is unimportant for this weak informal test, the only relevantinformation is what is available from the diagram, that is, whether the jumpprobabilities are zero or nonzero. The test is described in figure 3.2.Figure 3.2: The test for a circle with a center at the origin.If for each j, p1,j = 0 and p−1,j = 0, or p1,j 6= 0 and p−1,j 6= 0, then it ispossible that the contour is a circle with a center at the origin. Otherwiseit is impossible.What this means for usage is that arrows must come in pairs reflected acrossthe horizontal axis for L to be a circle. A similar result holds for the verticalaxis and M . If one test fails, then actually computing that vector equationis unnecessary.Equation (3.22) and the related test in figure 3.2 have results that co-incide with the determinant characterization in equation (3.1), as matriceswith linearly dependent rows or columns have a determinant of 0. This trickcan be extended somewhat to the circles with moved centers case. We havekp1,1 + y0p1,0p0,1 + y0(p0,0 − 1)p−1,1 + y0p−1,0 =p1,−1p0,−1p−1,−1 , (3.23)which has the equivalent test depicted in figure 3.3.The tests presented in figures 3.2 and 3.3 are only fast visual checksone can make to quickly obtain negative results, that is to say, to quickly303.3. A Fast Informal Test for Circular ContoursFigure 3.3: The test for a general circular contour.Formally, the test states ∀j ∈ {1,−1}, (p1,j = 0 ∧ p0,j = 0 ∧ p−1,j 6=0) =⇒ the contour is not a circle. A similar implication is derived by di-viding through by k.What this means for usage is to look for adjacent pairs of arrows not drawn(probability zero) along the edge. If such a pair exists without the third onein the line also not being drawn, then one of the contours cannot be a circle,so actual computation of the related vector equation is unnecessary.determine that a given random walk can not have a circular contour whenone approaches it using BVP methods.Please also note that the reflection of these tests along the diagonalalso gives a valid test: the tests pictured in figures 3.2 and 3.3 are basedon the equations (3.3) and (3.9), which have counterparts r =√c˜(y)a˜(y) andr2−x02 =c˜(y)+x0b˜(y)a˜(y) , which correspond to vector equations that correspondto the transposed tests.The cases where one wishes to determine if a contour is an ellipse arecomplicated enough that one is better off seeing if formulae (3.17) and (3.20)hold than attempting to use a test analogous to those that appear in thissection.31Chapter 4Homogeneous QueueingBVPsIn this chapter, we extend a known approach for obtaining the generat-ing function for the stationary probability to be compatible with a broadclass of F.E.s known as the homogeneous case, which is described shortly.The approach we take follows Avrachenkov’s in [1], wherein we attack theproblem with a Riemann-Hilbert BVP. Our extension to the mathematicsof the existing literature consists of an algorithmic method for determiningthe special branches of Y and a proof extending a required result in [1] tocover the arbitrary contours that may appear in the homogeneous case. Wealso hope to make it clear and explicit exactly what has to be done to fitthis approach to a given fundamental equation.The class of problems that this approach is designed to address is theclass of random walks in the quarter plane which satisfies the conditionsof section 2.2. From this class, we will take a representative with a F.E.constructed arbitrarily. Additionally, this approach also requires that inequation (2.23), q0(x, y) ≡ 0; This generating function being zero is why wecall this class of F.E.s the homogeneous case.Important noteThe construction of section 2.3 is effectively incompatible with the ho-mogeneous case that appears here. Other constructions, (like [1],) can resultin a homogeneous F.E., but observing closely the definitions applied in sec-tion 2.3, setting q0(x, y) ≡ 0 in the setup above implies that the Markovchain can never leave the empty state, which, as one would imagine, does324.1. General Theory for Homogeneous Fundamental Equationsnot require much labor for one to determine a stationary distribution.We will proceed under the assumption that the starting point of ouranalysis is the homogeneous F.E. itself, composed of generating functionsanalytic over the unit disc, that fits the form of equation (2.23), and proceedto make definitions as appear as in section 2.4.Chapter overviewOur goal is to manipulate the F.E. into a form amenable to treatmentby a Riemann-Hilbert BVP. We take measures to satisfy the requirementson the boundary condition, but this requires some setup. First, we willdetail the construction and provide an example construction of additionalbranches of Y (x) and X(y) labeled X0(y) and Y0(x). These are created firstto allow easy manipulation of the F.E., and second to have useful proper-ties that allow us to make analytic continuations up to the boundary, asrequired by the BVP setup. The details of these necessary properties aretreated immediately after the construction of the branches. Second, we dothe actual manipulation into the form of a Riemann-Hilbert BVP. Afterthe aforementioned analytic continuation, and modification to create theboundary condition, we eliminate one of the unknown piis and and coercethe F.E. into the BVP form in one step. Third, we meet the rest of theBVP requirements by verifying analyticity, and moving around roots, anduse the BVP as intended. Finally, we recover pi1 and pi2, then make a briefargument allowing us to recover pi(x, y), and we are done.4.1 General Theory for HomogeneousFundamental Equations4.1.1 Construction of New BranchesIn table 2.4.2, we documented how the existence of a different factoriza-tion at the beginning induced counterpart definitions, mirroring our work.These bodies of work are more than similarly constructed, in fact they areclosely linked through a relationship between the branches ofX(y) and Y (x).334.1. General Theory for Homogeneous Fundamental EquationsAccording to [4], there exists a choice of branches of these functions, callthem X0(y) and Y0(x) such that X0 ◦ Y0 is the identity map on GM\x1x2,and Y0◦X0 is the identity map on GL\y1y2. These functions are also analyticeverywhere except the cut between the first and second branch points, andthe cut between the third and fourth branch points. Additionally, recallequation (2.49), which says we can construct L and M using these branches.Constructing these new branches allows us to easily move from the x-planeto the y-plane and vice versa, as they also have several useful properties thatwill make our analysis easier later on.Constructing these branches explicitly has mainly been glossed over inthe literature, as their existence is well proven, but we will discuss it forthe sake of completeness here. We are primarily interested in the curveIm(D(x)) = 0, as this is where the choice of the branch of the square rootmatters. As long as we make changes of branch only here, we do not riskbreaking analyticity anywhere else, changes of branch here cannot worsenthings.We consider the partitions generated by looking at the signs of Re(D(x))and Im(D(x)) to construct possible branches for Y0(x). Once all of thesecases are enumerated, we can thin them out by looking at which ones areanalytic on the real line minus the cuts, e.g. (−∞, y1) ∪ (y2, y3) ∪ (y4,∞).We do the same with their image under D˜(y), and test to see which pairingresults in a pair of inverse functions, which will result in our choice of branchfor X0 and Y0. The ‘opposite’ choice of branch on each partition in eachcase will be meromorphic except on the cuts as well, and will be declared asX1 and Y1 respectively. (See [4].)Example Construction of Analytic BranchesThe explicit formulation of these branches comes directly from [8]. Theauthors also state that all of these branches are analytic in the cut plane,except X1, which has a pole at 0, and is thus meromorphic over the cutplane. These branches off of the branches Y +, Y −, X+, and X− as follows:344.1. General Theory for Homogeneous Fundamental EquationsX0(y) =X+(y) if y ∈{z∣∣∣∣∣Re(z) ≤ y2,Im(D˜(Re(z) + i|Im(z)|)) < 0},X+(y) if y ∈ (−∞, y1),X−(y) otherwise.(4.1)X1(y) =X−(y) if y ∈{z∣∣∣∣∣Re(z) ≤ y2,Im(D˜(Re(z) + i|Im(z)|)) < 0},X−(y) if y ∈ (−∞, y1),X+(y) otherwise.(4.2)Y0(x) =Y +(x) if x ∈{z∣∣∣∣∣Re(z) ≤ x2Im(D(Re(z) + i|Im(z)|)) < 0},Y +(x) if x ∈{z∣∣∣∣∣Re(z) ≥ x3,Im(D(Re(z) + i|Im(z)|)) > 0},Y +(x) if x ∈ (−∞, x1) ∪ (x4,∞),Y −(x) otherwise.(4.3)Y1(x) =Y −(x) if x ∈{z∣∣∣∣∣Re(z) ≤ x2Im(D(Re(z) + i|Im(z)|)) < 0},Y −(x) if x ∈{z∣∣∣∣∣Re(z) ≥ x3,Im(D(Re(z) + i|Im(z)|)) > 0},Y −(x) if x ∈ (−∞, x1) ∪ (x4,∞),Y +(x) otherwise.(4.4)4.1.2 Bounds for X0(y) and Y0(x)We will also require several facts about the behavior of these branchesfor our future use, bounds to help prove analyticity when the functions arecomposed. Here we introduce the notation Γ to represent the unit circle,and hence GΓ to represent the unit disc, abusing notation so that Γ and GΓ354.1. General Theory for Homogeneous Fundamental Equationsare part of the x-plane or y-plane as needed.We state our bounds-related results here:1. |X0(y)| ≤ 1 if |y| = 1, with equality only possible at 1.2. X0(y) ∈ GM ∪GMext ,∀y ∈ C\y1y2.3. GL\x1x2 and GM\y1y2 are homotopically equivalent via X0 and Y0.4. |X0(y)| ≤ 1 for y ∈ GL\GΓ.Although the first three points here have backing in the literature3 [4,pp 109, 118, 119], the general case of point number four that we need todayhas only been treated in specific examples [1, 3]. We present a proof of pointnumber four for our general case here.Proof. We want to show that |X0(y)| ≤ 1 for y ∈ GL\GΓ. The overarch-ing proof structure is as follows: We show that the desired property holdsover the image Y0(GΓ ∩ GM\x1x2), and then demonstrate that that imagecontains GL\GΓ. Effectively, we are illustrating how X0 and Y0 map someregions to others.− Lemma 1We demonstrate that |X0(s)| < 1,∀s ∈ Y0(GΓ ∩ GM\x1x2). First,we know that X0 ◦ Y0(t) = t,∀t ∈ GM\x1x2. Second, we take anelement e ∈ GΓ ∩ GM\x1x2. Then X0(Y0(e)) = e ∈ GΓ ∩ GM\x1x2.Note that since e ∈ GΓ, |X0(Y0(e))| ≤ 1. Finally, we identify s =Y0(e) ∈ Y0(GΓ ∩GM\x1x2). This results the range of possible s beingidentically Y0(GΓ ∩GM\x1x2) and so the lemma holds.− Lemma 2We show that the image Y0(GΓ ∩ GM ) contains GL\GΓ. First somebackground facts: Since we know that |Y0(x)| ≤ 1 if |x| = 1, withequality only possible at 1, the contour Y0(Γ) ⊂ GΓ. Also, according3Point three is not explicitly declared as it appears here, but it meets the definition.364.1. General Theory for Homogeneous Fundamental Equationsto the definition of L in equation (2.49), it is the contour coming fromthe limit of contours approaching the line segment x1x2 under Y0.Now we proceed by a topological argument. Our definitions for ho-mology groups and homotopic equivalences come from [14]. DenoteA = GΓ ∩ GM\y1y2, and B = GL\GY0(bdy GΓ∩GM ). Our immediatelypreceding points indicate that B contains GL\GΓ.The boundaries of A map to the boundaries of B, we show that theimage of A under Y0 must be exactly B:We need to work with closed sets for the time being, so we create eitheran open annulus or the open unit disk to be homotopically equivalentto A, and take the closure of whichever was needed to be A. We alsodenote the closures of Y0(A) and B as A′ and B respectively. Finally,we extend the homotopic equivalence of A and Y0(A) through Y0 tothese new sets, denoting the extended map Y, which maps boundariesto boundaries by taking limit points to limit points. A commutativediagram appears below depicting this construction, with the arrowsbetween Y0(A) and A′, and also between B and B denoting inclusion.Our current goal is to show that A′ = B.A ooYY−1// A′ ? BAOOooY0X0// Y0(A)OOBOOFor the case where A is an annulus, consider a path p from one bound-ary of A to the other entirely within the interior of A except at theendpoints. Y(p) leaving and returning to A′ would imply that theimage of a non-boundary point was the same as the image of a pointon the boundary, contradicting the fact that Y is injective. The casewhere A is a disc admits a similar argument and result for images ofpaths inside the disc mapping to paths entering and leaving the disc.So B ⊇ A′.374.1. General Theory for Homogeneous Fundamental EquationsThe homology groups of A are the same as those of A′ due to thehomotopic equivalence. Further, the boundaries B are part of A′.Therefore, it is impossible for any point in B to not be part of A′as this would put in holes different from the ones that come from theboundaries of A and change the first homology group. So B ⊇ A′ ⊇ B.Hence, int A′ = int B, and therefore Y0(A) = B ⊇ GL\GΓ.The combination of these two lemmas gives the proof as desired.4.1.3 Solution via BVPBoundary Value Problem SetupOur goal at this point is to manipulate our F.E. into a BVP, usingQ(x, y) = 0 (equivalently, y = Y (x), either branch,) to reduce the complex-ity and set the values of y on L simultaneously. This gives us our startingpoint, with restrictions based on the intersections of the regions of analyt-icity of the functions:q1(x, Y0(x))pi1(x) + q2(x, Y0(x))pi2(Y0(x)) = 0, x ∈ GΓ\x1x2. (4.5)If we rewrite equation (4.5) asq1(x, Y0(x))pi1(x) = −q2(x, Y0(x))pi2(Y0(x)), x ∈ GΓ\x1x2, (4.6)we find that we can analytically continue the left hand side out to M . Thereason that the RHS is analytic is that problems can only occur when thearguments of the functions pii leave the unit circle, as we defined them tobe generating functions (2.31), which only have guaranteed convergence onthe unit circle. To extend this, we note that Y0’s range is the unit disc giventhat we take arguments in GM\GΓ, by our work in section (4.1.1). Theterm that appears on the RHS is the composition pi2 ◦ Y0, and so pi1 can beanalytically continued out to M , leaving it with a region of analyticity ofGM .The interval [x1, x2], and its relevance to the definition of L is all that is384.1. General Theory for Homogeneous Fundamental Equationsof interest for our BVP setup, so we ignore most of the region of analyticityin favor of the following:q1(x, Y0(x))pi1(x) = −q2(x, Y0(x))pi2(Y0(x)), as x approaches [x1, x2].(4.7)We now apply our next bit of information from (4.1.1), that Y0 ◦X0 isthe identity function. This allows us to switch to the other way of trackingthat Q(x, y) = 0, by setting x = X0(y). This change of variable applied toequation (4.7) results in a condition over L:q1(X0(y), y)pi1(X0(y)) = −q2(X0(y), y)pi2(y), y ∈ L. (4.8)The next thing we do is observe the zeroes of q1(X0(y), y). We want toisolate pi1(X0(y)) via dividing through by the coefficient. pi1(X0(y)) doesn’thave any zeroes on L, so we are free to divide as we like, resulting in ournext step:q2(X0(y), y)q1(X0(y), y)pi2(y) = −pi1(X0(y)), y ∈ L. (4.9)When y ∈ L,X0(y) ∈ [x1, x2], so −pi1(X0(y)) ∈ R. pi1 here is a generatingfunction for a real data set, so putting real values into pi1 returns real values.We may now proceed to our next step by exploiting this fact, eliminatingpi1(X0(y)):Im(q2(X0(y), y)q1(X0(y), y)pi2(y))=12[q2(X0(y), y)q1(X0(y), y)pi2(y)−q2(X0(y), y)q1(X0(y), y)pi2(y)]= 0, y ∈ L.(4.10)A small amount of algebraic manipulation of (4.10) results in−q2(X0(y),y)q1(X0(y),y)(q2(X0(y),y)q1(X0(y),y))pi2(y) = pi2(y), y ∈ L. (4.11)Now, equation (4.11) would fit a BVP perfectly, if the coefficient of pi2(y)394.1. General Theory for Homogeneous Fundamental Equationswere nonvanishing and pi2 itself were analytic inside, and continuous up tothe boundary of GL. We will show that this is in fact the case.We handle the analyticity first. We make the argument via analyticcontinuation on equation (4.9), specifically the right hand side. We knowthat pi1 is analytic on the unit disc, and that X0 is analytic inside GL\x1x2,and continuous up to the boundary, and maps the same to the unit disc.Hence the composition is analytic inside GL\x1x2, as desired.The nonvanishing part takes a bit more work: Beginning with equation(4.9), we note that if q2(X0(y), y) had any zeroes, say wi, in GL, i ∈ I, (I ={1, 2, . . . , n}, ) wherein multiple roots correspond to multiple wi, then wecould define the following functions to move the roots into pi2.U∗(y) =q2(X0(y), y)q1(X0(y), y)∏i∈I(y − wi)P (y) = pi2(y)∏i∈I(y − wi). (4.12)We can now claim thatU∗(y)P (y) =q2(X0(y), y)q1(X0(y), y)pi2(y), (4.13)allowing us to claim that the BVP setup is now complete for the boundaryconditionU(y)P (y) = pi2(y), y ∈ L (4.14)where we defineU(y) = −U∗(y)(U∗(y)) . (4.15)Boundary Value Problem SolutionNow that we have a BVP, we can begin to specify a solution. First,however, we need to calculate the index of the BVP. Attentive readers mayhave wondered why we took all the zeroes in GL rather than just on L, asrequired for the BVP setup. Doing it the way we did allows us to greatlysimplify our index calculation. Moving all of the zeroes elsewhere, andtherefore not having any roots in GL ensured that the index is zero, by the404.1. General Theory for Homogeneous Fundamental EquationsCauchy argument principle.Thus the BVP solution, using the Riemann-Hilbert BVP setup, is asfollows:P (y) = exp(12pii∫Llog(U(t))dtt− y)(4.16)Equation (4.16) is a formal solution, it is rare that this integral can beevaluated analytically.We may now use this to solve for the generating functions that remainunknown. From the BVP solution in equation (4.16) we can easily extractpi2(y), aspi2(y) =1∏i∈I(y − yi)exp(12pii∫Llog(U(t))dtt− y), (4.17)and by a quick reference to table 2.4.2, we can determine pi1(x) by the sameapproach. Equation (4.17), with a quick reference to equation (2.23), allowsus to determine pi(x, y), aspi(x, y) =q1(x, y)pi1(x) + q2(x, y)pi2(y)Q(x, y). (4.18)At this point the pii,j may be extracted by taking derivatives, as is usual forgenerating functions.41Chapter 5ConclusionIt is our hope that our work in developing tests for identifying elegantlyreparameterizable contours will be useful in working towards analytic solu-tions of the associated Riemann-Hilbert BVPs for this subset of problems,and that our efforts in creating a formal general solution that extends ex-isting methods for approaching homogeneous F.E.s associated with randomwalks in the quarter plane will be of use for students and queueing theo-rists addressing those F.E.s in the future. Non-homogeneous BVPs are stillbest attacked with the methods found in chapters 4, 5, and 6 of the text byFayolle et al. on the subject [4].The bodies of our work in chapters 3 and 4 are developed independently.Worse yet, they are fundamentally incompatible, as chapter 3 relies on aconstruction that renders the work in chapter 4 trivial. Even if a solutionto the incompatibility were found, there would be a poor return on theinvestment of going through the work in chapter 3 in addition to that ofchapter 4, as the results of chapter 4 stop at the formal solution. Theaforementioned analytic solutions that we hope to use the results of chapter3 to find remain a goal rather than a reality at the time being. In any case,we would still be restricting the general class of F.E.s that chapter 4 treatsto a limited subset.There are several avenues that future work may take to extend whatappears here. Our approach for identifying special contours could be ex-tended to any locus-characterizable contour, but it is certain that more workneeds to be done before our approach is a usable tool for queueing theorists.Integration over these simpler contours should be easier than the naturalparameterization resulting from drawing values from [x1, x2] in equation(2.48), but a general solution for the integrals that occur in the solution of42Chapter 5. ConclusionRiemann-Hilbert BVPs is yet to be found. It has been suggested that look-ing at the non-BVP solution found in [5] alongside the formal BVP solutionfound here could give hints on how to compute this integral analytically,and others like it.43Bibliography[1] K. Avrachenkov, P. Nain, and U. Yechiali, A retrial system withtwo input streams and two orbit queues, Queueing Systems, (2013),pp. 1–31.[2] J. Cohen and O. Boxma, Boundary Value Problems in QueueingSystem Analysis, North-Holland Mathematics Studies, North-HollandPublishing Company, 1983.[3] G. Fayolle and R. Iasnogorodski, Two coupled processors: The re-duction to a riemann-hilbert problem, Zeitschrift fr Wahrscheinlichkeit-stheorie und Verwandte Gebiete, 47 (1979), pp. 325–351.[4] G. Fayolle, R. Iasnogorodski, and V. A. Malyshev, RandomWalks in the Quarter-Plane: Algebraic Methods, Boundary Value Prob-lems and Applications, Springer, 1999.[5] L. Flatto, Two parallel queues created by arrivals with two demandsII, SIAM Journal on Applied Mathematics, 45 (1985), pp. 861–878.[6] L. Flatto and S. Hahn, Two parallel queues created by arrivals withtwo demands i, SIAM Journal on Applied Mathematics, 44 (1984),pp. 1041–1053.[7] F. Guillemin, C. Knessl, and J. S. Van Leeuwaarden, Wirelessmultihop networks with stealing: Large buffer asymptotics via the raymethod, SIAM Journal on Applied Mathematics, 71 (2011), pp. 1220–1240.44Bibliography[8] F. Guillemin, C. Knessl, and J. S. van Leeuwaarden, Wirelessthree-hop networks with stealing II: exact solutions through boundaryvalue problems, Queueing Systems, 74 (2013), pp. 235–272.[9] F. Guillemin, C. Knessl, and J. S. H. van Leeuwaarden, Firstresponse to letter of g. fayolle and r. iasnogorodski, Queueing Systems,76 (2014), pp. 109–110.[10] H. Li and Y. Q. Zhao, Tail asymptotics for a generalized two-demandqueueing model - a kernel method, Queueing Systems, 69 (2011), pp. 77–100.[11] B. P. Palka, An Introduction to Complex Function Theory, Under-graduate Texts in Mathematics, Springer-Verlag, 1991.[12] J. Resing and L. O¨rmeci, A tandem queueing model with coupledprocessors, Operations Research Letters, 31 (2003), pp. 383–389.[13] S. Ross, Introduction to Probability Models, Academic Press, 2010.[14] J. Vick, Homology Theory: an Introduction to Algebraic Topology,Graduate Texts in Mathematics, Springer-Verlag, 1994.[15] P. E. Wright, Two parallel processors with coupled inputs, Advancesin applied probability, (1992), pp. 986–1007.45


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items