# Performance Analysis of Production Systems with Correlated Demand via Diffusion Approximations

Authors: Yingdong Lu
Journal: International Journal of Stochastic Analysis
 (2012)

---

## Abstract

We investigate the performance of a production system with correlated demand through diffusion approximation. The key performance metric under consideration is the extreme points that this system can reach. This problem is mapped to a problem of characterizing the joint probability density of a two-dimensional Brownian motion and its coordinate running maximum. To achieve this goal, we obtain the stationary distribution of a reflected Brownian motion within the positive quarter-plane, which is of independent interest, through investigating a solution of an extended Helmhotz equation.


---

## Body

## 1. Introduction

There are extensive studies on some classic one-dimensional models in the field of operations research, most notably, the work on the behavior and scheduling of single-server queuing system, as well as the understanding of performance and management of single-item inventory system. Probabilistic tools and techniques such as random walks and integral transformations are traditionally used in analyzing them. More recently, concepts and methodologies in dynamical systems and diffusion processes are brought in through fluid and diffusion approximations, and they extend our understanding and capability of analyzing and control of these systems substantially.While highly desirable, extending these studies to their multidimensional counterparts is a rather difficult problem. Extensions of classic probabilistic methods and techniques, such as random walk and integral transform, introduce a new level of complication that requires deeper understanding in algebraic and complex geometry for their general treatment. For the approximations methods, their multidimensional counterparts, multidimensional dynamical systems, and multidimensional diffusion processes, also pose significant barriers for either qualitative understanding and quantitative computation. While a general and sophisticated methodology is lacking, some very interesting results have been obtained in establishing convergence results that lay the foundation to both fluid and diffusion approximation schemes. Meanwhile, the quantitative aspect of the problem seriously lagged behind, and results in approximation and control of multidimensional systems are very limited. This calls for more focus be put on computational efforts on such system, so that new tools and techniques can be added to the arsenal of attacking these problem.In this paper, we aim at extend some theoretical and computational understanding to a simple but representative two-dimensional operations research model. Not only can it serve as a building block for analyzing much more complicated systems, but also provides insights to the understanding of basic structure of general systems with correlated demand. This model is motivated by several typical applications in production systems and scheduling. In such a system, it is very common that one resource is demanded by multiple demand processes, similarly, one class of customer demand could contain a combination of different products. The correlation induced by this commonality poses difficulty even for the simplest version of the model. In this paper, after briefly introducing the mathematical model and stochastic processes that characterize the basic underlying relationship, we provide an diffusion approximation to the stochastic processes in the form of a reflected two-dimensional Brownian motion, which is just an abstraction of previous known results in various forms. Then our main effort will be concentrated on the analysis of some key quantities of this diffusion process. More specifically, we are interested at the joint probability density of the Brownian motion and its coordinate running maximum. To produce the computational result, we relate its computation to finding a proper solution to a classic elliptic partial differential equation, the extended Helmhotz equation.The rest of the paper will be organized as follows, in Section2, we will introduce the mathematical model and some preliminary results; in Section 3 we will present the a key Brownian bridge argument; then, in Section 4, we solve an extended Helmhotz equation related to the invariant measure; in Section 5, we present the final calculation results for the probability density, and we conclude the paper with come comments on related future works in Section 6.


## 2. Model

Suppose that there are two types of products, type1 and type 2. Each of them is supplied by an independent renewal process. Let us denote them as A1(t) and A2(t), respectively. Statistically, we only require that the interarrival time follows a distribution with finite second moments. Meanwhile, there are two classes of demand, a class-(1) demand requires one unit of type 1 product, a class-(2) demand requires one unit of type 1 product, and a class-(1,2) demand requires one unit of each product 1 and 2. Note that class-(1,2) demand can only be fulfilled when both two types of products are provided. If any of the product requirement can not be satisfied, the demand will be backlogged. Although the model is quite simplistic on the surface, it embodies the basic structure and trade-off for many widely considered systems, for example, assemble-to-order systems in inventory management and queuing systems with flexibility servers. Results on this problem can have many indications to the understanding of those more complicated systems.Diffusion approximation can be obtained for this system, see, for example [1, 2], the limiting process is a two-dimensional Brownian motion with correlation. One of the most important problems for characterizing the process is the probabilistic distribution of the coordinate running maximum. For many applications, the coordinate running maximums are closely related to the capacity provisioning since they represent how far the process will reach at different directions. For probability theory, it is also a fundamental problem. And this will be the focus of the paper.From now on, we will focus on the calculation of the probability density function calculation. More specifically, let(X1(t),X2(t)) be a two-dimensional Brownian motion with constant correlation ρ>0 and initial state (X1(0),X2(0))=(0,0). Therefore, we have E[Xi(t)]=0, E[Xi2(t)]=t, E[Xi(s)Xi(t)]=min{s,t} for i=1,2, and E[X1(s)X2(t)]=ρmin{s,t}. For detailed properties of such process, see, for example, [3]. A quantity that of general interest is, for z≥x and w≥y and any fixed time t0>0(2.1)P[X1(t0)∈dx,M1(t0)∈dz;X2(t0)∈dy,M2(t0)∈dw],
where Mi(t) are the running maximum for Xi(t) for i=1,2 respectively, that is Mi(t)=sup0≤s≤t{Xi(s)}. It is easy to see that we can just focus on the case of t0=1, general case can be easily obtained through scaling. To be more specific, we have, for any t0,
(2.2)P[X1(t0)∈dx,M1(t0)∈dz;X2(t0)∈dy,M2(t0)∈dw]=P[X1(1)∈dxt0,M1(1)∈dzt0;X2(1)∈dyt0,M2(1)∈dwt0],
because (X1(t),X2(t)) and (t0X1(t/t0),t0X2(t/t0)) have the same law due to the scaling invariant property of Brownian motion.In the case of one-dimensional Brownian motion, the density of the Brownian motion and its running maximum is a classical result obtained through reflection principle. In two dimension, the dependence between the two coordinates creates a barrier for generalizing the result easily. In this paper, first, we convert the calculation of the quantity (2.1) to the calculation of density function of the invariant measure for a reflected two-dimensional Brownian motion in the positive orthant with negative drifts through a Brownian bridge argument; then we obtain a Laplace transform of the invariant measure of a reflected two-dimensional Brownian motion in the positive orthant with negative drifts through solving a partial differential equation.


## 3. A Brownian Bridge Argument

There are many different approaches to the classic problem of the joint distribution of a standard Brownian motion and its running maximum, see, for example, [4]. Here we introduce an approach connecting probability laws of a conditional Brownian motion to a Brownian bridge. To the best of our knowledge, no similar approaches appeared before in the literature. A standard Brownian motion, X(t), which starts at 0, conditioning on the event that it ends at x at time 1, has the same law as
(3.1)W(t)+xt,
where W(t) is a standard Brownian bridge, that is, a Brownian motion satisfies W(0)=W(1)=0. Therefore, from the representation of Brownian bridge, we have, for t∈[0,1],
(3.2)W(t)+xt=dtB(1t-1)+xt,
where B(t) is a standard Brownian motion, =d means equal in distribution. Thus, for any z≥x, the event of
(3.3){tB(1t-1)+xt≤z,∀t∈[0,1]},
will have the same probability as,
(3.4){B(1t-1)≤zt-x,∀t∈[0,1]}.
Apply a variable transform s=1/t-1, we get,
(3.5){B(s)≤zs+(z-x),∀s≥0},
or equivalently,
(3.6){B(s)-zs≤z-x,∀s≥0}.
This probability of the above event is equivalent to the stationary distribution of B(s)-zs reflected at zero, which, in turn, equals
(3.7)1-exp[-2z(z-x)].
Take derivative with respect to z, we have,
(3.8)2(2z-x)exp[-2z(z-x)].
So,
(3.9)P[X(1)∈dz,M(1)∈da]=P[M(1)∈da∣X(1)∈dz]=2(2z-x)2πexp[-(2z-x)22].
This is consistent with the classic results obtained from other approaches such as the reflection principle.For the two-dimensional Brownian motion, applying the same argument, calculating(3.10)P[X1(1)∈dx,M1(1)≤z,X2(1)∈dy,M2(1)≤w],
for z≥x and w≥y can be reduced to calculating
(3.11)P[Y1(t)-zt≤z-x,Y2(t)-wt≤w-y],∀t≥0.
Apply the construction of reflected Brownian motion through maximum of Brownian motion with negative drift, see, for example, [5], as well as the two-dimensional version of P. Lévy theorem, see [3], we can conclude that (3.11) is the same as
(3.12)P[Y^1(t)≤z-x,Y^2(t)≤w-y],∀t≥0,
where (Y^1(t),Y^2(t)) is a two-dimensional Brownian motion with negative drifts that reflected at the two positive axes. To be more precise,
(3.13)Y^1(t)=-zt+W1(t)+L1(t),Y^2(t)=-wt+W2(t)+L2(t),
where (W1,W2) is a two-dimensional driftless Brownian motion with covariance matrix,
(3.14)Γ=(1ρρ1),
and Li(t),i=1,2 are the local time defined as
(3.15)Li(t)≜sup0≤s≤t-[θit+Wi(t)]-,
where θ1=-z, θ2=-w. Therefore, following the same logic, to obtain (3.10), we need to compute the stationary distribution of (Y^1(t),Y^2(t)).


## 4. Stationary Distribution

To calculate the stationary distribution of the reflected Brownian motionY^1(t),Y^2(t), we make use the following result from Harrison and Williams [6].Theorem 4.1.
Functionsπ(x,y)∈C2(ℝ+2) and σi(x)∈C2(ℝ+1) are the density of the stationary distribution of X(t) on the orthant and the boundary, if and only if the following basic adjoint relationship (BAR) is satisfied:(4.1)∫ℝ+2ℒf(z)π(x,y)dxdy+∑j=12∫FjDjf(z)σi(z)dvj(z)=0,
where(4.2)ℒ=DtΓD+(θ1,θ2)tD
is the infinitesimal generator and vj(z) is the surface measures defined on the boundary such that the normal vector is pointed inward.Remark 4.2.
Notice that the boundary in the problem we study consists of the two axis. So the surface measure is basically the Lebesgue measure onℝ1 with the orientation that guarantees that the normal vector points towards the inside of the orthant.The main idea in this paper is to find a functionu(x,y) such that ℒu=e-zx-wy for any pair (w,z), therefore, the first term of (4.1) becomes the Laplace transformation of the stationary distribution, and then to investigate the basic properties of u(x,y) to compute the second term. To find u(x,y), we need to solve the following equation:
(4.3)∂2u∂x2+2ρ∂2u∂x∂y+∂2u∂y2+θ1∂u∂x+θ2∂u∂y=e-zx-wy.
Let us first consider the homogeneous version of the equation,
(4.4)∂2u∂x2+2ρ∂2u∂x∂y+∂2u∂y2+θ1∂u∂x+θ2∂u∂y=0.
Apply transform,
(4.5)ξ=x,η=-ρ1-ρ2x+11-ρ2y,
then, u(ξ,η) satisfies
(4.6)Δξηu+θ~1uξ+θ~2uη=0,
where θ~1=θ1 and θ~2=(θ1-ρθ2)/1-ρ2. Let v(ξ,η)=eθ~1ξ+θ~2ηu(ξ,η), then v(ξ,η) is the solution of the following problem,
(4.7)Δξηv-θ~1θ~2v=0.
This is so called extended Helmholtz equation. Separation of variables can be applied to obtain the general forms of the solution. Hence, we know that the general solution is
(4.8)v(ξ,η)=[c1e-kξ+c2ekξ][c3e-θ~1θ~2-kη+c4eθ~1θ~2-kη],
for any constant c1,c2,c3,c4, and k≤θ~12/4+θ~22/4. Therefore,
(4.9)u(x,y)=c1c3exp{[-(k+θ~1)+ρ(θ~1θ~2-k+θ~2)1-ρ2]x-θ~1θ~2-k+θ~21-ρ2y}+c1c4exp{[-(k+θ~1)-ρ(θ~1θ~2-k+θ~2)1-ρ2]x+θ~1θ~2-k-θ~21-ρ2y}+c2c3exp{[(k-θ~1)+ρ(θ~1θ~2-k+θ~2)1-ρ2]x-θ~1θ~2-k+θ~21-ρ2y}+c2c4exp{[(k-θ~1)-ρ(θ~1θ~2-k-θ~2)1-ρ2]x+θ~1θ~2-k-θ~21-ρ2y}.
For the inhomogeneous version, it is easy to see that u-(1/(z2+2ρzw+w2-θz-θw))e-zx-wy will be the solution.Now plugu(x,y) into the BAR (4.1), we have(4.10)π^(z,w)+c1c3{[-(k+θ~1)+ρ(θ~1θ~2-k+θ~2)1-ρ2]σ^(θ~1θ~2-k+θ~2/21-ρ2)-θ~1θ~2-k+θ~21-ρ2σ^(-(k+θ~1)+ρ(θ~1θ~2-k+θ~2)1-ρ2)}+c1c4{[-(k+θ~1)-ρ(θ~1θ~2-k-θ~2)1-ρ2]σ^(θ~1θ~2-k-θ~21-ρ2)-θ~1θ~2-k-θ~21-ρ2σ^(-(k+θ~1)-ρ(θ~1θ~2-k-θ~2)1-ρ2)}+c2c3{[(k-θ~1)+ρ(θ~1θ~2-k+θ~2)1-ρ2]σ^(θ~1θ~2-k+θ~21-ρ2)-θ~1θ~2-k+θ~21-ρ2σ^((k-θ~1)+ρ(θ~1θ~2-k+θ~2)1-ρ2)}+c2c4{[(k-θ~1)-ρ(θ~1θ~2-k-θ~2)1-ρ2]σ^(θ~1θ~2-k-θ~21-ρ2)-θ~1θ~2-k-θ~21-ρ2σ^((k-θ~1)-ρ(θ~1θ~2-k-θ~2)1-ρ2)}-σ^(z)-σ^(w)z2+2ρzw+w2-θ1z-θ2w=0,
where π^ and σ^ denote the Laplace transform of π(x,y) and σ(x). Given that ci, i=1,…,4 and k can be arbitrary, it is necessary that
(4.11)[-(k+θ~1)+ρ(θ~1θ~2-k+θ~2)1-ρ2]σ^(θ~1θ~2-k+θ~21-ρ2)-θ~1θ~2-k+θ~21-ρ2σ^(-(k+θ~1)+ρ(θ~1θ~2-k+θ~2)1-ρ2)=0,[-(k+θ~1)-ρ(θ~1θ~2-k-θ~2)1-ρ2]σ^(θ~1θ~2-k-θ~21-ρ2)-θ~1θ~2-k-θ~21-ρ2σ^(-(k+θ~1)-ρ(θ~1θ~2-k-θ~2)1-ρ2)=0,[(k-θ~1)+ρ(θ~1θ~2-k+θ~2)1-ρ2]σ^(θ~1θ~2-k+θ~21-ρ2)-θ~1θ~2-k+θ~21-ρ2σ^((k-θ~1)+ρ(θ~1θ~2-k+θ~2)1-ρ2)=0,[(k-θ~1)-ρ(θ~1θ~2-k-θ~2)1-ρ2]σ^(θ~1θ~2-k-θ~21-ρ2)-θ~1θ~2-k-θ~21-ρ2σ^((k-θ~1)-ρ(θ~1θ~2-k-θ~2)1-ρ2)=0.
In fact, they are all in the form of
(4.12)Aσ^(B)-Bσ^(A)=0,
and A and B satisfy the quadratic equation
(4.13)(A+ρB+θ~1)2+(1-ρ2B+θ~2)2=θ~12+θ~224,
which is equivalent to
(4.14)A2+2ρAB+B2+2θ1A+2θ2B+3(θ12+θ22-2ρθ1θ2)4(1-ρ2)=0.
For any x∈ℝ, denote η(x) as the positive conjugate point of x with respect to the curve defined by (4.14). Therefore,
(4.15)η(x)=4(ρ2-1)x2+(8ρθ2-8θ1)x+4θ22+3(θ12+θ22-2ρθ1θ2)/4(1-ρ2)-(2θ2+2ρx)2.Assume thatσ(x)=η(x)-1. It is easy to see that the equation has positive solution only when θi<0. Hence, we have
(4.16)σ^(z)=24(ρ2-1)z2+(8ρθ2-8θ1)z+4θ22+3(θ12+θ22-2ρθ1θ2)/4(1-ρ2)-(2θ2+2ρz).
Therefore, from (4.10), we obtain the following.Theorem 4.3.
The Laplace transform of the stationaryπ(x,y) bears the following form:(4.17)π^(z,w)=2z2+2ρzw+w2-θ1z-θ2w×{w4(ρ2-1)w2+(8ρθ2-8θ1)w+4θ22+3(θ12+θ22-2ρθ1θ2)/4(1-ρ2)-(2θ2+2ρw)-z4(ρ2-1)z2+(8ρθ2-8θ1)z+4θ22+3(θ12+θ22-2ρθ1θ2)/4(1-ρ2)-(2θ2+2ρz)}.


## 5. Final Calculations

Recall that our goal is to get an expression for(5.1)P[X1(t)∈dx,M1(t)∈dz;X2(t)∈dy,M2(t)∈dw],
which is, of course,
(5.2)∂2∂z∂wP[X1(t)∈dx,M1(t)≤dz;X2(t)∈dy,M2(t)≤dw].
The results in the previous section enable us to obtain the Laplace transform of
(5.3)P[X1(t)∈dx,M1(t)≤dz;X2(t)∈dy,M2(t)≤dw].
Apply the Laplace transform with respect to x,y, we get
(5.4)∬-∞∞esx+tyP[X1(1)∈dx,M1(1)≤z;X2(1)∈dy,M2(1)≤w]=∬-∞∞esx+tyP[M1(1)≤z,M2(1)≤w∣X1(1)∈dx,X2(1)∈dy]∬-∞∞×P[X1(1)∈dx,X2(1)∈dy]=∬-∞∞esx+tyP[Y1≤z-x,Y2(1)≤w-y]=∬-∞∞esx+ty(∫-∞z-x∫-∞w-yπz,w(u,v)dudv)P[X1(1)∈dx,X2(1)∈dy]dxdy=1ste-sz+twπ^(-s,-t;z,w)*ϕ^(s,t),
where ϕ^(s,t) is the Laplace transform of P[X1(1)∈dx,X2(1)∈dy]. Hence, we obtain the following.Theorem 5.1.
The laplace transform of(5.5)P[X1(1)∈dx,M1(1)≤dz;X2(1)∈dy,M2(1)≤dw]
with respect (x,y) has the form,
(5.6)1ste-sz+twπ^(-s,-t;z,w)*ϕ(s,t).


## 6. Conclusions

We studied a simple but representative operations research model for performance analysis and capacity planning. Through diffusion approximation, we identify the problem of characterizing the joint distribution of a two-dimensional Brownian motion and its coordinate running maximum, which is of independent interests in the theory of probability, as well as other applications. Using some probabilistic techniques, we are able to reduce the problem to the solution of an extended Helmholtz equation, and a careful exam of the solution to this well-known partial differential equation leads to the calculation of the Laplace transform of the desired joint probability distribution.Standard inversion methods exist for estimating the probabilities from the Laplace transform, in the ongoing research, we are exploring the special structure of our solution, and tailor the inversion method to our special needs, thus provide a very efficient and accurate way of estimating these fundamentally important quantities for both probability theory and operations research.



---

*Source: 109417-2012-12-31.xml*