Limiting the impact of supply chain disruptions in the face of distributional uncertainty in demand

Abstract

Service‐level requirements play a crucial role in eliminating stock‐outs in a production pipeline. However, delivering a specific service level can become an unattainable goal given the various uncertainties influencing both the production pipeline and customer demand and causing the manufacturer to adapt the initial strategy in response to disruptions. Such deviations from optimality frequently result in unexpected (and potentially very high) costs and are complex to manage. On the one hand, the manufacturer can use a robust or distributionally robust approach to prepare for the worst‐case disruption, ensuring that the realized cost will be lower than the estimated cost with high probability. On the other hand, this solution may lead to overly conservative production schedules. In this article, we take a different approach and develop a bilevel stochastic optimization model with chance constraints, which allows us to make the production more predictable in the event of disruptions by driving costs and optimal schedules closer to the benchmark for each scenario considered. We introduce doubly probabilistic service‐level requirements to account for two interdependent layers of uncertainty, that is, production disruptions and distributional uncertainties in customer demand. This allows us to make high‐quality production decisions with only a limited understanding of the demand pattern. Approximating the problem for a numerical solution, we guarantee tight optimality gaps for high service levels and propose an efficient solution scheme, combining robust scenario reduction with a customized Benders decomposition procedure. In the managerial section, we use the Omega× Swatch MoonSwatches example to demonstrate that a desirable doubly probabilistic service level can be attained for disruptions with a drop in demand. In the case of disruptions followed by a peak in demand, one can tighten the optimality gaps if the service level is reduced.

Keywords

Benders decomposition bilevel optimization production disruptions service level

INTRODUCTION AND LITERATURE REVIEW

Analytical models that incorporate internal uncertainty in the form of a random supply variable into the production optimization process were first studied by Arrow et al. (1958). Manufacturing decisions resilient to production disruptions lead to mitigation strategies with optimized safety stock, risk mitigation inventory, or alternatives for dual sourcing (Chopra & Sodhi, 2004; Lücker et al., 2020; Snyder et al., 2016; Tomlin, 2006). However, one of the main complexities in determining optimal quantities stems from the fact that production disruptions can affect demand (Hendricks et al., 2019; Timonina‐Farkas et al., 2020), driving it either down (e.g., in the case of damage to a manufacturer's reputation) or up (e.g., in the case of unprecedented events such as COVID‐19 or natural disasters). Thus, the production uncertainties and demand perturbations cannot be studied in isolation of each other, particularly in the presence of service‐level constraints (Craig et al., 2016; Katok et al., 2008). Considered as extremes, internal production disruptions cause manufacturers to reduce production capacity (Dong & Tomlin, 2012; Snyder et al., 2016), which, in turn, may influence demand (Craig et al., 2016).

The work of Ciarallo et al. (1994) took an important step forward by combining demand and capacity uncertainties into a single model, in which the demand follows a given probability distribution. However, during a disruption scenario, both demand and its distribution may remain unknown. Manufacturers are frequently unaware of the exact reason for a disruption before it occurs, which causes the results of scenario‐specific customer research to be biased and dependent on distributional assumptions. Nourelfath (2011) extended the literature to the case of robust service‐level requirements, where the probability of meeting a service level was evaluated using Wiener processes and a first passage time theory. From the retailer's perspective, the interdependence between demand and the optimal product assortment was studied in the work of Timonina‐Farkas et al. (2020). In our article, we develop a stochastic cost‐minimization problem with uncertain service‐level requirements extending the literature to the case where demand distributions depend on the realized disruption scenarios. Furthermore, via the notion of a regret function in our model, we control the cost difference between the manufacturer's strategy at the beginning of the planning horizon and its optimal adjustment in the event of a disruption. This also allows to bound the opportunity cost of knowing the disruption scenario a priori. To the best of our knowledge, our model is the first in the literature to combine these aspects and to account for underlying interdependencies, including those implied by random changes in demand distributions.

We propose a bilevel optimization model for determining the production schedule over a planning horizon and introduce doubly probabilistic service‐level constraints that account for two layers of uncertainty. From the manufacturer's perspective, the first layer of uncertainty is caused by the variability in demand and, moreover, the variability in demand distribution. The second layer of uncertainty arises from production disruptions, implying that there is a probability that the service level cannot be maintained. The two types of uncertainty cannot be viewed as independent of one another because production disruptions can originate outside of the firm (i.e., externally) for reasons that are directly linked to demand distribution perturbations (e.g., as during COVID‐19). We account for both types of uncertainty by utilizing chance constraints.

In our bilevel optimization model, the manufacturer addresses the upper‐level cost‐minimization problem and must evaluate the loss incurred under different disruption scenarios. The evaluation of this loss builds on the lower‐level optimization problem, which is dependent on an actual disruption. We introduce the loss incurred under a disruption scenario by means of a regret function that measures the cost difference between a reactive and an anticipative stochastic optimization problem. The anticipative problem optimizes the production strategy of a manufacturer who prepares for a specific production disruption at the beginning of the planning horizon and differentiates between full and partial anticipation. In particular, the manufacturer does not necessarily have the option of conducting customer research to assess the potential demand consequences of the disruption scenario because there could be multiple reasons for the same breakdown. In contrast, the reactive plan assumes that the decision maker only adapts production when the disruption starts. Because the cause of the breakdown is known at the time of the disruption, the manufacturer can estimate possible changes in the demand distribution. Importantly, our regret function contributes to the opportunity cost of knowing about the disruption beforehand. Thus, the benefit of the regret function is its ability to measure the cost impact of different decisions made by the manufacturer, for example, the decision to wait for the disruption before adapting the service level, or preparing for a specific disruption with or without additional customer research.

The resulting bilevel problem is a convex nonsmooth minimum cost network flow (MCF) problem (e.g., Prékopa & Boros, 1991; Zheng et al., 2015). Taking the regret and double probabilistic service‐level requirements into account, we propose an efficient numerical solution scheme that combines a robust scenario reduction approach with a customized Benders decomposition procedure. The solution approach is motivated by the work of Khang and Fujiwara (1993), which is devoted to dynamic minimum cost flow problems without disruptions, while the scenario reduction scheme is inspired by the works of Timonina et al. (2015) and Hochrainer‐Stigler et al. (2019), in which a complex network structure of European river basins is reduced to an ordered list of nodes. Introducing a severity measure and ranking our production scenarios from most to least distant ones according to the robustness principle, we propose an ordered solution scheme with a warm start, which speeds up the Benders decomposition procedure. By approximating the bilevel problem for a numerical solution, we guarantee tight optimality gaps between the initial and approximate problems for high service levels even when the true scenario‐specific demand distribution remains unknown. The guarantees for optimal solution proximity are derived from the results developed by Granot and Veinott Jr (1985), who provide a bound on the distance between optimal solutions of two convex MCF problems. In the managerial section of this article, we demonstrate that the best approximation quality with a desirable doubly probabilistic service level can be attained for disruptions with a drop in demand. In the case of disruptions with a peak in demand, tighter optimality gaps can be guaranteed if the service level is reduced. Our model encompasses a variety of practical applications due to the fact that our solution method is gradient‐free and, thus, a zero‐order method.

This article is organized as follows: The theoretical framework is developed in Section 2. First, the initial model and its approximation are presented for the case without disruptions, followed by a description of the extension to production disruptions with uncertain demand distributions. In Section 3, we develop the robust scenario reduction scheme and present a customized Benders decomposition with novel feasibility cuts and warm‐start techniques. In Section 4, we demonstrate the managerial implications of doubly probabilistic service‐level requirements in the optimization problem with production disruptions.

MATHEMATICAL MODEL

Consider a manufacturer deciding on the production quantity

P_{t}

in a multistage setup with discrete time

t \in {1, …, T}

. Let

κ_{t}

be the production capacity in period t, and

C_{t}^{prod} (\cdot)

denote a convex but not necessarily differentiable function determining the cost of producing

P_{t}

units. The manufacturer faces uncertainty in demand, which is a random variable denoted by

D_{t}

\forall t

. The inventory built up to period t and available in the beginning of period

t + 1

is equal to

I_{t} = \sum_{u = 1}^{t} P_{u} - \sum_{u = 1}^{t} D_{u}

. Some linear inventory cost

h_{t}

is to be paid if

I_{t}

is nonnegative. Conversely, a negative

I_{t}

represents a backlog with a linear cost

b_{t}

. We set

h_{T} = b_{T} = 0

and denote the starting inventory by I₀.

Next, we let

P_{t}^{c}

be the cumulative production and

D_{t}^{c}

the cumulative demand up to period t, that is,

P_{t}^{c} = \sum_{u = 1}^{t} P_{u}

and

D_{t}^{c} = \sum_{u = 1}^{t} D_{u}

. Working in a purely distributional setup, we denote the probability density and the cumulative distribution functions of the cumulative demand by

f_{t}

and

F_{t}

, respectively. If demands

(D_{1}, D_{2}, …, D_{T})

are independent and normally distributed random variables, the distribution of their sum is Gaussian. If the demands are independent but not necessarily identically distributed, one can use independent convolution techniques to estimate their joint density. Furthermore, coupling methods can be employed to estimate the density of the sum of dependent random variables (Hochrainer‐Stigler et al., 2019; Timonina‐Farkas et al., 2020). Such methods should be used in the presence of seasonality issues or product life cycle patterns that imply that the independence condition should be relaxed. To account for service‐level conditions, we specify a type 1 service‐level requirement in each period t. To do this, we bound the maximal allowed probability of stock‐out by

α_{t}

. In our setup, the decision maker aims to determine optimal production amounts

P_{t}

that minimize the total cost while not violating the probability of a stock‐out

α_{t}, \forall t

. The total cost in each period consists of the production cost

C_{t}^{prod} (P_{t})

, expected inventory cost

h_{t} E_{F_{t}} {(P_{t}^{c} - D_{t}^{c})}^{+}

, and expected backlogging cost

b_{t} E_{F_{t}} {(D_{t}^{c} - P_{t}^{c})}^{+}

, where

{(x)}^{+} = \max {0, x}

. The problem with the stage‐wise cost function

C_{t}^{orig} = C_{t}^{prod} (P_{t}) + h_{t} E_{F_{t}} {(P_{t}^{c} - D_{t}^{c})}^{+} + b_{t} E_{F_{t}} {(D_{t}^{c} - P_{t}^{c})}^{+}

, thus, yields:

\begin{matrix} v^{orig} = \min_{P, P^{c}} & \sum_{t = 1}^{T} C_{t}^{o r i g} \\ s . t . & \sum_{u = 1}^{t} P_{u} = P_{t}^{c}, 0 \leq P_{t} \leq κ_{t}, \forall t \in 1, …, T \\ P (P_{t}^{c} \leq D_{t}^{c}) \leq α_{t}, \forall t \in 1, …, T, \end{matrix}

where the sequences

P = (P_{1}, …, P_{T})

and

P^{c} = (P_{1}^{c}, …, P_{T}^{c})

determine the production schedules, and chance constraints are used to ensure that the probability of stock‐out does not exceed

α_{t}

. Problem (1) has similarities with the model proposed in Khang and Fujiwara (1993) and generalizes its production costs to arbitrary convex functions accounting for nonzero backlogging costs. The problem cannot be addressed in a closed form due to its variational structure. As an approximation, consider optimization problem (2) with deterministic demands, where the mean

μ_{t} = E_{F_{t}} (D_{t}^{c})

is used instead of random variables

D_{t}^{c}

in the stage‐wise cost function

C_{t}^{approx} = C_{t}^{prod} (P_{t}) + h_{t} {(P_{t}^{c} - μ_{t})}^{+} + b_{t} {(μ_{t} - P_{t}^{c})}^{+}

\begin{matrix} v^{approx} = \min_{P, P^{c}} & \sum_{t = 1}^{T} C_{t}^{a p p r o x} \\ s . t . & \sum_{u = 1}^{t} P_{u} = P_{t}^{c}, 0 \leq P_{t} \leq κ_{t}, \forall t \in 1, …, T \\ P (P_{t}^{c} \leq D_{t}^{c}) \leq α_{t}, \forall t \in 1, …, T . \end{matrix}

Note that problem (2) without chance constraints would be equivalent to an MCF problem (Bertsekas, 1998). Furthermore, the problem including chance constraints remains an MCF problem because the chance constraint of the form

P (P_{t}^{c} \leq D_{t}^{c}) \leq α_{t}

can be represented as a lower bound constraint

P_{t}^{c} \geq F_{t}^{- 1} (1 - α_{t})

where

F_{t}^{- 1} (1 - α_{t})

, also denoted by

l_{t}

, is the

(1 - α_{t})

‐th quantile of the distribution

F_{t}

and can be computed before solving the optimization problem (e.g., Vajda, 2014). We refer to this set of constraints as probabilistic service‐level requirements. The following result holds for optimization problems (1) and (2): Theorem 1

Let

(P^{*}, P^{c *})

be the optimal flow of problem (2). Then, the feasible regions of problems (1) and (2) coincide, and

(Lower Bound) Problem (2) imposes a lower bound on problem (1), that is,

v^{orig} \geq v^{approx}

;

(Upper Bounds) The following two upper bounds hold:

\begin{matrix} v^{orig} \leq v^{approx} + \sum_{t = 1}^{T} (h_{t} + b_{t}) \int_{P_{t}^{c *}}^{+ \infty} D_{t}^{c} f_{t} (D_{t}^{c}) d D_{t}^{c}, \end{matrix}

\begin{matrix} v^{orig} \leq v^{approx} + \sum_{t = 1}^{T} (h_{t} + b_{t}) \int_{0}^{P_{t}^{c *}} F_{t} (D_{t}^{c}) d D_{t}^{c} . \end{matrix}

Moreover, the optimality gap

(v^{o r i g} - v^{a p p r o x}) \geq 0

can be bounded by Δ, where

\begin{matrix} Δ = \sum_{t = 1}^{T} (h_{t} + b_{t}) \min {\int_{P_{t}^{c *}}^{+ \infty} D_{t}^{c} f_{t} (D_{t}^{c}) d D_{t}^{c}, \int_{0}^{P_{t}^{c *}} F_{t} (D_{t}^{c}) d D_{t}^{c}} . \end{matrix}

(Relative Error) The relative error of the approximation follows from statement 2, that is,

\begin{matrix} \frac{v^{orig} - v^{approx}}{v^{approx}} \leq \frac{Δ}{\sum_{t = 1}^{T} (C_{t}^{prod} (P_{t}^{*}) + h_{t} {(P_{t}^{c *} - μ_{t})}^{+} + b_{t} {(μ_{t} - P_{t}^{c *})}^{+})} . \end{matrix}

Proof

To prove this theorem, we note that the feasible regions of problems (1) and (2) coincide. Therefore, the first statement of the theorem is a direct consequence of Jensen's inequality for convex functions, and it also results from two lower bounds stated in Appendix A.1 in the Supporting Information. To derive the upper bounds, we consider an optimal flow

(P^{*}, P^{c *})

of problem (2). By evaluating the objective function of (1) at this point, we obtain a value larger or equal to

v^{orig}

and conclude the second statement of the theorem. See Appendix A.1 in the Supporting Information for the detailed proof.

□

The upper bound (3) is tighter than the upper bound (4) if

P_{t}^{c *} F_{t} (P_{t}^{c *}) > μ_{t} \forall t

. The optimal production

P_{t}^{c *}

grows if the service level

1 - α_{t}

increases due to the constraint

P_{t}^{c} \geq F_{t}^{- 1} (1 - α_{t})

. Therefore, the upper bound (3) is tight if the service‐level requirements are high. Conversely, the upper bound (4) is tight if the service‐level requirements are low, that is,

\begin{matrix} Δ = \{\begin{matrix} \sum_{t = 1}^{T} (h_{t} + b_{t}) \int_{P_{t}^{c *}}^{+ \infty} D_{t}^{c} f_{t} (D_{t}^{c}) d D_{t}^{c}, & if P_{t}^{c *} F_{t} (P_{t}^{c *}) \geq μ_{t} \forall t \\ \sum_{t = 1}^{T} (h_{t} + b_{t}) \int_{0}^{P_{t}^{c *}} F_{t} (D_{t}^{c}) d D_{t}^{c}, & if P_{t}^{c *} F_{t} (P_{t}^{c *}) < μ_{t} \forall t \\ Equation (5), & otherwise. \end{matrix} \end{matrix}

Figure 1 shows the bound Δ and the area to which the optimality gap belongs. The estimate is provided for Gaussian and lognormal demands with

μ = 5

σ = 1

h_{t} = b_{t} = 0.5

t = T = 1

FIGURE 1

Bounds for the optimality gap

As expected, the optimality gap is bounded tightly for cases with high or low service‐level requirements. The bound is least tight at point

P_{t}^{c *} F_{t} (P_{t}^{c *}) = μ_{t} \forall t

. In this article, we are primarily interested in the case

P_{t}^{c *} F_{t} (P_{t}^{c *}) > > μ_{t} \forall t

, which represents high service levels. Nevertheless, in the managerial Section 4, we also demonstrate the implications of a modeling error in the demand distribution

F_{t}

leading to a misfit between the realized and desired service levels.

Production disruptions

In this section, we extend the initial model (1) and its approximation (2) to the case of production disruptions and distributional perturbations in demand. The discrete set of all production disruptions is denoted by Ω, while each particular scenario is denoted by

ω \in Ω

. Introducing a change in the convex cost function

C_{t}^{prod} (P_{t})

to a convex function

{\bar{C}}_{t, ω}^{prod} (P_{t}, ω)

incurred during time period t of a disruption ω, we define a scenario ω as a set of periodic cost functions

{\bar{C}}_{t, ω}^{prod} (P_{t}, ω)

, which may differ from

C_{t}^{prod} (P_{t})

between periods

t = t_{ω}^{start}

and

t = t_{ω}^{end} \geq t_{ω}^{start}

. By this, the set Ω includes scenarios with multiple consecutive breaks, that is, those ω for which

\exists t : {\bar{C}}_{t, ω}^{prod} (P_{t}, ω) = C_{t}^{prod} (P_{t})

and

t_{ω}^{start} < t < t_{ω}^{end}

. Also, the set Ω includes a trivial scenario with no disruptions, that is,

ω_{0} \in Ω

such that

{\bar{C}}_{t, ω}^{prod} (P_{t}, ω) = C_{t}^{prod} (P_{t}), \forall t

. The probability of a scenario ω is denoted by

p_{ω} \in [0, 1]

with

\sum_{ω \in Ω} p_{ω} = 1

. The expected cost

E_{Ω} ({\bar{C}}_{t, ω}^{prod} (P_{t}, ω))

is denoted by

{\bar{C}}_{t}^{prod} (P_{t})

. Further, we suppose that the holding and backlogging costs do not change during the course of a disruption. This approach covers a wide range of potential types of disruptions, including cases where production becomes completely unavailable and cases where only maximum production capacity is reduced. The latter is directly linked to the literature on the implications of uncertain capacities on optimal production schedules (e.g., Ciarallo et al., 1994; Nourelfath, 2011). Nevertheless, in contrast to this literature stream, we assume that the disruption ω (

ω \neq ω_{0}

) can influence the demand by driving its distribution parameters up (e.g., as happened for particular products during COVID‐19 and as took place after the release of Omega× Swatch MoonSwatches) or down (e.g., in case of reputation damage), influencing the distribution

F_{t}

of a scenario with no disruptions and changing it to the distribution

G_{t, ω}

with density

g_{t, ω}

and mean

{\bar{μ}}_{t, ω}

. We differentiate between situations in which the manufacturer has full, limited, or no ability to estimate true distributions

G_{t, ω}

. Clearly, a special case in which the true distribution

G_{t, ω}

can be evaluated for a disruption scenario ω before the production starts may only serve as an idealistic benchmark. Nevertheless, the manufacturer may be able to estimate scenario‐specific distributions

F_{t, ω}

deviating from

G_{t, ω}

. The case with no ability to estimate scenario‐specific distributions can be defined as

F_{t, ω} = F_{t}, \forall ω

When faced with the risk of disruptions in the production pipeline, the manufacturer wishes to keep the production economically resilient. On the one hand, one can rely solely on a single‐level minimization of the expected total cost, taking disruption scenarios into account. In this case, if a disruption occurs, the manufacturer must adapt its strategy and react to the disruption in order to meet demand. By doing so, the manufacturer deviates from the expected optimum and observes the realized (and possibly very high) cost dependent on the disruption scenario. On the other hand, the manufacturer can use a robust or distributionally robust approach to prepare for the worst‐case disruption, ensuring that the realized cost is lower than estimated with a high probability (Ben‐Tal & Nemirovski, 1988, Nourelfath, 2011). Nevertheless, this solution may be overly conservative. In addition to these methods, there is another option that has not yet been well represented in the literature. In particular, the manufacturer can constrain the difference in optimal costs between (i) the strategy that can be implemented immediately in response to a disruption and (ii) a benchmark, which is an idealistic, but not realistically implementable, optimal production schedule. Under conditions of Lipschitz‐continuity of the cost functions, this approach bounds the difference in optimal production schedules between the above‐mentioned strategies, making the production more predictable by driving the schedule closer to the benchmark for each disruption scenario. Moreover, as we demonstrate further in this article, the cost difference between the production schedule implemented in normal times and the reactive strategy is also bounded in this case.

Overall, we control pair‐wise differences in costs between the following strategies:

Implemented strategy: We consider the manufacturer, whose upper‐level goal is to minimize the expected total cost of production and account for both production disruptions and demand perturbations. In this article, we aim for a resilient production schedule, which guarantees that the difference in costs and, importantly, optimal solutions, are bounded in the event that the decision needs to be adjusted in response to a specific disruption ω. We construct the bilevel optimization problem (i.e., problem (11)), which can guarantee this behavior. The resulting problem depends on the reactive (implementable in the case of a particular scenario ω) and the anticipative (i.e., benchmark) plans described next. These plans support the decision maker in constructing the solution with a bounded distance to the benchmark. Note that the implemented strategy corresponds to the upper‐level production plan, while the reactive strategy provides the lower‐level solution for each disruption scenario.

Reactive strategy: The reactive plan (8) is adopted starting at time

t_{ω}^{start}

in the event that a specific scenario ω occurs. The production schedule used prior to time

t_{ω}^{start}

equals

(P_{t}, P_{t}^{c}), \forall t = 1, …, t_{ω}^{start} - 1

, and corresponds to the implemented strategy, which accounts for all possible disruptions in probability and which the manufacturer aims to compute on the upper level. As the decision maker only adjusts the production when the disruption begins, we assume that the manufacturer can estimate likely changes in the demand distribution functions after time

t_{ω}^{start}

because the cause of the breakdown is known at time

t_{ω}^{start}

. The following reactive problem can be stated in mathematical terms, where

{\bar{C}}_{t, G_{t, ω}}^{orig} = {\bar{C}}_{t, ω}^{prod} (Q_{t}, ω) + h_{t} E_{G_{t, ω}} {(Q_{t}^{c} - D_{t}^{c})}^{+} + b_{t} E_{G_{t, ω}} {(D_{t}^{c} - Q_{t}^{c})}^{+}, \forall t

are stage‐wise cost functions taken in expectation over the true demand distribution

G_{t, ω}

\begin{matrix} v_{G_{ω}}^{react} (P, P^{c}, ω) = \min_{Q, Q^{c}} \sum_{t = 1}^{T} {\bar{C}}_{t, G_{t, ω}}^{orig} \\ s . t . & \sum_{u = 1}^{t} Q_{u} = Q_{t}^{c}, 0 \leq Q_{t} \leq κ_{t}, \forall t \in 1, …, T \\ Q_{t} = P_{t}, Q_{t}^{c} = P_{t}^{c}, \forall t \in 1, …, t_{ω}^{start} - 1 \\ Q_{t}^{c} \geq G_{t, ω}^{- 1} (1 - α_{t}), \forall t \in t_{ω}^{start}, …, T . \end{matrix}

Note that the service‐level constraints are automatically satisfied until

t_{ω}^{start} - 1

in problem (8) due to the adoption of the upper‐level optimal solution

(P_{t}, P_{t}^{c})

at times

t = 1, …, t_{ω}^{start} - 1

Partly anticipative strategy: In contrast to the reactive plan, anticipative strategies require the provision of a disruption scenario ω prior to the start of production. If this would be the case (e.g., the scenario without disruptions with

p_{ω_{0}} = 1

is often implemented in practice), the best possible production schedule would correspond to the anticipative plan. However, in the uncertain environment with the probability space Ω, the scenario ω can only be predicted with probability

p_{ω} ≪ 1

, which can be very low for high‐impact disruptions. Thus, scenario‐specific anticipative strategies may only serve as a low‐cost benchmark for both the implemented and the reactive production plans. The anticipative problems (9) and (10) optimize the strategy of a manufacturer preparing for a specific disruption ω at the beginning of the planning horizon. In a partly anticipative case (i.e., problem (9)), the manufacturer lacks the ability to conduct customer research to evaluate the true demand distribution of the scenario ω and can only estimate the distribution

F_{t, ω}

instead of

G_{t, ω}

. By this, the optimization problem with the stage‐wise cost function

{\bar{C}}_{t, F_{t, ω}}^{orig} = {\bar{C}}_{t, ω}^{prod} (Q_{t}, ω) + h_{t} E_{F_{t, ω}} {(Q_{t}^{c} - D_{t}^{c})}^{+} + b_{t} E_{F_{t, ω}} {(D_{t}^{c} - Q_{t}^{c})}^{+}

yields:

\begin{matrix} v_{F_{ω}}^{anticipate} (ω) = \min_{Q, Q^{c}} \sum_{t = 1}^{T} {\bar{C}}_{t, F_{t, ω}}^{orig} \\ s . t . & \sum_{u = 1}^{t} Q_{u} = Q_{t}^{c}, 0 \leq Q_{t} \leq κ_{t}, \forall t \in 1, …, T \\ Q_{t}^{c} \geq F_{t, ω}^{- 1} (1 - α_{t}), \forall t \in 1, …, T . \end{matrix}

Fully anticipative strategy: In contrast to the case outlined in problem (9), the manufacturer is able to estimate the true demand distribution

G_{t, ω}

for each time period t of the scenario ω in a fully anticipative case (i.e., problem (10)). This idealistic scenario with the stage‐wise cost function

{\bar{C}}_{t, G_{t, ω}}^{orig} = {\bar{C}}_{t, ω}^{prod} (Q_{t}, ω) + h_{t} E_{G_{t, ω}} {(Q_{t}^{c} - D_{t}^{c})}^{+} + b_{t} E_{G_{t, ω}} {(D_{t}^{c} - Q_{t}^{c})}^{+}

yields:

\begin{matrix} v_{G_{ω}}^{anticipate} (ω) = \min_{Q, Q^{c}} \sum_{t = 1}^{T} {\bar{C}}_{t, G_{t, ω}}^{orig} \\ s . t . & \sum_{u = 1}^{t} Q_{u} = Q_{t}^{c}, 0 \leq Q_{t} \leq κ_{t}, \forall t \in 1, …, T \\ Q_{t}^{c} \geq G_{t, ω}^{- 1} (1 - α_{t}), \forall t \in 1, …, T . \end{matrix}

Given the reactive and anticipative setups (8), (9), and (10), we state the bilevel cost‐minimization problem (11) with the stage‐wise objective

{\bar{C}}_{t}^{orig} = E_{Ω} [{\bar{C}}_{t, ω}^{prod} (P_{t}, ω) + h_{t} {(P_{t}^{c} - D_{t}^{c})}^{+} + b_{t} {(D_{t}^{c} - P_{t}^{c})}^{+}]

. Problem (11) takes production disruptions and demand perturbations explicitly into account and yields the implemented strategy described above. Importantly, the solutions to anticipative problems (9) and (10) can be computed independently of the optimization (11). In contrast, the reactive problem (8) is a lower‐level optimization problem, the solution of which depends on the optimal solution of the problem (11).

\begin{matrix} {\bar{v}}_{\bar{G}}^{orig} & = \min_{P, P^{c}} \sum_{t = 1}^{T} {\bar{C}}_{t}^{o r i g} \\ s.t. \sum_{u = 1}^{t} P_{u} = P_{t}^{c}, 0 \leq P_{t} \leq κ_{t}, \forall t \in 1, …, T \\ P_{Ω} (P_{t}^{c} \geq G_{t, ω}^{- 1} (1 - α_{t})) \geq 1 - β, \forall t \in 1, …, T (service level constr.) \\ P_{Ω} (v_{G_{ω}}^{react} (P, P^{c}, ω) - v_{F_{ω}}^{anticipate} (ω) \leq γ) \geq 1 - θ (market value constr.) \\ P_{Ω} (v_{G_{ω}}^{react} (P, P^{c}, ω) - v_{G_{ω}}^{anticipate} (ω) \leq γ) \geq 1 - θ (time value constr.) \\ \forall ω \in Ω : v_{G_{ω}}^{react} (P, P^{c}, ω) = \min_{Q, Q^{c}} \sum_{t = 1}^{T} {\bar{C}}_{t, G_{t, ω}}^{orig} \\ \sum_{u = 1}^{t} Q_{u} = Q_{t}^{c}, 0 \leq Q_{t} \leq κ_{t}, \forall t \in 1, …, T \\ Q_{t} = P_{t}, Q_{t}^{c} = P_{t}^{c}, \forall t \in 1, …, t_{ω}^{start} - 1 \\ Q_{t}^{c} \geq G_{t, ω}^{- 1} (1 - α_{t}), \forall t \in t_{ω}^{start}, …, T . \end{matrix}

We note that the objective function of problem (11) is equivalent to the function

\sum_{t = 1}^{T} [{\bar{C}}_{t}^{orig} = = {\bar{C}}_{t}^{prod} (P_{t}) + h_{t} E_{{\bar{G}}_{t}} {(P_{t}^{c} - D_{t}^{c})}^{+} + b_{t} E_{{\bar{G}}_{t}} {(D_{t}^{c} - P_{t}^{c})}^{+}]

, where

{\bar{G}}_{t} (D_{t}^{c}) = \sum_{ω \in Ω} p_{ω} G_{t, ω} (D_{t}^{c})

is a mixture distribution with mean

{\bar{μ}}_{t} = \sum_{ω \in Ω} p_{ω} {\bar{μ}}_{t, ω}

{\bar{C}}_{t}^{prod} (P_{t}) = E_{Ω} ({\bar{C}}_{t, ω}^{prod} (P_{t}, ω))

and

{\bar{C}}_{t}^{orig}

is the stage‐wise expected cost function. Importantly, due to the constraints of the bilevel problem (11), the difference in optimal costs

v_{G_{ω}}^{react} (P, P^{c}, ω) - {\bar{v}}_{\bar{G}}^{orig}

can obviously be bounded from above in line with the following implicit inequality:

v_{G_{ω}}^{react} (P, P^{c}, ω) - {\bar{v}}_{\bar{G}}^{orig} \leq \leq γ + \min {v_{F_{ω}}^{anticipate} (ω) - {\bar{v}}_{\bar{G}}^{orig}; v_{G_{ω}}^{anticipate} (ω) - {\bar{v}}_{\bar{G}}^{orig}}, \forall ω \in Ω

. Thus, for scenarios ω, which have low anticipative costs (i.e.,

v_{G_{ω}}^{anticipate} (ω) \leq {\bar{v}}_{\bar{G}}^{orig}

v_{F_{ω}}^{anticipate} (ω) \leq {\bar{v}}_{\bar{G}}^{orig}

), the difference in optimal values between reactive and implemented strategies can be bounded by γ. For high anticipative costs, the difference can be bounded explicitly as

v_{G_{ω}}^{react} (P, P^{c}, ω) - - {\bar{v}}_{\bar{G}}^{orig} \leq γ + \min {v_{F_{ω}}^{anticipate} (ω); v_{G_{ω}}^{anticipate} (ω)}, \forall ω

, that is,

\begin{matrix} v_{G_{ω}}^{react} (P, P^{c}, ω) - {\bar{v}}_{\bar{G}}^{orig} \leq \\ \leq \{\begin{matrix} γ, if v_{G_{ω}}^{anticipate} (ω) \leq {\bar{v}}_{\bar{G}}^{orig} or v_{F_{ω}}^{anticipate} (ω) \leq {\bar{v}}_{\bar{G}}^{orig} \\ γ + \min {v_{F_{ω}}^{anticipate} (ω); v_{G_{ω}}^{anticipate} (ω)}, otherwise. \end{matrix} \end{matrix}

Furthermore, the solution of the reactive problem for scenario ω is feasible in the upper‐level problem (11). Thus, the difference is nonnegative in expectation, that is,

E_{Ω} (v_{G_{ω}}^{react} (P, P^{c}, ω) - {\bar{v}}_{\bar{G}}^{orig}) \geq 0

Overall, the scenario‐specific differences in costs between anticipative and reactive strategies can, therefore, be bounded as schematically shown in Figure 2, where we introduce the notion of regret functions

R_{F_{ω}}^{G_{ω}} (P, P^{c}, ω) : = v_{G_{ω}}^{react} (P, P^{c}, ω) - v_{F_{ω}}^{anticipate} (ω)

and

R_{G_{ω}}^{G_{ω}} (P, P^{c}, ω) : = v_{G_{ω}}^{react} (P, P^{c}, ω) - v_{G_{ω}}^{anticipate} (ω)

, which measure the scenario‐specific cost surplus (or, possibly, cost reduction) of the reactive strategy with respect to the anticipative plans.

FIGURE 2

Bounded regrets

Importantly, the regret

R_{G_{ω}}^{G_{ω}} (P, P^{c}, ω)

measures the time value of the information in the disruption scenario ω (the course of the disruption is known at the start of the planning horizon in the fully anticipative plan, and at time

t_{ω}^{start}

in the reactive plan). The constraint on this regret is referred to as the time value constraint in problem (11). Further, the regret

R_{F_{ω}}^{G_{ω}} (P, P^{c}, ω)

measures the value of knowing the demand in the reactive plan in comparison to preparing for a specific production disruption scenario ω without conducting market research (i.e., as in the partly anticipative plan). The constraint on this regret is referred to as the market value constraint in problem (11). Both constraints are chance constraints requiring that the regrets do not exceed some threshold γ with probability

1 - θ

, where

θ \in [0, 1]

is a given risk budget. Note that the regret

R_{G_{ω}}^{G_{ω}} (P, P^{c}, ω)

takes only nonnegative values. In contrast, the regret

R_{F_{ω}}^{G_{ω}} (P, P^{c}, ω)

can take negative values for specific relationships between distributions

F_{t, ω}

and

G_{t, ω}

(see Section 2.3, Appendix A.2 in the Supporting Information).

High service‐level requirements and the regret reduction

Next, we demonstrate that high service levels play a crucial role in reducing regrets beyond the threshold γ. For this, we introduce the cost misfit implied by the difference in holding and backlogging amounts due to distributions

G_{t, ω}

and

F_{t, ω}

with respective means

{\bar{μ}}_{t, ω}

and

μ_{t, ω}

, that is,

\begin{matrix} δ_{F_{ω}}^{G_{ω}} (Q) & = & \sum_{t = 1}^{T} [h_{t} \int_{0}^{Q} (Q - D_{t}^{c}) (g_{t, ω} (D_{t}^{c}) - f_{t, ω} (D_{t}^{c})) d D_{t}^{c} + \\ + b_{t} \int_{Q}^{\infty} (D_{t}^{c} - Q) (g_{t, ω} (D_{t}^{c}) - f_{t, ω} (D_{t}^{c})) d D_{t}^{c}] . \end{matrix}

Using the above and following Figure A5 in Appendix A.2 in the Supporting Information, one can claim that the time value constraint is not active for high service‐level requirements and a right shift of distributions

G_{t, ω}

w.r.t.

F_{t, ω}, \forall t \in [t_{ω}^{start}, t_{ω}^{end}]

. This is due to the fact that the misfit is negative under the anticipative solution

Q_{G_{t, ω}}^{*, anticipate}

of problem (10) in this case. Similarly, the market value constraint is not active for high service levels when the distribution

G_{t, ω}

w.r.t.

F_{t, ω}

shifts to the left because the difference in holding and backlogging costs evaluated at the optimal solution

Q_{F_{t, ω}}^{*, anticipate}

of the anticipative problem (9) is positive. Based on this, one can conclude that high service‐level requirements (i.e.,

α_{t}

is high enough) play a very important role in reducing regrets between different strategies beyond the threshold γ; for example,

α_{t} \to \infty

leads to

R_{G_{ω}}^{G_{ω}} (P, P^{c}, ω) \leq γ - \sum_{t = 1}^{T} h_{t} ({\bar{μ}}_{t, ω} - μ_{t, ω})

{\bar{μ}}_{t, ω} \geq μ_{t, ω}

, as well as to

R_{F_{ω}}^{G_{ω}} (P, P^{c}, ω) < γ - \sum_{t = 1}^{T} h_{t} (μ_{t, ω} - {\bar{μ}}_{t, ω})

μ_{t, ω} > {\bar{μ}}_{t, ω}

. Further, we refer to the set of constraints

P_{Ω} (P_{t}^{c} \geq G_{t, ω}^{- 1} (1 - α_{t})) \geq 1 - β, \forall t

as doubly probabilistic service‐level requirements, because they account both for the demand uncertainty and production disruptions. As one can observe, we allow the service constraint to be violated with probability not exceeding the risk budget β.

Expected opportunity costs (EOCs) and the regret reduction

Further, before we proceed with the approximation of problem (11), we analyze the potential loss incurred by a decision maker whenever a disruption scenario ω occurs. For this, we account for expected opportunity costs (EOCs) (i.e., forgone benefits) of anticipative strategies (9) and (10) with respect to the reactive strategy (8). In particular, the EOC of knowing about the production disruption is defined as the difference in the manufacturer's profits between the anticipative (i.e., foregone) and reactive strategies. Though problem (11) is a cost‐minimization problem, its optimal solution coincides with the decision in the corresponding profit‐maximization problem if the total revenue from sales depends on the uncertain demand rather than production quantity (see Appendix A.2 in the Supporting Information). As a result, the EOCs of fully and partly anticipative strategies (denoted by

{EOC}_{G_{ω}}^{G_{ω}}

and

{EOC}_{F_{ω}}^{G_{ω}}

, respectively) can be bounded by respective regret functions for high service levels under conditions described in Appendix A.2 in the Supporting Information:

\begin{matrix} \begin{matrix} {EOC}_{F_{ω}}^{G_{ω}} \leq v_{G_{ω}}^{react} - v_{F_{ω}}^{anticipate} = R_{F_{ω}}^{G_{ω}} (P, P^{c}, ω), \\ {EOC}_{G_{ω}}^{G_{ω}} = v_{G_{ω}}^{react} - v_{G_{ω}}^{anticipate} = R_{G_{ω}}^{G_{ω}} (P, P^{c}, ω) . \end{matrix} \end{matrix}

Based on these bounds, the manufacturer is able to reduce scenario‐specific opportunity costs requiring high service‐level constraints while imposing an upper bound on the regret functions.

It should be noted that, similar to the opportunity cost, the regret function

R_{F_{ω}}^{G_{ω}} (P, P^{c}, ω)

can take negative values. This is especially true if the disruption scenario ω results in a left shift in the demand distribution

G_{t, ω}

w.r.t.

F_{t, ω}

. In this case, the manufacturer can reduce production while maintaining service levels in response to the disruption caused by a likely drop in demand (Figure 3a). Otherwise, if the demand distribution is not affected in scenario ω, the regret function,

R_{F_{ω}}^{G_{ω}} (P, P^{c}, ω)

, and the expected opportunity cost,

{EOC}_{F_{ω}}^{G_{ω}}

, become nonnegative (see Figure 3b).

FIGURE 3

Possible optimal values of the anticipative and reactive problems

It is also worth noting that the probability of negative regrets should not be constrained because they reduce EOCs and appreciate the ability to respond to an upcoming disruption at a lower cost than preparation in advance would allow. This is also why the mirror constraint

P_{Ω} (- R_{F_{t, ω}}^{G_{t, ω}} (P, P^{c}, ω) \leq γ) \geq 1 - θ

is not included in the optimization problem, despite the fact that such a generalization is straightforward.

Approximation for numerical solution

The optimal solution of (11) is called a resilient plan for problem (1). Our goal is to accurately and efficiently develop such a plan in the absence of complete knowledge about distributions

G_{t, ω} \forall t, ω

. For this, we construct the approximation (15), where

{\bar{R}}_{F_{ω}}^{F_{ω}} (P, P^{c}, ω) = {\bar{v}}_{F_{ω}}^{react} (P, P^{c}, ω) - - {\bar{v}}_{F_{ω}}^{anticipate} (ω)

denotes the regret dependent on the optimal values

{\bar{v}}_{F_{ω}}^{react} (P, P^{c}, ω)

and

{\bar{v}}_{F_{ω}}^{anticipate} (ω)

of reactive and anticipative problems with the modified objective function (see (15)) and allows time and market value constraints to be combined in one approximating chance constraint. The modified problem (15) with the stage‐wise cost

{\bar{C}}_{t}^{approx} = {\bar{C}}_{t}^{prod} (P_{t}) + h_{t} {(P_{t}^{c} - μ_{t})}^{+} + b_{t} {(μ_{t} - P_{t}^{c})}^{+}

yields:

\begin{matrix} {\bar{v}}_{F}^{approx} = \min_{P, P^{c}} \sum_{t = 1}^{T} {\bar{C}}_{t}^{approx} \\ s.t. \sum_{u = 1}^{t} P_{u} = P_{t}^{c}, 0 \leq P_{t} \leq κ_{t}, \forall t \in 1, …, T \\ P_{Ω} (P_{t}^{c} \geq F_{t, ω}^{- 1} (1 - {\bar{α}}_{t})) \geq 1 - β, \forall t \in 1, …, T (service-level constr.) \\ P_{Ω} ({\bar{R}}_{F_{ω}}^{F_{ω}} (P, P^{c}, ω) \leq γ) \geq 1 - θ (time and market value constr.) \\ \forall ω \in Ω : {\bar{v}}_{F_{ω}}^{react} (P, P^{c}, ω) = \\ = \min_{Q, Q^{c}} \sum_{t = 1}^{T} ({\bar{C}}_{t, ω}^{prod} (Q_{t}, ω) + h_{t} {(Q_{t}^{c} - μ_{t, ω})}^{+} + b_{t} {(μ_{t, ω} - Q_{t}^{c})}^{+}) \\ \sum_{u = 1}^{t} Q_{u} = Q_{t}^{c}, 0 \leq Q_{t} \leq κ_{t}, \forall t \in 1, …, T \\ Q_{t} = P_{t}, Q_{t}^{c} = P_{t}^{c}, \forall t \in 1, …, t_{ω}^{start} - 1 \\ Q_{t}^{c} \geq F_{t, ω}^{- 1} (1 - {\bar{α}}_{t}), \forall t \in t_{ω}^{start}, …, T . \\ \forall ω \in Ω : {\bar{v}}_{F_{ω}}^{anticipate} (ω) = \\ = \min_{Q_{t}, Q_{t}^{c}} \sum_{t = 1}^{T} ({\bar{C}}_{t}^{prod} (Q_{t}, ω) + h_{t} {(Q_{t}^{c} - μ_{t, ω})}^{+} + b_{t} {(μ_{t, ω} - Q_{t}^{c})}^{+}) \\ \sum_{u = 1}^{t} Q_{u} = Q_{t}^{c}, 0 \leq Q_{t} \leq κ_{t}, \forall t \in 1, …, T \\ Q_{t}^{c} \geq F_{t, ω}^{- 1} (1 - {\bar{α}}_{t}), \forall t \in 1, …, T . \end{matrix}

Note that the approximate problem (15) relies on biased distribution functions

F_{t}

and

F_{t, ω}

with means

μ_{t}

and

μ_{t, ω}

at each period t instead of the true distributions

G_{t}

and

G_{t, ω}

because the manufacturer only has limited information about customer demand. Furthermore, we use the service level

1 - {\bar{α}}_{t}

, which is different from the service level

1 - α_{t}

. We state the condition (16) to guarantee a fine approximation of (11).

Optimization problem (15) constitutes a bilevel decision‐making program because the regret function form is dependent on the reactive solution. In contrast, the anticipative problem is a single‐level optimization problem that can be solved efficiently for all scenarios. To constitute a fine approximation of problem (11), the optimal solution of (15) must (1) be feasible in problem (11) (see Proposition 1) and (2) provide a tight optimality gap on the optimal value of problem (11) (see Theorem 2). Under the condition of Lipschitz‐continuity, the latter also implies a small difference in optimal solutions of problems (11) and (15). Proposition 1 Feasibility

With a small abuse of notations, let

(P^{*}, P^{c *})

be the optimal flow of problem (15). If the following conditions on the fraction of unsatisfied customers hold:

\begin{matrix} {\bar{α}}_{t} \leq 1 - F_{t, ω} (G_{t, ω}^{- 1} (1 - α_{t})), \forall t, ω, \end{matrix}

and if

\exists \tilde{Ω} \in Ω

such that

R_{G_{ω}}^{G_{ω}} (P^{*}, P^{c *}, ω) \leq γ

R_{F_{ω}}^{G_{ω}} (P^{*}, P^{c *}, ω) \leq γ

\forall ω \in \tilde{Ω}

and

\sum_{ω \in \tilde{Ω}} p_{ω} \geq 1 - θ

, the optimal flow

(P^{*}, P^{c *})

is feasible in the problem (11).

Proof

See Appendix A.3 in the Supporting Information for the complete proof and note that by increasing the set of scenarios Ω, the manufacturer increases the probability of finding its subset

\tilde{Ω}

, for which

R_{G_{ω}}^{G_{ω}} (P^{*}, P^{c *}, ω) \leq γ

R_{F_{ω}}^{G_{ω}} (P^{*}, P^{c *}, ω) \leq γ, \forall ω \in \tilde{Ω}

and

\sum_{ω \in \tilde{Ω}} p_{ω} \geq 1 - θ

□

Theorem 2

Let the conditions of Proposition 1 be satisfied. Denote by

{\bar{v}}_{G}^{approx}

the optimal value of the problem (15) if found with the true perturbed distributions

G_{t, ω} \forall t, ω

. Thus, the following bounds hold:

(Lower Bound) The lower bound on the original problem (11) yields

{\bar{v}}_{G}^{orig} \geq {\bar{v}}_{G}^{approx}

(Upper Bounds) The following upper bounds (denoted u₁ and u₂) hold:

\begin{matrix} {\bar{v}}_{G}^{orig} & - & {\bar{v}}_{F}^{approx} \leq \sum_{t = 1}^{T} (h_{t} + b_{t}) \int_{P_{t}^{c *}}^{+ \infty} D_{t}^{c} {\bar{g}}_{t} (D_{t}^{c}) d D_{t}^{c} + \\ + & \sum_{t = 1}^{T} h_{t} (μ_{t} - {\bar{μ}}_{t}) = u_{1}, \end{matrix}

\begin{matrix} {\bar{v}}_{G}^{orig} & - & {\bar{v}}_{F}^{approx} \leq \sum_{t = 1}^{T} (h_{t} + b_{t}) \int_{0}^{P_{t}^{c *}} {\bar{G}}_{t} (D_{t}^{c}) d D_{t}^{c} - \\ - & \sum_{t = 1}^{T} b_{t} (μ_{t} - {\bar{μ}}_{t}) = u_{2}, \end{matrix}

where

{\bar{g}}_{t} (D_{t}^{c}) = \sum_{ω \in Ω} p_{ω} g_{t, ω} (D_{t}^{c})

is a density of a mixture distribution

{\bar{G}}_{t} (D_{t}^{c})

The optimality gap

∣ {\bar{v}}_{G}^{o r i g} - {\bar{v}}_{F}^{a p p r o x} ∣

can be bounded by

Δ_{(G, F)}

, where

\begin{matrix} Δ_{(G, F)} = \max [{\bar{v}}_{F}^{approx} - {\bar{v}}_{G}^{approx}, \min {u_{1}, u_{2}}] . \end{matrix}

Moreover,

Δ_{(G, F)} \to Δ

{\bar{μ}}_{t} \to μ_{t}

{\bar{v}}_{G}^{approx} \to {\bar{v}}_{F}^{approx}

and if the service‐level requirements of the original problem (11) approach those of the approximate problem (15).

(Relative Error) The relative error of the approximation follows from statement 2, that is,

\begin{matrix} \frac{∣ {\bar{v}}_{G}^{orig} - {\bar{v}}_{F}^{approx} ∣}{{\bar{v}}_{F}^{approx}} \leq \frac{Δ_{(G, F)}}{\sum_{t = 1}^{T} ({\bar{C}}_{t}^{prod} (P_{t}^{*}) + h_{t} {(P_{t}^{c *} - μ_{t})}^{+} + b_{t} {(μ_{t} - P_{t}^{c *})}^{+})} . \end{matrix}

Proof

See Appendix A.4 in the Supporting Information for the complete proof, which is based on the results of Proposition 1 and proceeds from stating lower bounds for the problem's objective function to its upper bounds.

□

When Theorems 1 and 2 are compared, one can see that the upper bounds and optimality gap in the case of production disruptions have additional terms due to the difference in demand distributions. These terms describe the deviation of production quantiles in the perturbed distribution.

Further in the article, we consider the following two cases: Case 1:

The service level

1 - {\bar{α}}_{t}

is high enough to guarantee condition (16).

Case 2:

The chosen service level

1 - {\bar{α}}_{t}

violates condition (16) for some

t, ω

Case 1 is described in Theorem 2 and guarantees the existence of the upper bound for the optimality gap given the set of scenarios

\tilde{Ω}

that satisfies the chance constraint on the regret function (see Proposition 1). If Case 2 holds, the optimal point

(P_{t}^{*}, P_{t}^{* c})

of the approximation problem (15) does not necessarily satisfy the conditions of Theorem 2. Thus, the difference between Cases 1 and 2 is that the producer can maintain the desired service level in Case 1 by solving the optimization problem (15), whereas this is not always possible in Case 2. Note that Case 1 holds true when a production disruption leads to a drop in demand, whereas Case 2 corresponds to scenarios in which the mean demand and variance tend to increase.

SOLUTION METHOD

We introduce the inventory level

K_{t} = {(P_{t}^{c} - μ_{t})}^{+}

and the shortage

M_{t} = {(μ_{t} - P_{t}^{c})}^{+}

for the solution of problem (15). This allows us to simplify the objective function using the set of constraints (21c) and to state the following upper‐level manufacturer's optimization problem:

\min \sum_{t = 1}^{T} ({\bar{C}}_{t}^{prod} (P_{t}) + h_{t} K_{t} + b_{t} M_{t})

21a

\begin{matrix} s . t . & \sum_{u = 1}^{t - 1} P_{u} + P_{t} = P_{t}^{c}, 0 \leq P_{t} \leq κ_{t}, \forall t \in 1, …, T, \end{matrix}

21b

K_{t} - M_{t} = P_{t}^{c} - μ_{t}, K_{t} \geq 0, M_{t} \geq 0, \forall t \in 1, …, T,

21c

\begin{matrix} P_{t}^{c} \geq F_{t, ω}^{- 1} (1 - {\bar{α}}_{t}) - H y_{ω}, \sum_{ω \in Ω} p_{ω} y_{ω} \leq β, \\ y_{ω} \in {0, 1}, \forall ω \in Ω, t \in 1, …, T, \end{matrix}

21d

\begin{matrix} {\bar{R}}_{F_{ω}}^{F_{ω}} (P, P^{c}, ω) \leq H z_{ω} + γ, \sum_{ω \in Ω} p_{ω} z_{ω} \leq θ, \\ z_{ω} \in {0, 1}, \forall ω \in Ω, \end{matrix}

21e

where H is a sufficiently high number and

y_{ω}, z_{ω}

are binary decision variables, which determine scenarios violating thresholds β and γ.

Here, demand distributions

F_{t, ω}

may vary for different disruption scenarios, that is, the manufacturer can attempt to estimate scenario‐specific demands. Nevertheless, if

F_{t, ω} = F_{t} \forall t, ω

, that is, if the manufacturer assumes equal demand distributions in all disruption scenarios, the service‐level requirements

P_{t}^{c} \geq F_{t}^{- 1} (1 - {\bar{α}}_{t})

are satisfied in all time periods in the upper‐level problem, as long as

\exists ω : y_{ω} = 0

in the constraint (21d). This is in line with the fact that the decision maker aims to deliver the desired service level and to recollect all the scenarios for which it will be violated (i.e., those ω for which

y_{ω} = 1

Next, the lower‐level optimization problem (22) arises because the reactive setup in the regret function

{\bar{R}}_{F_{ω}}^{F_{ω}} (P, P^{c}, ω) = {\bar{v}}_{F_{ω}}^{react} (P, P^{c}, ω) - {\bar{v}}_{F_{ω}}^{anticipate} (ω)

depends on the optimal solution of the manufacturer's cost‐minimization problem (15). Furthermore, the regret

{\bar{R}}_{F_{ω}}^{F_{ω}} (P, P^{c}, ω)

uses the solution of the anticipative problem (23), which, however, does not depend on the production

(P, P^{c})

and, thus, can be solved prior to problem (15)

\forall ω

, that is,

\begin{matrix} {\bar{v}}_{F_{ω}}^{react} (P, P^{c}, ω) = & \min_{Q, Q^{c}} & \sum_{t = 1}^{T} ({\bar{C}}_{t}^{prod} (Q_{t}, ω) + h_{t} K_{t} + b_{t} M_{t}) \\ s . t . & \sum_{u = 1}^{t} Q_{u} = Q_{t}^{c}, 0 \leq Q_{t} \leq κ_{t}, \forall t \in 1, …, T, \\ Q_{t} = P_{t}, Q_{t}^{c} = P_{t}^{c}, \forall t \in 1, …, t_{ω}^{start} - 1, \\ K_{t} - M_{t} = Q_{t}^{c} - μ_{t}, K_{t} \geq 0, M_{t} \geq 0, \forall t \in 1, …, T, \\ Q_{t}^{c} \geq F_{t, ω}^{- 1} (1 - {\bar{α}}_{t}), \forall t \in 1, …, T . \end{matrix}

\begin{matrix} {\bar{v}}_{F_{ω}}^{anticipate} (ω) = & \min_{Q, Q^{c}} & \sum_{t = 1}^{T} ({\bar{C}}_{t}^{prod} (Q_{t}, ω) + h_{t} K_{t} + b_{t} M_{t}) \\ s . t . & \sum_{u = 1}^{t} Q_{u} = Q_{t}^{c}, 0 \leq Q_{t} \leq κ_{t}, \forall t \in 1, …, T, \\ K_{t} - M_{t} = Q_{t}^{c} - μ_{t}, K_{t} \geq 0, M_{t} \geq 0, \forall t \in 1, …, T, \\ Q_{t}^{c} \geq F_{t, ω}^{- 1} (1 - {\bar{α}}_{t}), \forall t \in 1, …, T . \end{matrix}

The presented bilevel optimization problem is a mixed‐integer nonlinear program that can be addressed using solvers, such as BARON in small‐dimensional cases. Additionally, if production cost functions

{\bar{C}}_{t}^{prod}

are linear or piecewise linear functions, the problem can be solved by standard solvers such as CPLEX or Gurobi. However, the complexity increases depending on the length of the planning horizon T. In this article, we develop customized feasibility cuts and combine Benders decomposition with network speed‐up methods. Overall, we propose the following procedure for the solution of the problem:

Step 1

(Scenario generation) : Define the form of convex cost functions

{\bar{C}}_{t}^{prod}

and generate the scenario set Ω using discrete‐time Markov processes.

Step 2

(Scenario reduction) : Employ scenario reduction techniques aiming for high‐impact scenarios with high cumulative probability.

Step 3

(Anticipative solution) : Solve the anticipative problem (23) for each

ω \in Ω

using network speed‐up techniques.

Step 4

(Bilevel optimization) : Solve the bilevel optimization problem using a Bender's decomposition procedure with customized feasibility cuts.

Scenario generation

We consider the following three uncertainty types influencing the production pipeline:

Type 1: Production disruptions that set the capacity

κ_{t}

to zero in disrupted periods. To generate such scenarios, we employ discrete‐time Markov decision processes and fix a disruption and a recovery probability at each time stage (Schmitt et al., 2010).

We let

ζ_{t}

denote the transition probability from a state with no disruption to a disrupted state at time period t (i.e., the disruption probability). Also, we let

χ_{t}

be the transition probability from a disrupted state to a state with no disruption at time period t (i.e., the recovery probability). We assume that these probabilities are independent of each other. By varying

ζ_{t}

and

χ_{t}

, we are able to account for different types of disruptions as demonstrated in Figure 4. This includes those disruptions that are more probable during particular time periods, for example, during vacations when a large proportion of a company's employees are out of the office (see Figure 4c).

FIGURE 4

Probability of production disruptions with different disruption lengths (

T = 25

)

The steady‐state probabilities of this Markov chain are unique due to the well‐known Perron–Frobenius Theorem for nonnegative irreducible transition matrices. We let ξ₀ be the steady‐state probability of no disruption and

ξ_{s}

be the steady‐state probability of

s \geq 1

disrupted periods. If

ζ_{t} = ζ

and

χ_{t} = χ \forall t

stay constant over time, the steady‐state probabilities can be written similarly to those computed by Schmitt et al. (2010), that is,

ξ_{0} = \frac{χ}{ζ + χ} and ξ_{s} = \frac{ζ χ}{ζ + χ} {(1 - χ)}^{s - 1}, s \geq 1,

with the following expected duration and variance:

\begin{matrix} E (s) & = & 0 \cdot ξ_{0} + \sum_{s = 1}^{\infty} s ξ_{s} = \frac{ζ}{(ζ + χ) χ}, \\ V a r (s) & = & E^{2} (s) (1 + 2 \frac{χ}{ζ}) - E (s) . \end{matrix}

Note that the expected duration and variance of finite‐state Markov processes with varying

ζ_{t}

and

χ_{t}

can be computed using the probabilities of scenarios with different disruption lengths.

Given

T = 20

periods, Figure 5a demonstrates the cases when the expected duration increases or decreases, while Figure 5b shows the corresponding variance. Similar to the behavior outlined in Equation (25) for infinite‐state Markov chains, the expected duration grows as the stage‐wise probabilities of disruptions increase. Furthermore, the expected duration drops if the recovery probabilities

χ_{t}

decrease. However, if the early and late recovery probabilities change in opposite directions, resulting in a bell‐shaped function, the expected duration enlarges. This behavior is not accounted for in Equation (25) and emphasizes the importance of avoiding underestimation of an average disruption duration, especially in the case of probability model errors.

FIGURE 5

Expected duration and variance of finite horizon disruptions

Type 2: Production disruptions that reduce the available production capacity to a nonzero value. To extend the scenario generation method to production disruptions with a nonzero remaining capacity level, we combine Markov decision processes with an optimal quantization of a random variable ξ describing the loss in the capacity level

κ_{t}

κ_{t} ξ

and the resting capacity as

κ_{t} (1 - ξ)

(see Figure 6a for the truncated Pareto distribution). Similarly, we could increase the cost

{\bar{C}}_{t}^{prod}

{\bar{C}}_{t}^{prod} (1 + ξ)

by using a different random variable ξ (see Figure 6b for the lognormal case). The optimal quantization is used in the well‐known sense of the minimal Kantorovich–Wasserstein distance between a continuous distribution of ξ and its discrete approximation. The distance is introduced by Kantorovich (2006) and is based on the transportation plan of a given total volume to optimal locations (Villani, 2008).

FIGURE 6

Optimal quantization in the sense of the minimal Kantorovich–Wasserstein distance

Type 3: A combination of production disruptions with uncertain demand. To account for the uncertainty in demand, we impose service‐level constraints, which must hold for production disruption scenarios

ω \in Ω

with a probability of at least β. Importantly, in the case of high service‐level requirements, the knowledge of the true demand distribution function

G_{t, ω}

is not required according to Theorem 2. This is a crucial contribution of our research, allowing us to make high‐quality production decisions with only a limited understanding of the demand pattern (i.e., using distributions

F_{t, ω}

instead of

G_{t, ω}, \forall ω

The total number of disruption scenarios generated is equal to

S = 2^{T} N

, where T is the number of time periods and N is the number of optimal quantizers of the random variable ξ. The value grows exponentially in the number of periods T and is linear in the number of quantizers N. This motivates the use of scenario reduction techniques to improve the efficiency of the algorithms. In particular, we use a three‐step scenario reduction procedure (see Appendix A.5 in the Supporting Information) to select disruptions for the numerical solution of optimization problems (21a) to (23). The procedure is based on (1) measuring, (2) ranking, and (3) reducing the scenario set, and it also allows for the solution of the anticipative problem (23) to be accelerated through the use of warm‐start techniques. Indeed, we solve the anticipative problem (23) for each scenario prior to the solution of the upper‐level and reactive problems (21a) and (22) because the anticipative solution does not depend on problem (21a). Proceeding with the reversed ranking order of Appendix A.5 in the Supporting Information, we set the solution for the previous scenario as a starting point to the next problem. This warm‐start approach follows the theoretical framework developed by Granot and Veinott Jr (1985) and Veinott Jr (2005) for two parametric convex MCF problems with parameter sets l,

l^{'}

and optimal solutions

\hat{x} (l)

\hat{x} (l^{'})

. The result is known as the smoothing theorem and states that

∥ \hat{x} (l^{'}) - \hat{x} (l) ∥_{\infty} \leq {∥ l^{'} - l ∥}_{1}

(see Appendix A.6 in the Supporting Information). Figure 7a–c compares the well‐known interior‐point, sequential quadratic programming (SQP), and active‐set methods with and without warm‐start techniques and demonstrates a respective reduction in computation time by approximately 75%, 67%, and 50% for warm‐start cases.

FIGURE 7

The anticipative solution with the warm‐start and a quadratic production cost function

Bilevel optimization using Benders decomposition

To solve the optimization problem (21a), we apply a generalized Benders decomposition introduced in Geoffrion (1972). A feasible point

(P, P^{c}, y, z)

satisfying constraints (21b), (21c), and (21d) is the optimal solution of problem (21a) if

(P, P^{c})

minimizes the objective function and the inequality

{\bar{R}}_{F_{ω}}^{F_{ω}} (P, P^{c}, ω) \leq H z_{ω} + γ

holds for every scenario ω. In order to evaluate the regret function, we compute the objective value of problem (22) for each scenario ω such that

z_{ω} = 0

and check whether it is smaller than

{\bar{v}}_{F_{ω}}^{anticipate} (ω) + γ

. To do this, we iterate over binary solutions z, customizing the Benders decomposition procedure by efficient feasibility cuts added to problem (21a), which we refer to as the relaxed master problem. At each iteration k, the Benders procedure obtains a solution

(P^{k}, P^{c, k}, z^{k})

of problem (21a) under additional cut constraints and verifies whether the feasibility of

{\bar{R}}_{F_{ω}}^{F_{ω}} (P^{k}, P^{c, k}, ω) \leq γ

holds for all ω such that

z_{ω}^{k} = 0

. For each scenario that violates the constraint, we add a feasibility cut to the optimization by solving the Benders feasibility problem:

\min_{Δ_{ω}} | Δ_{ω} | s.t. {\bar{R}}_{F_{ω}}^{F_{ω}} (P^{k}, P^{c, k}, ω) - γ \leq Δ_{ω},

where the optimal value

Δ_{ω} = Δ_{Feas}^{k} (P^{k}, P^{c, k}, ω)

is the amount by which the constraint is violated. Due to the form of the regret function, the optimization problem (26) can be written as a convex MCF problem. The dual multipliers

ρ_{t}^{+}

ρ_{t}^{-}

ρ_{t}^{c +}

, and

ρ_{t}^{c -}

\forall t \in 1, …, t_{ω}^{start} - 1

corresponding to its network flow equilibrium constraints can be derived as described in Bertsekas (1998). Using these variables, we are able to specify feasibility cuts for the relaxed master problem:

\begin{matrix} H z_{ω} & \geq & Δ_{Feas}^{k} (P^{k}, P^{c, k}, ω) + \sum_{t \in 1, …, t_{ω}^{start} - 1} (ρ_{t}^{+} + ρ_{t}^{-}) (P_{t} - P_{t}^{k}) \\ + \sum_{t \in 1, …, t_{ω}^{start} - 1} (ρ_{t}^{c +} + ρ_{t}^{c -}) (P_{t}^{c} - P_{t}^{c, k}) . \end{matrix}

Note that the left‐hand side of the cut controls whether the regret constraint is active for scenario ω in the iteration k of the Benders procedure (Algorithm 1). A similar approach was applied in Fischetti et al. (2016) to construct cuts for fixing constraints.

ALGORITHM 1

Benders decomposition with customized feasibility cuts.

Set $k = 0$ and
$z_{ω}^{k} = 0, \forall ω \in Ω$
. Solve the relaxed master problem (21) obtaining its solution
$(P^{k}, P^{c, k}, y_{ω}^{k})$
.

while
$\exists ω : z_{ω}^{k} = 0 and {\bar{R}}_{F_{ω}}^{F_{ω}} (P^{k}, P^{c, k}, ω) > γ$
do
$\begin{matrix} \min \sum_{t = 1}^{T} ({\bar{C}}_{t}^{prod} (P_{t}) + h_{t} K_{t} + b_{t} M_{t}) \\ s . t . \sum_{u = 1}^{t - 1} P_{u}^{c} + P_{t} = P_{t}^{c}, 0 \leq P_{t} \leq κ_{t}, \\ K_{t} - M_{t} = P_{t}^{c} - μ_{t}, K_{t} \geq 0, M_{t} \geq 0, \forall t \in 1, …, T \\ P_{t}^{c} \geq F_{t, ω}^{- 1} (1 - {\bar{α}}_{t}) - H y_{ω}, \sum_{ω \in Ω} p_{ω} y_{ω} \leq β, \\ y_{ω} \in {0, 1}, \forall ω \in Ω, t \in 1, …, T \\ H z_{ω} \geq Δ_{Feas}^{k} (P^{k}, P^{c, k}, ω) + \sum_{t \in 1, …, t_{ω}^{start} - 1} (ρ_{t}^{+} + ρ_{t}^{-}) (P_{t} - P_{t}^{k}) \\ + \sum_{t \in 1, …, t_{ω}^{start} - 1} (ρ_{t}^{c +} + ρ_{t}^{c -}) (P_{t}^{c} - P_{t}^{c, k}) \sum_{ω \in Ω} p_{ω} z_{ω} \leq θ, z_{ω} \in {0, 1}, \\ \forall ω \in Ω . \end{matrix}$
computing an improved solution
$(P^{k + 1}, P^{c, k + 1}, y_{ω}^{k + 1}, z_{ω}^{k + 1})$
.

else
$(P^{k}, P^{c, k}, y_{ω}^{k}, z_{ω}^{k})$
is the desired optimal solution. endwhile

In the next section, we solve problem (21a) using Algorithm 1 and emphasize the following insights: 1.
Without knowledge of the true scenario‐specific demand distributions
$G_{t, ω}$
, one can guarantee the accurate solution of the optimization problem (11) using high service‐level requirements and the unperturbed demand distribution
$F_{t}$
. The accuracy is especially high when the realized demand
$G_{t, ω}^{- 1}$
is lower than the predicted demand
$F_{t}^{- 1}$
(Case 1).
2.
Underestimating the demand and, consequently, requiring a lower service level increases the optimality gap (see Case 2 and the case study on Omega× Swatch MoonSwatches in Section 4).
3.
In the setup with production disruptions, service‐level requirements can be kept high by utilizing safety stock. Nevertheless, the difficulty associated with this is that the safety stock grows exponentially for service levels above about 95% (Simchi‐Levi, 2010). The solution to this issue is to produce more during the periods with lower disruption probabilities, or during the periods with disruption scenarios included in the risk budget (i.e., neglected).

MANAGERIAL IMPLICATIONS

In this section, we emphasize the managerial importance of maintaining high service‐level requirements, the fulfillment of which results in accurate production decisions even when estimated demand distributions
$F_{t, ω}$
are offset with respect to realized functions
$G_{t, ω}, \forall ω$
. This also holds true if the manufacturer assumes the distributions
$F_{t, ω}$
to be equal for all disruption scenarios, that is,
$F_{t} = F_{t, ω} \forall ω$
as in the case when no scenario‐specific customer research can be performed.

Considering a quadratic production function
${\bar{C}}_{t}^{prod} (P_{t})$
, we naturally distinguish between the average cost of producing smaller and larger quantities. Overall, manufacturers often face different unit production costs depending on the amount produced. For example, this is consistent with the law of diminishing marginal returns for larger quantities, which is especially true in the presence of production disruptions, inefficient labor, and/or inappropriate or: misaligned technology selection (e.g., using 3D‐printing technologies can increase the cost of production). We account for both inventory and production capacity constraints, which allow the manufacturer to hold and produce only a limited number of products. However, the decision maker may keep a safety stock on hand to prevent stock‐outs up to the scenario‐specific service level. Furthermore, unsatisfied cumulative demand beyond the specified service level may result in a backlog.

We start by analyzing different types of disruption scenarios and their effects on the optimal production plan. The occurrence and duration of a disruption at time
$t_{ω}^{start}$
depend on the choice of disruption and recovery probabilities
$ζ_{t}$
and
$χ_{t} \forall t$
in the Markov process. In the simplest case, we distinguish disruptions that are more probable early, in the middle, or late in the planning horizon by varying probabilities
$ζ_{t} = ζ$
and
$χ_{t} = χ, \forall t$
: (i)
Early disruptions (
$ζ_{t} = ζ = 0.2$
and
$χ_{t} = χ = 0.2, \forall t$
): scenarios in which the disruption tends to start early and persist;
(ii)
Middle‐start disruptions (
$ζ_{t} = ζ = 0.1$
and
$χ_{t} = χ = 0.2, \forall t$
): scenarios with persisting disruptions tending to start at the middle of the horizon;
(iii)
Late disruptions (
$ζ_{t} = ζ$
and
$χ_{t} = χ = 0.2, \forall t$
): scenarios in which production disruptions tend to start late in the planning horizon.
Furthermore, we distinguish disruptions by their expected duration
$E (s)$
, which can be computed in line with Equation (25) for an infinite state Markov process with constant
$ζ_{t} = ζ$
and
$χ_{t} = χ, \forall t$
.

Figure 8 demonstrates the influence of the starting period and the disruption's expected duration on the optimal value of the optimization problem (21a). Figure 8a shows that disruptions with higher probability to start earlier tend to be more costly. This is also consistent with the fact that multiple disruptions are possible over the planning horizon. Furthermore, as demonstrated in Figure 8b, longer disruptions result in higher production costs. Naturally, longer disruptions provide shorter periods of production recovery, putting the manufacturer at a disadvantage.

FIGURE 8
Optimal values dependent on the expected start and duration of disruptions

Further, Figures 9 and 10 show that the optimal mitigation strategy for scenarios with a high periodic disruption probability
$ζ_{t}$
and a long expected duration
$E (s)$
necessitates stockpiling inventory to reduce the loss at time T and the overall backlog. Indeed, because the time to prepare for such disruptions is limited, inventory amounts for early and persistent disruptions obtained via problems (22) and (23) are close to each other. Furthermore, the backlog cost is high because the disruptions are long. Thus, the regret contributing to the EOC tends to be high. Also, the cumulative probability of such scenarios is high.

FIGURE 9
Mitigation inventories dependent on the expected start and duration of disruptions

FIGURE 10
Optimal backlog dependent on the expected start and duration of disruptions

As a result, mitigating for such scenarios creates a significant peak in inventory level and necessitates the keeping of a safety stock. Due to the availability of such inventory, the backlog tends to be lower for early and longer disruptions in the second half of Figure 10's planning horizon.

Next, we assess the impact of service‐level requirements on the level of inventory, backlog, and safety stock amounts (Figures 11 and 12). The manufacturer holds more inventory earlier in the planning horizon in order to maintain the same service level as in the absence of disruptions (Figure 11). This matches the behavior of the inventory policy early in the planning horizon in Figure 9. Note that the inventory is at the limit of available capacity (60 units in Figure 11) for service levels approaching perfect customer satisfaction (i.e., beyond 96%). Also, because the safety stock is designed to prevent the majority of stock‐outs early in the planning horizon, its level is increasing in service level given the high probability for early disruptions (Figure 12a). When the service level reaches approximately 95%, the safety stock starts growing exponentially (Simchi‐Levi, 2010).

FIGURE 11
Comparison of optimal inventories with and without disruptions

FIGURE 12
Influence of service‐level requirements on optimal safety stock and backlog amounts

Case study on the Omega × Swatch MoonSwatches

Importantly, we assess the impact of doubly probabilistic service‐level requirements on the optimal cost in problem (21a) using Omega× Swatch MoonSwatches as an example. Eleven new watches devoted to bodies in our solar system (Moon, Earth, Sun, Mars, Mercury, Neptune, Pluto, Uranus, Venus, Jupiter, and Saturn) were released on March 26, 2022, in selected Swatch stores and were followed by a massive and unexpected demand, to the point where Swatch needed to provide clarifications on social media. There could be a variety of reasons for such demand, including the Omega branding on the watches, as well as a price of only US$260 per piece instead of about US$6,600 for the original Omega Speedmaster. Considering the time horizon of
$T = 12$
days in a selected Swatch store, and assuming the cost of production to be lower than US$260, we model the influence of realized demand on the optimal total cost in problem (21a).

We specifically consider two layers of uncertainty in the problem: (1) demand variability and (2) production disruptions (e.g., due to a lack of specific color pigments used to produce some of the MoonSwatches), implying that the service level cannot be maintained for some scenarios with probability β. The second layer of uncertainty impacts the service‐level constraint and takes the variability in the unperturbed demand distributions
$F_{t} = F_{t, ω} \forall ω$
into account in line with Theorem 2. In Figures 13 and 14, we conduct a sensitivity analysis, allowing the demand distribution quantile to change as
$(1 + λ ξ) F_{t}^{- 1} (1 - {\bar{α}}_{t})$
. Here, ξ is the random variable influencing production capacity as described in Section 3.1 and
$λ \in [- 1, 1]$
is a multiplier, which defines the direction of the change in demand distribution.

FIGURE 13
Comparison of optimal values for optimization problems with and without a perturbation in demand

FIGURE 14
Optimal value and average backlog in the optimization problem (21a) given different demand offsets and risk thresholds

As shown in Figure 13a, a negative multiplier λ, implying that a lower production should suffice for the demand, results in a tight optimality gap for the problem (21a) with a fixed probability threshold β. This is in line with Case 1 described in Section 2. Conversely, a positive multiplier λ results in a larger optimality gap and, thus, an underestimation of total costs for a given service level. This is the case for Omega× Swatch MoonSwatches, for which the demand was initially underestimated and described by Swatch as “phenomenal” when observed. Note that the size of the optimality gap and, thus, the cost increase is dependent on the quantile offset and increases on average for higher service levels (see Case 2 and Figure 14a, where the optimal cost reaches about US$24,500 for the average demand of 60 pieces per store and a service level of 99%). Moreover, as demonstrated in Figures 13b and 14, larger‐than‐expected backlogs resulting from demand underestimation and a lack of safety stock may remain unresolved until the end of the planning horizon if no additional actions, such as changing the pricing strategy or expanding the production, are taken.

In conclusion, we numerically test the sensitivity of the results to changes in the demand distribution mean. We observe that the growing mean demand naturally increases the necessary level of safety stock (Figure 15a), backlog in later time periods (Figure 15b), and, thus, the lost sales at time T. In the Omega× Swatch case and in the absence of sufficient safety stock, optimal backlog plays an important role due to the fact that a large portion of customers are willing to return and purchase the watch for US$260. Nevertheless, Swatch will have to reevaluate their long‐term decision, which could include limiting the edition or increasing the price to control demand.

FIGURE 15
Influence of mean demand on the optimal strategy

CONCLUSION AND OUTLOOK

Optimization models accounting for production uncertainties and capable of determining resilient strategies are crucial for the mitigation of supply chain disruptions. In this article, we account for two layers of uncertainty that affect the production pipeline. The first layer results from the variability in demand and demand distribution. The second layer is triggered by production disruptions, when the service level cannot be maintained with some probability. We develop a bilevel optimization model that takes both demand uncertainty and production disruptions into account by using chance constraints and doubly probabilistic service‐level requirements that incorporate the dependence of demand distribution and production disruptions. Our approach is particularly advantageous for manufacturers who wish to build resilient production strategies that do not require significant adjustments in the event of a disruption. Even though the manufacturer has the option of focusing solely on lowering the expected total cost, the corresponding optimal strategy would need to be adjusted in the event of a disruption. By doing so, one deviates from the expected optimum and observes the realized (and, potentially, extremely high) cost dependent on the disruption scenario. Clearly, the manufacturer can use a robust framework and prepare for the worst‐case disruption guaranteeing that the realized cost is lower than estimated with a high probability. Nevertheless, this solution may be overly conservative. Aside from these options, there is one more possibility that we address in this article. Naturally, the manufacturer can constrain the difference in optimal costs between a strategy that can be implemented quickly in the event of a disruption (i.e., reactive setup) and a benchmark, which is an idealistic, but unattainable production schedule (i.e., anticipative setup). Under conditions of Lipschitz‐continuity, our approach bounds the difference in production schedules between the above‐mentioned strategies, making the production more predictable. The cost difference between the production schedule implemented at the start of the planning horizon and the reactive strategy is also bounded in this case.

Approximating our problem for numerical purposes, we demonstrate theoretically that the optimality gap is tight for high service‐level requirements when disruptions reduce demand. In order to solve the resulting problem numerically, we customize the Benders decomposition procedure by introducing novel feasibility cuts and using warm‐start techniques for computational speed‐ups. In the article's managerial section and by using Omega× Swatch MoonSwatches as an example, we demonstrate that the best approximation quality with a desirable doubly probabilistic service level can be achieved for disruptions with a drop in demand.

Footnotes
ORCID
Anna Timonina‐Farkas
René Y. Glogg
Ralf W. Seifert
References
1.
Arrow K. J. Karlin S. Scarf H. (1958). Studies in the mathematical theory of inventory and production. Stanford University Press.
2.
Ben‐Tal A. Nemirovski A. (1988). Robust convex optimization. Mathematics of Operations Research, 23(4), 769–805.
3.
Bertsekas D. P. (1998). Network optimization: Continuous and discrete models. Athena Scientific.
4.
Chopra S. Sodhi M. (2004). Supply‐chain breakdown. MIT Sloan Management Review, 46(1), 53–61.
5.
Ciarallo F. W. Akella R. Morton T. E. (1994). A periodic review, production planning model with uncertain capacity and uncertain demand—Optimality of extended myopic policies. Management Science, 40(3), 320–332.
6.
Craig N. DeHoratius N. Raman A. (2016). The impact of supplier inventory service level on retailer demand. Manufacturing & Service Operations Management , 18(4), 461–474.
7.
Dong L. Tomlin B. (2012). Managing disruption risk: The interplay between operations and insurance. Management Science, 58(10), 1898–1915.
8.
Fischetti M. Ljubić I. Sinnl M. (2016). Benders decomposition without separability: A computational study for capacitated facility location problems. European Journal of Operational Research, 253(3), 557–569.
9.
Geoffrion A. M. (1972). Generalized Benders decomposition. Journal of Optimization Theory and Applications, 10(4), 237–260.
10.
Granot F. Veinott A. F. Jr (1985). Substitutes, complements and ripples in network flows. Mathematics of Operations Research, 10(3), 471–497.
11.
Hendricks K. B. Jacobs B. W. Singhal V. R. (2019). Stock market reaction to supply chain disruptions from the 2011 Great East Japan Earthquake. Manufacturing & Service Operations Management , 22(4), 683–699.
12.
Hochrainer‐Stigler S. Timonina‐Farkas A. Silm K. Balkovič J. (2019). Large scale extreme risk assessment using copulas: An application to drought events under climate change for Austria. Computational Management Science, 16(4), 651–669.
13.
Kantorovich L. V. (2006). On the translocation of masses. Journal of Mathematical Sciences, 133(4), 1381–1382.
14.
Katok E. Thomas D. Davis A. (2008). Inventory service‐level agreements as coordination mechanisms: The effect of review periods. Manufacturing & Service Operations Management , 10(4), 609–624.
15.
Khang D. B. Fujiwara O. (1993). Multi‐period stochastic network flow problems with service level requirements: Approximation and error estimates. IIE Transactions, 25(2), 104–110.
16.
Lücker F. Chopra S. Seifert R. (2020). Mitigating product shortages due to disruptions in multi‐stage supply chains. Production and Operations Management, 30(4), 941–964.
17.
Nourelfath M. (2011). Service level robustness in stochastic production planning under random machine breakdowns. European Journal of Operational Research, 212(1), 81–88.
18.
Prékopa A. Boros E. (1991). On the existence of a feasible flow in a stochastic transportation network. Operations Research, 39(1), 119–129.
19.
Schmitt A. J. Snyder L. V. Shen Z.‐J. M. (2010). Inventory systems with stochastic demand and supply: Properties and approximations. European Journal of Operational Research, 206(2), 313–328.
20.
Simchi‐Levi D. (2010). Operations rules: Delivering customer value through flexible operations. MIT Press.
21.
Snyder L. V. Atan Z. Peng P. Rong Y. Schmitt A. J. Sinsoysal B. (2016). OR/MS models for supply chain disruptions: A review. IIE Transactions, 48(2), 89–109.
22.
Timonina A. Hochrainer‐Stigler S. Pflug G. Jongman B. Rojas R. (2015). Structured coupling of probability loss distributions: Assessing joint flood risk in multiple river basins. Risk Analysis , 35(11), 2102–2119.
23.
Timonina‐Farkas A. Katsifou A. Seifert R. W. (2020). Product assortment and space allocation strategies to attract loyal and non‐loyal customers. European Journal of Operational Research , 285(3), 1058–1076.
24.
Tomlin B. (2006). On the value of mitigation and contingency strategies for managing supply chain disruption risks. Management Science, 52(5), 639–657.
25.
Vajda S. (2014). Probabilistic programming. Academic Press.
26.
Veinott A. F. Jr (2005). Lectures in supply‐chain optimization. Management Science and Engineering , 361.
27.
Villani C. (2008). Optimal transport: Old and new. Vol. 338. Berlin: springer, 2009.
28.
Zheng Q. P. Shen S. Shi Y. (2015). Loss‐constrained minimum cost flow under arc failure uncertainty with applications in risk‐aware kidney exchange. IIE Transactions, 47(9), 961–977.