Uncertainty quantification with risk measures in production planning

TL;DR

This paper develops a simulation-based optimization framework for production planning under uncertainty, using risk measures inspired by finance to evaluate system performance amid random capacity fluctuations and workforce availability.

Contribution

It introduces a novel application of finance-inspired risk measures to stochastic production and workforce planning problems, integrating them into a simulation-based optimization approach.

Findings

Effective risk measures for production capacity variability

Insights into workforce availability impacts

Framework adaptable to various uncertain environments

Abstract

This paper is concerned with a simulation study for a stochastic production network model, where the capacities of machines may change randomly. We introduce performance measures motivated by risk measures from finance leading to a simulation based optimization framework for the production planning. The same measures are used to investigate the scenario when capacities are related to workers that are randomly not available. This corresponds to the study of a workforce planning problem in an uncertain environment.

Tables

Table 1. Table 1: Comparison of the different profit evaluations

	$(0.4,1)$	$(0.9,0)$	$(0.5,0.7)$	$(0.6,0.9)$	$(0.5,0.8)$
$\overline{\Pi \left(T\right)}$	6.2087	$1.9018$	$5.9841$	$6.1532$	$6.1820$
$\sigma \left(\Pi \left(T\right)\right)$	$2.2415$	1.9261	$2.1142$	$2.1117$	$2.0602$
	$0.0210$	$0.1890$	0.0120	$0.0140$	0.0120
${V@R}_{0.1}\left(\Pi \left(T\right)\right)$	$-3.0895$	$0.8949$	$-2.9776$	-3.2336	$-3.1459$
${\mathrm{AV}@R}_{0.1}\left(\Pi \left(T\right)\right)$	$-1.3008$	$1.6910$	$-1.4979$	$-1.5509$	-1.6523

Table 2. Table 2: Comparison of the different profit evaluations for ( N 1 , N 2 ) superscript 𝑁 1 superscript 𝑁 2 (N^{1},N^{2})

	$(10,12)$	$(15,15)$	$(10,13)$	$(7,9)$	$(8,10)$
$\overline{\Pi }\left[{10}^3\right]$	2.3535	$-4.0108$	$1.7406$	$1.4955$	$1.8488$
$\sigma \left(\Pi \right)\left[{10}^3\right]$	$1.6026$	0.6733	$1.2503$	$1.1219$	$1.2804$
	0.0804	$1.0000$	$0.0883$	$0.0966$	$0.0831$
${V@R}_{0.1}\left(\Pi \right)\left[{10}^3\right]$	-0.2283	$4.8542$	$-0.0995$	$-0.0271$	$-0.1818$
${\mathrm{AV}@R}_{0.1}\left(\Pi \right)\left[{10}^3\right]$	$0.7212$	$5.6145$	$0.6568$	$0.6682$	0.6255
${V@R}_{0.01}\left(\Pi \right)\left[{10}^3\right]$	$1.8642$	$6.5899$	1.5683	$1.5774$	$1.6160$
${\mathrm{AV}@R}_{0.01}\left(\Pi \right)\left[{10}^3\right]$	$2.6565$	$7.2945$	$2.1790$	2.0632	$2.3147$

Equations

$${\partial }_t\rho (x,t)+{\partial }_x\min \{v\rho (x,t),\mu \}=0,$$

$${\partial }_{t}q(t)=g_{\text{in}}(t)-g_{\text{out}}(t),$$

$$g_{\text{out}}(t)=\{\begin{array}{ll}\min \{G_{\text{in}}(t),\mu \},& \text{ if }q(t)=0,\\ \mu ,& \text{ if }q(t)>0.\end{array}$$

$$g_{\text{in}}^e(t)$$

$$g_{\text{out}}^e(t)$$

$${\Phi }_{st}:E$$

$$(r_0,q_0,{\rho }_0)$$

$$E=\{1,\dots ,C^1\}\times \cdots \times \{1,\dots ,C^N\}\times {\mathbb{R}}_{\ge 0}^N\times L^1((a_1,b_1))\times \cdots \times L^1((a_N,b_N)),$$

$$r(t)=r_0,(q(t),\rho (t))=S_{st}^{\mu (r_0)}(q_0,{\rho }_0).$$

$${\lambda }_{ii}^e=\sum _{j\ne i}^{C^e}{\lambda }_{ij}^e.$$

$$\psi (t,y)$$

$$\eta (t,y,B)$$

$${\mu }^e(t)=\sum _{i=1}^{N^e}{X}_i^e(t),$$

$$\mathcal{Q}=\left(\begin{array}{cc}-{\lambda }_0& {\lambda }_0\\ {\lambda }_1& -{\lambda }_1\end{array}\right)$$

$$X(t)=\sum _{i=1}^N{X}_i(t)$$

$$q_{jk}=\{\begin{array}{ll}-(j{\lambda }_1+(N-j){\lambda }_0)& \text{ if }k=j\text{,}\\ (N-j){\lambda }_0& \text{ if }k=j+1\text{,}\\ j{\lambda }_1& \text{ if }k=j-1\text{,}\\ 0& \text{ else}\end{array}$$

$$P(X(t)=j)=\left(\genfrac{}{}{0ex}{}{N}{j}\right)P(X_1(t)=1{)}^{j}P(X_1(t)=0{)}^{N-j}\text{.}$$

$$X(t)\sim \mathrm{Bin}(N,P(X_1(t)=1)).$$

$$\underset{t\to \infty }{\lim }P(X(t)=j)=\left(\genfrac{}{}{0ex}{}{N}{j}\right){\left(\frac{{\lambda }_0}{{\lambda }_0+{\lambda }_1}\right)}^j{\left(\frac{{\lambda }_1}{{\lambda }_0+{\lambda }_1}\right)}^{N-j},$$

$$\underset{t\to \infty }{\lim }P(X_1(t)=1)=\frac{{\lambda }_0}{{\lambda }_0+{\lambda }_1}.$$

$${\lambda }_{ij}^e(t,(q^e,{\rho }^e))=\{\begin{array}{ll}j{\lambda }_1^e+(N^e-j){\lambda }_0^e& \text{ if }j=i\text{,}\\ (N^e-j){\lambda }_0^e& \text{ if }j=i+1\text{,}\\ j{\lambda }_1^e& \text{ if }j=i-1\text{,}\\ 0& \text{ else}\end{array}$$

$$G_{\text{out}}^{\text{net}}(t)$$

$$q^{\text{net}}(t)$$

$$\mathbb{E}[G_{\text{out}}^{\text{net}}(t)],\mathbb{E}[q^{\text{net}}(t)],$$

$${\sigma }^2(G_{\text{out}}^{\text{net}}(t)),{\sigma }^2(q^{\text{net}}(t))$$

$$P(\Pi (t)<0).$$

$${V@R}_{\lambda }(X)=\inf \{m\in \mathbb{R}:P(X+m<0)\le \lambda \}$$

$${V@R}_{\lambda }(X)=-q_X^{+}(\lambda ).$$

$$q_X^{+}(t)=\sup \{x\in \mathbb{R}:P(X<x)\le t\}.$$

$${AV@R}_{\lambda }(X)=\frac{1}{\lambda }{\int }_0^{\lambda }{V@R}_{\gamma }(X)d\gamma .$$

$$\Pi (t)={\int }_0^t\left(\sum _{v\in V_{\text{out}}}\sum _{e\in {\delta }_v^{-}}\min \{v^e{\rho }^e(b^e,s),{\mu }^e(r^e(t))\}\cdot p(s)-\sum _{e\in \mathcal{A}}{q}^e(s)\cdot C_q^e(s)\right)ds.$$

Full text

AMS Classification: 90B30, 90B15, 62P30

Keywords: uncertainty quantification, production networks, coupled PDE-ODE system

1 Introduction

Uncertainty quantification is currently an active research topic including a wide field of applications. In this work, we focus on the numerical evaluation of performance measures for a production network model whose dynamics is stochastic in the sense that machine failures or capacity drops can randomly occur. Measuring the performance of this stochastic production is related to the quantification of uncertainty or risk. In the production context, we can directly relate risk to the profit of an enterprise and therefore focus on the comparison of classical evaluation methods with monetary risk measures from finance.

Having tailored performance measures at hand, production planning considerations and questions of control can be studied. For example, production planning and control of continuous production models have been studied in [1, 18] and production planning of discrete production models with randomness in [15]. Obviously, there is a high interest in the definition of performance measures that support planning decisions and also allow for optimization purposes, see [16, 17].

The stochastic production network model under consideration has been originally introduced in [11] and is based on the deterministic production network model from [5]. There have been previous stochastic extensions in [10, 12] of the deterministic model but in [11] the dependence of the machine failures on the actual workload of the machine has been introduced leading to a more complex dynamics. In contrast to agent-based models, the continuous equations govern the evolution of aggregated quantities (such as the density of goods) and are valid in the case of homogenous mass production. By using this model, we can predict future outcomes of production and also machine breakdowns or capacity drops. A simulation based optimization approach is applied to study the distribution (or routing) of goods within the network, or, how the capacity should be chosen in advance for a fixed time period. To do so, we nee performance measures, i.e. measures that translate the stochastic outcome into reasonable quantities, to validate decisions.

Considering monetary quantities, appropriate performance measures have been introduced in the literature of finance [8] and are based on so-called risk measures. Famous examples are the expectation, Value at Risk and Average Value at Risk (also called conditional Value at Risk). They are originally introduced in the finance area to quantify financial risk. Since risk measures are fairly general, they can be also used in other contexts. In work [4], risk measures have been introduced for the optimization of oil production. Therein, the focus is the combination and comparison of risk measures in optimization formulations. However, our major goal is to focus on risk measures, their evaluation and impact for the stochastic production network model. Therefore, we examine how the distribution rates for an optimal routing should be chosen on the base of the introduced performance measures in a simulation based optimization. Furthermore, we introduce a setting, where the capacity is determined as a sum of individual capacities, which can be on or off, respectively. This allows for a numerical analysis of the maximal possible capacity, e.g. number of workers, with respect to the performance measures.

This paper is organized as follows: in section 2 we introduce the stochastic production network model and its extensions. Section 3 is devoted to the performance measures. Numerical simulation results are studied in section 4, where mainly two cases are considered: the optimal routing and the workforce planning.

2 Stochastic model equations

The core model is a production network model consisting of a coupled system of partial and ordinary differential equations. We start with the introduction of the deterministic model and present then two possible stochastic extensions.

2.1 The deterministic model

We briefly recall the production network model from [5, 9] and its stochastic extension to a load-dependent model from [11]. To focus on the main ideas, we restrict on the case of a production network consisting of one queue-processor unit first. This means, we consider a processor with an unbounded queue in front, which represents the storage. We assume that the processor is defined on an interval $(a,b)\subset\mathbb{R}$, i.e., with length $L=b-a$, and use $\rho(x,t)$ as the density of goods at $x\in(a,b)$ and time $t\geq 0$. Here, the interpretation of $x$ can be the spatial position of the goods within the processor or the so-called degree of completion. The dynamics of the density is given by the following hyperbolic partial differential equation

[TABLE]

where $\mu\geq 0$ is the maximal capacity and $v>0$ the production velocity. If we prescribe a processing time $T^{\text{prod}}$, we have the relation $v=\frac{L}{T^{\text{prod}}}$. The queue is placed in front of the processor and its length $q$ is modeled by the following ordinary differential equation

[TABLE]

i.e. as the balance between inflow into and outflow. Here, $g_{\text{in}}(t)$ is the inflow into the queue, which can be an externally given inflow $G_{\text{in}}(t)$ in the case if there is no predecessor. If the queue has predecessors, the inflow into the queue is the weighted sum of the outflow out of the incoming processors. The outflow of the queue can be described as follows: if the queue is non-empty, the processor takes its maximal capacity $\mu$, and if the queue is empty, the processor can take the inflow into the queue, which is bounded by the maximal capacity $\mu$. Summarizing, this reads as

[TABLE]

The coupling of the processor and the corresponding queue is prescribed by a boundary condition $\rho(a,t)=\frac{g_{\text{out}}(t)}{v}$ and initial conditions $\rho(x,0)=\rho_{0}(x)\in L^{1}((a,b))$, $q(0)=q_{0}\in\mathbb{R}_{\geq 0}$ are given.

Figure 1 summarizes the modeling equations graphically, where the gray boxes at $a$ represent the queue load following the ODE and the black solid line the density on $(a,b)$ with dynamics described by the hyperbolic PDE.

We now describe the extension to the network case, see [5]. Suppose that equation (1) holds for a density $\rho^{e}$ on interval $(a^{e},b^{e})$ with production velocity $v^{e}>0$ and capacity $\mu^{e}\geq 0$ for every arc $e\in\mathcal{A}$ in a directed network $\mathcal{G}=(\mathcal{V},\mathcal{A})$, where $\mathcal{G}$ describes the production network topology with nodes $\mathcal{V}$ and edges $\mathcal{A}=\{1,\dots,N\}$. We denote by $\delta_{v}^{-}$ and $\delta_{v}^{+}$ the set of all ingoing and outgoing arcs for every vertex $v\in\mathcal{V}$. At vertices without any predecessor $v\in V_{\text{in}}=\{v\in\mathcal{V}\colon\delta_{v}^{-}=\emptyset\}$, we prescribe a time-dependent inflow function $G^{v}_{\text{in}}(t)$ and for every $v\in\mathcal{V}$ with $|\delta_{v}^{+}|>0$ we assume distribution rates $A^{v,e}\in[0,1],e\in\delta_{v}^{+}$, i.e. how much of the flow is distributed the subsequent processors. They satisfy $\sum_{e\in\delta_{v}^{+}}A^{v,e}(t)=1$.

An example of a production network is shown in figure 2, where, for example, at node $2$ the distribution rates are given as $A^{2,2}=\alpha_{1}$ and $A^{2,3}=1-\alpha_{1}$.

We adapt the in- and outflow in and out of the queue of processor $e$ for the network case as follows:

[TABLE]

and

[TABLE]

where $s(e)$ is the starting node of edge $e$. This deterministic model is well-defined if the initial densities and the inflow functions are of bounded total variation, see [5]. Moreover, we can define a solution if the inflow and initial densities are in $L^{1}$ by using an extended solution operator $S^{\mu}$, i.e. $S^{\mu}_{st}u$ is the solution starting at time $s$ with initial condition $u=(q_{0}^{1},\dots,q_{0}^{N},\rho_{0}^{1},\dots,\rho_{0}^{N})$ at time $t\geq s$ and with capacities $\mu=(\mu^{1},\dots,\mu^{N})$, see [3, 5]. This deterministic flow $S^{\mu}$ will determine the deterministic evolution between the machine failures or capacity drops.

2.2 Stochastic extension: load-dependent capacities

In the following, we introduce the stochastic production network model from [11]. We assume that the capacities $\mu^{e}$ can take finitely many non-negative values $\mu^{e}(i)$ for $i\in\{1,\dots,C^{e}\}$, $C^{e}\in\mathbb{N}$ and we introduce the variable $r(t)=(r^{1}(t),\dots,r^{N}(t))\in\{1,\dots,C^{1}\}\times\cdots\times\{1,\dots,C^{N}\}$ determining the capacities used in the production network at time $t$. Combining all ingredients leads to deterministic dynamics

[TABLE]

with

[TABLE]

where

[TABLE]

To incorporate random capacity drops in the production network, we follow the theory of piecewise deterministic Markov processes, see e.g. [6, 14]. Since we have the deterministic evolution between the jump times given by $\Phi$, we only need to specify the intensity $\psi$ at which jumps occur and the distribution of the jumps $\eta$. To do so, we use rate functions $\lambda_{ij}^{e}$ describing the rate that processor $e$ has a capacity change from $\mu^{e}(i)$ to $\mu^{e}(j)$ and assume

[TABLE]

That means, we assume for all $y=(r,q,\rho)\in E$ and $B\in\sigma(E)$

[TABLE]

where $\varepsilon_{x}$ is the Dirac measure with unit mass in $x$ and $\sigma(E)$ is the $\sigma$-algebra on $E$. If the rate functions are continuous with respect to $(t,y)$ and uniformly bounded, then there exists a stochastic process $Y=(Y(t),t\in[0,T])\subset E$ on some probability space $(\Omega,\mathcal{F},P)$, which is piecewise deterministic between the jumps and follows the deterministic production network equations, see [11].

To construct sample paths of the stochastic process $Y$, we use a numerical scheme for the deterministic evolution $\Phi$ and a thinning algorithm for the jump times simultaneously. The approximation of the deterministic evolution consists of a left-sided Upwind scheme for the densities $\rho^{e}$ and the forward Euler method for the queue-length evolution $q^{e}$. Let $T_{n}\geq 0$ be the time of the $n$th jump to the value of $Y_{n}\in E$, then a thinning algorithm produces the next jump time $T_{n+1}$ and post-jump location $Y_{n+1}$ as it is shown schematically in figure 3.

For the description of the procedure we assume a uniform bound $\bar{\lambda}$ on $\psi$. Starting from $T_{n}$ with the value $Y_{n}$, we take an exponentially distributed time $\xi_{1}$ with mean $\bar{\lambda}^{-1}$ and use the deterministic evolution to obtain the value $\Phi_{T_{n},T_{n}+\xi_{1}}Y_{n}$ at time $T_{n}+\xi_{1}$. With an acceptance rejection method, we decide whether a jump is accepted with probability $\tfrac{\Psi(T_{n}+\xi_{1},\Phi_{T_{n},T_{n}+\xi_{1}}Y_{n})}{\bar{\lambda}^{-1}}$. If a jump is accepted, the new state of the system $Y_{n+1}$ is produced by the kernel $\eta$. This procedure is repeated until the final time horizon is reached.

2.3 Stochastic extension: capacities as clusters

In the following, we study the scenario where the capacity $(\mu^{e}(t),t\geq 0)$ at processor $e$ consists of $N^{e}\in\mathbb{N}$ individual capacities. To be more precisely, we assume that $\mu^{e}(t)$ can be written in the form

[TABLE]

where $(X^{e}_{i}(t),t\geq 0)$ is the capacity of cluster part $i$. This is important if the production capacity depends on $N^{e}$ workers that are randomly not available. We assume that every part of the cluster can be on or off, i.e. $X^{e}_{i}\in\{0,1\}$, and leads to capacity drops of the capacity process $\mu^{e}$. In the context of individuals, $X_{i}^{e}(t)$ represents whether worker $i$ is available for work or not.

Figure 4 shows possible realized capacities in the case of the diamond network, cf. figure 2, at time $t_{1}$ and a later time $t_{2}$. We note that the total capacities are not conserved, i.e. workers are either available or not.

The next question is how we can incorporate these ideas into the setting of the model presented in subsection 2.1. The mathematical idea is to interpret every $(X_{i}^{e}(t),t\geq 0)$ as a continuous time Markov Chain (CTMC) on a common probability space $(\Omega,\mathcal{F},P)$ and that they are independent of each other. The following lemma 2.1 provides the main tool to numerically evaluate the workforce planning problem in sub section 4.2. To simplify the notation, we neglect the index $e$ in the following.

Lemma 2.1**.**

Let, for some $N\in\mathbb{N}$, a family $(X_{1}(t),\dots,X_{N}(t),t\geq 0)$ of independent CTMCs on a common probability space $(\Omega,\mathcal{F},P)$ with values in $\{0,1\}$ and $\mathcal{Q}$-matrix

[TABLE]

for $\lambda_{0},\lambda_{1}>0$ be given. The stochastic process defined by

[TABLE]

is a CTMC with $\mathcal{Q}$-matrix satisfying

[TABLE]

for $j,k=0,\dots,N$. Furthermore, we have

[TABLE]

The proof can be found in the appendix A. Note that the proof consists of basically two steps: first proving the Markov property and second using combinatorics to compute the generator $\mathcal{Q}$ of the process.

We also have the following remark.

Remark 2.2*.*

Equation (3) shows that $X(t)$ is binomially distributed, i.e.

[TABLE]

The steady-state distribution is consequently

[TABLE]

because of

[TABLE]

The latter allows for a simple availability analysis of the capacities because the parameters $\lambda_{0}$ and $\lambda_{1}$ are known from estimations.

Since the entries $q_{ij}$ of the generator $\mathcal{Q}$ describe the transition rate from state $i$ to $j$ with $i\neq j$, we have

[TABLE]

for $i,j\in\{1,\dots,N^{e}\}$.

This choice embeds this load-independent model into the load-dependent model by the special choice of the rate functions $\lambda_{ij}^{e}(t,(q^{e},\rho^{e}))$ independent of $t,q^{e}$ and $\rho^{e}$.

3 Performance measures

As we have seen, the stochastic production network model introduced allows for random capacities capturing load-dependent failure rates. Since the model is driven by a stochastic process, there is a need for tailored evaluation tools, so-called performance measures, if we have given a set of sample paths. This section is devoted to classical performance measures for production models and performance measures based on risk measures motivated from the finance area.

3.1 Classical performance measures

Considering sample paths of the production network at every point in time is too detailed in most cases and aggregated quantities are of high interest. According to [10, 12], the aggregated outflow until time $t\geq 0$ of the complete network can be computed as

[TABLE]

where $V_{\text{out}}=\{v\in\mathcal{V}\colon\delta^{+}_{v}=\emptyset\}$ are the nodes without a subsequent processor. We can also define

[TABLE]

as the cumulative sum of all queue-loads up to time $t\geq 0$. Both quantities above are real-valued random variables. Classical performance measures are for example the expectation

[TABLE]

and variance

[TABLE]

of these quantities.

3.2 Risk measures as performance measures

If we consider the profit until time $t$ as a random variable $\Pi(t)$ (a functional of $Y$), where $Y$ is a stochastic production network model, then we can include monetary aspects. We could also use classical measures as the expectation or higher order moments for the random variable $\Pi(t)$ to describe the performance of the production but it turns out that they are not the best choice especially in the context of optimization, see section 4. The reason is that a high expected profit can incorporate a high risk that we intend to measure in an appropriate manner. One possibility is the quantification of probability of a bankruptcy, i.e.,

[TABLE]

One has to keep in mind that this probability does not include any information about the needed surplus to capture “bad” events. There, we can help us with the so-called Value at Risk and the Average Value at Risk, see e.g. [8], which have been introduced in the context of finance and insurance. They correspond to the class of monetary risk measures, which we introduce in the following definition 3.1 taken from [8].

Definition 3.1** ((Coherent) risk measure).**

Let $\mathcal{H}$ be a linear space of bounded, real-valued functions containing constants. The mapping $\varrho\colon\mathcal{H}\to\mathbb{R}$ is called a monetary risk measure if

for every $X,\tilde{X}\in\mathcal{H}$ with $X\leq\tilde{X}$, we have $\varrho(X)\geq\varrho(\tilde{X})$ (Monotonicity) 2. 2.

for every $X\in\mathcal{H}$, $m\in\mathbb{R}$, it holds that $\rho(X+m)=\rho(X)-m$ (Cash invariance)

and it is called a coherent monetary risk measure if additionally it holds that

for every $\lambda\geq 0$, we have $\varrho(\lambda X)=\lambda\varrho(X)$ (Positive homogeneity) 2. 4.

for every $X,\tilde{X}\in\mathcal{H}$, it follows that $\varrho(X+\tilde{X})\leq\varrho(X)+\varrho(\tilde{X})$ (Subadditivity).

Property 1. of definition 3.1 states that, if $X$ is interpreted as the profit, the risk of a less profitable company is higher. Assume that we have a fixed surplus of $m$, then the risk is reduced by $m$, see property 2. Property 3. induces a normalization, i.e. $\varrho(0)=0$, and property 4. describes risk diversification effect.

One very common risk measure is the Value at Risk ($\operatorname{V@R}_{\lambda}(X)$). It is defined as

[TABLE]

for some level $\lambda\in(0,1)$ and a real-valued random variable $X$ on some probability space $(\Omega,\mathcal{F},P)$; see [8]. The Value at Risk is simply a quantile and cam be rewritten as

[TABLE]

with the upper quantile function

[TABLE]

One can show that $\operatorname{V@R}_{\lambda}$ is a monetary risk measure, which is positive homogeneous but not subadditive. Since the profit of a production network will contain sums of individual costs of machines, we can not guarantee that the Value at Risk incorporates risk diversification effects. One can easily construct a coherent risk measure from the Value at Risk, which is the Average Value at Risk, and it is defined by

[TABLE]

From the computational point of view, we can estimate the Value at Risk (5) from a sample of profit realizations by an estimation of the quantile function such as quantile111Documentation: https://de.mathworks.com/help/stats/quantile.html in Matlab. This can again be used to compute the Average Value at Risk with, e.g.  the Matlab function integral222Documentation: https://de.mathworks.com/help/matlab/ref/integral.html.

4 Computational results

In this section, we numerically investigate the performance measures, their similarities and differences. First, we investigate the optimal routing problem and comment on different combinations of distribution rates in a diamond network with respect to the performance measures. Second, we study the impact of cluster sizes, i.e. available workforce, on the performance measures.

4.1 Distribution parameter planning in the load-dependent case

Let the topology of a production network given as a diamond network, see figure 2, where $\alpha_{1},\alpha_{2}\in[0,1]$ are the two distribution parameters, i.e., a percentage of $\alpha_{1}$ is fed from processor one into queue two ($A^{1,2}(t)=\alpha_{1}$), $1-\alpha_{1}$ from one to three ($A^{1,3}=1-\alpha_{1}$) and the same for $\alpha_{2}$ from processor two to queue five and $1-\alpha_{2}$ to queue four.

As before in section 4.2, we analyze the profit but here for different distribution rates $\alpha_{1}$ and $\alpha_{2}$. To do so, we adopt the profit functional, which is now given as

[TABLE]

The price of the product is $p(s)\geq 0$ and $C_{q}^{e}(s)\geq 0$ is the storage cost at time $s$ for storage $e\in\mathcal{A}$.

We assume $p(s)=1$ and $C_{q}^{e}(s)=0.1$ for every $e=1,\dots,7$ and $s\in[0,T]$. The queue-processor units all have a production velocity $v^{e}=1$ and a length of one, and we start with an empty system at full capacity. The capacities are given by (ordered by states) $\mu^{1}\in\{0,3\}$, $\mu^{2},\mu^{3}\in\{0,1,2\}$, $\mu^{4}\in\{0,1\}$, $\mu^{5}\in\{0,2\}$, $\mu^{6}\in\{0,1,3\}$, and $\mu^{7}\in\{0,2,3\}$. To include the load-dependency we follow [11] and define the Utilization Ratio $\operatorname{UR}$ by

[TABLE]

and the Ratio of Work In Progress and the maximal amount of goods in the machine $\operatorname{RWIP}$ as

[TABLE]

The rate functions of processors $e\in\{1,4,5\}$ are then given by

[TABLE]

with $\lambda^{\text{rep,max,e}}=10$, $\lambda^{\text{rep,min,e}}=4$ and $\lambda^{\text{down},e}=1$.

In the case of three states, i.e., processors 2, 3, 6 and 7, the rate functions read as

[TABLE]

with $\lambda^{\text{rep,max,e}}=10$, $\lambda^{\text{rep,min,e}}=4$ and $\lambda^{\text{down},e}=2$. This means we consider two different types of processors with two or three capacity states, identical repair rates but different breakdown rates. We can not state an explicit stationary expected capacity as in section 4.2 and simulations are needed to determine which combination $(\alpha_{1},\alpha_{2})$ performs in a “optimal” way.

In figure 6 (a)-(b), the sample mean and standard deviation of the profit depending on different choices of $\alpha_{1}$ and $\alpha_{2}$ are shown, where $\alpha_{1},\alpha_{2}\in\{0,0.1,\dots,0.9,1\}$ are evaluated. In figures 6 (c) and 5 (a)-(b), the bankruptcy probability and both the Value at Risk and Average Value at Risk, respectively, are drawn for different combinations of $\alpha_{1}$ and $\alpha_{2}$. The inflow into the network is constant and given by $1.5$, and the time horizon is chosen to be $T=10$.

We see a strong influence of the distribution parameters on the profit. The vertical black line describes the allocation of $(\alpha_{1},\alpha_{2})$, which is the best choice for the corresponding evaluation of the profit. In the case of the sample mean, the choice $(\alpha_{1},\alpha_{2})=(0.4,1)$ leads to the highest value of $6.2087$; see table 1. In contrast, the sample standard deviation is small for the choice $(\alpha_{1},\alpha_{2})=(0.9,0)$, with a value of $1.9261$. The lowest bankruptcy probability follows from the choice $(0.5,0.7)$, which is close to the best choices of the Value at Risk $(0.6,0.9)$ and the Average Value at Risk $(0.5,0.8)$ in figure 5. Regarding the different evaluations of the profit sample, the choice of combination $(\alpha_{1},\alpha_{2})$ strongly depends on the choice of the evaluation criterion, but there is a tendency to equally distribute at node two and more into processor five than in processor four at node three. Feeding more into processor five can also be motivated by the higher capacity compared to the capacity of processor four. The values of the mean capacity, standard deviation, bankruptcy probability and (Average) Value at Risk for all the best choices are listed in table 1. There, the Average Value at Risk performs well because all other evaluation criteria are not worse in this case, specifically, because the Average Value at Risk is a coherent risk measure from the theory of risk measures and is thus a reasonable result.

Obviously, the bankruptcy probability is a worse measure in optimizing $(\alpha_{1},\alpha_{2})$ for this example since the relevant region is very flat and disturbed by the errors from the Monte Carlo simulation. In high contrast, both, the Value at Risk and Average Value at Risk form a nice shaped measure and neglecting the Monte Carlo errors it seems to be convex as well.

4.2 Study of cluster sizes and workforce planning

We analyze the stochastic production network model with two arcs, see figure 7, where capacities are given as in subsection 2.3 and we define the profit $\Pi(t)$ as the cumulative revenue reduced by storage and cluster-size costs until time $t\geq 0$. In formulas, this reads as

[TABLE]

where we denote with $p(s)\geq 0$ the price of the product, $C_{q}^{e}(s)\geq 0$ and $C_{N}^{e}(s)\geq 0$ are the storage and cluster size costs, respectively, at time $s$. Typical examples for the cluster size cost are maintenance costs or salaries. If $\varrho$ is a monetary risk measure, we can interpret $\varrho(\Pi(T))$ as the risk of the company of the aggregated loss and gains until time $T$.

In the following, we study an example and consider a production network model in the form of a chain of two processors. The capacities of the processors satisfy $\mu^{e}(t)\in\{0,\dots,N^{e}\}$ for given cluster sizes $N^{1},N^{2}\in\mathbb{N}$. We further assume that the times between failures and the repair times of each part of a cluster are exponentially distributed with mean $\operatorname{MTBF}^{1}=80$, $\operatorname{MTBF}^{2}=50$ and $\operatorname{MRT}^{1}=10$, $\operatorname{MRT}^{2}=20$. This results in the rates

[TABLE]

We can compute the expected capacities in steady state as

[TABLE]

which implies

[TABLE]

and means that the expected capacity of the second processor for $N^{1}=N^{2}$ is lower than the capacity of the first one. We assume a processing velocity $v^{e}=1$ for both processors, a time step $\Delta t=1$ and a time horizon of $T=365$. The length of each processor is one, and to satisfy the Courant-Friedrichs-Lewy (CFL) stability condition, we choose a spatial step size of $\Delta x=1$. Given that $G_{\text{in}}^{1}(t)=10$ is a constant inflow of ten parts per unit time, we start with an empty production. The only missing parameters are the cluster sizes $N^{1}$ and $N^{2}$ and our goal is to analyze how the profit $\Pi(T)$ changes if we choose different combinations of cluster sizes. To evaluate the profit, we assume storage costs $C^{1}_{q}=C^{2}_{q}=0.01$, cluster size costs $C_{N}^{1}=4$, $C_{N}^{2}=6$ and a price $p=10.02$.

The following simulation results are based on $M=10^{4}$ Monte Carlo samples. Figure 8 contains the sampled values of the expected profit (a), the standard deviation (b) and bankruptcy probability (c).

In this example, the expected profit strongly depends on the cluster sizes, i.e. there are few combinations that lead to a positive expected profit. The highest expected profit is achieved with the choice $N^{1}=10$ and $N^{2}=12$ with a profit of $2.3535\cdot 10^{3}$, as table 2 states.

With respect to the standard deviation, we obtain the combination $N^{1}=N^{2}=15$ having the lowest standard deviation of $0.6733\cdot 10^{3}$; see table 2. This is reasonable because in this case, the production is less affected by capacity drops due to the large cluster sizes.

The bankruptcy probability in figure 8 (c) shows a tight region in which the probability of going bankrupt is low. The combination with $N^{1}=10$ and $N^{2}=12$ leads to a bankruptcy probability of $0.0804$ and the same combination as the expected profit.

For a level of $\lambda=0.1$, figure 9 contains the Value and Average Value at Risk for this level. The cluster sizes $N^{1}=10$ and $N^{2}=12$ lead to a Value at Risk of $-0.2283\cdot 10^{3}$. This can be interpreted as follows: even if we have debts of $-0.2283\cdot 10^{3}$, the probability to face bankruptcy at $T=365$ is lower than $0.1$. The more pessimistic risk measure, Average Value at Risk, leads to a best choice $N^{1}=8$ and $N^{2}=10$ with a $\operatorname{AV@R}_{0.1}(\Pi(T))$ of $0.6255\cdot 10^{3}$, i.e., we should have a surplus of $0.6255\cdot 10^{3}$ for these cluster sizes.

In the case of a level of $\lambda=0.01$, the combination $N^{1}=10$ and $N^{2}=13$ yields a $\operatorname{V@R}_{0.01}(\Pi(T))$ of $1.5683\cdot 10^{3}$; see figure 10 (a) and table 2. The combination $N^{1}=7$ and $N^{2}=9$ implies an $\operatorname{AV@R}_{0.01}(\Pi(T))$ of $2.0632\cdot 10^{3}$; see figure 10 (b).

All introduced performance measure evaluations in table 2 are given for reasonable cluster size choices. Obviously, the best choice for the standard deviation leads to bad performance in all the other performance measures. The allocation $N^{1}=10$ and $N^{2}=12$ implies the best solution in terms of expectation, bankruptcy probability and Value at Risk with level $0.1$. The Value and Average Value at Risk seem to work quite well here and give information about the surplus one has to hold to capture bad events and avoid a bankruptcy with a high probability (given by $1-\lambda$). One has to be careful with the interpretation here. The profit is defined as the cumulative difference of earnings and costs and does not imply detailed information at every time between zero and $T$.

From the point of optimization, the shape of the Value at Risk and Average Value at Risk are the best, since starting from any combination we are directed into the “valley” of good combinations. In the case of the bankruptcy probability we have two disadvantages, namely totally flat regions and high jumps in the values.

5 Conclusions

We have analyzed performance measures for a stochastic production network model. It turned out that different performance measures might lead to totally different results and the choice of the appropriate measure should be a crucial point in the development of an optimization framework. In our examples, the measures Value at Risk and Average Value at Risk always lead to reasonable results and a suitable curvature in the set of feasible solutions.

Future work will deal with the formulation of rigorous optimization problems for the stochastic production network model with respect to relevant performance measures. We will study the quality of solutions and the speed of convergence to an optimum.

Acknowledgments

This work was financially supported by BMBF project ENets (05M18VMA) and the DFG grant No. GO 1920/7-1.

Appendix A Proof of lemma 2.1

Proof.

The proof consists of two steps. First, we prove the Markov property of $(X(t),t\geq 0)$, and we compute the $\mathcal{Q}$-matrix in the second part. Clearly, the stochastic process $(X(t),t\geq 0)$ takes only values in the finite, discrete space $E^{X}:=\{0,\dots,N\}$, which together with the $\sigma$-algebra $\mathcal{E}^{X}=\mathcal{P}(E^{X})$ builds a measurable space. We prove the Markov property with Dynkin’s criterion [7, theorem 10.13] and follow the ideas of [2, problem 4.14]. We have $N$ CTMCs with state space $\{0,1\}$ given, and we construct the $N$-dimensional Markov process $(Y(t),t\geq 0)$ with $Y(t)=(X_{1}(t),\dots,X_{N}(t))$ on the product space $(\Omega^{N},\mathcal{A}^{N},P^{N})$. This process takes values in the measurable space $(E^{Y},\mathcal{E}^{Y})=(\{0,1\}^{N},\mathcal{P}(\{0,1\}^{N}))$, and we obtain the process $(X(t),t\geq 0)$ defined by the measurable surjective mapping

[TABLE]

as $X(t)=\psi(Y(t))$. Let $(U_{t},t\geq 0)$ be the semigroup of Markovian kernels given by $(Y(t),t\geq 0)$, i.e.,

[TABLE]

for every $x\in E^{Y}$, $s,t\geq 0$ and $B\in\mathcal{E}^{Y}$. If we can show that for every $x,\tilde{x}\in E^{Y}$ with $\psi(x)=\psi(\tilde{x})$ the equality

[TABLE]

for every $t\geq 0$ and $\Gamma\in\mathcal{E}^{X}$ holds (Dynkin’s criterion), then we know that the stochastic process $(X(t),t\geq 0)$ is a Markov process. From $x\in\psi^{-1}(\{m\})$, it follows that $x_{\sigma}\in\psi^{-1}(\{m\})$ for every $m\in E^{X}$ and every permutation $\sigma$ because of the structure of $\psi$. Therefore, we can compute

[TABLE]

for a permutation such that $x_{\sigma}=y$ holds. Hence, we have the Markov property of $(X(t),t\geq 0)$ by Dynkin’s criterion, and it is a CTMC. We calculate the transition probabilities of the stochastic process $(X(t),t\geq 0)$ to obtain the $\mathcal{Q}$-matrix. First, we know for every $i=1,\dots,N$ and for every $l,m\in\{0,1\}$ that the transition probability fulfills

[TABLE]

in the limit $\Delta t\to 0$; see [13, pp. 28].

We choose $k,j\in\{0,\dots,N\}$, and we set $C=E^{N}$ in the following. Then, we can write

[TABLE]

using the independence of the stochastic processes and the multi-index notation

[TABLE]

At this point, we can use the transition probability (7) and merge all terms of order $\mathrm{o}(\Delta t)$. This yields

[TABLE]

To simplify the computation, we analyze both summands in the last equation separately, and we start with the first one (8). The product of the indicator function implies $\alpha=\beta$, and we put higher orders of $\Delta t$ into $\mathrm{o}(\Delta t)$ again. By doing this and observing that $k$ has to equal $j$, we see

[TABLE]

The last equality follows from

[TABLE]

since $j$ entries of $\beta$ are one and all the others are zero. Now, we analyze the second summand (9), which we simplify by merging terms to $\mathrm{o}(\Delta t)$. We distinguish two cases, i.e., the case $k=j+1$ and the case $k=j-1$, since all remaining cases are impossible. In detail, if, e.g., $k=j+2$, then $\alpha$ has two entries more with value one, which implies that $\alpha$ and $\beta$ are different in at least two entries. Hence, the product in the second summand is zero. We have

[TABLE]

by using that $(X_{i}(t),i=1,\dots,N)$ are iid. We easily count the remaining sums with combinatorics, i.e.,

[TABLE]

and observe that

[TABLE]

Summarizing all computations yields the transition probability

[TABLE]

as $\Delta t\to 0$ and we conclude the generator of $(X(t),t\geq 0)$ as the matrix $\mathcal{Q}$ from the statement of this lemma. ∎