Uncertainty quantification with risk measures in production planning

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.

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. According to [21], this scenario fits into the context of parametric variability which is one possible source of uncertainty. In order to measure the performance of a production regarding different production plans, risk-adjusted performance measures might be used. They are formally given as a mapping from the expected outcome (e.g. profit) and its corresponding risk to some value, see [12, 26]. Hence, measuring the performance of a stochastic production model is related to the quantification of risk or uncertainty since the risk is caused by equipment failure being one part of business risk.

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, 25] and production planning of discrete production models with randomness in [20]. Obviously, there is a high interest in the definition of performance measures that support planning decisions and also allow for optimization purposes, see [23, 24]. However, we remark that the model under consideration is only academic so far. If production data is available, continuous production models can be adapted to real production purposes using data-fitted approaches as for instance done in [11, 14].

The stochastic production network model under consideration has been originally introduced in [16] and is based on the deterministic production network model from [6]. There have been previous stochastic extensions in [15, 17] of the deterministic model but in [16] 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 need performance measures, i.e. measures that translate the stochastic outcome into reasonable quantities, to validate decisions. The extended numerical study can be also used as a tool to support decision making processes.

Considering monetary quantities, appropriate performance measures have been introduced in the literature of finance [2, 10, 27] 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 [5], risk measures have been introduced for the optimization of oil production. There, 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 Sect. 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 Sect. 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 (PDE) and ordinary (ODE) differential equations. Different to the existing literature on flow models, continuous production models (PDE and ODE) allow for the description of transient production processes. 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 [6, 13] and its stochastic extension to a load-dependent model from [16]. 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 \geq0$. 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

$$ \partial_{t} \rho(x,t) + \partial_{x} \min\bigl\{ v \rho(x,t),\mu\bigr\} =0, $$

(1)

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

$$ \partial_{t} q(t) = g_{\mathrm{in}}(t)-g_{\mathrm{out}}(t), $$

(2)

i.e. as the balance between inflow into and outflow. Here, $g_{\mathrm{in}}(t)$ is the inflow into the queue, which can be an externally given inflow $G_{\mathrm{in}}(t)$ in case 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 μ, and if the queue is empty, the processor can take the inflow into the queue, which is bounded by the maximal capacity μ. Summarizing, this reads as

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

The coupling of the processor and the corresponding queue is described by a boundary condition $\rho(a,t) = \frac{g_{\mathrm{out}}(t)}{v}$ and initial conditions $\rho(x,0) = \rho_{0}(x) \in L^{1}((a,b))$, $q(0) = q_{0} \in \mathbb {R}_{\geq0}$ 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.

Figure 1

Production dynamics on one edge

We now describe the extension to the network case, see [6]. 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} \geq0$ for every arc (or processor) $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, \ldots,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_{\mathrm{in}} = \{v \in \mathcal {V}\colon\delta_{v}^{-} = \emptyset\}$, we define a time-dependent inflow function $G^{v}_{\mathrm{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 to the subsequent processors. They satisfy $\sum_{e \in\delta_{v}^{+}} A ^{v,e} (t) = 1$.

An example of a production network is shown in Fig. 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}$.

Figure 2

Diamond network with seven processors

Letting $s(e)$ denote the starting node for a given arc e, we adapt the flow into and out of the queue of processor e as follows:

$$ g_{\mathrm{in}}^{e}(t) = \textstyle\begin{cases} A^{s(e),e}(t)\sum_{\tilde{e} \in\delta_{s(e)}^{-}} \min\{v^{ \tilde{e}} \rho^{\tilde{e}}(b^{\tilde{e}},t),\mu^{\tilde{e}}\} &\text{if } s(e) \notin V_{\mathrm{in}}, \\ G_{\mathrm{in}}^{s(e)}(t) &\text{if } s(e) \in V_{\mathrm{in}}, \end{cases} $$

and

$$ g_{\mathrm{out}}^{e}(t) = \textstyle\begin{cases} \min\{g_{\mathrm{in}}^{e}(t),\mu^{e}\} &\text{if } q^{e}(t) = 0, \\ \mu^{e} &\text{if } q^{e}(t)>0. \end{cases} $$

This deterministic model is well-defined if the initial densities and the inflow functions are of bounded total variation, see [6]. 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},\ldots,q _{0}^{N},\rho_{0}^{1},\ldots,\rho_{0}^{N})$ at time $t\geq s$ and with capacities $\mu= (\mu^{1},\ldots,\mu^{N})$, see [4, 6]. This deterministic flow $S^{\mu}$ will determine the deterministic evolution between the machine failures or capacity drops. As already mentioned in the introduction, due to the missing data, we focus on the mathematical description of the stochastic model in the following section.

2.2 Stochastic extension: load-dependent capacities

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

$$\begin{aligned}& \varPhi_{st} \colon E \to E, \\& (r_{0},q_{0},\rho_{0}) \mapsto \bigl(r(t),q(t),\rho(t)\bigr), \end{aligned}$$

with

$$ E = \bigl\{ 1,\ldots,C^{1}\bigr\} \times\cdots\times\bigl\{ 1, \ldots,C^{N}\bigr\} \times \mathbb {R}_{\geq0}^{N} \times L^{1}\bigl((a_{1},b_{1})\bigr)\times\cdots \times L ^{1}\bigl((a_{N},b_{N})\bigr), $$

where

$$ r(t) = r_{0}, \qquad \bigl(q(t),\rho(t)\bigr) = S^{\mu(r_{0})}_{st}(q_{0}, \rho _{0}). $$

To incorporate random capacity drops in the production network, we follow the theory of piecewise deterministic Markov processes, see e.g. [7, 19]. Since we have the deterministic evolution between the jump times given by Φ, we only need to specify the intensity ψ at which jumps occur and the distribution of the jumps η. 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

$$ \lambda _{ii}^{e} = \sum_{j\neq i}^{C^{e}} \lambda _{ij}^{e}. $$

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

$$\begin{aligned}& \psi(t,y) = \sum_{e=1}^{N} \lambda ^{e}_{r_{e}r_{e}}\bigl(t,(q_{e}, \rho_{e})\bigr), \\& \eta(t,y,B) = \sum_{e=1}^{N}\sum _{l\neq r_{e}}^{C^{e}} \frac{\lambda ^{e}_{r_{e} l}(t,(q^{e},\rho^{e}))}{\psi(t,y)} \varepsilon_{(r_{1},\ldots,r_{e-1},l,r_{e+1},\ldots,r_{N},q,\rho)}, \end{aligned}$$

where $\varepsilon_{x}$ is the Dirac measure with unit mass in x and

$$\begin{aligned} \sigma(E) =& \mathcal {P}\bigl(\bigl\{ 1,\ldots,C^{1}\bigr\} \bigr)\otimes \mathcal {P}\bigl( \bigl\{ 1,\ldots,C^{N}\bigr\} \bigr)\otimes \mathcal {B}\bigl(\mathbb {R}^{N}_{\geq0} \bigr) \\ &{}\otimes\sigma \bigl(L^{1}\bigl((a_{1},b_{1}) \bigr)\bigr)\otimes\sigma\bigl(L^{1}\bigl((a_{N},b_{N}) \bigr)\bigr) \end{aligned}$$

is the σ-algebra on E. Here, $\mathcal {P}$ denotes the power set, $\mathcal {B}(\mathbb {R}^{N}_{\geq0})$ the Borel σ-algebra on $\mathbb {R}^{N}_{\geq 0}$ and $\sigma(L^{1}((a_{k},b_{k})))$ the Borel σ-algebra generated by the open sets induced by the $L^{1}$-norm. 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 $(\varOmega,\mathcal{F},P)$, which is piecewise deterministic between the jumps and follows the deterministic production network equations, see [16].

To construct sample paths of the stochastic process Y, we use a numerical scheme for the deterministic evolution Φ 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} \geq0$ be the time of the nth 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 Fig. 3.

Figure 3

Thinning algorithm

For the description of the procedure we assume a uniform bound λ̄ on ψ. 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 $\varPhi_{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 $\frac{\varPsi(T_{n}+\xi_{1},\varPhi_{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 η. 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 \geq0)$ 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

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

where $(X^{e}_{i}(t), t\geq0)$ 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. Fig. 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.

Figure 4

Worker allocation

The next question is how we can incorporate these ideas into the setting of the model presented in Sect. 2.1. The mathematical idea is to interpret every $(X_{i}^{e}(t), t \geq0)$ as a continuous time Markov Chain (CTMC) on a common probability space $(\varOmega,\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 Sect. 4.2. To simplify the notation, we neglect the index e in the following.

Lemma 2.1

Fix$N \in \mathbb {N}$and suppose we are given a family$(X_{1}(t),\ldots ,X_{N}(t),t\geq0)$of independently identically distributed CTMCs on the probability space$(\varOmega,\mathcal{F},P)$taking values in$\{0,1\}$. In addition, suppose we are given a$\mathcal {Q}$-matrix

$$ \mathcal {Q}= \begin{pmatrix} -\lambda _{0} & \lambda _{0} \\ \lambda _{1} & -\lambda _{1} \end{pmatrix} $$

for$\lambda _{0}, \lambda _{1}>0$. The stochastic process defined by

$$ X(t)= \sum_{i=1}^{N} X_{i}(t) $$

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

$$ q_{jk} = \textstyle\begin{cases} -(j\lambda _{1}+(N-j)\lambda _{0}) &\textit{if } k=j, \\ (N-j)\lambda _{0} &\textit{if } k = j+1, \\ j\lambda _{1} &\textit{if } k = j-1, \\ 0 &\textit{else} \end{cases} $$

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

$$ P\bigl(X(t) = j\bigr)= \binom{N}{j}P\bigl(X_{1}(t) = 1 \bigr)^{j}P\bigl(X_{1}(t) = 0\bigr)^{N-j} . $$

(3)

The proof can be found in the Appendix. 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.

$$ X(t) \sim\operatorname{Bin}\bigl(N,P\bigl(X_{1}(t)= 1\bigr)\bigr). $$

The steady-state distribution is consequently

$$ \lim_{t \to\infty} P\bigl(X(t) = j\bigr) = \binom{N}{j} \biggl(\frac{\lambda _{0}}{\lambda _{0}+\lambda _{1}} \biggr)^{j} \biggl(\frac{\lambda _{1}}{\lambda _{0}+\lambda _{1}} \biggr) ^{N-j}, $$

(4)

because of

$$ \lim_{t \to\infty} P\bigl(X_{1}(t) = 1\bigr) = \frac{\lambda _{0}}{\lambda _{0}+\lambda _{1}}. $$

The latter allows for a simple availability analysis of the capacities because the parameters $\lambda _{0}$ and $\lambda _{1}$ are known from estimations. Equation (4) describes the long term probabilities, which give an easy representation of the capacity one can expect in this production unit. The expected capacity is $\mathbb {E}[X(\infty)] = N \frac{\lambda _{0}}{\lambda _{0}+\lambda _{1}}$, which is useful in practical applications to identify bottlenecks in the production.

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

$$ \lambda _{ij}^{e}\bigl(t,\bigl(q^{e}, \rho^{e}\bigr)\bigr) = \textstyle\begin{cases} j\lambda _{1}^{e}+(N^{e}-j)\lambda _{0}^{e} &\text{if } j=i, \\ (N^{e}-j)\lambda _{0}^{e} &\text{if } j = i+1, \\ j\lambda _{1}^{e} &\text{if } j = i-1, \\ 0 &\text{else} \end{cases} $$

for $i,j \in\{1,\ldots,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 are 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. Performance and in particular risk measures are statistical measures that can be used as predictors of investment risk and volatility in many applications. Even though such measures are mostly applied in finance, they can be also used to evaluate production systems.

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 [15, 17], the aggregated outflow until time $t\geq0$ of the complete network can be computed as

$$ G^{\mathrm{net}}_{\mathrm{out}}(t) = \int_{0}^{t} \sum_{v \in V_{ \mathrm{out}}} \sum_{e \in\delta^{-}_{v}} \min\bigl\{ v^{ e} \rho^{ e}\bigl(b ^{ e},s\bigr),\mu^{e} \bigl(r^{e}(s)\bigr)\bigr\} \,ds, $$

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

$$ q^{\mathrm{net}}(t) = \sum_{e\in \mathcal {A}} \int_{0}^{t} q^{e}(s) \,ds, $$

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

$$ \mathbb {E}\bigl[G^{\mathrm{net}}_{\mathrm{out}}(t)\bigr], \qquad \mathbb {E}\bigl[q^{\mathrm{net}}(t) \bigr], $$

and variance

$$ \sigma^{2}\bigl(G^{\mathrm{net}}_{\mathrm{out}}(t)\bigr),\qquad \sigma^{2}\bigl(q^{ \mathrm{net}}(t)\bigr) $$

of these quantities. They only measure the outflow or the network queue loads and do not include any information on the profit of the company.

3.2 Risk measures in performance measures

If we consider the profit until time t as a random variable $\varPi(t)$ (a functional of Y), where Y is a stochastic production network model, then we can include monetary aspects. We could also use not risk-adjusted performance measures, as e.g. the expectation for the random variable $\varPi(t)$, to describe the performance of the production. However, it turns out that they are not the best choice especially in the context of optimization, see Sect. 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.,

$$ P\bigl(\varPi(t)< 0\bigr), $$

which has different unit (probability) than the profit (monetary unit). 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. [2, 10, 27], 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 [10].

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

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

• 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

• 3. for every $\lambda \geq0$, we have $\varrho(\lambda X) = \lambda \varrho(X)$    (Positive homogeneity)

• 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 risk-free surplus of m, then the risk is reduced by m, see property 2. Property 3 induces a normalization, i.e. $\varrho(0) = 0$ and ensures an easy scaling of different currencies. Property 4 describes how aggregation leads to a reduction of risk.

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

$$ \operatorname{V@R}_{\lambda }(X) = \inf\bigl\{ m \in \mathbb {R}\colon P(X+m< 0)\leq \lambda \bigr\} $$

(5)

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

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

(6)

with the upper quantile function

$$ q^{+}_{X}(t) = \sup\bigl\{ x \in \mathbb {R}\colon P(X< x) \leq t \bigr\} . $$

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

$$ \operatorname{AV@R}_{\lambda }(X) = \frac{1}{\lambda } \int_{0}^{\lambda }\operatorname{V@R}_{\gamma}(X) \,d\gamma. $$

From the computational point of view, we can estimate the Value at Risk (6) from a sample of profit realizations by an estimation of the quantile function such as quantile¹ in Matlab. This can again be used to compute the Average Value at Risk with, e.g. the Matlab function integral².

Remark 3.2

Following [10], AV@R can also be expressed as

$$ \operatorname{AV@R}_{\lambda }(X) = \frac{1}{\lambda } \inf _{t \in \mathbb {R}} \bigl(\mathbb {E}\bigl[(t-X)^{+}\bigr]-\lambda t \bigr) = \frac{1}{\lambda } \mathbb {E}\bigl[(q-X)^{+}\bigr]-q, $$

where $(z)^{+} = \max\{0,z\}$ and q is a λ-quantile (e.g. $q = q^{+}_{X}(\lambda )$) of X. This yields for continuous distributed X an alternative approach from a computational point of view since solving $\frac{1}{\lambda } \inf_{t \in \mathbb {R}} (\mathbb {E}[(t-X)^{+}]-\lambda t )$ directly leads to the $\operatorname{V@R}_{\lambda }(X)$ and the $\operatorname{AV@R}_{\lambda }(X)$.

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. The simulation results should help decision makers to understand risks, uncertainties and complexities. Therefore, all results are described and interpreted in detail.

4.1 Distribution parameter planning in the load-dependent case

Let the topology of a production network be given as a diamond network, see Fig. 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 Sect. 3.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

$$ \varPi(t) = \int_{0}^{t} \biggl(\sum _{v \in V_{\mathrm{out}}} \sum_{e \in\delta^{-}_{v}} \min\bigl\{ v^{e}\rho^{e}\bigl(b^{e},s\bigr), \mu^{e}\bigl(r ^{e}(t)\bigr)\bigr\} \cdot p(s) -\sum _{e\in \mathcal {A}} q^{e}(s)\cdot C_{q}^{e}(s) \biggr)\,ds. $$

(7)

The price of the product is $p(s) \geq0$ and $C_{q}^{e}(s) \geq0$ 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,\ldots ,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 [16] and define the Utilization Ratio UR by

$$ \operatorname{UR}^{e}\bigl(r^{e},q^{e}, \rho^{e}\bigr) = \frac{1}{\max_{i}\{ \mu ^{e}(i)\}(b^{e}-a^{e})} \int_{a^{e}}^{b^{e}}\min\bigl\{ \mu^{e} \bigl(r^{e}(t)\bigr),v ^{e}\rho^{e}(x,t)\bigr\} \,dx, $$

and the Ratio of Work In Progress and the maximal amount of goods in the machine RWIP as

$$ \operatorname{RWIP}^{e} \bigl(r^{e},q^{e}, \rho^{e}\bigr) = \frac{v^{e}}{\max_{i}\{\mu^{e}(i)\}(b^{e}-a^{e})} \int_{a^{e}}^{b^{e}} \rho^{e}(x,t)\,dx. $$

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

$$\begin{aligned}& \lambda _{12}^{e}\bigl(t,q^{e}, \rho^{e}\bigr) = \lambda ^{\text{rep,max,e}}-\bigl(\lambda ^{ \text{rep,max,e}}- \lambda ^{\text{rep,min,e}}\bigr)\operatorname{RWIP}^{e}\bigl(1,q ^{e},\rho^{e}\bigr), \\& \lambda _{21}^{e}\bigl(t,q^{e}, \rho^{e}\bigr) = \lambda ^{\text{down},e} \operatorname{UR}^{e} \bigl(2,q^{e},\rho^{e}\bigr), \end{aligned}$$

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

$$\begin{aligned}& \begin{aligned} \lambda _{13}^{e}\bigl(t,q^{e}, \rho^{e}\bigr) &= \lambda _{23}^{e} \bigl(t,q^{e},\rho^{e}\bigr) \\ &= \lambda ^{\text{rep,max,e}}-\bigl(\lambda ^{\text{rep,max,e}}-\lambda ^{ \text{rep,min,e}}\bigr) \operatorname{RWIP}^{e}\bigl(1,q^{e},\rho^{e} \bigr), \end{aligned} \\& \lambda _{21}^{e}\bigl(t,q^{e}, \rho^{e}\bigr) = \lambda ^{\text{down},e} \operatorname{UR} \bigl(2,q^{e},\rho^{e}\bigr), \\& \lambda _{31}^{e}\bigl(t,q^{e}, \rho^{e}\bigr) = \lambda _{32}^{e} \bigl(t,q^{e},\rho^{e}\bigr)=\lambda ^{\text{down},e} \operatorname{UR}^{e}\bigl(3,q^{e},\rho^{e} \bigr) \end{aligned}$$

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 cannot state an explicit stationary expected capacity as in Sect. 4.2 and simulations are needed to determine which combination $(\alpha_{1},\alpha_{2})$ performs in a “optimal” way.

Figure 5

Value and Average Value at Risk with level $\lambda= 0.1$ for 1000 samples

In Fig. 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,\ldots,0.9,1\}$ are evaluated. In Figs. 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$.

Figure 6

Sample mean, standard deviation and bankruptcy probability of profit for 1000 samples

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 Fig. 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.

Table 1 Comparison of the different profit evaluations

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 Fig. 7, where capacities are given as in Sect. 2.3 and we define the profit $\varPi(t)$ as the cumulative revenue reduced by storage and cluster-size costs up to time $t \geq0$.

Figure 7

Serial network with two processors

In formulas, this reads as

$$\begin{aligned} \varPi(t) =& \int_{0}^{t} \biggl(\sum _{v \in V_{\mathrm{out}}} \sum_{e \in\delta^{-}_{v}} \min\bigl\{ v^{e}\rho^{e}\bigl(b^{e},s\bigr), \mu^{e}\bigl(r ^{e}(t)\bigr)\bigr\} \cdot p(s) \\ &{}-\sum _{e\in \mathcal {A}}\bigl(q^{e}(s)\cdot C_{q}^{e}(s)+N ^{e} \cdot C_{N}^{e}(s)\bigr) \biggr) \,ds, \end{aligned}$$

where we denote with $p(s) \geq0$ the price of the product, $C_{q}^{e}(s) \geq0$ and $C_{N}^{e}(s) \geq0$ are the storage and cluster size costs, respectively, at time s. Typical examples for the cluster size cost are maintenance costs or salaries. If ϱ is a monetary measure of the risk, we can interpret $\varrho(\varPi(T))$ as the risk of the company of the aggregated loss and gains up to 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,\ldots,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

$$ \lambda _{0}^{1} = \frac{1}{10},\qquad \lambda _{1}^{1} = \frac{1}{80}, \qquad \lambda _{0}^{2} = \frac{1}{20}, \qquad \lambda _{1}^{2} = \frac{1}{50}. $$

We can compute the expected capacities in steady state as

$$ \mathbb {E}\bigl[\mu^{e}(\infty)\bigr] = \sum_{n=0}^{N^{e}}n \binom{N^{e}}{n} \biggl(\frac{\lambda _{0}^{e}}{\lambda _{0}^{e}+\lambda _{1}^{e}} \biggr)^{n} \biggl(\frac{\lambda _{1} ^{e}}{\lambda _{0}^{e}+\lambda _{1}^{e}} \biggr)^{N^{e}-n} = N^{e} \frac{\lambda _{0}^{e}}{\lambda _{0}^{e}+\lambda _{1}^{e}}, $$

which implies

$$ \mathbb {E}\bigl[\mu^{1}(\infty)\bigr] = N^{1} \frac{8}{9}, \qquad \mathbb {E}\bigl[\mu^{2}(\infty)\bigr] = N^{2} \frac{5}{7}, $$

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_{\mathrm{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 $\varPi(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).

Figure 8

Profit evaluation for different cluster sizes

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 \cdot10^{3}$, as Table 2 states.

Table 2 Comparison of the different profit evaluations for $(N^{1},N^{2})$

With respect to the standard deviation, we obtain the combination $N^{1} = N^{2} = 15$ having the lowest standard deviation of $0.6733\cdot10^{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 Fig. 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$, Fig. 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 \cdot10^{3}$. This can be interpreted as follows: even if we have debts of $-0.2283 \cdot10^{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}(\varPi(T))$ of $0.6255 \cdot10^{3}$, i.e., we should have a surplus of $0.6255 \cdot10^{3}$ for these cluster sizes.

Figure 9

Value and Average Value at Risk for different cluster sizes and level $\lambda = 0.1$

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}(\varPi(T))$ of $1.5683 \cdot10^{3}$; see Fig. 10(a) and Table 2. The combination $N^{1} = 7$ and $N^{2} = 9$ implies an $\operatorname{AV@R}_{0.01}(\varPi(T))$ of $2.0632 \cdot10^{3}$; see Fig. 10(b).

Figure 10

Value and Average Value at Risk for different cluster sizes and level $\lambda = 0.01$

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. In particular, from a practical point of view, a major benefit of our approach is a comprehensive numerical analysis of a fully time-dependent stochastic production model that keeps track of random capacity reductions. The combined study with risk measures also allows for a performance evaluation in a simulation environment. So far, the model is limited to the consideration of homogeneous products but could be extended to multicommodity flows in a straightforward way. Furthermore, the well-known buffer allocation problem [8, 22, 28] used for the design of manufacturing systems can be analyzed within this framework to find optimal buffer sizes.

Therefore, future work will deal with the formulation of rigorous optimization problems for the stochastic production network model with respect to relevant performance measures. Mathematically, a focus will be the study of the quality of solutions and the speed of convergence to an optimum. We also hope to get more insight in the decision making for more realistic flow lines.

Appendix:  Proof of Lemma 2.1

Proof

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

$$\begin{aligned} \psi\colon E^{Y}\to E^{X},\qquad x \mapsto \sum_{i=1}^{N} x _{i} \end{aligned}$$

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

$$ U_{t}(x,B) = P^{N}\bigl(Y(t+s)\in B|Y(s)=x\bigr) $$

for every $x \in E^{Y}$, $s,t\geq0$ 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

$$ U_{t}\bigl(x,\psi^{-1}(\varGamma)\bigr) = U_{t}\bigl(\tilde{x},\psi^{-1}(\varGamma)\bigr) $$

for every $t\geq0$ and $\varGamma\in \mathcal {E}^{X}$ holds (Dynkin’s criterion), then we know that the stochastic process $(X(t),t\geq0)$ 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 σ because of the structure of ψ. Therefore, we can compute

$$\begin{aligned}& U_{t}\bigl(x,\psi^{-1}(\varGamma)\bigr) \\& \quad = \sum_{z \in\psi^{-1}(\varGamma)}P^{N}\bigl(Y(t) = z|Y(0) = x\bigr) \\& \quad = \sum_{z \in\psi^{-1}(\varGamma)}\prod_{i=1}^{N} P\bigl(X_{i}(t) = z _{i}|X_{i}(0) = x_{i}\bigr) \\& \quad = \sum_{z \in\psi^{-1}(\varGamma)}\prod_{i=1}^{N} P\bigl(X_{1}(t) = z _{i}|X_{1}(0) = x_{i}\bigr) \\& \quad = \sum_{z \in\psi^{-1}(\varGamma)}\prod_{i=1}^{N} P\bigl(X_{1}(t) = z _{\sigma(i)}|X_{1}(0) = x_{\sigma(i)}\bigr) \\& \quad = \sum_{z_{\sigma} \in\psi^{-1}(\varGamma)}\prod_{i=1}^{N} P\bigl(X _{i}(t) = z_{\sigma(i)}|X_{i}(0) = y_{i}\bigr) \\& \quad = \sum_{z_{\sigma} \in\psi^{-1}(\varGamma)}P^{N}\bigl(Y(t) = z_{\sigma }|Y(0) = y\bigr) \\& \quad = U_{t}\bigl(y,\psi^{-1}(\varGamma)\bigr) \end{aligned}$$

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

$$P(X_i(t+\mathrm{\Delta }t)=m|X_i(t)=l)={\mathbb{1}}_l(m)(1-\mathrm{\Delta }t{\lambda }_l)+{\mathbb{1}}_{1-l}(m)\mathrm{\Delta }t{\lambda }_l+\mathrm{o}(\mathrm{\Delta }t)$$

(8)

in the limit $\Delta t\to0$; see [18, pp. 28].

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

$$\begin{aligned}& P\bigl(X(t+\Delta t) = k,X(t) = j\bigr) \\& \quad = \sum_{\substack{\alpha \in C\\ \vert \alpha \vert = k}} \sum_{\substack{\beta \in C\\ \vert \beta \vert = j}}P \bigl(X_{1}(t+\Delta t) = \alpha _{1},\ldots,X_{N}(t+ \Delta t)=\alpha _{N}, X_{1}(t) = \beta _{1},\ldots,X _{N}(t) = \beta _{N} \bigr) \\& \quad = \sum_{\substack{\alpha \in C\\ \vert \alpha \vert = k}} \sum_{\substack{\beta \in C\\ \vert \beta \vert = j}} \prod_{i=1}^{N} P\bigl(X_{i}(t+ \Delta t) = \alpha _{i}|X_{i}(t) = \beta _{i} \bigr)P\bigl(X_{i}(t) = \beta _{i}\bigr) \end{aligned}$$

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

$$ \vert \alpha \vert = \sum_{i=1}^{N} \alpha _{i}. $$

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

$$\begin{array}{rl}& P(X(t+\mathrm{\Delta }t)=k,X(t)=j)\\ & =\sum _{\begin{array}{c}\alpha \in C\\ |\alpha |=k\end{array}}\sum _{\begin{array}{c}\beta \in C\\ |\beta |=j\end{array}}(\prod _{i=1}^N{\mathbb{1}}_{{\beta }_i}({\alpha }_i)(1-\mathrm{\Delta }t{\lambda }_{{\beta }_i})P(X_i(t)={\beta }_i))\end{array}$$

(9)

$$\begin{array}{rl}& +\sum _{\begin{array}{c}\alpha \in C\\ |\alpha |=k\end{array}}\sum _{\begin{array}{c}\beta \in C\\ |\beta |=j\end{array}}(\sum _{l=1}^N{\mathbb{1}}_{1-{\beta }_l}({\alpha }_l)\mathrm{\Delta }t{\lambda }_{{\beta }_l}P(X_l(t)={\beta }_l)\prod _{\begin{array}{c}i=1\\ i\ne l\end{array}}^N{\mathbb{1}}_{{\beta }_i}({\alpha }_i)(1-\mathrm{\Delta }t{\lambda }_{{\beta }_i})P(X_i(t)={\beta }_i))\\ & +\mathrm{o}(\mathrm{\Delta }t).\end{array}$$

(10)

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

$$\begin{array}{c}\sum _{\begin{array}{c}\alpha \in C\\ |\alpha |=k\end{array}}\sum _{\begin{array}{c}\beta \in C\\ |\beta |=j\end{array}}(\prod _{i=1}^N{\mathbb{1}}_{{\beta }_i}({\alpha }_i)(1-\mathrm{\Delta }t{\lambda }_{{\beta }_i})P(X_i(t)={\beta }_i))\hfill \\ ={\mathbb{1}}_j(k)\left(\sum _{\begin{array}{c}\beta \in C\\ |\beta |=j\end{array}}P(X_i(t)={\beta }_i)(1-\mathrm{\Delta }t\sum _{l=1}^N{\lambda }_{{\beta }_l})\right)+\mathrm{o}(\mathrm{\Delta }t)\hfill \\ ={\mathbb{1}}_j(k)(1-\mathrm{\Delta }t(j{\lambda }_1+(N-j){\lambda }_0))P(X(t)=j)+\mathrm{o}(\mathrm{\Delta }t).\hfill \end{array}$$

The last equality follows from

$$ \sum_{l=1}^{N} \lambda _{\beta _{l}} = \bigl(j\lambda _{1}+(N-j)\lambda _{0}\bigr) $$

since j entries of β are one and all the others are zero. Now, we analyze the second summand (10), 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 α has two entries more with value one, which implies that α and β are different in at least two entries. Hence, the product in the second summand is zero. We have

$$\begin{array}{c}\sum _{\begin{array}{c}\alpha \in C\\ |\alpha |=k\end{array}}\sum _{\begin{array}{c}\beta \in C\\ |\beta |=j\end{array}}\sum _{l=1}^N({\mathbb{1}}_{1-{\beta }_l}({\alpha }_l)\mathrm{\Delta }t{\lambda }_{{\beta }_l}P(X_l(t)={\beta }_l)\prod _{\begin{array}{c}i=1\\ i\ne l\end{array}}^N{\mathbb{1}}_{{\beta }_i}({\alpha }_i)(1-\mathrm{\Delta }t{\lambda }_{{\beta }_i})P(X_i(t)={\beta }_i))\hfill \\ =\mathrm{\Delta }t\sum _{l=1}^N\sum _{\begin{array}{c}\alpha \in C\\ |\alpha |=k\end{array}}\sum _{\begin{array}{c}\beta \in C\\ |\beta |=j\end{array}}({\mathbb{1}}_{1-{\beta }_l}({\alpha }_l){\lambda }_{{\beta }_l}P(X_l(t)={\beta }_l)\prod _{\begin{array}{c}i=1\\ i\ne l\end{array}}^N{\mathbb{1}}_{{\beta }_i}({\alpha }_i)P(X_i(t)={\beta }_i))+\mathrm{o}(\mathrm{\Delta }t)\hfill \\ ={\mathbb{1}}_{j+1}(k)\mathrm{\Delta }t{\lambda }_0\sum _{l=1}^N\sum _{\begin{array}{c}\beta \in C\\ |\beta |=j\\ {\beta }_l=0\end{array}}P{(X_1(t)=1)}^{j}P{(X_1(t)=0)}^{N-j}\hfill \\ +{\mathbb{1}}_{j-1}(k)\mathrm{\Delta }t{\lambda }_1\sum _{l=1}^N\sum _{\begin{array}{c}\beta \in C\\ |\beta |=j\\ {\beta }_l=1\end{array}}P{(X_1(t)=1)}^{j}P{(X_1(t)=0)}^{N-j}+\mathrm{o}(\mathrm{\Delta }t)\hfill \end{array}$$

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

$$\begin{aligned}& \sum_{l=1}^{N}\sum _{\substack{\beta \in C\\|\beta | = j\\\beta _{l} = 0}} 1 = N\cdot\binom{N-1}{j} = (N-j)\cdot \binom{N}{j}, \\& \sum_{l=1}^{N}\sum _{\substack{\beta \in C\\|\beta | = j\\\beta _{l} = 1}} 1 = N\cdot\binom{N-1}{j-1} = j\cdot \binom{N}{j} \end{aligned}$$

and observe that

$$ P\bigl(X(t) = j\bigr)= \binom{N}{j}P\bigl(X_{1}(t) = 1 \bigr)^{j}P\bigl(X_{1}(t) = 0\bigr)^{N-j}. $$

Summarizing all computations yields the transition probability

$$\begin{array}{rl}P(X(t+\mathrm{\Delta }t)=k|X(t)=j)=& {\mathbb{1}}_j(k)(1-\mathrm{\Delta }t(j{\lambda }_1+(N-j){\lambda }_0))\\ & +{\mathbb{1}}_{j+1}(k)\mathrm{\Delta }t(N-j){\lambda }_0\\ & +{\mathbb{1}}_{j-1}(k)\mathrm{\Delta }tj{\lambda }_1+\mathrm{o}(\mathrm{\Delta }t)\end{array}$$

as $\Delta t \to0$ and we conclude the generator of $(X(t),t\geq0)$ as the matrix $\mathcal {Q}$ from the statement of this lemma. □