Performance Analysis of Production Systems with Correlated Demand via Diffusion Approximations

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.


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 International Journal of Stochastic Analysis 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 Section 2, 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.

Model
Suppose that there are two types of products, type 1 and type 2. Each of them is supplied by an independent renewal process.Let us denote them as A 1 t and A 2 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 X 1 t , X 2 t be a two-dimensional Brownian motion with constant correlation ρ > 0 and initial state X 1 0 , X 2 0 0, 0 .Therefore, we have , and E X 1 s X 2 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 t 0 > 0 where M i t are the running maximum for X i t for i 1, 2 respectively, that is M i t sup 0≤s≤t {X i s }.It is easy to see that we can just focus on the case of t 0 1, general case can be easily obtained through scaling.To be more specific, we have, for any t 0 , because X 1 t , X 2 t and √ t 0 X 1 t/t 0 , √ t 0 X 2 t/t 0 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 twodimensional Brownian motion in the positive orthant with negative drifts through solving a partial differential equation.

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 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 , where B t is a standard Brownian motion, d means equal in distribution.Thus, for any z ≥ x, the event of will have the same probability as, Apply a variable transform s 1/t − 1, we get, {B s ≤ zs z − x , ∀s ≥ 0}, 3.5 or equivalently, {B s − zs ≤ z − x, ∀s ≥ 0}.

3.6
This probability of the above event is equivalent to the stationary distribution of B s − zs reflected at zero, which, in turn, equals 1 − exp −2z z − x .

3.7
Take derivative with respect to z, we have,

3.9
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 for z ≥ x and w ≥ y can be reduced to calculating 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 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
where W 1 , W 2 is a two-dimensional driftless Brownian motion with covariance matrix, Γ 1 ρ ρ 1 , 3.14 and L i t , i 1, 2 are the local time defined as 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 .

Stationary Distribution
To calculate the stationary distribution of the reflected Brownian motion Y 1 t , Y 2 t , we make use the following result from Harrison and Williams 6 .
is the infinitesimal generator and v j 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 R 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 function u x, y such that Lu 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: Let us first consider the homogeneous version of the equation, Apply transform, then, u ξ, η satisfies where

4.9
For the inhomogeneous version, it is easy to see that u − 1/ z 2 2ρzw w 2 − θz − θw e −zx−wy will be the solution.Now plug u x, y into the BAR 4.1 , we have where π and σ denote the Laplace transform of π x, y and σ x .Given that c i , i 1, . . ., 4 and k can be arbitrary, it is necessary that

Theorem 4 . 1 . 2 1 6International
Functions π x, y ∈ C 2 R 2 and σ i x ∈ C 2 R 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: R Lf z π x, y dx dy 2 j 1 F j D j f z σ i z dv j z 0, 4.Journal of Stochastic Analysis where