Robust condition‐based production and maintenance planning for degradation management

Abstract

We study the robust production and maintenance control for a production system subject to degradation. A periodic maintenance scheme is considered, and the system production rate can be dynamically adjusted before maintenance, serving as a proactive way of degradation management. Optimal control of the degradation rate aims to strike a balance between the risk of failure and the production profit. We first consider the scenario in which the degradation rate increases linearly with the production rate. Different from the existing literature that posits a parametric stochastic degradation process, we suppose that the degradation increment during a period lies in an uncertainty set, and our objective is to minimize the maintenance cost in the worst case. The resulting model is a robust mixed‐integer linear program. We derive its robust counterpart and establish structural properties of the optimal production plan. These properties are then used for real‐time condition‐based control of the production rate through reoptimization. The model is further generalized to the nonlinear production–degradation relation. Based on a real production–degradation dataset from an extruder system, we conduct comprehensive numerical experiments to illustrate the application of the model. Numerical results show that our model significantly outperforms existing methods in terms of the mean and variance of cost rate when degradation model misspecification is presented.

INTRODUCTION

Background and motivation

Most production systems are subject to degradation that ultimately leads to system failures. Specifically, degradation refers to a cumulative change in a system's performance characteristics over time. The system is deemed failed when its degradation level exceeds a failure threshold (Elwany et al., 2011; Chen et al., 2022). To mitigate the failure risk, a common practice is to perform condition‐based maintenance (CBM) that preemptively makes a system overhaul based on the in situ system degradation level (Kouvelis et al., 2005). CBM relies on monitoring the degradation and predicting the system residual useful life by assuming stochastic system deterioration driven by extrinsic environments or random usage. Preventive maintenance is then called for when the degradation level is higher than a threshold. However, the CBM paradigm ignores the fact that future usage of the system that drives system degradation might be controllable. Moreover, last‐minute maintenance scheduling can be infeasible because of staffing issues and the preparation of tools and parts.

Instead of passively waiting for the degradation to hit the threshold, the failure risk can also be mitigated if we actively adjust the workload to the system to manage the system degradation before a prefixed replacement epoch. This strategy is attractive by virtue of its fixed replacement time and is inherently feasible for production systems because we can dynamically control the production or manufacturing rate, for example, adjusting the speed of a stamping machine (Hao et al., 2015) and the filtration rate of a rapid gravity filter in a waterworks (Sun et al., 2019). Recent developments of the Internet‐of‐Things (IoT) technology further facilitate a remote adjustment of the production rate (uit het Broek et al., 2020). Generally, the system degradation rate increases in the production rate, giving rise to a higher failure risk but a higher production profit. Hence, the workload adjustment essentially strikes a balance between production revenues and failure costs.

Despite the benefits of dynamic production rate adjustment, the joint planning of production and maintenance is generally challenging due to the stochastic nature of degradation. To capture the uncertainty in degradation, a common approach is to model the degradation as a stochastic process such as Wiener, gamma, or inverse Gaussian processes (Panagiotidou & Tagaras, 2010; Li & Ryan, 2011; Ye & Xie, 2015; Yildirim et al., 2016; Hu, Sun, Ye, & Ling, 2021; Zhang et al., 2022), based on which dynamic production planning can be numerically solved through the Markov decision process (MDP). When the degradation follows a process different from the presumed model, however, our simulation results in Section 6 reveal that such a model misspecification can increase the maintenance cost.

To address the above issues, we resolve the production and maintenance planning problem through robust optimization. The proposed model accounts for the inherent uncertainty of degradation, and the objective is to minimize the worst‐case maintenance cost. We establish structural properties of the model, which enable us to determine the optimal replacement time when optimizing the production system at time zero. Through reoptimization, these properties further allow online production control when in situ degradation signals are regularly monitored.

Related literature

A large amount of research effort has been dedicated to maintenance optimization problems. The early literature mainly focuses on planning maintenance based on system's failure rate (Dayanik & Gürler, 2002; Dogramaci & Fraiman, 2004). Due to the development of sensing technology, in situ degradation data are further considered in maintenance scheduling, and a substantial amount of work has been done that enriches the arsenal of CBM models. (See Arts et al., 2016; Khojandi et al., 2018; Abbou and Makis, 2019; Bensoussan et al., 2020; Zhu et al., 2021; Li and Tomlin, 2022; Wang, 2022; Zhao et al., 2022; Mookerjee and Samuel, 2023; Zhang and Zhang, 2023, for instance, and de Jonge and Scarf, 2020, for a review.) Several maintenance models are further developed for production systems by considering restricted time windows for production, potential short of raw materials, random production waits, or finite inventory capacity for finished goods (Iravani & Duenyas, 2002; Drozdowski et al., 2017; Hu, Sun, & Ye, 2021). However, the above works treat the production plan as predetermined and do not allow dynamic adjustment of the production rate based on the system degradation level. To tackle the problem, Sloan and Shanthikumar (2000), Batun and Maillart (2012), and Kouvelis et al. (2023) investigate the joint optimization of a production and maintenance program, but they all implicitly assume that the system degradation evolves independently of the production tasks.

To overcome these potential disadvantages, recent studies assume the degradation process could be adjusted by the controllable production rate. Sun et al. (2019) analyze a degradation dataset from a waterworks, where workloads of rapid gravity filters are dynamically adjusted. Basciftci et al. (2020) investigate a data‐driven maintenance and operations scheduling problem with a load‐dependent degradation process. Similar studies also include Yang et al. (2007), Hao et al. (2015), and Wiebe et al. (2018). All these studies show that the dynamic adjustment of the production rate can significantly reduce the operational cost and failure risk. Because of their complicated formulations, however, metaheuristic or approximation methods have to be adopted to deliver high‐quality solutions.

uit het Broek et al. (2020) and Drent et al. (2022) are the two most closely related works to ours. uit het Broek et al. (2020) initiate the condition‐based production planning and derive the optimal production policy when the degradation is a deterministic process. They further propose two heuristic policies when degradation is stochastic, but both lack theoretical support and require a brute‐force search for optimization. Drent et al. (2022) formulate the condition‐based production problem as a continuous‐time MDP for a discrete‐state system and derive insightful structural results about the optimal production policy. Compared with the two studies, our model deals with continuous degradation signals using a different modeling paradigm, that is, robust optimization, and our closed‐form solution for the optimal production rate allows us to further optimize the replacement time. The closed form makes our model easy to implement, and our model from the robust optimization perspective enriches the arsenal of the condition‐based production and maintenance planning problem that has become popular recently. The structural results derived from Sections 3–5 complement the theoretical results of condition‐based production from the perspective of robust optimization.

In summary, most existing studies on production planning assume a known relation between the workload and the degradation rate, which could be difficult to accurately determine or estimate in practice. Drent et al. (2022) pioneer in using Bayesian learning to solve this problem with known priors for model parameters. As degradation data are often scarce, degradation model misspecification can be common in practice. To tackle the challenges, robust optimization is a promising tool to protect against the potential misspecification of the degradation model and the production–degradation relation. Robust optimization has been widely used in operations management problems such as inventory control (Bertsimas & Thiele, 2006), sharing economy (He et al., 2020), portfolio selection (Lim et al., 2012), manufacturing (Bandi et al., 2015), etc. There is relatively scant literature on using robust optimization for reliability problems. Wang et al. (2020) formulate a distributionally robust optimization model to design a serial‐parallel system to meet a required reliability constraint. Their focus is to formulate the optimization problem such that the model can be efficiently solved by existing solvers such as CPLEX and Gurobi. Kim (2016) formulates a robust MDP to solve a maintenance optimization problem of a deteriorating system in the presence of model misspecification, and the structure of the optimal policy is established. The penalty approach in Kim (2016), if applied to our problem, leads to a dynamic program where analysis of the corresponding value function is intractable.

Overview and contributions

The objective of this study is to resolve the production control problem from the perspective of robust optimization, with a primary focus on establishing structural properties of the optimal solution. We discretize the time horizon and systematically investigate the optimization of the dynamic workload adjustment at the beginning of each time interval under both linear and nonlinear production–degradation relations. Our contributions can be summarized as follows.

Compared with the prior literature that assumes a deterministic degradation path, we allow the degradation to be uncertain, which is more realistic and challenging. As suggested in Section 6, our proposed robust optimization solution is more robust against model misspecification, even when the working degradation model to determine the uncertainty set is misspecified.

Under the robust optimization formulation, we successfully derive the structural optimality of the production planning problem given our uncertainty on the future degradation when the production–degradation relation is linear or convex. The structural results enable us to solve the optimization problem in real time without relying on existing solvers. Thanks to this feature, the optimal solution can also be updated in real time once a new degradation measurement is taken.

When the production–degradation relation is concave, the robust production‐planning problem becomes nonconvex (see Section 5.3), which is difficult to solve in general. Motivated by the structural results under the deterministic degradation setting, we propose an algorithm to approximately solve the optimization problem. Comprehensive numerical experiments reveal that the proposed algorithm beats the state‐of‐the‐art interior point algorithm.

We demonstrate our methods using a degradation dataset on an extruder system. Through comparisons with several existing methods, we show that the robust model outperforms these methods in terms of a lower mean, lower worst‐case value, and lower variation of the production cost rates when the stochastic degradation model is misspecified. The significantly smaller variability in the cost rate makes the robust approach attractive to risk‐averse decision‐makers without sacrificing much in terms of the average cost rate.

The remainder of this article is organized as follows. Section 2 states our production planning problem, where we take parameter uncertainty into consideration. Section 3 formulates a static robust optimization model and derives its robust counterpart under a linear production–degradation relation, which sets the stage for analysis. Section 4 then proposes a rolling‐horizon formulation and derives the structural properties of the optimal production solution. Section 5 further generalizes our model to the case of nonlinear production–degradation relations. Section 6 uses comprehensive numerical experiments to illustrate the advantages of our model. Section 7 concludes this study and discusses future research. All technical proofs are provided in Appendix A in the Supporting Information.

PROBLEM STATEMENT

Consider a production system subject to production‐induced degradation and periodic condition monitoring. The system has an initial degradation level x ₀ and fails when the degradation level exceeds a failure threshold L. The threshold L is assumed to be exogenous and normally determined by industrial standards (Li & Tomlin, 2022). The replacement interval is fixed at T. To strike a balance between the production profit and system reliability, we adopt the following joint maintenance and production control scheme from uit het Broek et al. (2020). The system is new at time 0, and its degradation is periodically inspected at time epochs

$\delta ,2\delta ,\text{…},N\delta$

, where δ is a fixed interinspection interval, and

$N\delta =T$

. Let

$x_n$

be the observed degradation level at

$n\delta$

,

$n=0,\text{…},N$

. If

$x_n\le L$

for all n, then the system is preventively replaced with a cost

$c_{\mathrm{r}}>0$

at time T. Otherwise, the system fails before T. As with uit het Broek et al. (2020), we assume the system is not correctively replaced until T, with a corrective replacement cost

$c_{\mathrm{r}}+c_{\mathrm{f}}$

, where

$c_{\mathrm{f}}>0$

is a failure penalty cost. A fixed replacement time is often desired in practice because it facilitates budget allocations of decision‐makers (Sun et al., 2019). The time T is optimized in Section 3.2. We call the time interval

$[(n-1)\delta ,n\delta )$

period n for

$n\in [N]\triangleq \{1,\text{…},N\}$

and use a constant production rate

$u_n$

during this period. We assume

$0\le u_n\le 1$

without loss of generality, and the production rates

$u_n$

are adjusted according to the following rules:

• At the beginning of period n, we select a feasible production rate

  $u_n\in [0,1]$
  and let the system produce at this rate throughout the whole period with a production profit
  $\pi u_n$
  per time unit, provided that the system does not fail during period n.

• If the system fails during period n, the production profit at this period is lost. Notationally, this amounts to enforcing

  $u_n=0$
  when any failure occurs during period n.

The second assumption above is for tractability in a discrete‐time setting. In practice, δ is typically small due to advances in the IoT technology, making it reasonable to impose the assumption.

At period

$n\in [N]$

, we model the system degradation rate as a power function of

$u_n$

:

$$h_n(u_n)=b_n{u}_n^{\eta }+a_n,$$

1where

$a_n,b_n,\eta \ge 0$

. The constant degradation rate is a reasonable approximation to the true degradation when the time interval δ is short enough. The production rate is commonly known as a covariate in the degradation literature, and the power law is a popular link function to relate a covariate to the degradation process (Hao et al., 2015; Chen & Ye, 2018). The shape parameter η is usually determined by the degradation mechanism and assumed to be a known constant (Lawless & Crowder, 2004). The values of

$a_n$

and

$b_n$

depend on some unobserved dynamic environmental factors, such as the temperature and humidity within the machinery system, variations of material properties in this period, and operators; thus, they are subject to uncertainty.

We now quantify the uncertainty of

$a_n$

and

$b_n$

to account for the uncertainty of the degradation process. Instead of specifying a distribution for

$(a_n,b_n)$

, we assume

$a_n$

and

$b_n$

can only take values within the intervals

$[{\mathtt{a}}_n-{\stackrel{\sim }{a}}_n,{\mathtt{a}}_n+{\stackrel{\sim }{a}}_n]$

and

$[{\mathtt{b}}_n-{\stackrel{\sim }{b}}_n,{\mathtt{b}}_n+{\stackrel{\sim }{b}}_n]$

, respectively (Bertsimas & Sim, 2004). Here,

${\mathtt{a}}_n$

and

${\mathtt{b}}_n$

represent the nominal values and

${\stackrel{\sim }{a}}_n$

and

${\stackrel{\sim }{b}}_n$

specify the largest possible deviation. We normalize

$a_n$

and

$b_n$

as

$y_n\triangleq (a_n-{\mathtt{a}}_n)/{\stackrel{\sim }{a}}_n$

and

$z_n\triangleq (b_n-{\mathtt{b}}_n)/{\stackrel{\sim }{b}}_n$

so

$y_n,z_n\in [-1,1]$

for all

$n\in [N]$

. Then it is sufficient to quantify the uncertainty of

$(a_n,b_n)$

through

$(y_n,z_n)$

. It is unlikely that all the

$(y_n,z_n)$

are close to 1 or −1. Motivated by this observation, we follow Bertsimas and Thiele (2006) and introduce a sequence of uncertainty budgets

${({\mathrm{\Gamma }}_n)}_{n\in [N]}$

to bound the possible realizations of

$(y_n,z_n)$

. We consider the following feasible uncertainty set for

$({\mathit{y}}_n,{\mathit{z}}_n)\triangleq (y_1,\text{…},y_n,z_1,\text{…},z_n)$

:

$$\begin{array}{c}\hfill {\mathcal{P}}_n\triangleq \left\{({\mathit{y}}_n,{\mathit{z}}_n)||y_j|\le 1,|z_j|\le 1,\forall j\in [n],\sum _{j=1}^n(|y_j|+|z_j|)\le {\mathrm{\Gamma }}_n\right\}.\end{array}$$

2Note that

$|y_n|+|z_n|\le 2$

for all

$n\in [N]$

by their definitions. We then assume that

${\mathrm{\Gamma }}_n$

increases in

$n\in [N]$

and

${\mathrm{\Gamma }}_n-{\mathrm{\Gamma }}_{n-1}\le 2$

for all

$n=2,\text{…},N$

. These assumptions imply that uncertainty accumulates over time, but the growth rate of the budgets of uncertainty

${\mathrm{\Gamma }}_n$

is bounded. Tractability is the main reason for adopting a polyhedral uncertainty set as in (2) rather than an ellipsoidal set. As will be shown in Section 4, the linear structure of (2) enables us to derive closed‐form expressions of the optimal robust production solution. Moreover, Section 6 demonstrates that reasonable

${\mathtt{a}}_n,{\mathtt{b}}_n,{\stackrel{\sim }{a}}_n,{\stackrel{\sim }{b}}_n,{\mathrm{\Gamma }}_n$

,

$n\in [N]$

, can be readily determined based on historical degradation data, making our model easy to implement in practice.

Under the above formulation, we aim to determine the optimal replacement time

$T^{\ast }\triangleq N^{\ast }\delta$

to minimize the worst‐case system operational cost rate (maximizing the worst‐case system revenue rate) at time 0, as discussed in Section 3. The worst‐case cost rate at time 0 for any fixed production schedule is the maximum possible cost rate during the planning horizon computed immediately before production, which only depends on N, x ₀, and

${({\mathrm{\Gamma }}_n)}_{n\in [N]}$

. After

$N^{\ast }$

is chosen and the system is put into production, our setting assumes that the system degradation

$x_1,\text{…},x_{N^{\ast }}$

is monitored at the beginning of each period. Based on these in situ degradation measurements, the corresponding worst‐case cost rates estimated at

$\delta ,\text{…},N^{\ast }\delta$

will change over time. As such, we design an efficient reoptimization algorithm to dynamically adjust the production rate based on in situ degradation data, and the optimal production rate sequence

${\mathit{u}}_{N^{\ast }}^{\ast }\triangleq {(u_n^{\ast })}_{n\in [N^{\ast }]}$

is obtained in Section 4.

STATIC ROBUST OPTIMIZATION PROBLEM

This section confines to a linear relation between the production and degradation rate, that is,

$\eta =1$

in (1). In Section 3.1, we fix N and formulate a static robust optimization problem to minimize a system's worst‐case operational cost at time 0 where only x ₀ is known and

$x_n$

,

$n\in [N]$

, are not observed yet. The static model will be used to determine the optimal N in Section 3.2, and it sets the stage for reoptimization in Section 4.

Robust counterpart

We consider a fixed N and thus fixed

$T=N\delta$

. At time 0, let

$C_N({\mathit{u}}_N)$

be the worst‐case total operational cost during

$[0,N\delta ]$

given N and production control variables

${\mathit{u}}_N$

. The optimal

${\mathit{u}}_N$

that minimizes

$C_N({\mathit{u}}_N)$

can be obtained by solving the following mixed‐integer program at time 0:

3a

3b

$$\begin{array}{c}\hfill u_n\le 1-g_n,\forall n\in [N],\end{array}$$

3c

$$\begin{array}{c}\hfill 0\le u_n\le 1,\forall n\in [N],\end{array}$$

3d

$$\begin{array}{c}\hfill g_n\in \{0,1\},\forall n\in [N],\end{array}$$

3e

where

${\mathit{g}}_N\triangleq {(g_n)}_{n\in [N]}$

and M is a sufficiently large constant. The objective (3a) is the total operational cost, which is the difference between the maintenance cost and the production revenue. Constraints (3b) are based on the worst case of the following system dynamics:

$$x_n=x_{n-1}+\delta (b_n{u}_n+a_n)=x_0+\delta \sum _{j=1}^n(b_j{u}_j+a_j),n\in [N],$$

for any

$a_n$

and

$b_n$

in the uncertainty set. With the help of a sufficiently large constant M, the auxiliary binary decision variable is forced to take

$g_n=1$

if the worst

$x_n$

exceeds the failure threshold L. Constraints (3c) enforce the production rate

$u_n$

to be zero if the system fails before or during period n. Constraints (3d) and (3e) specify the support of the decision variables. The uncertainties in the degradation process affect Constraints (3b). By strong duality of linear programming, we then derive the robust counterpart of the above optimization problem as follows. Lemma 1

Problem (3a)–(3e) has the following robust counterpart:

4a

4b

$$\begin{array}{c}\hfill r_{nj}+q_n\ge {\stackrel{\sim }{b}}_j{u}_j,\forall j\in [n],n\in [N],\end{array}$$

4c

$$\begin{array}{c}\hfill s_{nj}+q_n\ge {\stackrel{\sim }{a}}_j,\forall j\in [n],n\in [N],\end{array}$$

4d

$$\begin{array}{c}\hfill u_n\le 1-g_n,\forall n\in [N],\end{array}$$

4e

$$\begin{array}{c}\hfill 0\le u_n\le 1,\forall n\in [N],\end{array}$$

4f

$$\begin{array}{c}\hfill g_n\in \{0,1\},\forall n\in [N],\end{array}$$

4g

4h where

${\mathit{q}}_N\triangleq {(q_n)}_{n\in [N]}$

,

${\mathit{s}}_N\triangleq {(s_{nj})}_{j\in [n],n\in [N]}$

, and

${\mathit{r}}_N\triangleq {(r_{nj})}_{j\in [n],n\in [N]}$

.

Problem (4a)–(4h) is a mixed integer program that can be solved by off‐the‐shelf solvers. We can further strengthen the integer programming formulation based on the following proposition. Proposition 1

At optimality of Problem (4a)–(4h), the variables

${(g_n)}_{n\in [N]}$

satisfy

$$g_n\le g_{n+1},\forall n\in [N-1].$$

5

The proposition is intuitive in the sense that the optimal

$g_n$

is an indicator of system failure. If a system fails before or during period n (i.e.,

$g_n=1$

), then it cannot be functional during the remaining periods without maintenance (i.e.,

$g_k=1$

for

$k=n+1,\text{…},N$

). Our numerical experiments show that including Constraints (5) in the root node of the branch‐and‐cut algorithm improves its numerical performance.

The robust counterpart formulated in (4a)–(4h) introduces auxiliary decision variables

${\mathit{q}}_N$

,

${\mathit{s}}_N$

, and

${\mathit{r}}_N$

to account for parameter uncertainty. To derive analytical insights into this formulation, we show that our robust optimization problem is equivalent to a production planning problem for an (artificial) deterministic system. To this end, let

${(q_n^{\ast })}_{n\in [N]}$

,

${(s_{nj}^{\ast })}_{j\in [n],n\in [N]}$

and

${(r_{nj}^{\ast })}_{j\in [n],n\in [N]}$

be the optimal solutions to the mixed‐integer program in Lemma 1. Define

$A_0\triangleq 0$

and

$$A_n\triangleq \delta \left[\sum _{j=1}^n\left(s_{nj}^{\ast }+r_{nj}^{\ast }\right)+q_n^{\ast }{\mathrm{\Gamma }}_n\right],n\in [N],$$

which can be understood as the worst‐case deviation of the degradation level at time

$n\delta$

from its nominal values by strong duality. Denote

6as the modified degradation level upon all decision epochs. Moreover, define

$${\overline{u}}_n\triangleq \min \left\{1,\underset{i=n,\text{…},N}{\min }\left\{(r_{in}^{\ast }+q_i^{\ast })/{\stackrel{\sim }{b}}_n\right\}\right\},n\in [N],$$

as the modified upper bound of the production rate. Then we can interpret Problem (4a)–(4h) by considering an artificial production system, such that (i) its production rate

$u_n$

takes values in the range

$[0,{\overline{u}}_n]$

during period

$n\in [N]$

and (ii) when the production rate is fixed at

$u_n$

during period n, the system degradation level evolves deterministically according to the modified degradation level in Equation (6). Then, the optimal production rate sequence for the above deterministic system is the same as that for the robust counterpart in (4a)–(4h). This observation suggests that we can derive the optimal production plan for our static robust optimization problem by analyzing a deterministic system. By the definition of

${\overline{u}}_n$

, the modified upper bound of the production rate is not greater than the original upper bound of value 1. This implies that decision‐makers should produce at a lower rate to protect system reliability against uncertainty. For comparison, consider the degradation levels under nominal parameters:

$${\mathtt{x}}_0=x_0\text{and}{\mathtt{x}}_n={\mathtt{x}}_{n-1}+\delta ({\mathtt{b}}_n{u}_n+{\mathtt{a}}_n).$$

7Routine algebra then shows that

$x_n^{\prime }={\mathtt{x}}_n+A_n\ge {\mathtt{x}}_n$

, where the inequality holds since

$A_n\ge 0$

for all

$n\in [N]$

. This suggests that a production plan that uses the nominal degradation level

${\mathtt{x}}_n$

would underestimate the underlying risk of the system failure when there is inherent uncertainty.

Interestingly, the equivalence revealed above is closely related to the two heuristic policies in uit het Broek et al. (2020). Both heuristics simplify the problem by assuming the degradation evolves deterministically with nominal parameters

${\mathtt{a}}_n$

and

${\mathtt{b}}_n$

as in (7). Then at time

$(n-1)\delta$

:

• The first heuristic uses the current degradation level

  $x_{n-1}$
  as an input and obtains the optimal production rate
  ${\mathtt{u}}_n^{\ast }$
  for the deterministic system with nominal values
  $({\mathtt{a}}_k,{\mathtt{b}}_k)$
  ,
  $n\le k\le N$
  . The production rate is then set to be
  ${\mathtt{u}}_n^{\ast }-{\alpha }_1$
  for some tuning parameter
  ${\alpha }_1\ge 0$
  .

• The second heuristic modifies the current degradation level

  $x_{n-1}$
  to
  $x_{n-1}+{\alpha }_2$
  for some tuning parameter
  ${\alpha }_2\ge 0$
  . The production rate is then set to be the optimal value for the deterministic system with nominal values
  $({\mathtt{a}}_k,{\mathtt{b}}_k)$
  ,
  $n\le k\le N$
  , and initial degradation level
  $x_{n-1}+{\alpha }_2$
  .

Basically, the first heuristic directly decreases the production rate by an amount of α₁, which resembles the (smaller) modified upper bound

${\overline{u}}_n$

of the production rate. The second heuristic adds a safety margin of value α₂ to the observed degradation level to prevent unexpected failure under uncertainty. This coincides with the (higher) modified degradation level

$x_n^{\prime }$

in (6). A brute‐force search has to be implemented to find the optimal tuning parameters

${\alpha }_1^{\ast }$

and

${\alpha }_2^{\ast }$

for the heuristics. For the proposed model, the values of

${\overline{u}}_n$

and

$A_n$

can be readily acquired from solving the robust counterpart in Problem (4a)–(4h), for example, by Gurobi or CPLEX. Under some mild conditions, these variables even admit closed‐form expressions (see Section 4.2). Therefore, our model is more computationally efficient and lends theoretical support to the two heuristics in uit het Broek et al. (2020) from the perspective of robust optimization. Moreover, the tuning parameters α₁ and α₂ for the heuristics in uit het Broek et al. (2020) are constant across periods. By contrast, our model is more flexible in the sense that both

$A_n$

and

${\overline{u}}_n$

are allowed to be dependent on n.

Optimal time to replacement

The robust counterpart in Section 3.1 is derived for a fixed number of periods N and thus fixed replacement interval

$T=N\delta$

. It can be optimized to minimize the long‐term cost rate, so the replacement time T would neither be too long nor too short in practice. Given the optimal

${\mathit{u}}_N^{\ast }$

to Problem (4a)–(4h) for a fixed N, the worst‐case cost rate is

$C_N({\mathit{u}}_N^{\ast })/(N\delta )$

. This function is generally nonconvex in

$N\in {\mathbb{N}}_{+}$

, and enumeration can be used to find the optimal

$N^{\ast }$

that minimizes

$C_N({\mathit{u}}_N^{\ast })/(N\delta )$

since N can only take positive integer values. To do so, we need to specify a range for the above one‐dimensional search. Under some mild assumptions, we are able to identify a bounded range that contains

$N^{\ast }$

. We list these assumptions as follows, which will also be used to derive the optimal dynamic production solution in Section 4.2. Assumption 1

The nominal and the maximum deviation parameters are constant over n, that is,

$${\mathtt{a}}_n\equiv \mathtt{a},{\mathtt{b}}_n\equiv \mathtt{b},{\stackrel{\sim }{a}}_n\equiv \stackrel{\sim }{a},{\stackrel{\sim }{b}}_n\equiv \stackrel{\sim }{b},\text{for all}n\in \mathbb{N}.$$

8

Assumption 2

The optimal production rates

${\mathit{u}}_{N^{\ast }}^{\ast }$

under the optimal

$N^{\ast }$

satisfies

$$\frac{C_{N^{\ast }}({\mathit{u}}_{N^{\ast }}^{\ast })}{N^{\ast }\delta }<0.$$

Assumption 1 requires constant nominal values and maximum deviations during the entire planning horizon. The assumption is commonly adopted in maintenance optimization (Chen & Ye, 2018; uit het Broek et al., 2020) and mild when the degradation process has stationary increments under a fixed production rate, which is reasonable for degradation resulting from wear, for example, filters in a waterworks and car tires (Sun et al., 2020). When the degradation process is nonlinear, taking the logarithm of the degradation data often yields a linear path, for example, rolling bearings studied in Elwany et al. (2011). Assumption 2 requires positive (negative) optimal revenue (cost). This assumption holds for reasonable parameter settings; otherwise, the optimal strategy is simply choosing to halt production. Based on Assumptions 1 and 2, we can derive the following results. Lemma 2

Under Assumptions 1 and 2, the optimal solution to Problem (4a)–(4h) satisfies

$g_{N^{\ast }}^{\ast }=0$

for the optimal

$N^{\ast }$

if the production–degradation relation follows (1) with

$\eta =1$

for all

$n\in \mathbb{N}$

.

Since

$g_{N^{\ast }}^{\ast }=0$

indicates that the system is functional upon replacement even in the worst case, Lemma 2 tallies with the common practice that decision‐makers would not plan a system failure at the beginning of the production period. Under Assumption 1, let

where

$x_{+}\triangleq \max \{x,0\}$

. The variable

$\underline{N}$

can be understood as the maximum possible number of periods such that the worst‐case degradation level is still lower than the failure threshold L even if the system operates under the maximum production rate throughout

$[0,\underline{N}\delta ]$

. We can then derive the following proposition based on Lemma 2. Proposition 2

Under Assumptions 1 and 2, the optimal

$N^{\ast }$

satisfies

$N^{\ast }\ge \underline{N}$

. Further we assume that there exists

$N^{†}$

such that

${\mathrm{\Gamma }}_n<n$

for all

$n>N^{†}$

. Then we have

$N^{\ast }\le \overline{N}$

for some

$\overline{N}\in {\mathbb{N}}_{+}$

satisfying

$\overline{N}\ge N^{†}$

, and the optimal solution to Problem (4a)–(4h) for

$N=\overline{N}$

satisfies either of the following two conditions: (i)

$g_{\overline{N}}^{\ast }=1$

; (ii)

$g_{\overline{N}}^{\ast }=0$

, the corresponding modified degradation level is

$x_{\overline{N}}^{\prime }=L$

, and the optimal production rates obtained at time 0 satisfy

$u_n^{\ast }\le \stackrel{\sim }{a}/\stackrel{\sim }{b}$

for all

$n\in [\overline{N}]$

.

Proposition 2 yields a range

$[\underline{N},\overline{N}]$

for the optimal number of planning periods

$N^{\ast }$

. Motivated by the central limit theorem,

${\mathrm{\Gamma }}_n$

is usually chosen to satisfy

${\mathrm{\Gamma }}_n=\mathcal{O}(\sqrt{n})$

, which ensures the existence of such an

$N^{†}$

as required in the proposition. Moreover, when

$\mathtt{a}>0$

, it is easy to see that the condition

$g_N^{\ast }=1$

holds for sufficiently large N. Thus, we can safely use the conditions in Proposition 2 to terminate the enumerative search when we find

$\overline{N}$

that satisfies the condition in Proposition 2.

PRODUCTION CONTROL IN A DYNAMIC ENVIRONMENT

We continue with the linear production–degradation with

$\eta =1$

in (1). With

$N^{\ast }$

from Section 3, we fix the optimal replacement time

$T^{\ast }=N^{\ast }\delta$

and use reoptimization for production control in a dynamic environment based on in situ degradation signals. To ease notation, this section uses N to denote the total number of periods instead of

$N^{\ast }$

. Section 4.1 formulates the reoptimization problem. Section 4.2 then derives structural properties of the reoptimization problem, which enable us to find the optimal solution for production control without using existing solvers.

Reoptimization problem in a rolling horizon

A common approach for the optimization problem in a dynamic environment is to use reoptimization whenever new information, which in our case is the newly inspected degradation level, becomes available. After the in situ degradation measurement

$x_n$

is taken at the beginning of period

$n+1$

,

$n\in [N-1]$

, the reoptimization is run. We feed

$x_n$

into the reoptimization problem to obtain the optimal production rates

$u_{n+1}^{\ast },\text{…},u_N^{\ast }$

, and set

$u_{n+1}^{\ast }$

as the production rate for period

$n+1$

. This process is repeated until period N. This so‐called rolling‐horizon implementation is widely used in dynamic operations management, for example, call center routing (Chan et al., 2014), repositioning of vehicle sharing systems (He et al., 2020), and inventory dual sourcing (Sun & Van Mieghem, 2019). All numerical experiments therein show a good performance of the reoptimization.

Theoretically, this reoptimization method is a type of static policy in multistage robust optimization problems, which often have good performance in practice (Delage & Iancu, 2015). Similar to Bertsimas and Thiele (2006) and Mamani et al. (2017), the main purpose of our model formulation in (3) is to derive closed‐form solutions that can be easily understood by practitioners. Alternative methods such as affinely adjustable robust optimization and distributionally robust optimization (e.g., Jackson et al., 2019; He et al., 2020) may not admit such closed‐form solutions, which hinders the real‐time implementation of the model for IoT‐enabled production systems.

To reoptimize our production planning problem, the budget of uncertainty at every period needs to be updated. Let

$\{{\mathrm{\Gamma }}_j^{(n)}:j>n\}$

be the sequence of uncertainty budget for the reoptimization at time

$n\delta$

. As with Bertsimas and Thiele (2006), we adopt the following updating scheme. First, we let

${\mathrm{\Gamma }}_n^{(0)}={\mathrm{\Gamma }}_n$

,

$n\in [N]$

, for Problem (4a)–(4h) to highlight that the problem is solved at time 0. At time

$n\delta$

, we then update

${\mathrm{\Gamma }}_j^{(n)}={\mathrm{\Gamma }}_j^{(n-1)}-{\mathrm{\Gamma }}_n$

,

$j=n+1,\text{…},N$

, for the remaining

$N-n$

periods to go. Routine algebra then shows

$0\le {\mathrm{\Gamma }}_j^{(n)}\le 2(j-n)$

. Such an update scheme maintains the uncertainty budget allocated to each period (Bertsimas & Thiele, 2006). However, the structural properties derived in Section 4.2 do not depend on any specific update scheme on

${\mathrm{\Gamma }}_n$

.

After updating

$\{{\mathrm{\Gamma }}_j^{(n)}:j>n\}$

, the reoptimization problem at time

$n\delta$

,

$n\in [N-1]$

, is given by

$$\begin{array}{c}\hfill \min c_{\mathrm{r}}+c_{\mathrm{f}}{g}_N-\pi \delta \sum _{j=n+1}^N{u}_j\end{array}$$

9a

9b

9c

9d

$$\begin{array}{c}\hfill u_j\le 1-g_j,\forall j=n+1,\text{…},N,\end{array}$$

9e

$$\begin{array}{c}\hfill 0\le u_j\le 1,\forall j=n+1,\text{…},N,\end{array}$$

9f

$$\begin{array}{c}\hfill g_j\le g_{j+1},\forall j=n+1,\text{…},N-1,\end{array}$$

9g

$$\begin{array}{c}\hfill g_j\in \{0,1\},\forall j=n+1,\text{…},N,\end{array}$$

9h

9i

In the same spirit as Problem (4a)–(4h), Problem (9) is a mixed‐integer linear program, which is generally time‐consuming to solve. Nevertheless, the next section derives some structural properties to this reoptimization problem, based on which closed‐form solutions are available. This feature makes the implementation of our model as easy as the heuristics in uit het Broek et al. (2020).

Structural properties

Solving (9) directly is generally hard. Nevertheless, we have

$g_{n+1}=0$

implies

$g_n=0$

for all

$n\in [N-1]$

by the cuts derived in Proposition 1. Hence, the feasible region of the binary indicators

${(g_n)}_{n\in [N]}$

is small, which enables us to enumerate all values of

${(g_n)}_{n\in [N]}$

for optimization. With fixed

${(g_n)}_{n\in [N]}$

, Problem (9) reduces to a linear program (LP). The following theorem formalizes the above idea that decomposes Problem (9) into several LPs. Theorem 1

At time

$n\delta$

, let

$C_{n^{\prime }}$

,

$n^{\prime }=n+1,\text{…},N$

, be the optimal value of the following LP

10Then the optimal objective value of Problem (9) at time

$n\delta$

equals

${\min }_{n^{\prime }=n+1,\text{…},N}{C}_{n^{\prime }}$

.

By convention, we set

$C_{n^{\prime }}=+\infty$

if the above LP is infeasible for any

$n^{\prime }=n+1,\text{…},N$

. Essentially, each LP aims to prevent the system failure (when

$n^{\prime }=N$

) or deliberately induce a failure at period

$n^{\prime }$

(when

$n^{\prime }<N$

) by adjusting the production rates. When

$n^{\prime }=N$

, the worst‐case degradation level is still lower than the failure threshold L upon replacement, so the failure penalty cost c _f is excluded when computing the operational cost. On the other hand, the production rates when

$n^{\prime }<N$

are generally higher, which may be able to compensate for the failure penalty cost. Therefore, we need to solve

$N-n$

LPs in total at time

$n\delta$

. For the LP that achieves the minimal production cost, the corresponding optimal solution

$u_{n+1}^{\ast }$

is adopted for production during the next period. Theorem 1 can also be used to solve Problem (4a)–(4h) by setting

$n=0$

to avoid the big‐M formulation. The computational burden to solve these

$N-n$

LPs can be remarkably reduced under Assumptions 1 and 2, with which we can derive a closed‐form expression for the optimal production rate as follows. Theorem 2

Consider period

$n\in [N]$

and the production–degradation relation (1) with

$\eta =1$

. Under Assumptions 1 and 2, if

$u_j^{\ast }=1$

for all

$j=n+1,\text{…},n^{\prime }$

is feasible to Problem (10) for any fixed

$n^{\prime }$

, then it is the optimal solution. Otherwise, the optimal solution

$u_{n+1}^{\ast }$

to (10) has the closed‐form expressions:

• If

  ${\mathrm{\Gamma }}_{n^{\prime }}^{(n)}\le n^{\prime }-n$
  , and either
  $\stackrel{\sim }{b}<\stackrel{\sim }{a}$
  or
  $$x_n+\delta \left[(n^{\prime }-n)\left(\mathtt{b}\frac{\stackrel{\sim }{a}}{\stackrel{\sim }{b}}+\mathtt{a}\right)+{\mathrm{\Gamma }}_{n^{\prime }}^{(n)}\stackrel{\sim }{a}\right]>L,$$
  11then the optimal production rate of the next period is
  $u_{n+1}^{\ast }=\frac{L-x_n-\delta [(n^{\prime }-n)\mathtt{a}+{\mathrm{\Gamma }}_{n^{\prime }}^{(n)}\stackrel{\sim }{a}]}{\delta (n^{\prime }-n)\mathtt{b}}$
  .

• If

  ${\mathrm{\Gamma }}_{n^{\prime }}^{(n)}\le n^{\prime }-n$
  ,
  $\stackrel{\sim }{b}\ge \stackrel{\sim }{a}$
  , and (11) is violated, then the optimal production rate of the next period is
  $u_{n+1}^{\ast }=\frac{L-x_n-(n^{\prime }-n)\delta \mathtt{a}}{\delta [(n^{\prime }-n)\mathtt{b}+{\mathrm{\Gamma }}_{n^{\prime }}^{(n)}\stackrel{\sim }{b}]}$
  .

• If

  $n^{\prime }-n<{\mathrm{\Gamma }}_{n^{\prime }}^{(n)}\le 2(n^{\prime }-n)$
  , and either
  $\stackrel{\sim }{b}<\stackrel{\sim }{a}$
  or (11) holds, then the optimal production rate of the next period is
  $u_{n+1}^{\ast }=\frac{L-x_n-(n^{\prime }-n)\delta (\mathtt{a}+\stackrel{\sim }{a})}{\delta [(n^{\prime }-n)\mathtt{b}+({\mathrm{\Gamma }}_{n^{\prime }}^{(n)}-n^{\prime }+n)\stackrel{\sim }{b}]}$
  .

• If

  $n^{\prime }-n<{\mathrm{\Gamma }}_{n^{\prime }}^{(n)}\le 2(n^{\prime }-n)$
  ,
  $\stackrel{\sim }{b}\ge \stackrel{\sim }{a}$
  , and (11) is violated, then the optimal production rate of the next period is
  $u_{n+1}^{\ast }=\frac{L-x_n-\delta [(n^{\prime }-n)\mathtt{a}+({\mathrm{\Gamma }}_{n^{\prime }}^{(n)}-n^{\prime }+n)\stackrel{\sim }{a}]}{\delta (n^{\prime }-n)(\mathtt{b}+\stackrel{\sim }{b})}$
  .

The key to prove Theorem 2 is to construct a pair of primal‐dual feasible solutions to Problem (10) and show their optimality by strong duality. This construction is largely enabled by Assumption 1. If we can find similar stationary assumptions for other systems, our analysis may also be useful to derive closed‐form solutions for the corresponding control problem. When parameter uncertainty is neglected, that is,

$\stackrel{\sim }{a}=0$

,

$\stackrel{\sim }{b}=0$

, and

${\mathrm{\Gamma }}_{n^{\prime }}^{(n)}\equiv 0$

for all

$n^{\prime }=n+1,\text{…},N$

, the closed‐form expression of

$u_{n+1}^{\ast }$

in Theorem 2 degenerates to

$u_{n+1}^{\ast }=\frac{L-x_n-(n^{\prime }-n)\delta \mathtt{a}}{\delta (n^{\prime }-n)\mathtt{b}}$

. Compared with this solution, the additional terms involving

$\stackrel{\sim }{a}$

,

$\stackrel{\sim }{b}$

, and

${\mathrm{\Gamma }}_{n^{\prime }}^{(n)}$

in our solution correspond to the safety margin that commonly appears when characterizing the optimal solution to robust optimization problems (see Bertsimas and Thiele, 2006; Mak et al., 2015; Zhang et al., 2021, to name a few). With the closed‐form expression for any fixed

$n^{\prime }=n+1,\text{…},N$

, searching for the optimal

$n^{\prime }$

is simple. The closed‐form solutions are crucial for IoT‐enabled systems that entail real‐time production control. It is also worth mentioning that existing literature has successfully derived closed‐form optimal robust solutions for newsvendor problems (Bertsimas & Thiele, 2006; Mamani et al., 2017). The worst‐case scenario for newsvendor problems is generally independent of the decision variables, that is, the order quantities. In contrast, the worst case of our condition‐based production problem is decision‐dependent, making the analysis of our problem more complicated. Consequently, we need to consider four cases as shown in Theorem 2, and our analytical results shed light on how to deal with decision‐dependent uncertainty in a robust optimization setting.

Recall that our model reveals a relationship between our robust optimization problem and the heuristics proposed in uit het Broek et al. (2020). Specifically,

${(A_n,{\overline{u}}_n)}_{n\in [N]}$

plays a similar role to the tuning parameters α₁ and α₂ in uit het Broek et al. (2020). As a by‐product of Theorem 2, we can also derive closed‐form expressions of

${(A_n,{\overline{u}}_n)}_{n\in [N]}$

(see Appendix A.7 in the Supporting Information). This enables us to control the production process in real time. In contrast, the two tuning parameters α₁ and α₂ for the heuristics in uit het Broek et al. (2020) have to be determined through a brute‐force search. This again shows the computational efficiency of the proposed model.

NONLINEAR PRODUCTION–DEGRADATION RELATIONS

Robust counterpart

We consider a general η in (1). With this production–degradation relation, the static production planning problem is the same as Problem (3a)–(3e), except that Constraint (3b) is replaced with

$$M{g}_n\ge x_0+\delta \underset{({\mathit{y}}_n,{\mathit{z}}_n)\in {\mathcal{P}}_n}{\max }\sum _{j=1}^n(b_j{u}_j^{\eta }+a_j)-L,n\in [N].$$

12Recall that we assume a known value of η since the shape parameter η is usually determined by the degradation mechanism. When the value of η cannot be precisely known, we can assume that η lies in a box uncertainty set

$[{\eta }_{\min },{\eta }_{\max }]$

. Since

$0\le u_j\le 1$

, the right‐hand side of (12) is nonincreasing in η for any fixed

$a_j,b_j\ge 0$

. It then follows from the monotonicity in η that we only need to consider (12) for

$\eta ={\eta }_{\min }$

in the worst case. Parallel to Lemma 1, we can reformulate the robust counterpart in the nonlinear case as Problem (EC.11) in Corollary EC.1 in Appendix A.8 in the Supporting Information.

The variables

${(A_n)}_{n\in [N]}$

defined in Section 3 can still be used to compute the modified degradation level at each decision epoch. By contrast, the modified upper bound of production rate in the nonlinear case is expressed as

${\overline{u}}_n\triangleq \min \{1,{\min }_{i=n,\text{…},N}\{{[(r_{in}^{\ast }+q_i^{\ast })/{\stackrel{\sim }{b}}_n]}^{1/\eta }\}\}$

for

$n\in [N]$

. However, it is challenging to directly solve the robust counterpart with existing solvers. The optimization problem with given

${(g_n)}_{n\in [N]}$

is even nonconvex when

$0<\eta <1$

. Nevertheless, under Assumptions 1 and 2, it is possible to obtain an optimal or approximately optimal solution for production planning based on some structural properties. To this end, we discuss the optimal control policy in the dynamic environment under the convex production–degradation relation (

$\eta >1$

) and the concave production–degradation relation (

$0<\eta <1$

), respectively. Lemma 2 can be readily generalized to nonlinear production–degradation relations, as shown in Corollary EC.2 in Appendix A.9 in the Supporting Information.

Based on Corollary EC.2, we follow the same procedure as in Sections 3 and 4 to derive the optimal solution for the dynamic production planning problem under Assumptions 1 and 2. We first solve Problem (EC.11) at time 0 for any fixed N and search for the optimal

$N^{\ast }$

. We then solve reoptimization problems in a rolling horizon. The reoptimization problem is the same as Problem (9) except that we replace all

$u_j$

in Constraints (9b) and (9c) with

$u_j^{\eta }$

. We derive the optimal solution to the production planning problem in Section 5.2 when

$\eta >1$

and propose an algorithm to approximately solve the optimization problem in Section 5.3 when

$0<\eta <1$

.

Convex production–degradation relation

Directly solving the robust counterpart is time‐consuming due to the integer decision variables

${(g_j)}_{j\in [N]}$

. Similar to the linear case, the robust counterpart becomes a convex program if

$\eta >1$

and the values of

${(g_j)}_{j\in [N]}$

are fixed. Parallel to Theorem 1, we can decompose this integer program into N convex programs (see Problem (EC.13) in Lemma EC.1 in Appendix A.10 in the Supporting Information).

Lemma EC.1 also reveals a structural property of the optimal solution to each convex program corresponding to the robust counterpart: for the optimal production plan, if the system does not produce at its maximum or minimum rate during some periods, then these production rates are equal. A similar result is obtained in uit het Broek et al. (2020, Lemma 2) when the degradation is deterministic. Under parameter uncertainty, our proof to Lemma EC.1 is significantly more complex due to the introduction of auxiliary variables

$q_j$

,

$r_{ji}$

, and

$s_{ji}$

. This structural property in Lemma EC.1 enables us to solve the optimization problem without resorting to solvers, and such results are summarized in Theorem EC.1 in Appendix A.11 in the Supporting Information.

Based on Theorem EC.1, one can readily solve the static planning problem in (EC.11) for every fixed N. Then we can follow the same procedure in Section 3.2 to search the optimal

$N^{\ast }$

for periodic maintenance. By replacing

$u_n^{\ast }\le \stackrel{\sim }{a}/\stackrel{\sim }{b}$

in condition (ii) of Proposition 2 with

$u_n^{\ast }\le {(\stackrel{\sim }{a}/\stackrel{\sim }{b})}^{\frac{1}{\eta }}$

, the proposition can still be used to terminate the search for convex production–degradation relations.

Next, we fix

$N=N^{\ast }$

and focus on the production control with in situ degradation signals. As illustrated in Section 5.1, the reoptimization problem for control in a rolling horizon can be formulated by replacing all

$u_j$

in Constraints (9b) and (9c) with

$u_j^{\eta }$

. The resulting mixed‐integer program is essentially the same as the static planning problem (EC.11). Indeed, we observe the in situ degradation level

$x_n$

and update the budget of uncertainty to

$\{{\mathrm{\Gamma }}_j^{(n)}:j>n\}$

at time

$n\delta$

. If we replace x ₀ and

${\mathrm{\Gamma }}_j$

in Problem (EC.11) with

$x_n$

and

${\mathrm{\Gamma }}_j^{(n)}$

, respectively, and keep decision variables that are only relevant to periods

$n+1,\text{…},N$

, then we can transform the static optimization problem at time 0 to the dynamic reoptimization problem at time

$n\delta$

. Based on this observation, we only need to slightly modify the closed‐form solution derived in Theorem EC.1 for control in a dynamic environment. The detailed closed‐form expressions are relegated to Appendix A.11 in the Supporting Information.

Concave production–degradation relation

Similar to the previous cases, we first use (5) to eliminate the integer decision variables

${(g_j)}_{j\in [N]}$

in (EC.11) by fixing a value

$n^{\prime }\in [N]$

so that

$g_j=0$

for

$j=1,\text{…},n^{\prime }$

and

$g_j=1$

for

$j=n^{\prime }+1,\text{…},N$

. Parallel to Theorems 1 and EC.1, we can enumerate

$n^{\prime }\in [N]$

and formulate an optimization problem for each fixed

$n^{\prime }$

, so that the optimal objective value of these N optimization problems coincides with that of Problem (EC.11). For any fixed

$n^{\prime }\in [N]$

, the corresponding optimization problem has the same form as (EC.13). Since

$0<\eta <1$

, however, this optimization problem is nonconvex. Due to the auxiliary decision variables

$q_j$

,

$s_{ji}$

, and

$r_{ji}$

that factor in uncertainty, it is challenging to solve the problem to optimality. If there is no parameter uncertainty for the degradation model (

${\mathrm{\Gamma }}_j\equiv 0$

for all j), then all of these auxiliary variables vanish. In this case, the analysis in uit het Broek et al. (2020) can be used to exactly obtain the optimal solution to Problem (EC.13) when

$\eta \in (0,1)$

. Based on this observation, we develop an algorithm that approximately solves the production planning problem for the concave production–degradation relation.

Specifically, we confine to

$\eta \in (0,1)$

and a fixed

$n^{\prime }\in [N]$

. When

${\mathrm{\Gamma }}_j\equiv 0$

for all j (i.e.,

$a_j=\mathtt{a}$

and

$b_j=\mathtt{b}$

for all j), Lemmas 2 and 5 in uit het Broek et al. (2020) indicate that there exists

$m\in [n^{\prime }]$

such that the optimal solution to Problem (EC.13) satisfies

13a

13b

13cHere, (13a) and (13b) imply that the optimal solution produces at the maximum rate during the first

$m-1$

periods, then produces at an intermediate rate

$u_m$

during period m, and finally shuts down the production for the remaining periods. Moreover,

$x_0+\delta [\mathtt{b}(m-1+u_m^{\eta })+n^{\prime }\mathtt{a}]$

in (13c) can be understood as the nominal degradation level at time

$n^{\prime }\delta$

if the degradation evolves deterministically according to (7). By carefully setting m and

$u_m$

, we can make this degradation level exactly equal to the failure threshold L.

Motivated by the structural properties of the optimal solution when

${\mathrm{\Gamma }}_j\equiv 0$

for all j, we propose approximately solving Problem (EC.13) when

$\eta \in (0,1)$

by adding Constraints (13a) and (13b) to cut down its feasible region. Moreover, we modify Constraint (13c) to

$$x_0+\delta \underset{({\mathit{y}}_{n^{\prime }},{\mathit{z}}_{n^{\prime }})\in {\mathcal{P}}_{n^{\prime }}}{\max }\left[\sum _{j=1}^{m-1}{b}_j+b_m{u}_m^{\eta }+\sum _{j=1}^{n^{\prime }}{a}_j\right]=L$$

14and also add this constraint to Problem (EC.13). Constraint (14) considers uncertainty in the parameters

$(a_j,b_j)$

through

$(y_j,z_j)$

and the uncertainty set in (2). By further optimizing the integer m, we arrive at the following optimization problem:

15The optimal value of the above program is an upper bound of

$C_{n^{\prime }}$

for Problem (EC.13). To minimize this upper bound, we need to deal with the maximization in (14) as well as the one‐dimensional search for the optimal m. Lemma 3 shows that the maximization in (14) can be solved analytically for any fixed value of m. Proposition 3 further facilitates the search for m. Lemma 3

Under Assumption 1, if

${(u_j)}_{j\in [n^{\prime }]}$

satisfies (13a) and (13b), then Constraint (14) can be written in closed form as follows:

• When

  $\stackrel{\sim }{a}\ge \stackrel{\sim }{b}$
  , Constraint (14) becomes
  
  16

• When

  $\stackrel{\sim }{a}<\stackrel{\sim }{b}$
  , Constraint (14) becomes
  
  17

Proposition 3

When

$\eta \in (0,1)$

, the following statements hold if there is at least one feasible solution

$(m,u_m)$

to (14).

• (i)

  When

  $\stackrel{\sim }{a}\ge \stackrel{\sim }{b}$
  , there is only one feasible solution
  $(m,u_m)$
  that satisfies (16). If, furthermore, we have
  ${\mathrm{\Gamma }}_{n^{\prime }}\le n^{\prime }$
  , then the optimal objective values of Problems (EC.13) and (15) coincide.

• (ii)

  When

  $\stackrel{\sim }{a}<\stackrel{\sim }{b}$
  , there are at most two feasible solutions that satisfy (17). If there are two feasible solutions
  $(m_1,u_{m_1})$
  and
  $(m_2,u_{m_2})$
  , we have
  $m_2=m_1+1$
  .

Statement (i) in Proposition 3 suggests that adding (13a), (13b), and (14) to Problem (EC.13) does not affect its optimal solution when

$\stackrel{\sim }{a}\ge \stackrel{\sim }{b}$

and

${\mathrm{\Gamma }}_N\le N$

, which may justify the proposed solution approach. For statement (ii), if there are two feasible solutions of

$(m,u_m)$

to (17) when

$\stackrel{\sim }{a}<\stackrel{\sim }{b}$

, it is easy to see that the one with a larger m yields a higher production revenue and hence is adopted as the approximately optimal solution to Problem (EC.13). Based on the above observation, we can then use a binary search that requires only

$\mathcal{O}(\log (N))$

steps to obtain the optimal m. The details are deferred to Algorithm EC.1 in Appendix B.1 in the Supporting Information. Comprehensive numerical experiments are conducted in Section 6.3 that empirically show the good performance of our algorithm.

We next optimize

$n^{\prime }\in [N]$

. For the linear and convex production–degradation relations, it entails enumerating all

$n^{\prime }\in [N]$

and hence solve N optimization problems to find the optimal

$n^{\prime }$

(see Theorems 1 and EC.1). By contrast, we only need to compare two values of

$n^{\prime }$

for Problem (15) to approximately solve Problem (EC.11). Theorem 3

Under Assumption 1, let

${\overline{C}}_{n^{\prime }}$

be the optimal objective value of Problem (15) when

$0<\eta <1$

. Then

$\arg {\min }_{n^{\prime }\in [N]}{\overline{C}}_{n^{\prime }}$

equals either

$n^{†}$

or N, where

$n^{†}\triangleq \min \{n^{\prime }\in [N]:x_0+\delta [n^{\prime }(\mathtt{a}+\mathtt{b})+\min \{{\mathrm{\Gamma }}_{n^{\prime }},n^{\prime }\}\max \{\stackrel{\sim }{a},\stackrel{\sim }{b}\}+{({\mathrm{\Gamma }}_{n^{\prime }}-n^{\prime })}_{+}\min \{\stackrel{\sim }{a},\stackrel{\sim }{b}\}]>L\}$

.

Based on Theorem 3, we can readily compute

$\arg {\min }_{n^{\prime }\in [N]}{\overline{C}}_{n^{\prime }}$

to obtain an approximately optimal solution to Problem (EC.11) for any fixed N. The remaining procedure is the same as that for the convex production–degradation relation in Section 5.2.

In practice, we can further input our approximately optimal solution as an initial solution of an interior point algorithm for a warm start, which can possibly find a better production plan with the cost of a longer computational time. Nevertheless, our numerical experiments in Section 6.3 show that the improvement brought by this strategy is insignificant. This shows that the approximately optimal production plan obtained by our algorithm has performed sufficiently well.

APPLICATIONS TO EXTRUDER PRODUCTION

We apply the proposed model to control the production rate of an extruder system, and a graphical illustration of this system is displayed in Figure 1. An extruder produces bottle caps by using a rotating screw to move the melted plastic raw material through a cylinder. Due to wear of the screw, the melted material may leak within the cylinder, and a thermal resistor is then activated for cooling. Hence, we use the mean daily usage frequency of the thermal resistor as a surrogate degradation signal for the screw wear. Section 6.1 analyzes a real‐world extruder degradation dataset, based on which we determine our simulation setting. Section 6.2 uses simulation to compare the performance of our robust model with several existing control methods for the convex production–degradation relation. Section 6.3 further assesses the performance of our model for the concave production–degradation relation. All optimization algorithms are implemented in Python, and all computational times are obtained from a computer with an Intel Core i7 CPU 2.2 GHz processor.

FIGURE 1

A graphical illustration of the extruder production system in this study. The screw and cylinder are subject to the degradation‐induced failure.

Data analysis

We collect degradation data of an extruder in an Italian factory from the commission date, that is, July 2017 to July 2019. During the study period, the rotational speed

$u_{\mathrm{o},i}$

(normalized from 0 to 1) of the screw and the degradation indicator

$x_{\mathrm{o},i}$

(usage of the thermal resistor) are recorded every month, so there are in total

$n=24$

observations. The initial degradation level is

$x_0=0$

. Denote the observed data as

$\mathbb{D}=\{(\mathrm{\Delta }{x}_{\mathrm{o},i},u_{\mathrm{o},i}),i=1,\text{…},n\}$

, where

$\mathrm{\Delta }{x}_{\mathrm{o},i}=x_{\mathrm{o},i}-x_{\mathrm{o},i-1}$

for

$i\in [n]$

with

$x_{\mathrm{o},0}\triangleq 0$

.

We adopt two popular degradation models to fit the data

$\mathbb{D}$

, that is, the stationary Wiener and gamma processes (Hao et al., 2015; Chen & Ye, 2018). The dynamic covariate u (i.e., the rotational speed) is incorporated through an accelerated failure time model with link function

$h(u)=a+b{u}^{\eta }$

, a commonly used model for handling covariates in reliability. The resulting models are shown in Table 1. By assuming a constant production rate in each period, the log‐likelihoods (up to a constant) for the Wiener and gamma processes (with the shape‐scale parameterization) can be specified as

where

$\mathrm{\Gamma }(·)$

denotes the gamma function,

$h_i\triangleq b{u}_{\mathrm{o},i}^{\eta }+a$

, and

$\mathrm{\Delta }t=4.5$

weeks. Through numerical maximization, the maximum likelihood estimates (MLEs) for parameters in the Wiener process are

$(\widehat{a},\widehat{b},\widehat{\eta },\widehat{\sigma })=(0.65,1.34,2.38,2.64)$

, while those for the gamma process are

$(\widehat{a},\widehat{b},\widehat{\eta },\widehat{\theta })=(0.55,1.45,2.41,7.49)$

. Under some mild conditions, the estimator converges weakly to a zero‐mean Gaussian distribution as detailed in Appendix A.15 in the Supporting Information.

TABLE 1

Distributions of the degradation increment at a period under different stochastic process models.

Model	Distribution	Mean	Variance
Wiener process	$\mathcal{N}(h(u)\delta ,{\sigma }^{2}h(u)\delta )$	$h(u)\delta$	${\sigma }^{2}h(u)\delta$
Gamma process	$\mathsf{Gamma}(h(u)\delta /\theta ,\theta )$	$h(u)\delta$	$\theta h(u)\delta$

Note: The distributions are displayed under a fixed production rate u. Here,

$\mathcal{N}(\mu ,{\sigma }^2)$

and

$\mathsf{Gamma}(k,\theta )$

denote the normal distribution with mean μ and variance σ² and gamma distribution with shape k and scale θ, respectively.

We can use the asymptotic normality results in Appendix A.15 to determine the parameter values in the proposed robust model. We first identify a “working degradation model,” for example, either a Wiener or gamma process, and obtain MLEs of its parameters. Here, we demonstrate using the Wiener process as the working model. The procedure is almost the same for the gamma process. In the robust model, we set

$(\mathtt{a},\mathtt{b})=(\widehat{a},\widehat{b})$

, and

$\stackrel{\sim }{a}$

and

$\stackrel{\sim }{b}$

are set to be 1.96 times of the asymptotic standard deviation of

$\widehat{a}$

and

$\widehat{b}$

, respectively. With the inferential results for the Wiener process above, we have

$(\mathtt{a},\mathtt{b},\stackrel{\sim }{a},\stackrel{\sim }{b})=(0.65,1.34,0.52,1.34)$

. To determine the budgets of uncertainties

${\mathrm{\Gamma }}_n$

, note that if we produce at a nominal rate

$\mathtt{u}$

throughout

$[0,n\delta ]$

, the estimated standard deviation of the Wiener degradation process at time

$n\delta$

is

$\widehat{\sigma }\sqrt{\widehat{h}(\mathtt{u})n\delta }$

. While by Theorem EC.1, the worst‐case deviation of

$x_n$

from its nominal value is

${\mathrm{\Gamma }}_n\delta \max \{\stackrel{\sim }{b}{\mathtt{u}}^{\eta },\stackrel{\sim }{a}\}$

when

${\mathrm{\Gamma }}_n\le n$

. Similar to Bertsimas and Thiele (2006), we match this worst‐case deviation with 1.96 times the estimated standard deviation of the degradation level of the working model, which gives

$${\mathrm{\Gamma }}_n=\frac{1.96\widehat{\sigma }}{\max \{\stackrel{\sim }{b}{\mathtt{u}}^{\eta },\stackrel{\sim }{a}\}}\sqrt{\frac{\widehat{h}(\mathtt{u})n}{\delta }}.$$

18With

$\mathtt{u}$

equal to the average of the historical production rates, the above display equals

${\mathrm{\Gamma }}_n=\phi \sqrt{n}$

with

$\phi =9.66$

for our extruder example. In addition, the setting of our robust model ensures that

${\mathrm{\Gamma }}_n\le 2n$

. Therefore, we let

${\mathrm{\Gamma }}_n=\min \{9.66\sqrt{n},2n\}$

for

$n\in {\mathbb{N}}_{+}$

. Sensitivity analysis in Section 6.2 shows that our optimization result is insensitive to the choice of φ.

Experiments on convex production–degradation relations

We first assess the performance of the proposed robust model under convex production–degradation relations. In the simulation, we set the workload adjustment interval as

$\delta =1$

and the failure threshold as

$L=500$

. The replacement cost is set to be

$c_{\mathrm{r}}=0.2$

, and a system failure incurs an additional penalty cost

$c_{\mathrm{f}}=19.8$

. When the system produces at the maximum rate, the production revenue is

$\pi =0.03$

per time unit. We then compare the long‐run production costs per time unit obtained by the following control methods:

• the proposed robust model with parameters determined as in Section 6.1 (referenced as “RO”) using the Wiener process as the working model;

• the MDP using the correctly specified degradation model (referenced as “MDP”);

• the MDP using a misspecified degradation model (referenced as “MDP‐MIS”);

• the deterministic degradation process (referenced as “Deterministic”).

Using the above setting, we consider two scenarios based on the data‐fitting results above. In the first scenario, we use the stationary Wiener process with parameter

$(a,b,\eta ,\sigma )=(0.65,1.34,2.38,2.64)$

as the true model that governs the degradation of the extruder, and the stationary gamma process with parameters

$(a,b,\eta ,\theta )=(0.55,1.45,2.41,7.49)$

as the misspecified model. In the second scenario, we use the stationary gamma process with

$(a,b,\eta ,\theta )=(0.55,1.45,2.41,7.49)$

as the true model, and the stationary Wiener process with

$(a,b,\eta ,\sigma )=(0.65,1.34,2.38,2.64)$

as the misspecified model. In both scenarios, the deterministic degradation process assumes that the degradation rate is

$h(u)=\mathtt{a}+\mathtt{b}{u}^{\eta }$

under production rate u, with

$(\mathtt{a},\mathtt{b},\eta )=(0.65,1.34,2.38)$

. The number of periods N for replacement is separately optimized for each control method.

For each of the two settings above, we replicate the simulation 1000 times. Table 2 summarizes the computational time to solve the static planning problem at time 0 for all control methods by averaging over the 1000 replications. As with Bertsimas and Thiele (2006) and Mamani et al. (2017), we then assess the performance of each control method by examining the mean, upper semivariance (USV), and p‐quantiles (

$p\in \{0.75,0.85,0.95,0.99\}$

) of the production cost rate, where

$USV(X)\triangleq E{(X-E[X])}_{+}^2$

measures the variability of production costs above the mean. Unsurprisingly, when we are able to correctly specify the degradation model, the MDP yields the lowest mean cost rate as shown in Table 3. We also observe that the robust model outperforms the deterministic algorithm and the MDP using a misspecified model, and it achieves a significantly lower USV than other control methods (even when the working degradation model to determine the uncertainty sets is misspecified). This is because a system can hardly fail under the robust control method. As such, the variation in the resulting production cost is quite small. When the actual degradation follows a gamma process, the MDP using a misspecified Wiener process has the risk of high worst‐case production cost (see the 0.99 quantile) and large USV. In summary, the proposed model enjoys both computational efficiency and robustness under model misspecification compared with all benchmarks.

TABLE 2

Computational times for optimization of the replacement time and static production planning.

	RO‐CF	RO‐Gurobi	MDP	Deterministic
Time (in min)	≈0	—	104.98	3.67

Note: “RO‐CF” and “RO‐Gurobi” refer to using the closed‐form expressions and Gurobi for optimization, respectively. The symbol “—” means that we cannot find the optimal solution within a 300‐min time limit for the joint optimization problem under our mixed‐integer second‐order cone program (MISOCP) formulation. (See Appendix B.2 for implementation details of all control methods and Appendix B.3 for the detailed MISOCP formulation in the Supporting Information.)

TABLE 3

Relative performances between robust and benchmark solutions under different degradation models.

		MDP	MDP‐MIS	Deterministic
	Mean	−1.15	0.048	77.10
	USV	1.10 × 104	38.80	6.82 × 105
True model:	0.99 quantile	4.95	0.61	119.49
Wiener process	0.95 quantile	3.93	0.41	95.75
	0.85 quantile	−1.04	0.24	95.75
	0.75 quantile	−1.29	0.15	89.62
	Mean	−0.69	4.95	66.75
	USV	159.28	3.93 × 103	3.78 × 104
True model:	0.99 quantile	−0.58	17.76	109.91
gamma process	0.95 quantile	−0.12	16.58	84.98
	0.85 quantile	−0.45	11.14	84.98
	0.75 quantile	−0.59	6.15	79.33

Note: All performance measures are in terms of the long‐run production cost per time unit and empirically computed from 1000 simulation replications. For any performance measure

$\mathsf{PM}$

, the relative difference is computed as

$({\mathsf{PM}}_{\text{Benchmark}}-{\mathsf{PM}}_{\text{Robust}})/|{\mathsf{PM}}_{\text{Robust}}|\times 100\%$

(the subscript indicates the model) and displayed in percent. A positive value indicates that the robust solution performs better than the corresponding benchmark solution with regard to the particular measure.

Abbreviations: Deterministic, the deterministic degradation process; MDP, the MDP using the correctly specified degradation model; MDP‐MIS, the MDP using a misspecified degradation model; USV, upper semivariance.

The above simulation uses

${\mathrm{\Gamma }}_n=\min \{\phi \sqrt{n},2n\}$

with

$\phi =9.66$

for the robust model based on the Wiener working model. As can be seen from (18), φ is a function of the estimated parameters. We further use the delta method to construct its 90% confidence interval as [7.24, 12.08], vary φ within this range, and run simulations again to see the impacts of varying

${\mathrm{\Gamma }}_n$

on the optimization results. Figure 2 shows the empirical mean and 0.99 quantile of the production cost rates of our robust model as a function of φ when the actual degradation follows the gamma process. For all values of φ, our robust solutions perform better than the MDP solutions with model misspecification. Generally, the optimization result of the proposed robust model is not very sensitive to φ.

FIGURE 2

Sensitivity analysis on φ when the actual degradation is the gamma process. Dotted and dash horizontal lines indicate the performance measures obtained from the MDP solution with and without model misspecification, respectively. Red and black separately indicate the empirical mean and 0.99 quantile.

Historical data to determine

$\stackrel{\sim }{a}$

and

$\stackrel{\sim }{b}$

may not be available in some cases. Hence, we also assess the performance of the robust model when

$\stackrel{\sim }{a}$

and

$\stackrel{\sim }{b}$

take different values. Due to the space limit, the detailed setting and findings are provided in Appendix C in the Supporting Information.

Experiments on concave production–degradation relations

Section 5.3 proposes an algorithm to approximately solve Problem (EC.13) under concave production–degradation relations. We compare the performance of the approximately optimal robust production plan with the MDPs with and without model misspecification. For simplicity, we assume that the actual degradation follows a gamma process; the parameter values used in Section 6.2 are adopted here, except

$\eta =0.8$

for all models. Figure 3 shows the comparison result from 1000 simulation replications. Similar to Section 6.2, the MDP without misspecification yields the lowest mean cost, but our robust solution is comparable to the optimal MDP solution. On the other hand, the robust model outperforms the MDP with model misspecification. Moreover, the variations of the optimal cost rates obtained by the robust model are much smaller than the MDP solutions. These observations tally with the results in Section 6.2.

FIGURE 3

The box plot of the production cost rates obtained from 1000 simulation replications under our robust model and the MDPs with and without model misspecification. MDP, the MDP using the correctly specified degradation model; MDP‐MIS, the MDP using a misspecified degradation model.

Finally, we compare the proposed algorithm in Section 5.3 with the state‐of‐the‐art interior point algorithm. Based on the simulation setting above, we generate a problem instance through randomly sampling the parameter values for the robust model:

$\mathtt{a}\sim \mathsf{U}(0.5,0.8)$

,

$\mathtt{b}\sim \mathsf{U}(1.2,1.8)$

,

$\stackrel{\sim }{a}\sim \mathsf{U}(0.8\mathtt{a},\mathtt{a})$

,

$\stackrel{\sim }{b}\sim \mathsf{U}(0.8\mathtt{b},\mathtt{b})$

, and then

$\eta \sim \mathsf{U}(0.7,0.9)$

, where

$\mathsf{U}(c,d)$

is the uniform distribution over

$[c,d]$

. We repeat the above random sampling to generate 500 instances. For each instance, the values of L, δ, and

${\mathrm{\Gamma }}_n$

are set to be the same as in the previous sections, and we fix

$N=300$

. We use the proposed algorithm and the Pyomo interior point IPOPT solver (version 3.14.4; released in September 2021) to solve each problem instance. With the default starting point of the IPOPT solver, we observe that the interior point algorithm usually converges to a local optimum that is much worse than our approximately optimal solution. Hence, we input our approximately optimal solution as an initial solution for the IPOPT solver as a warm start. Summary statistics of the relative differences are presented in Table 4. Out of the 500 problem instances, the IPOPT solver with a warm start only finds a better solution with a lower cost than 89 problem instances. This shows the satisfactory performance of our proposed algorithm.

TABLE 4

Summary statistics of relative differences (in %) between our algorithm and the interior point algorithm.

	Min	First quartile	Median	Third quartile	Max	Mean
Default starting point	−29.54	−17.04	−10.74	−2.97	34.25	−9.01
Warm start	0	0	0	0	57.8	1.94

Note: The relative difference of every pair of solutions is computed as

$(C_{\text{prop}}-C_{\text{solver}})/|C_{\text{solver}}|\times 100\%$

, where

$C_{\text{solver}}$

and

$C_{\text{prop}}$

are the optimal values of Problem (EC.13) computed by the solver and the proposed algorithm, respectively. A negative value indicates that the proposed model finds a better solution with a lower cost than the IPOPT solver.

Sensitivity analysis

This section further examines the sensitivity of the optimization result with respect to changes in model parameters. We set the gamma process as the actual degradation model as in Section 6.3. For each parameter, we vary the value from 25% to twice of its base value in Sections 6.1 and 6.2 and keep other parameters unchanged. We determine

$(\mathtt{a},\mathtt{b},\stackrel{\sim }{a},\stackrel{\sim }{b},{({\mathrm{\Gamma }}_n)}_{n\in {\mathbb{N}}_{+}})$

for our robust model as in Appendix B.2 in the Supporting Information. Then we rerun the simulation in Section 6.2 under each parameter setting and compare the results between our model and the benchmarks.

Figure 4 compares the long‐run revenue rates (obtained by 1000 simulation replications) under our model and the benchmarks when the cost parameters π, c _r, and c _f change, respectively. We can see that our robust solution is comparable to the optimal MDP solution with correct model specification under all parameter settings, while the MDP under model misspecification yields a significantly lower revenue in some cases. The long‐run average revenue rate as a function of π, c _r, and c _f, respectively, exhibits different patterns. The costs c _r and c _f are not very significant in affecting the gaps between the robust model and the optimal MDP policy. However, the MDP with model misspecification is likely to significantly increase the operational cost rate when these costs are large. If we zoom in Figure 4c (see Table EC.1 for the detailed values in the Supporting Information), the gap between our method and the benchmarks is sensitive to π, and magnitudes of the gaps decease in π. This is because when the production profit decreases, an unexpected failure would greatly affect the long‐run revenue rate. Since our robust model generally has fewer failures in simulations, the gap between our model and the MDP with model misspecification is large when

$\pi =0.0075$

.

FIGURE 4

The long‐run production revenues of the robust model and the MDPs with and without model misspecification (obtained from 1000 simulation replications) under different cost parameters. MDP, the MDP using the correctly specified degradation model; MDP‐MIS, the MDP using a misspecified degradation model.

Figure 5 compares the long‐run revenue rates (obtained by 1000 simulation replications) under the robust model and benchmarks when the degradation parameters a, b, θ, and L change, respectively. Note that the parameters

$(\mathtt{a},\mathtt{b},\stackrel{\sim }{a},\stackrel{\sim }{b},{({\mathrm{\Gamma }}_n)}_{n\in {\mathbb{N}}_{+}})$

used for the robust model change with a, b, and θ. The robust model performs significantly better than the MDP with model misspecification when a is not too large. Under different values of b, the gap between the robust model and the optimal MDP policy is small, and our model performs significantly better than the MDP with model misspecification. The production revenue rates of all models are less sensitive to the degradation volatility θ in the sense that all the revenue rates do not vary significantly with respect to θ. When the failure threshold L is small, it is easier to observe a system failure under a nonoptimal production decision. Consequently, both the robust model and MDP with model misspecification yield a significantly lower revenue rate than the optimal MDP solution. In contrast, the robust model performs comparably to the optimal MDP policy when L is large, while the MDP with model misspecification still suffers from a low revenue rate. Generally, the performances of the robust model and the optimal MDP are close when the system failure risk during a period is relatively low. This happens when a, b, θ are large or L is small. This finding supports the application of the IoT technology in practice for real‐time control, where the duration of a period is relatively short.

FIGURE 5

The long‐run production revenues of the robust model and the MDPs with and without model misspecification (obtained from 1000 simulation replications) under different degradation parameters. MDP, the MDP using the correctly specified degradation model; MDP‐MIS, the MDP using a misspecified degradation model.

CONCLUSION

This study develops a robust optimization framework for the optimal production and maintenance control problem for both linear and nonlinear production–degradation relations. The framework circumvents the difficulty of predetermining a stochastic degradation model for optimization. As seen in Section 6, the parameters for the proposed model can be easily determined based on historical degradation data. Based on the derived structural properties, we obtain closed‐form solutions to the optimal (in the case of linear or convex production–degradation relations) or approximately optimal (in the case of concave production–degradation relations) dynamic production rates. Our approximate algorithm for concave production–degradation relations showcases how to use structural results under the deterministic degradation setting to solve nonconvex optimization problems. Both the closed‐form solution and the approximate algorithm ensure real‐time implementation of the proposed model. This appealing feature could enable us to adopt a small enough interinspection interval δ for real‐time control of an IoT‐enabled production system. Moreover, our model lends theoretical support to the heuristics proposed in uit het Broek et al. (2020) from robust optimization. Finally, the application of the extruder example has shown that the proposed model is robust against model misspecification, whereas the optimal cost rate produced from MDP can be 20% higher than the robust solution. The USV of the cost rates obtained from the robust model is also shown to be remarkably smaller than the MDP solutions.

Several topics could be further investigated in the future. First, we determine all model parameters before optimization. In the future, we would investigate the possibility of integrating robust optimization and learning of uncertain parameters (e.g., Lim et al., 2012) for our problem whenever a new degradation level is revealed at a decision epoch. Second, this study estimates the nominal parameters by a working degradation model as implemented in Section 6. We would further investigate the joint estimation and robustness optimization (Zhu et al., 2022) for the condition‐based production planning problem. Last, we could integrate stochastic (multiperiod) demands into the proposed model, so our production plan can be dynamically adjusted to the realized demands over time. The coordination among multiple production systems can also be considered in the presence of stochastic production demands.

Footnotes

ACKNOWLEDGMENTS

The authors would like to thank the editor, a senior editor, and two anonymous referees for their valuable comments and constructive suggestions, which have led to a significant improvement in the quality and presentation of this work. This work was supported in part by the National Science Foundation of China (72071071, 72101145, 72221001), the Singapore MOE AcRF Tier 2 Grant (A‐8001052‐00‐00), and the Shanghai Pujiang Program (21PJC071).

ORCID

Qiuzhuang Sun