Exact Solution for the Production Planning Problem with Several Regimes Switching over an Infinite Horizon Time

1. Introduction and proposal of the paper

We consider a factory producing $N\ge 1$ types of economic goods that stores them in an inventory-designated place. The model is described mathematically in the next.

Let $(\Omega ,\mathcal{F},\mathbb{F}$,$P)$ be a complete filtered probability space, where P is the historical probability and

$$\mathbb{F}=\{\left.{\mathcal{F}}_t\right|t\in \left[0,\infty \right)\},$$

is generated by an ${\mathbb{R}}^N$-valued Brownian motion denoted by $w=\left(w_1,...,w_N\right)$ with respect to the probability P.

In the production planning problem, the regime switching is captured by a continuous time homogeneous Markov chain $\epsilon \left(t\right)$ adapted to $\mathbb{F}$ that can take k different values, modelling k regimes which should be noted by $1,2,...,k$. The Markov chain’s rate matrix that denotes the strongly ireductible generator of $\epsilon$, is denoted by $G={\left[{\vartheta }_{ij}\right]}_{k\times k}$ where

$${\vartheta }_{ii}=-a_{ii}<0\mathrm{for}\mathrm{all}i,{\vartheta }_{ij}=a_{ij}\ge 0\mathrm{for}\mathrm{all}i\ne j,$$

and the diagonal elements ${\vartheta }_{ii}$ may be expressed as

$${\vartheta }_{ii}=-\underset{j\ne i}{\Sigma }{\vartheta }_{ij}.$$

(1)

In this case, if $P_t\left(t\right)=\mathbb{E}\left[\epsilon \left(t\right)\right]\in \mathbb{R}$, then

$$\frac{d{P}_t\left(t\right)}{dt}=G\epsilon \left(t\right).$$

(2)

Moreover, $\epsilon \left(t\right)$ it is explicitly described by the integral form

$$\epsilon \left(t\right)=\epsilon \left(0\right)+{\int }_0^{t}G\epsilon \left(u\right)du+M\left(t\right),$$

(3)

where $M\left(t\right)$ is a martingale with respect to $\mathbb{F}$. Here and hereafter, we use the notation from other papers to keep the applicative character of the problem,

$$p\left(t\right)=\left(p_1\left(t\right),...,p_N\left(t\right)\right),$$

represent the production rate at time t (control variable) adjusted for the demand rate.

These adjusted for demand inventory levels are modeled by the following system of stochastic differential equations

$$d{y}_i\left(t\right)=p_{i}dt+{\sigma }_{\epsilon \left(t\right)}d{w}_i,y_i\left(0\right)=y_i^0\mathrm{for}i=1,...,N,$$

(4)

where $y_i\left(t\right)$ is an Itô process in $\mathbb{R}$ (i.e., the inventory level of good i, at times $t,$adjusted for demand), $p_i$ is the deterministic part, ${\sigma }_{\epsilon \left(t\right)}$ is a random regime-dependent constant (non-zero) diffusion coefficient taking on the values ${\sigma }_1$, ${\sigma }_2$, ... , ${\sigma }_k$ and $y_i^0$ is the initial condition (i.e., initial inventory level of goods i).

The stochasticity here is due to demand adjustment, which is random and dependent on the regime. This is the most commonly used process when the demand is more volatile in some periods (e.g., some states of the Markov chain) and less volatile in other periods.

The performance over time of a demand-adjusted production

$$p\left(t\right)=\left(p_1\left(t\right),...,p_N\left(t\right)\right),$$

is measured by means of its cost. At this point, we introduce the cost functional which yields the cost

$$J\left(p_1,...,p_N\right):=E{\int }_0^{\infty }{\left(\right|p\left(t\right)|}^2+{\left|y\left(t\right)\right|}^2)e^{-{\alpha }_{\epsilon \left(t\right)}t}dt,y\left(t\right)=\left(y_1\left(t\right),...,y_N\left(t\right)\right),$$

(5)

which measures the quadratic loss.

We measure deviations from the demand, from what place the loss. Here ${\alpha }_{\epsilon \left(t\right)}$ is a regime dependent, taking on the values ${\alpha }_1\ge 0$, ${\alpha }_2\ge 0$, ... , ${\alpha }_k\ge 0$, constant psychological rate of time discount from what place the exponential discounting.

At the moment, we are ready to frame our objective, which is to minimize the cost functional, i.e.,

$$\underset{p_1,...,p_N}{inf}J\left(p_1,...,p_N\right),$$

(6)

subject to the Itô equation (4).

This model problem was proposed by Bensoussan, Sethi, Vickson and Derzko [1] in the context of no regime switching in the economy and for the case of a factory producing one type of economic goods. Later, many other authors are concerned with regime switching.

In production management, Cadenillas, Lakner and Pinedo [2] adapted the model problem in [1] to study the optimal production stochastic control planning problem of a company within an economy characterized by the two-state regime switching with limited/unlimited information. Later, Dong, Malikopoulos, Djouadi and Kuruganti [9] applied in the civil engineering the model described by [2] to the study of the optimal stochastic control problem for home energy systems with solar and energy storage devices when the demand is subject to Brownian motion; the two switching regimes are the peak and off peak energy demand.

A good deal of attention to this subject has been also devoted by Pirvu and Zhang [17] where the authors studied the effect of high versus low discount rates to a consumption-investment decision problem.

After that, there have been numerous applications of regime switching in many important problems in economics, operations research, actuarial science, finance, reinsurance, and other fields, see the works of: Capponi and Figueroa-López [3], Elliott and Hamada [11], Gharbi and Kenne [13], Yao, Zhang and Zhou [22] and Wang, Chang and Fang [23] for more details.

There are of course other research studies that may also serve to better explain the importance of regime switching in the real world.

In a precursor to this article, Covei and Pirvu [5], formulate and analyze the production-planning problem in the continuous time case, with no regime switching in the economy over an infinite time. In the paper [7], the author improved the results of [5], in the sense that the value function in the production model is given in the closed form. Related works that deal with no regime switching in the economy are Sheng-Zhu-Wang [20] and Qin-Bai-Ralescu [18].

Recently, Canepa, Covei and Pirvu [4], considered the production planning problem with regime switching in the economy over a finite horizon time. Here, the solution is obtained through numerical approaches. Although a closed form expression for the corresponding case of regime switching on a particular state space consisting of two regimes over an infinite horizon time is available in the paper of [6]. So, at least one question suggested by the paper of [14] has some nice features: can we obtain a closed form solution when the state space consists of several numbers of states? Our present paper fills the gap in the literature by proving a closed form solution to the stochastic production planning problem with regime switching in the economy over an infinite horizon in a general state space.

The technique presented in this paper makes a methodological contribution that is of independent interest in other considerable number of works on regime switching.

To conclude this introduction, our paper is structured as follows. In Section 2 we give the relationship of our model with a system of partial differential equations (PDE system). Section 3 presents a closed form solution and the uniqueness of solution for our production planning problem. A numerical approximation of the solution for the production planning problem is also given in Section 4. In Section 5 we present a verification result. We introduced in Section 6 the equilibrium production rates as the the subgame perfect production rates. They are the output of an interpersonal game between the present self and future selves. The equilibrium production rates are time consistent meaning there is no incentive to deviate from them. It turns out that in our setting the optimal production rates are among the equilibrium ones so they are time consistent. In Section 7, we give some applications. Finally, in Section 8 we want to discuss our strategy.

Having presented the model that we want to solve, now we provide our means to tackle it.

2. Reduction of the model to a PDE system

Our approach is based on the value function and dynamic programming, which leads to the Hamilton-Jacobi-Bellman (HJB) system of equations.

To characterize the value function, we apply the probabilistic approach. We search for functions $V\left(x,1\right)$, ..., $V\left(x,k\right)$ such that the stochastic process $S^p\left(t\right)$ defined below

$$S^p\left(t\right)=e^{-{\alpha }_{\epsilon \left(t\right)}t}V\left(y\left(t\right),\epsilon \left(t\right)\right)-{\int }_0^t{\left[\right|p\left(s\right)|}^2+{\left|y\left(s\right)\right|}^2]e^{-{\alpha }_{\epsilon \left(s\right)}s}ds,$$

(7)

is supermartingale for all

$$p\left(t\right)=\left(p_1\left(t\right),...,p_N(t\right)),$$

and martingale for the optimal control

$$p^{*}\left(t\right)=\left(p_1^{*}\left(t\right),...,p_N^{*}\left(t\right)\right).$$

As shown by [5], if this is achieved, with the following transversality condition

$$\underset{t\to \infty }{lim}E\left[e^{-{\alpha }_{\epsilon \left(t\right)}t}V\left(y\left(t\right),\epsilon \left(t\right)\right)\right]=0,$$

(8)

and some estimates on the value function yield that

$$-V(x,i)=infJ\left(p_1,...,p_N\right),$$

(9)

where $x=\left(x_1,...,x_N\right)\in {\mathbb{R}}^N$ assumes values

$$(y_1\left(0\right),...,y_N\left(0\right)).$$

Once such a function is found, it turns out that $(u_1,...,u_k)$ with

$$u_1\left(x\right)=-V\left(x,1\right),...,u_k\left(x\right)=-V\left(x,k\right),$$

is the value function. We search for $u_1,....,u_k$ the functions in $C^2\left[0,\infty \right)$, and the supermartingale/martingale requirement yields by using Itô’s Lemma for Markov modulated diffusion, the HJB system of equations, which characterizes the value function

$$-\left(\begin{array}{c}\frac{{\sigma }_1^2}{2}\Delta u_1\\ ...\\ \frac{{\sigma }_k^2}{2}\Delta u_k\end{array}\right)+G_{a,\alpha }\left(\begin{array}{c}{u}_1\\ ...\\ u_k\end{array}\right)-\left(\begin{array}{c}{\left|x\right|}^2\\ ...\\ {\left|x\right|}^2\end{array}\right)=\left(\begin{array}{c}\underset{p}{inf}\{p\nabla u_1+{\left|p\right|}^2\}\\ ...\\ \underset{p}{inf}\{p\nabla u_k+{\left|p\right|}^2\}\end{array}\right),$$

(10)

where

$$G_{a,\alpha }=\left(\begin{array}{cccc}{a}_{11}+{\alpha }_1& -a_{12}& ...& -a_{1k}\\ -a_{21}& a_{22}+{\alpha }_2& ...& -a_{2k}\\ ...& ...& ...& ...\\ -a_{k1}& -a_{k2}& ...& a_{kk}+{\alpha }_k\end{array}\right).$$

For the transformation of the HJB system, it is essential to observe that

$$\underset{p}{inf}\{p\nabla u_i+{\left|p\right|}^2\}=-\frac{1}{4}{\left|\nabla u_i\right|}^2,i=1,2,...,k.$$

(11)

Thus, the HJB system (10) can be written as a PDE system

$$\left\{\begin{array}{c}-\frac{{\sigma }_1^2}{2}\Delta u_1+\left(a_{11}+{\alpha }_1\right)u_1-\sum _{i=2}^k{a}_{1i}{u}_i-{\left|x\right|}^2=-\frac{1}{4}{\left|\nabla u_1\right|}^2,\\ ...\\ -\frac{{\sigma }_k^2}{2}\Delta u_k+\left(a_{kk}+{\alpha }_k\right)u_k-\sum _{i=1}^{k-1}{a}_{ki}{u}_i-{\left|x\right|}^2=-\frac{1}{4}{\left|\nabla u_k\right|}^2.\end{array}\right.$$

(12)

To perform the verification, i.e., show that the HJB system gives the solution to the optimization problem, one should write (12) with the following boundary condition

$$u_1\left(x\right)\to \infty ,...,u_k\left(x\right)\to \infty ,\mathrm{as}\left|x\right|\to \infty .$$

(13)

The value function will give us in turn the candidate optimal control. The first-order optimality conditions on the left-hand side of (11) are sufficient for optimality since we deal with a quadratic (convex) function and they produce the candidate optimal control as follows:

$$p_i^{*}\left(t\right)={\overline{p}}_i(y_1\left(t\right),\dots ,y_N\left(t\right),\epsilon \left(t\right)),i=1,...,N,$$

and

$${\overline{p}}_i(x_1,...,x_N,j)=-\frac{1}{2}\frac{\partial u_j}{\partial x_i}\left(x_1,...,x_N\right),\mathrm{for}i\in \{1,...,n\},j\in \{1,...,k\}.$$

(14)

The production rate ${\overline{p}}_i$ is allowed to be negative. A negative production rate would correspond to a write-off or disposal of inventory (for example, due to obsolescence or perishability).

Our next goal of this paper is to determine the candidate optimal control in closed form.

3. Closed form solution for the PDE system

In spite of their clear simplicity, the PDE system (12) with boundary conditions (13) presents a host of mathematical difficulties arising from the presence of nonlinear gradient terms ${\left|\nabla u_1\right|}^2$, ..., ${\left|\nabla u_k\right|}^2$, see for details [8].

The following result will be proved and is the main original element of the article.

Theorem 1. 

Assume that $G_{a,\alpha }$ is a positive definite matrix with all elements of $G_{a,\alpha }^{-1}$ positive. Then, the PDE system (12) with boundary condition (13) has a unique radially symmetric convex positive classical solution with quadratic growth.

Proof of Theorem 1

In the following, we construct the function

$$\left(u_1,...,u_k\right)\in C^2\left[0,\infty \right)\times ...\times C^2\left[0,\infty \right),$$

which satisfies (12) with boundary condition (13). One way of solving this partial differential equation is to show that there exists

$$\left(u_1\left(x\right),...,u_k\left(x\right)\right)=\left({\beta }_1{\left|x\right|}^2+{\eta }_1,...,{\beta }_k{\left|x\right|}^2+{\eta }_k\right),\mathrm{with}{\beta }_1,...,{\beta }_k,{\eta }_1,...,{\eta }_k\in \left(0,\infty \right),$$

(15)

that solves (1).

The main task for the proof of existence of (15) is performed by proving that there exists

$${\beta }_1,...,{\beta }_k,{\eta }_1,...,{\eta }_k\in \left(0,\infty \right),$$

such that

$$\left\{\begin{array}{c}-\frac{2{\beta }_{1}N{\sigma }_1^2}{2}+\left(a_{11}+{\alpha }_1\right)\left({\beta }_1{\left|x\right|}^2+{\eta }_1\right)-\sum _{i=2}^k{a}_{1i}\left({\beta }_i{\left|x\right|}^2+{\eta }_i\right)-{\left|x\right|}^2=-\frac{1}{4}{\left(2{\beta }_1\left|x\right|\right)}^2,\\ ...\\ -\frac{2{\beta }_{k}N{\sigma }_k^2}{2}+\left(a_{kk}+{\alpha }_k\right)\left({\beta }_k{\left|x\right|}^2+{\eta }_k\right)-\sum _{i=1}^{k-1}{a}_{ki}\left({\beta }_i{\left|x\right|}^2+{\eta }_i\right)-{\left|x\right|}^2=-\frac{1}{4}{\left(2{\beta }_k\left|x\right|\right)}^2,\end{array}\right.$$

or equivalently, after grouping the terms

$$\left\{\begin{array}{c}{\left|x\right|}^2\left[-\sum _{i=2}^k{a}_{1i}{\beta }_i+(a_{11}+{\alpha }_1){\beta }_1+{\beta }_1^2-1\right]-{\beta }_{1}N{\sigma }_1^2-\sum _{i=2}^k{a}_{1i}{\eta }_i+(a_{11}+{\alpha }_1){\eta }_1=0,\\ ...\\ {\left|x\right|}^2\left[-\sum _{i=1}^{k-1}{a}_{ki}{\beta }_i+(a_{kk}+{\alpha }_k){\beta }_k+{\beta }_k^2-1\right]-{\beta }_{k}N{\sigma }_k^2-\sum _{i=1}^{k-1}{a}_{ki}{\eta }_i+\left(a_{kk}+{\alpha }_k\right){\eta }_k=0.\end{array}\right.$$

Now, we consider the system of equations

$$\left\{\begin{array}{c}-\sum _{i=2}^k{a}_{1i}{\beta }_i+(a_{11}+{\alpha }_1){\beta }_1+{\beta }_1^2-1=0\hfill \\ ...\hfill \\ -\sum _{i=1}^{k-1}{a}_{ki}{\beta }_i+(a_{kk}+{\alpha }_k){\beta }_k+{\beta }_k^2-1=0\hfill \\ -{\beta }_{1}N{\sigma }_1^2-\sum _{i=2}^k{a}_{1i}{\eta }_i+(a_{11}+{\alpha }_1){\eta }_1=0\hfill \\ ...\hfill \\ -{\beta }_{k}N{\sigma }_k^2-\sum _{i=1}^{k-1}{a}_{ki}{\eta }_i+\left(a_{kk}+{\alpha }_k\right){\eta }_k=0.\hfill \end{array}\right.$$

(16)

To solve (16), we can rearrange those equations 1, ... , k such

$$\left(\begin{array}{ccc}{a}_{11}+{\alpha }_1& ...& -a_{1k}\\ ...& ...& ...\\ -a_{k1}& ...& a_{kk}+{\alpha }_k\end{array}\right)\left(\begin{array}{c}{\beta }_1\\ ...\\ {\beta }_k\end{array}\right)=\left(\begin{array}{c}1-{\beta }_1^2\\ ...\\ 1-{\beta }_k^2\end{array}\right).$$

(17)

The arguments in [15,16] say that the system (17) has a unique positive solution. Next, letting

$$\left({\beta }_1,...,{\beta }_k\right)\in \left(0,\infty \right)\times ...\times \left(0,\infty \right)$$

a unique solution of (17) we observe that the Equations $k+1$, ... , $2k$ of (16) can be written equivalently as

$$\left(\begin{array}{c}{\beta }_{1}N{\sigma }_1^2\\ ...\\ {\beta }_{k}N{\sigma }_k^2\end{array}\right)=\left(\begin{array}{ccc}{a}_{11}+{\alpha }_1& ...& -a_{1k}\\ ...& ...& ...\\ -a_{k1}& ...& a_{kk}+{\alpha }_k\end{array}\right)\left(\begin{array}{c}{\eta }_1\\ ...\\ {\eta }_k\end{array}\right),$$

(18)

from where using the fact that $G_{a,\alpha }^{-1}$ has all elements positive, we can see that there exist and are unique ${\eta }_1$,..., ${\eta }_k\in \left(0,\infty \right)$ that solve (16) and then

$$\left(u_1\left(x\right),...,u_k\left(x\right)\right),$$

solve (12). This finishes the proof of Theorem 1.

Because our solution depends on solving a non-linear algebraic system of equations the exact solution of the PDE system cannot be determined using a computer software. In order to be implemented the solution of the PDE system (12) in a software application in the next section it is necessary to give the numerical approximation of solution to (16), and therefore the arguments in [15,16] are used again.

4. Numerical solution of an algebraic nonlinear system in building the solution for the PDE system

We intend to approximate ${\beta }_1,...,{\beta }_k,{\eta }_1,...,{\eta }_k\in \left(0,\infty \right)$ in (15) by the Newton-Raphson method. To do this, we denote

$$\begin{array}{c}{h}_1\left({\beta }_1,..,{\beta }_k\right)=-\sum _{i=2}^k{a}_{1i}{\beta }_i+(a_{11}+{\alpha }_1){\beta }_1+{\beta }_1^2-1,\hfill \\ ...\hfill \\ h_k\left({\beta }_1,..,{\beta }_k\right)=-\sum _{i=1}^{k-1}{a}_{ki}{\beta }_i+(a_{kk}+{\alpha }_k){\beta }_k+{\beta }_k^2-1,\hfill \end{array}$$

(19)

and

$$J_{\left(h_1,...,h_k\right)}=\left(\begin{array}{ccc}{a}_{11}+{\alpha }_1+2{\beta }_1& ...& -a_{1k}\\ ...& ...& ...\\ -a_{k1}{\beta }_1& ...& a_{kk}+{\alpha }_k+2{\beta }_k\end{array}\right),$$

the Jacobian matrix of (19). For $n=1,2,...$ we find the approximate of the unique parameters

$$\left({\beta }_1,...,{\beta }_k\right)\in \left(0,\infty \right)\times ...\times \left(0,\infty \right),$$

in the following way

$$\left(\begin{array}{c}{\beta }_1^{n+1}\\ ...\\ {\beta }_k^{n+1}\end{array}\right)=\left(\begin{array}{c}{\beta }_1^n\\ ...\\ {\beta }_k^n\end{array}\right)-{\left(\begin{array}{ccc}{a}_{11}+{\alpha }_1+2{\beta }_1^n& ...& -a_{1k}\\ ...& ...& ...\\ -a_{k1}& ...& a_{kk}+{\alpha }_k+2{\beta }_k^n\end{array}\right)}^{-1}\left(\begin{array}{c}{h}_1\left({\beta }_1^n,..,{\beta }_k^n\right)\\ ...\\ h_k\left({\beta }_1^n,..,{\beta }_k^n\right)\end{array}\right),$$

with ${\beta }_1^0,..,{\beta }_k^0\in \left(0,\infty \right)$. Clearly ${\eta }_1$,..., ${\eta }_k\in \left(0,\infty \right)$ are easy determined from (18).

Now we are moving on to the verification result which is also inspired from [6].

5. Verification

Next, we show that the control of (14) obtained in our reduction strategy is indeed optimal. We apply the supermartingale and martingale approach.

Repeating the same argument in [4], as the first step we can show that the stochastic process $S^p\left(t\right)$ defined below

$$S^p\left(t\right)=e^{-{\alpha }_{\epsilon \left(t\right)}t}V\left(y\left(t\right),\epsilon \left(t\right)\right)-{\int }_0^t{\left[\right|p\left(s\right)|}^2+{\left|y\left(s\right)\right|}^2]e^{-{\alpha }_{\epsilon \left(s\right)}s}ds,$$

is supermartingale for all

$$p\left(t\right)=\left(p_1\left(t\right),...,p_N\left(t\right)\right),$$

and martingale for the optimal control

$$p^{*}\left(t\right)=\left(p_1^{*}\left(t\right),...,p_N^{*}\left(t\right)\right).$$

Owing to the well-known Itô Lemma for Markov modulated diffusion (see [22] for more on this) we have

$$\begin{array}{ccc}\hfill d{S}^p\left(s\right)& =& e^{-{\alpha }_{\epsilon \left(s\right)}s}[\frac{{\sigma }_{\epsilon \left(s\right)}^2}{2}\Delta V\left(y\left(s\right),\epsilon \left(s\right)\right)-{\left|y\left(s\right)\right|}^2+p\left(s\right)\nabla V\left(y\left(s\right),\epsilon \left(s\right)\right)\hfill \\ & & -{\left|p\left(s\right)\right|}^2-({\alpha }_{\epsilon \left(s\right)}+a_{\epsilon \left(s\right)\epsilon \left(s\right)})V\left(y\left(s\right),\epsilon \left(s\right)\right)\hfill \\ & & +\sum _{i=1,i\ne \epsilon \left(s\right)}^k{a}_{\epsilon \left(s\right)i}V\left(y\left(s\right),i\right)]ds+dZ\left(s\right),\hfill \end{array}$$

for some martingale $Z\left(s\right)$, and $Z\left(0\right)=0$. Therefore

$$\begin{array}{ccc}\hfill E{S}^p\left(t\right)& =& S^p\left(0\right)+E\left[{\int }_0^t{e}^{-{\alpha }_{\epsilon \left(s\right)}s}[\frac{{\sigma }_{\epsilon \left(s\right)}^2}{2}\Delta V\left(y\left(s\right),\epsilon \left(s\right)\right)-{\left|y\left(s\right)\right|}^2+p\left(s\right)\nabla V\left(y\left(s\right),\epsilon \left(s\right)\right)]ds\right]\hfill \\ & & +E\left[{\int }_0^t{e}^{-{\alpha }_{\epsilon \left(s\right)}s}[-{\left|p\left(s\right)\right|}^2-({\alpha }_{\epsilon \left(s\right)}+a_{\epsilon \left(s\right)\epsilon \left(s\right)})V\left(y\left(s\right),\epsilon \left(s\right)\right)]ds\right]\hfill \\ & & +E\left[{\int }_0^t{e}^{-{\alpha }_{\epsilon \left(s\right)}s}\left[\sum _{i=1,i\ne \epsilon \left(s\right)}^k{a}_{\epsilon \left(s\right)i}V\left(y\left(s\right),i\right)\right]ds\right].\hfill \end{array}$$

Then, the claim yields considering HJB equation (10) and (12) which says that $S^p\left(t\right)$ is martingale for the optimal control and supermartingale otherwise. This last fact combined with the transversality condition yields the claim.

In the second step, let us establish the optimality of $\left(p_1^{*},...,p_N^{*}\right)$. Considering the quadratic estimate on the value function

$$V\left(x,1\right)=-{\beta }_1{\left|x\right|}^2-{\eta }_1,...,V\left(x,k\right)=-{\beta }_k{\left|x\right|}^2-{\eta }_k,$$

(20)

where ${\beta }_i$, ${\eta }_i\in \left(0,\infty \right)$ are the solution of (16).

Let us provide a lower bound estimate for ${\alpha }_1,...,{\alpha }_k$so that the transversality condition (8) is met

$$\underset{t\to \infty }{lim}E[e^{-{\alpha }_{ϵ\left(t\right)}t}|y\left(t\right){|}^2]=0,$$

holds true. The SDE system (4) in this case becomes

$$d{y}_i\left(t\right)=-{\beta }_{ϵ\left(t\right)}{y}_i\left(t\right)dt+{\sigma }_{\epsilon \left(t\right)}d{W}^i\left(t\right),i=1,\dots N.$$

Using Itô’s Lemma, one gets

$$\begin{array}{ccc}\hfill d{\left(y_i\left(t\right)\right)}^2& =& 2y_i\left(t\right)d{y}_i\left(t\right)+d{y}_i\left(t\right)d{y}_i\left(t\right)\hfill \\ & =& [-2{\beta }_{ϵ\left(t\right)}{\left(y_i\left(t\right)\right)}^2+{\sigma }_{ϵ\left(t\right)}^2]dt+2y_i\left(t\right){\sigma }_{ϵ\left(t\right)}d{W}^i\left(t\right).\hfill \end{array}$$

We introduce

$$F_i\left(t\right)=E\left[{\left(y_i\left(t\right)\right)}^2\right].$$

By taking expectations in the above equation, we get

$$\begin{array}{ccc}\hfill F_i\left(t\right)& =& E\left[{\int }_0^t[-2{\beta }_{ϵ\left(s\right)}{\left(y_i\left(s\right)\right)}^2+{\sigma }_{ϵ\left(s\right)}^2]ds+\left[{\left(y_i\left(0\right)\right)}^2\right]\right]\hfill \\ & =& E\left[{\int }_0^t[-2{\beta }_{ϵ\left(s\right)}{\left(y_i\left(s\right)\right)}^2+{\sigma }_{ϵ\left(s\right)}^2]ds\right]+y_i^2\left(0)\right).\hfill \end{array}$$

Let

$$D_2=max\{{\sigma }_1^2,...,{\sigma }_k^2\},D_3=max(\left[{\left(y_1\left(0\right)\right)}^2\right],...,\left[{\left(y_k\left(0\right)\right)}^2\right]).$$

Then, in the light of the above equation, we get

$$F_i\left(t\right)\le {\int }_0^t{D}_{2}ds+D_3.$$

Hence, we have that

$$F_i\left(t\right)\le D_{2}t+D_3.$$

Therefore, one must to choose ${\alpha }_1,...,{\alpha }_k\in \left(0,\infty \right)$ for the transversality condition to hold true and the proof is completed. Finally, a simple system of nonlinear equations (16) remains to be solved.

6. The Equilibrium Production

For a production rate ${\left\{p_i\left(t\right)\right\}}_{t\ge 0}$ and its corresponding inventory level ${\left\{y_i\left(t\right)\right\}}_{t\ge 0}$ given by (4), we introduce equilibrium production as the subgame perfect production in the definition below (for more on this economic concept see [10]).

Definition 1. 

Let $F=(F_i,i=1,\dots N):\mathbb{R}\times \{1,2,\dots k\}\to {\mathbb{R}}^N$ be a vector map such that for any $x>0$ and $i\in \{1,2,\dots ,k\}$

$$lim\underset{ϵ\downarrow 0}{inf}\frac{J\left({\overline{p}}_i\right)-J\left(p_i^{ϵ}\right)}{ϵ}\le 0,$$

(21)

where the subgame perfect production

$${\overline{p}}_i\left(s\right):=F_i({\overline{y}}_i\left(s\right),ϵ\left(s\right)).$$

Here, the process ${\left\{{\overline{y}}_i\left(s\right)\right\}}_{s\ge 0}$ is the inventory level process corresponding to ${\left\{{\overline{p}}_i\left(s\right)\right\}}_{s\ge 0}$. The production rate ${\left\{{p^{ϵ}}_i\left(s\right)\right\}}_{s\ge 0}$ is defined by

$${p^{ϵ}}_i\left(s\right)=\left\{\begin{array}{c}{\overline{p}}_i\left(s\right),s\in [0,\infty ]\setminus E_{ϵ,0}\hfill \\ p_i\left(s\right),s\in E_{ϵ,0},\hfill \end{array}\right.$$

(22)

with $E_{ϵ,0}=[0,ϵ];$${\left\{p_i\left(s\right)\right\}}_{s\in E_{ϵ,0}}$ is any production rate. If (21) holds true, then ${\overline{p}}_i\left(s\right),$$i=1\dots N,$ is a subgame perfect production rate.

The equilibrium production are by design time consistent meaning that they will be implemented at a future date even if the optimization criterion is updated. In some situations the optimal production may be time inconsistent meaning that they will fail to be implemented in the future because they are not optimal anymore if the optimization criterion is updated; they will be implementable only in the presence of a commitment mechanism, that is why sometimes they are referred as pre commitment production. Let us remark that in our setting the optimal production rate

$${\overline{p}}_i,i\in \{1,...,N\},$$

(23)

is a subgame perfect production with

$$F_i(x,j):=-\frac{1}{2}\frac{\partial u_j}{\partial x_i}\left(x\right),$$

since

$$({\overline{p}}_i,i=1\dots N)=arg\underset{p_1,...,p_N}{min}J\left(p_1,...,p_N\right)$$

and thus (21) is automatically satisfied. Therefore the equilibrium production is time consistent.

7. Applications

We offer some applications, which also are inspired by the paper of Ghosh, Arapostathis, Marcus [14].

Application 1. Suppose there is one machine producing two products and let $\epsilon \left(t\right)$ the machine state that can take values in two regimes 1=good and 2=bad, i.e., for every $t\in \left[0,\infty \right)$ we have $\epsilon \left(t\right)\in \{1,2\}$. We consider $\epsilon \left(t\right)$ a continuous time Markov chain with generator

$$\left(\begin{array}{cc}-\frac{1}{2}& \frac{1}{2}\\ \frac{1}{2}& -\frac{1}{2}\end{array}\right),$$

and the inventory $y_i\left(t\right)$ which is governed by the Itô system of stochastic differential equations (4) with the diffusion ${\sigma }_1={\sigma }_2=\frac{1}{\sqrt{2}}$ and let ${\alpha }_1={\alpha }_2=\frac{1}{2}$ the discount factor. Under these assumptions, the system (17) becomes

$$\left(\begin{array}{cc}{a}_{11}+{\alpha }_1& -a_{11}\\ -a_{22}& a_{22}+{\alpha }_2\end{array}\right)\left(\begin{array}{c}{\beta }_1\\ {\beta }_2\end{array}\right)=\left(\begin{array}{c}1-{\beta }_1^2\\ 1-{\beta }_2^2\end{array}\right),$$

or, with our data

$$\left\{\begin{array}{c}{\beta }_1^2+{\beta }_1-\frac{1}{2}{\beta }_2-1=0\hfill \\ {\beta }_2^2-\frac{1}{2}{\beta }_1+{\beta }_2-1=0\hfill \end{array}\right.$$

which has a unique positive solution

$${\beta }_1=\frac{1}{4}\left(\sqrt{17}-1\right),{\beta }_2=\frac{1}{4}\left(\sqrt{17}-1\right).$$

On the other hand the system (18) becomes

$$\left(\begin{array}{c}{\beta }_{1}N{\sigma }_1^2\\ {\beta }_{2}N{\sigma }_2^2\end{array}\right)=\left(\begin{array}{cc}{a}_{11}+{\alpha }_1& -a_{11}\\ -a_{22}& a_{22}+{\alpha }_2\end{array}\right)\left(\begin{array}{c}{\eta }_1\\ {\eta }_2\end{array}\right),$$

or, with our data

$$\left(\begin{array}{c}{\beta }_1\\ {\beta }_2\end{array}\right)=\left(\begin{array}{cc}1& -\frac{1}{2}\\ -\frac{1}{2}& 1\end{array}\right)\left(\begin{array}{c}{\eta }_1\\ {\eta }_2\end{array}\right),$$

which has a unique positive solution

$${\eta }_1=\frac{4}{3}{\beta }_1+\frac{2}{3}{\beta }_2=\frac{1}{2}\left(\sqrt{17}-1\right),{\eta }_2=\frac{2}{3}{\beta }_1+\frac{4}{3}{\beta }_2=\frac{1}{2}\left(\sqrt{17}-1\right).$$

Then

$$V\left(\left(x_1,x_2\right),1\right)=V\left(\left(x_1,x_2\right),2\right)=-\frac{1}{4}\left(\sqrt{17}-1\right)\left(x_1^2+x_2^2\right)-\frac{1}{2}\left(\sqrt{17}-1\right)$$

and furthermore, the production rate is

$${\overline{p}}_i(x_1,x_2,j)=-\frac{1}{2}\left(\sqrt{17}-1\right)x_i,\mathrm{for}i\in \{1,2\},j\in \{1,2\}.$$

We also give the approximate of ${\beta }_1$, ${\beta }_2$, ${\eta }_1$, ${\eta }_2$ by using the Newton-Raphson Method. Denote

$$\begin{array}{c}{h}_1\left({\beta }_1,{\beta }_2\right)=-a_{12}{\beta }_2+(a_{11}+{\alpha }_1){\beta }_1+{\beta }_1^2-1\hfill \\ h_2\left({\beta }_1,{\beta }_2\right)=-a_{21}{\beta }_1+(a_{22}+{\alpha }_2){\beta }_2+{\beta }_2^2-1\hfill \end{array}$$

and

$$J_{\left(h_1,...,h_k\right)}=\left(\begin{array}{cc}2{\beta }_1+1& -\frac{1}{2}\\ -\frac{1}{2}& 2{\beta }_2+1\end{array}\right).$$

We construct

$$\left\{\begin{array}{c}\left(\begin{array}{c}{\beta }_1^{n+1}\\ {\beta }_2^{n+1}\end{array}\right)=\left(\begin{array}{c}{\beta }_1^n\\ {\beta }_2^n\end{array}\right)-{\left(\begin{array}{cc}{a}_{11}+{\alpha }_1+2{\beta }_1^n& -a_{1k}\\ -a_{k1}& a_{kk}+{\alpha }_k+2{\beta }_k^n\end{array}\right)}^{-1}\left(\begin{array}{c}{h}_1\left({\beta }_1^n,{\beta }_2^n\right)\\ h_2\left({\beta }_1^n,{\beta }_2^n\right)\end{array}\right)\hfill \\ {\beta }_1^0={\beta }_2^0=0.1.\hfill \end{array}\right.$$

Using the standard computation, approximations to four digits are

$$\begin{array}{cccc}n=1⟹& {\beta }_1^1=1.4429& \mathrm{and}& {\beta }_2^1=1.4429\\ n=2⟹& {\beta }_1^2=0.9102& \mathrm{and}& {\beta }_2^2=0.9102\\ n=3⟹& {\beta }_1^3=0.7808& \mathrm{and}& {\beta }_2^3=0.7808\\ n=4⟹& {\beta }_1^4=0.7808& \mathrm{and}& {\beta }_2^4=0.7808\end{array}$$

On the other hand

$${\beta }_1={\beta }_2=\frac{1}{4}\left(\sqrt{17}-1\right)\simeq 0.7807.$$

Clearly, the approximations for ${\eta }_1$ and ${\eta }_2$ are

$${\eta }_1={\eta }_2\simeq 1.5616.$$

Application 2. Suppose there is one machine producing three products and let $\epsilon \left(t\right)$ the machine state that can take values in three regimes 1, 2, 3, i.e., for every $t\in \left[0,\infty \right)$ we have $\epsilon \left(t\right)\in \{1,2,3\}$. We consider $\epsilon \left(t\right)$ a continuous time Markov chain with generator

$$\left(\begin{array}{ccc}-3& 3& 0\\ 4& -7& 3\\ 0& 4& -4\end{array}\right),$$

and the inventory $y_i\left(t\right)$ which is governed by (4) with ${\sigma }_1={\sigma }_2={\sigma }_3=\frac{1}{\sqrt{3}}$ and let ${\alpha }_1={\alpha }_2={\alpha }_3=1$ the discount factor. Under these assumptions, the system (17) becomes

$$\left(\begin{array}{ccc}{a}_{11}+1& -a_{11}& 0\\ -a_{22}& a_{22}+a_{11}+1& -a_{11}\\ 0& -a_{22}& a_{22}+1\end{array}\right)\left(\begin{array}{c}{\beta }_1\\ {\beta }_2\\ {\beta }_3\end{array}\right)=\left(\begin{array}{c}1-{\beta }_1^2\\ 1-{\beta }_2^2\\ 1-{\beta }_3^2\end{array}\right),$$

or, with our data

$$\left\{\begin{array}{c}{\beta }_1^2+4{\beta }_1-3{\beta }_2-1=0\hfill \\ {\beta }_2^2+8{\beta }_2-4{\beta }_1-3{\beta }_3-1=0\hfill \\ {\beta }_3^2+5{\beta }_3-4{\beta }_2-1=0\hfill \end{array}\right.$$

which has a unique positive solution

$${\beta }_1={\beta }_2={\beta }_3=\frac{1}{2}\left(\sqrt{5}-1\right).$$

On the other hand, the system (18) becomes

$$\left(\begin{array}{c}{\beta }_{1}N{\sigma }_1^2\\ {\beta }_{2}N{\sigma }_2^2\\ {\beta }_{3}N{\sigma }_3^2\end{array}\right)=\left(\begin{array}{ccc}{a}_{11}+1& -a_{11}& 0\\ -a_{22}& a_{22}+a_{11}+1& -a_{11}\\ 0& -a_{22}& a_{22}+1\end{array}\right)\left(\begin{array}{c}{\eta }_1\\ {\eta }_2\\ {\eta }_3\end{array}\right),$$

or, with our data

$$\left(\begin{array}{c}{\beta }_1\\ {\beta }_2\\ {\beta }_3\end{array}\right)=\left(\begin{array}{ccc}3+1& -3& 0\\ -4& 4+3+1& -3\\ 0& -4& 4+1\end{array}\right)\left(\begin{array}{c}{\eta }_1\\ {\eta }_2\\ {\eta }_3\end{array}\right),$$

from where

$$\left(\begin{array}{c}{\eta }_1\\ {\eta }_2\\ {\eta }_3\end{array}\right)=\left(\begin{array}{ccc}\frac{7}{13}& \frac{15}{52}& \frac{9}{52}\\ \frac{5}{13}& \frac{5}{13}& \frac{3}{13}\\ \frac{4}{13}& \frac{4}{13}& \frac{5}{13}\end{array}\right)\left(\begin{array}{c}{\beta }_1\\ {\beta }_2\\ {\beta }_3\end{array}\right),$$

has a unique positive solution

$${\eta }_1={\eta }_2={\eta }_3=\frac{1}{2}\sqrt{5}-\frac{1}{2}.$$

Then

$$\begin{array}{ccc}\hfill V\left(\left(x_1,x_2,x_3\right),1\right)& =& V\left(\left(x_1,x_2,x_3\right),2\right)=V\left(\left(x_1,x_2,x_3\right),3\right)\hfill \\ & =& -\frac{1}{2}\left(\sqrt{5}-1\right)\left(x_1^2+x_2^2+x_3^2+1\right)\hfill \end{array}$$

and furthermore, the production rate is

$${\overline{p}}_i(x_1,x_2,x_3,j)=-\frac{1}{2}\left(\sqrt{5}-1\right)x_i,\mathrm{for}i\in \{1,2,3\},j\in \{1,2,3\}.$$

We also point out that the numerical approximations for ${\beta }_1$, ${\beta }_2$, ${\beta }_3$, using Newton-Raphson Method described, are

$$\begin{array}{cccc}n=1⟹& {\beta }_1^1=0.8418& {\beta }_2^1=1.017& {\beta }_3^1=1.2789\\ n=2⟹& {\beta }_1^2=0.6575& {\beta }_2^2=0.6761& {\beta }_3^2=0.7066\\ n=3⟹& {\beta }_1^3=0.6192& {\beta }_2^3=0.6196& {\beta }_3^3=0.6202\\ n=4⟹& {\beta }_1^4=0.618& {\beta }_2^4=0.618& {\beta }_3^4=0.618\end{array}$$

when ${\beta }_1^0=1$, ${\beta }_2^0=2$ and ${\beta }_3^0=3$. Clearly $\frac{1}{2}\left(\sqrt{5}-1\right)\simeq 0.618$.

8. Final Remark and Conclusion

When $w_i$ are correlated with correlation $\rho ,$ the HJB system (10) becomes

$$-\left(\begin{array}{c}\frac{{\sigma }_1^2}{2}\Delta u_1\\ ...\\ \frac{{\sigma }_k^2}{2}\Delta u_k\end{array}\right)+G_{a,\alpha }\left(\begin{array}{c}{u}_1\\ ...\\ u_k\end{array}\right)-\frac{\rho }{2}\left(\begin{array}{c}{\sigma }_1^2\sum _{i\ne j}\frac{{\partial }^2{u}_1}{\partial x_i\partial x_j}\\ ...\\ {\sigma }_k^2\sum _{i\ne j}\frac{{\partial }^2{u}_k}{\partial x_i\partial x_j}\end{array}\right)-\left(\begin{array}{c}{\left|x\right|}^2\\ ...\\ {\left|x\right|}^2\end{array}\right)=\left(\begin{array}{c}\underset{p}{inf}\{p\nabla u_1+{\left|p\right|}^2\}\\ ...\\ \underset{p}{inf}\{p\nabla u_k+{\left|p\right|}^2\}\end{array}\right),$$

which has the same solution as (10), due to the mixed derivative terms (see [8] for details).

In summary, we have reduced the stochastic production-planning problem with several regime switching in the economy to demonstrate that there is an exact solution for the PDE system which models the stochastic production problem.