The stochastic opportunistic replacement problem, part III: improved bounding procedures

We consider the problem to find a schedule for component replacement in a multi-component system, whose components possess stochastic lives and economic dependencies, such that the expected costs for maintenance during a pre-defined time period are minimized. The problem was considered in Patriksson et al. (Ann Oper Res 224:51–75, 2015b), in which a two-stage approximation of the problem was optimized through decomposition (denoted the optimization policy). The current paper improves the effectiveness of the decomposition approach by establishing a tighter bound on the value of the recourse function (i.e., the second stage in the approximation). A general lower bound on the expected maintenance cost is also established. Numerical experiments with 100 simulation scenarios for each of four test instances show that the tighter bound yields a decomposition generating fewer optimality cuts. They also illustrate the quality of the lower bound. Contrary to results presented earlier, an age-based policy performs on par with the optimization policy, although most simple policies perform worse than the optimization policy.


Introduction
This article studies a generalization of a multi-component maintenance scheduling problem with deterministic component lives (the so-called opportunistic replacement problem, ORP;  Almgren et al. 2012), to the case when component lives are considered to be stochastic variables. The resulting stochastic opportunistic replacement problem (SORP) was introduced and analyzed in Patriksson et al. (2015a, b).
We consider a stochastic programming approach for the minimization of the expected cost of maintenance over the remaining planning horizon, assuming that the failure risk functions are non-decreasing with increasing component ages. Our problem formulation leads to a two-stage approximation of a multi-stage stochastic programming problem, which can be decomposed, akin to integer Benders' decomposition (Benders 1962). The decomposition requires lower bounds on the value of the recourse function, which represents the expected future maintenance cost of a specific replacement decision at the current time, given the current ages of the components. We present a means to compute such lower bounds; in particular, we improve upon the corresponding bounds from Patriksson et al. (2015b), resulting in a more efficient decomposition method.
Two identical systems with identical component ages were considered in Patriksson et al. (2015b, Prop. 2); it claims that the value of the recourse function of either system is upper bounded by that of the other, plus a non-negative value defined by the (current) replacement decisions for the two systems. Based on our improved bounds on the recourse function, we state and verify a tighter bound than the one claimed.
Beside the improved efficiency of the decomposition, we establish a lower bound on the expected cost of the solution to the SORP; it is useful since our method does not in general provide an optimal solution to the SORP.
A review over current research is given in Sect. 2. A formal definition of the problem studied is given in Sect. 3, and in Sect. 4 it is formalized how time is discretized. Section 5 contains a lower bound on the value of the SORP. Sections 6, 7, and 8 contain (the new) bounds on the value of the recourse function, a brief description of the ORPIL (a deterministic optimization problem), and a description of a decomposition method, respectively. Together, these three sections form a presentation of the optimization policy. Section 9 contains some numerical tests comparing the new optimization policy with the old one. Section 10 concludes the paper.

Literature survey
While the literature abounds with articles on aspects of reliability and maintenance of multicomponent systems with deteriorating and/or stochastically failing components, we provide below a summary of the most interesting and (mostly) recent related articles in the area, sorted along a time-line.
van Noortwijk and Frangopol (2004) describe the difference between proactive and reactive preventive maintenance (PM) actions, applied before and after, respectively, an indication of a deterioration. As the rationality behind these principles have been questioned, an optimal PM strategy based on lifetime reliability and life-cycle costs is established. The former model was applied by the Netherlands Ministry of Transport, Public Works and Water Management (Rijkswaterstaat); the latter contributed to the further development of a bridge management methodology set up by the UK Highways Agency. Deloux et al. (2009) develop a predictive maintenance policy for a continuously deteriorating system subject to stress. The system considered possesses two failure mechanisms: an excessive deterioration level and a shock. To optimize the maintenance policy, a combined statistical process control (SPC) and condition-based maintenance (CBM) approach is proposed. The CBM is used to inspect and replace the system according to the deterioration level. SPC is used to monitor the stress covariate. In order to assess the performance of the proposed policy and to minimize the long-run expected maintenance cost per time unit, a mathematical model for the system cost is derived. An analysis of the numerical results highlights the properties of the proposed policy w.r.t. the different maintenance parameters. Laggoune et al. (2009) propose a PM planning approach for a multi-component series system subjected to random failures, the cost rate being minimized under a general life distribution. The expected total cost comprises corrective and preventive costs (w.r.t. the components) as well as common costs (w.r.t. the production loss during system shutdown). When the system is down (correctively or preventively) the opportunity to preventively replace nonfailed components is considered. A solution procedure combining Monte-Carlo simulations and a heuristic search is proposed and applied to component replacement in the hydrogen compressor in an oil refinery. Wang (2012) recognizes that spare parts demands are difficult to predict based on historical data of parts usages, wherefore an optimal inventory control policy may be difficult to obtain. A joint optimisation of the inventory control of spare parts and the PM inspection interval is proposed. Stochastic cost models for spare parts inventory and maintenance are derived and enumeration is employed to find optimal solutions over a finite time horizon. A delaytime concept is developed for inspection modelling and used to derive the probabilities of numbers of failures and defective items identified at a PM epoch. The inventory model follows a periodic review policy, the demand being governed by the need for spare parts due to maintenance. Carlos et al. (2013) argue that the maintenance frequency of wind farms is partly determined by wind velocity variations, which partly determine the degradation progress. Wind measurements and Monte-Carlo sampling are used to estimate the power generated, with a cost associated to production loss during maintenance. The bi-objective goal of the wind farm maintenance strategy is to minimize the total cost and maximize the energy produced, maintenance frequency being a decision variable. Jha et al. (2013) develop a maintenance scheduling optimization model for bridge infrastructure-utilizing a prediction model to account for stochastic aspects of deterioration-as well as an experimental procedure for examining bridge deterioration. The importance of regarding structural deterioration is discussed, and the work intend to imply environmentally sustainable structures. Specifically, the use of composite materials to replace steel rebarswhich are due to corrosion-within bridge decks is presented as a sustainable alternative, to reduce the need for repair and maintenance.
Stochastic control methods have a long history in risk management and life-cycle cost procedures. Papakonstantinou and Shinozuka (2014) combine stochastic control methods and Bayesian principles into partially observable Markov decision processes (POMDPs) that expand available policy options compared to some state-of-the-art methods. POMDPs enable optimum decisions, based on the best possible knowledge of a decision-maker at each time. The problem of finding optimal policies for the maintenance and management of aging structures through a POMDP framework with large state spaces is modelled and solved; the framework is formed using stochastic, physically based models. An example of a corroded existing structure is presented; it is based on non-stationary POMDPs, for an infinite and a finite horizon case with 332 and 14,009 states, respectively. Comparea et al. (2015) propose and compare maintenance optimization techniques based on genetic algorithms (GA), the parameters of the maintenance model being affected by uncertainty and the fitness values represented by cumulative distributions. A method to rank the uncertain fitness values and a novel Pareto dominance concept are developed. The GA-based methods are applied to a practical case concerning the setting of a CBM policy on the degrading nozzles of a gas turbine operating in an energy production plant. Do et al. (2015) consider a proactive CBM technique, with perfect and/or imperfect maintenance actions, in a deteriorating system. Perfect maintenance actions restore the system to an'as good as new' state-often at a high cost. A positive impact of imperfect maintenance actions is the relatively low cost; negative impacts are that (i) the system is restored to a state between 'as good as new' and 'as bad as old', and (ii) each action may accelerate the system's deterioration speed. An adaptive maintenance policy is proposed which, at each inspection time needed, selects optimal maintenance actions, the time interval between successive inspection times being determined w.r.t. a remaining useful life (RUL); the use of the policy is illustrated by a numerical example. Gunn and Diallo (2015) consider systems operating in critical environments, where operators and/or regulators often specify-for major components-replacement intervals, within which failure cannot occur. A preventive replacement of a component yields an opportunity to replace other components that are within their replacement intervals-to avoid repeating the large teardown cost in short term. For a fixed time horizon, the opportunistic indirect grouping of periodic events problem is represented by a tree of possible combinations of replacements of major component. A depth-first shortest path algorithm is developed; for moderate numbers of components, out of millions of nodes created only a small portion need to be examined. For larger numbers of components and longer time horizons, the depth-first search still rapidly finds improving solutions and serves as a good heuristic approach. Aghezzaf et al. (2016) consider the integration of production and maintenance planning in a failure-prone manufacturing system. The system's operating state (i.e., age) is stochastically predictable and it can receive PM during preplanned periods. PM is imperfect: the system is restored to a state between 'as bad as old' and 'as good as new', the latter reachable only by an overhaul. An integrated production and PM planning model-accounting for the system's manufacture capacity and operational reliability state-is formulated using mixed-integer non-linear optimization, and then reformulated as an extended mixed-integer linear program. Since the reformulation is computationally demanding, a 'fix-and-optimize' procedure is proposed, that utilizes some properties of the original model. The resulting procedure is tested-with good results-on instances adapted from the literature. Atashgar and Abdollahzadeh (2016) formulate a joint redundancy and imperfect block opportunistic maintenance optimization model, to determine the redundancy level and a maintenance strategy that simultaneously minimizes the wind farm loss-of-load probability and life-cycle cost. Reliability thresholds are introduced for imperfect maintenance of failed and working turbines, and preventive dispatch of maintenance teams. A simulation method developed evaluates the performance measures of a wind farm system w.r.t. different types of wind turbines, maintenance activation delays and durations, and a limited number of maintenance teams. The influence on the performance of the wind farm, of the assumptions and the parameters of the simulation model, is discussed based on a sensitivity analysis. Pareto optimal solutions are derived using a multi-objective particle swarm optimization algorithm. A comparative study with a commonly used maintenance policy shows that the proposed strategy significantly reduces maintenance costs and loss-of-load probability. Cherkaoui et al. (2016) assess the economic performance of CBM strategies through their long-run expected maintenance cost rate criterion as well as robustness. A cost model is developed to quantify the economic performance and robustness to assess CBM strategies. Two representative strategies-periodic and quantile-based-for inspection and replacement are compared, revealing the factors affecting their performance and robustness the most. Rasmekomen and Parlikad (2016) describe an approach to optimize CBM of multicomponent systems, where the state of certain components may affect the rate of degradation of other components. A real case is presented: an industrial cold box in a petrochemical plant; data collected on the fouling of its tubes show that the extent of fouling of one tube affects the rate of fouling of other tubes due to overloading. Shi and Zeng (2016) describe an opportunistic PM strategy and analyze multi-component systems with stochastic dependence (similar to that by Rasmekomen and Parlikad 2016). Assuming that measures are available of the RULs of components as well as of the impact of a component's degradation on the RULs of other components, filtering theory is used for predictions. An optimization model determines timings and groupings, based on an opportunistic maintenance principle: at failure maintenance is performed also on other components if time lies in an "opportunity zone". An optimal balance is sought between the cost of performing maintenance too early on some components, and the gain from opportunistic coordination. Tests are limited to three-component systems. Alaswad and Xiang (2017) review CBM optimization models for stochastically deteriorating single-and multi-unit systems. A CBM strategy collects and assesses real-time information, and suggests maintenance decisions based on the system's current condition. In recent decades, research on CBM has grown rapidly due to a rapid development of computer-based monitoring technologies. Research studies have shown that CBM-if planned properly-can effectively improve equipment reliability at reduced costs. The review is emphasised on mathematical modeling and optimization approaches. Focus lies on optimization criteria, inspection frequency, maintenance degree, and solution methodology. Since the modeling choice for the stochastic deterioration process greatly influences CBM strategy decisions, the CBM models are classified based on the underlying deterioration processes, namely discrete-and continuous-state deterioration, and the proportional hazard model. Let V be the set of system states, where each state (s, ξ ξ ξ, a) ∈ V is composed by the time s ∈ R + , the vector ξ ξ ξ ∈ B N of component states [ξ n = 0 (ξ n = 1) if component n ∈ N is functioning (in a failed state)] at time s, and the vector a ∈ R N + of ages at time s (a n denotes the age of component n ∈ N ).

Problem definition
We denote the set of feasible maintenance decisions with respect to the state (s, ξ ξ ξ, a) ∈ V, by X ξ ξ ξ = { x ∈ B N | x ≥ ξ ξ ξ }. A decision policy may be used to make a maintenance decision for a system in a given state.
Definition 1 (Maintenance policies) A deterministic maintenance policy is a function p : V → B N such that p(s, ξ ξ ξ, a) ∈ X ξ ξ ξ . Let B N denote the set of random vectors with outcome space B N . A stochastic maintenance policy is a function p : V → B N such that Pr(p(s, ξ ξ ξ, a) ∈ X ξ ξ ξ ) = 1.
While, in Patriksson et al. (2015a, b), the term maintenance policy denotes a deterministic maintenance policy (stochastic maintenance policies are not employed in Patriksson et al. 2015a, b), it is here used with either of the attributes deterministic or stochastic. Stochastic maintenance policies are employed in the proof of Proposition 1(i), (ii).
Costs corresponding to the different decisions and/or policies must be defined. We assume that the replacement of any component n ∈ N incurs a replacement cost c n , and that a startup cost d is incurred each time at least one component is replaced. The current cost of the maintenance decision x ∈ X ξ ξ ξ is then defined as Letting S > 0 denote the planning horizon, the period of time during which the system must work is defined as the interval [0, S).
Definition 2 (Stochastic opportunistic replacement problem, SORP) Given a planning horizon S, a start-up cost d, component replacement costs c n , and probability distributions of the components' lives, for n ∈ N , find a deterministic maintenance policy, such that the expected sum of current costs during the planning period [0, S) is at minimum. A deterministic maintenance policy which finds an optimal solution to the SORP is called an optimal maintenance policy.
Definition 3 (Current problem) Given a planning horizon S, a start-up cost d, component replacement costs c n and probability distributions of the components' lives, for n ∈ N , the system state (s, ξ ξ ξ, a) ∈ V, and assuming that an optimal maintenance policy will be used for all maintenance decisions during the time span (s, S), find a maintenance decision x ∈ X ξ ξ ξ such that the expected sum of current costs during the planning period [s, S) is at minimum. If x * (t, ξ ξ ξ, a) is an optimal solution to the current problem at state (t, ξ ξ ξ, a) ∈ V, and p * (s, ξ ξ ξ, a) = x * (s, ξ ξ ξ, a) for all (s, ξ ξ ξ, a) ∈ V, then p * is an optimal maintenance policy. Thus, to maintain a system optimally, it is sufficient to solve each current problem that occurs.

Discretization of time in the current problem and the life distributions
The current problem at state (s, ξ ξ ξ, a) ∈ V is defined w.r.t. the continuous time interval [s, S). We define the time unit δ := S−s T +1 and consider the set T := {0, . . . , T } representing T + 1 discrete time points. For each t ∈ T , any individual of any component which fails during the time interval [s + tδ, s + tδ + δ) is treated as it would fail already at the time s + tδ (or that the individual is observed to be so close to failure that a replacement is necessary before the time s + tδ + δ). 1 Hence, the discretized version of the current problem is defined for the discrete set {s, s + δ, . . . , s + T δ} of time points-hereafter represented by the set T . Since decisions made at time t ∈ T depend on realizations of the system's states at the times {1, . . . , t} and of the decisions made at the times {0, . . . , t − 1}, the discretized current problem is a multi-stage stochastic programming problem (Birge and Louveaux 1997, Ch. III.7).
We define the N × N -matrix A := diag(a) and the N -vector 1 := (1, . . . , 1) . Then, A(1 − x) ≡ a − Ax denotes the components' ages after the replacement decision x. We let the function f p t : B N × R N → R represent the expected cost incurred by following the stochastic maintenance policy p over the time steps {t, . . . , T }, defined as where the function Q a t : and η η η ∈ B N is a random vector of component states and whose distribution is given by Pr(η n = 1) = 1 − Pr(η n = 0) = G n (a n − a n x n + 1), where G n is computed as in (3), below. The recourse function Q a t represents the expected future costs at time t. In the special case of p being a deterministic maintenance policy, the right-hand-side of (1a) can be simplified to For an optimal maintenance policy p * , the expected sum of current costs over the times {t, . . . , T } may be expressed as subject to x ≥ ξ ξ ξ, (2b) The current problem is to find f p * 0 (ξ ξ ξ, a) [i.e., solve the program (2) for the state (0, ξ ξ ξ, a)] and the corresponding optimal solution, typically without explicit knowledge of the optimal policy p * [i.e., without minimizing (2a) for each (t, ξ ξ ξ, a)].
The lives of the individuals i ∈ Z + of component n ∈ N are represented by the continuous 2 and independent identically distributed (i.i.d.) stochastic variablesT ni . Let τ n : R + → [0, 1] denote the probability density function ofT n1 , so that its cumulative probability function is given by The failure risk function G n : R + → [0, 1], defined as G n (a) := Pr(T n1 < a + 1 :T n1 ≥ a), denotes the probability that component n fails during the next time step, provided that it is in a non-failed state and of age a ≥ 0; it is calculated as We propose a lower bound on the optimal value of the SORP, as given by Definition 2. The bound is restricted to instances of the SORP for which each component n ∈ N has a stochastic life following the distribution F n , such that the failure risk function G n , defined in (3), is non-decreasing. 3 The bound is determined by minimizing the expected replacement costs for each component individually, as well as minimizing the expected start-up costs for replacements individually, followed by a summation of those costs. In a one component system (i.e., for n = N = 1), component n is replaced only when it fails; we denote by φ n the expected number of failures of such a system. In a systemfor which it holds that N > 1, the expected number of times that component n ∈ N must be replaced is at least φ n . Hence, the expected cost for replacing component n (excluding start-up costs, d) is at least c n φ n . To calculate the value of φ n , let the stochastic variable Y The expected value of the number of replacements of component n, i.e., Z n , may be found by determining the probability function for Z n , or as the sum over m of the probability that the mth individual of component n must be replaced; 4 Finding an analytic expression for φ n is possible only for special cases of the function F n . UnlessT n1 lack expected value or variance, normal approximations may be used for large enough values of m, since Y 3 When time is treated as continuous rather than discrete, the failure risk function G n is replaced by the hazard rate function λ n : R + → R + ; it is for our lower bound assumed to be non-decreasing and defined by λ n (a) The equivalences in (4) hold due to the relations We denote by Φ the least expected number of times when at least one component must be replaced; it is given by the case when all components are replaced whenever any component fails; this follows since the failure risk functions (or, hazard rates) are non-decreasing, implying that the expected time until the next system failure may only decrease if not all components are replaced. Hence, the expected start-up cost cannot be less than dΦ. To calculate the value of Φ, letT i denote the time between system failures i − 1 and i, and assume that all components are replaced at each system failure. Then, The distribution function ofT 1 is given by Analogously to the derivations for individual components, let Y (m) := m i=1T i denote the time for the mth system failure, and Z := max { m ∈ Z | Y (m) < S } the number of times replacements must be made. It then holds that [cf. (4)] We conclude that the sum of the least possible expected start-up costs and the least possible expected replacement costs for all components forms a lower bound on the optimal value of the SORP, namely dΦ + n∈N c n φ n .
The bound is based on the minimization of the replacement costs for each component, without regards for the start-up costs then incurred, and on minimization of the start-up costs without regards for the replacement costs for the components then incurred. It is only in special cases that both component replacement costs and start-up costs can be minimized simultaneously, so the bound is not tight in the general case.

Bounds on the value of the recourse function
Consider the expected minimum total cost resulting from the employment of an optimal maintenance policy as expressed in (2), along with the definition (1b)-(1c) of the recourse function. We will in this section establish bounds on the value of the recourse function which will prove useful in a decomposition method used for solving the current problem. In particular, we reconsider (Patriksson et al. 2015b, Prop. 1), for which we present a new and clearer proof of part (a) [i.e., Proposition 1(i)], and an improved version of part (b) [i.e., Proposition 1(ii)]. 5 We let G(a) denote the probability that a functioning system with component ages a ∈ R N + fails within the next time step, i.e., G(a) := 1 − n∈N 1 − G n (a n ) .
For failure distributions with non-decreasing failure risk functions (see Sect. 4), the next proposition states that the expected future cost is non-decreasing with (non-decreasing) component age and provides an upper bound on the increase of the expected future cost with increasing age.
For any two vectors a,â ∈ R N + we define the setN (â, a) := { n ∈ N |â n > a n }. Then, for each t ∈ T , we define the function K t : Proposition 1 Assume that each component n ∈ N has a failure distribution with a nondecreasing failure risk function G n , let p * be an optimal maintenance policy, let the age vectors a,â ∈ R N + be such that a ≤â, let the functions Q a t and K t , t ∈ T , be defined by (1b)-(1c) and (5), respectively. Then, for all t ∈ T the following hold: Proof (i) The inequality Q a t (0, p * ) ≤ Q a t (0, p) holds for all policies p, since by definition, an optimal policy p * provides the lowest possible expected future cost. Our idea is to find a stochastic policy p which, as applied to a system with component agesâ at time t, simulates an optimal policy p * applied to a system with component ages a at time t, i.e., such that Q a t (0, p) = Qâ t (0, p * ) holds. This means that in the policy p, each component n ∈N (â, a) fails artificially with some probability, at each time step until it is replaced. The probability that such a component n fails artificially at time t + k should equal the probability that component n-for component agesâ-fails at time t + k, even though-for component ages a-it does not fail at time t + k, conditioned on that it is not replaced at any of the times t + 1, . . . , t + k − 1. The sought probability is given by H n (a n ,â n , k) := G n (â n +k)−G n (a n +k) 1−G n (a n +k) , if G n (a n + k) < 1, Since the failure risk functions are non-decreasing, i.e., for any 0 ≤ a ≤â, the inequalities 0 ≤ G n (a) ≤ G n (â) ≤ 1 hold, the relations in (6) imply that G n (â n + k) = G n (a n + k) + 1 − G n (a n + k) H n (a n ,â n , k).
It follows that the inequalities 0 ≤ H n (a n ,â n , k) ≤ 1 hold. Therefore, H n (a n ,â n , k) is indeed a probability whenever k ≥ 0 and 0 ≤ a n ≤â n . To construct the stochastic policy p needed, let t, a, andâ be fixed, and define the mutually independent stochastic variables X (k) n ∈ Bin(1, H n (a n ,â n , k)) for n ∈ N and k ∈ Z >0 . Further, let α α α ∈ Z N >0 and define, for all n ∈ N and k ∈ Z >0 , , if α n = a n + k and a n <â n , ξ n , otherwise, andα (k) n := â n + k, if α n = a n + k and a n <â n , α n , otherwise.
The stochastic policy p is then defined as Consider two systems which, at time t, function and possess component ages a andâ, respectively. The policy p applied to the system with ages a, acts as the policy p * applied to the system with agesâ. 7 Each component n is treated according to its actual age and state, unless it was too young at the start (i.e., if a n <â n ) and has not been replaced since then (i.e., if α n = a n + k). It follows from the construction of p that the correct artificial ages and failure probabilities are given otherwise.
Hence, the equality Q a t (0, p) = Qâ t (0, p * ) holds and the proposition follows. (ii) To construct the stochastic policy p needed, let t, a, andâ be fixed and define the mutually independent stochastic variables 8 and define Let where p n and p * n denotes the nth component of the policy p and p * , respectively, and n ∈N (â, a), a+k1 and n ∈N (â, a), and n ∈N (â, a), 0, otherwise, For k ≤ 0, without loss of generality we let p(t + k, ξ ξ ξ, α α α ) = p * (t + k, ξ ξ ξ, α α α ). 9 Let t + K be the first replacement time after time t, where K ≥ 1 is a stochastic variable, and let "case a" and "caseâ" denote the cases with component ages at time t being a andâ, respectively. Due to the definition (8), the distributions of K are equivalent in "case a" and in "caseâ". In particular, it holds that where the second equality can be established by induction over k. The independence between failure times of different components means that given K = k, at time t + k the distributions of component states for components n / ∈N (â, a) are equivalent for the two cases. Further, given K = k and { ξ n } n∈N \N (â,a) , due to the definitions (7) and (9), the probabilities of the possible maintenance decisions-and hence also the corresponding expected costs-at time t + k are equivalent in the two cases. Combined with the fact that all components in N (â, a) are replaced at the first replacement time after t, the distribution of component states after the first replacement after time t will be equivalent in the two cases. It follows that Qâ t (0, p) = Q a t (0, p). We now have that Qâ t (0, p * ) ≤ Qâ t (0, p) = Q a t (0, p). The time for the first replacement after time t is t + K , whence, in case a, policy p may perform replacements not performed by policy p * , causing extra costs. After these extra replacements, the policy p leaves a system in which no component individual is older than its equivalent in the system left by the policy p * . Hence, according to Proposition 1(i), N (â, a)) holds, where R t (a,N (â, a)) denotes the extra cost incurred by replacing the components n ∈N (â, a) at time t + K given case a. To calculate the extra costs in case a, we partition the event "K = k" in case a into the two events "max n∈N { ξ n } = 1" (i.e., the system fails at time t + k) and "max n∈N { ξ n } = 0 and X (k) = 1" (i.e. the system does not fail at time t + k): If the system fails at time t + k, no extra start-up cost must be paid, but for each component n ∈N (â, a) an extra cost c n must be paid, unless it is to be replaced according to policy p * . 10 The probability that component n has failed at time t + k, given a system failure at time t + k, is G n (a n + k)/G(a + k1) (the event of system failure is implied by the failure of component n). The extra cost incurred by the event "K = k" due to system failure at time t + k, multiplied by its probability, is at least ⎛ If K = k even though there is no system failure at time t + k, then an extra cost of d must be paid, as well as the cost c n for each n ∈N (â, a). The extra cost incurred for K = k and a system failure at time t + k, multiplied by the probability of this event, is ⎛ Hence, the expected value of the extra cost (from the component ages beingâ instead of a) incurred by system failures over the times t + 1, . . . , T , is not less than the sum, over k = 1, . . . , T − t, of (10) and (11), which equals K t (â, a), defined in (5). The proposition follows since the inequality K t (â, a) ≥ R t (a,N (â, a)) holds.
The result in Patriksson et al. (2015b, Prop. 2) relates the values of the recourse function for two systems without restrictions on component ages; its proof contains, however, an error (a missing superscript). Based on our improvement of the bounds in Proposition 1, we next present a (corrected and) stronger version of Patriksson et al. (2015b, Prop. 2).
We define the component-wise maximum as max { x, y } : Proposition 2 Assume that each component n ∈ N possesses a failure distribution with nondecreasing failure risk function G n , let p * be an optimal maintenance policy, and let K t be defined as in (5). Then, for all x, y ∈ B N and all t ∈ T it holds that Proof Since the equivalence G n (a n − a n max { x n , y n } + k) = G n (k) holds for any k ≥ 0 and all n ∈N (a − Ay, a − A max { x, y }), it follows that where the first and second inequalities are obtained by application of Proposition 1(ii) and (i), respectively. The proposition follows.

The ORPIL
The For full details of the ORPIL, see Patriksson et al. (2015a).

Solving the discretized current problem
Two approaches to solve a two-stage approximation of the discretized current problem were presented in Patriksson et al. (2015b): (i) a two-stage deterministic equivalent, and (ii) a twostage decomposition method, the latter of which was shown to be the computationally most efficient. Hence, we use approach (ii), as presented in Sect. 8.2. Both approaches require, in each iteration, the calculation of an approximate value of the recourse function, as presented in Sect. 8.1.

The two-stage sample average approximation
Let Ω denote the set of all possible scenarios ω for the discretized current problem (2). A scenario ω ∈ Ω is defined by a life T ni (ω) ∈ Z >0 for each individual i ∈ I n ∪Ĩ n of each component n ∈ N , and is realized with probability Pr(ω).
Each scenario, ω, defines an instance of ORPIL (see Sect. 7) with q n = |I n | = T + 1, n ∈ N ; its solution is a maintenance schedule, which is optimal provided that scenario ω is realized (see Def. 4). Such large values of q n make ORPIL computationally demanding, which is problematic when a large number of scenarios are to be realized. Since it is unlikely that T + 1 non-identical life individuals are needed for all components, it may suffice to let 11 In Patriksson et al. (2015a), all components were given the same number of non-identical life individuals, i.e., q n = q m for all n, m ∈ N . Such a restriction is not necessary, and may be computationally unfavorable when the ORPIL is used in a solution scheme for the SORP. 12 Since the life of each individual is at least δ, corresponding to one time step, at most T + 1 individuals of each component are needed. 13 A maintenance schedule is defined as a set of decisions regarding which components to replace at what time steps. The schedule is feasible if and only if each individual i ∈ I n (i ∈Ĩ n ) of each component n ∈ N spends no more than T ni (T n ) time units in the system. q n equal the number of individuals of component n ∈ N that are likely to be used prior to the planning horizon; we anticipate a computational advantage from this construction.
Let q = (q 1 , . . . , q N ), and let the sampled multisetΩ q approximate the set Ω, witĥ where E[T n ] denotes a discretization of the expected life (calculated based on continuous life) of an individual of component n. A scenario ω ∈Ω q is sampled by sampling T n1 (ω) from the life distribution of component n, given T n1 (ω) ≥ a n + 1; sampling T ni (ω) for i ∈ I n \ {1} from the life distribution of component n, given T ni (ω) ≥ 1; and letting T ni (ω) := E[T n ] for i ∈Ĩ n . We let Pr(ω) = |Ω q | −1 for all ω ∈Ω q . For a first-stage decisionx n ≡v 0 n1 , n ∈ N , define the constraints Further, define the function F ORPIL : B N ×Ω q → R such that F ORPIL (x, ω) equals the optimal objective value of the ORPIL for scenario ω, with the constraints (12) added to ensure that the first-stage decision is x. A two-stage sample average approximation of Q a 0 (x, p * ) is given byQ Given a decisionx, the expected sum of current costs during the planning period 14 [s, S) is approximated by computing the values F ORPIL (x, ω), ω ∈Ω q . For the optimization in (13) to be tractable, the number |Ω q | of scenarios must be sufficiently small, but too few scenarios would make the two-stage sample average approximation unstable. Effects of varying |Ω q | was studied by Patriksson et al. (2015b, Section 7.2), who concluded that |Ω q | = 100 is sufficiently large for the instances studied.

The decomposition method
We employ integer L-shaped decomposition; see Laporte and Louveaux (1993), to which our description below refers. The so-called deterministic equivalent optimization problem of our two-stage stochastic optimization problem is to where the recourse function Q : The recourse function is iteratively approximated using a set of optimality cuts; for our problem setting, the L-shaped decomposition consists of solving a sequence, for = 0, 1, . . ., of current problems to minimize x,z 0 ,θ c x + dz 0 + θ, subject to x ≥ ξ ξ ξ, (15b) For = 0 all the optimality cuts (15d) are relaxed. Let (x,z 0 ,θ) denote the solution to (15); if it holds thatθ ≥ Q(x,z 0 ), then (x,z 0 ,θ) is optimal in (14); otherwise let := + 1 and add an optimality cut to next current problem. Feasible maintenance at time s corresponds to the constraints (15b)-(15c); if feasible replacements are made at time s, the set of feasible future replacements is non-empty (e.g., replace all non-functional components at times s + 1, s + 2, . . .). The optimality cuts provide lower bounds (at least one of which is tight) on the recourse function Q(·), and the bounds from Prop. 2 are defined by highly non-linear functions. Further, for computational efficiency the exact evaluation of Q a t (x, p * ) is replaced by the approximation (13).
We assume that ξ n = 1 for at least one n ∈ N (i.e., at least one component is not functioning; recall Sect. 3), since otherwise no maintenance would be performed; hence, the start-up cost d will always be included in the current problem (Definition 3). Let L denote an index set corresponding to the optimality cuts generated so far, let y ∈ X ξ ξ ξ denote a maintenance decision, ∈ L (i.e., the optimal value for x in (15) when − 1 optimality cuts (15d) have been generated), let the function K 0 be defined as in (5). The current problem is then iteratively approximated as Proposition 3 Define the recourse function Q a t as in (1), assume that the failure risk functions G n , n ∈ N , are non-decreasing, and let the function f p * 0,L be defined as in (16). Then, it holds that Since the recourse function is non-negative, i.e., Q a 0 (x, p * ) ≥ 0, it then follows that Letx ∈ X ξ ξ ξ be optimal in (16). The inequality then holds for all x ∈ X ξ ξ ξ and the proposition follows.
If f p * 0,L (ξ ξ ξ, a) = c y + d + Q a 0 (y , p * ) holds for some ∈ L, then y is optimal in (2). ξ ξ, a) holds, i.e., (16) yields a lower bound on the optimal value of (2), and the recourse function for the best lower bound found is evaluated. The method is formalized in Algoritm 1, in whichQ a s (y ) approximates Q a 0 (y , p * ).

Numerical experiments
We have performed numerical experiments to show the relevance and effects of the theoretical results developed. The test instances and the policies employed are presented in Sects. 9.1 and 9.2, respectively. Section 9.3 reports on the goodness (i.e., tightness) of the lower bound developed in Sect. 5, while Sect. 9.4 presents the effect of the improved bound on the recourse function [Proposition 1(ii)] on the optimization process.
For more elaborate numerical tests of the optimization policy, including computation times, see Patriksson et al. (2015b), which investigates the effect of the number of nonidentical life individuals as well as of the number of scenarios used.

The test instances
For each of the four test instances, T1-T4, each component n has a Weibull-distributed life with scale parameter α n and shape parameter β n . The discrete time points for maintenance decisions are given by T := 0, . . . , max { (S − s)/δ , 3 } . The data are shown in Table 1.
• Instance T1 (Patriksson et al. 2015b) was constructed to illustrate a poor performance of (general) age-based policies (see Sect. 9.2). (Our use of a different age-based policy than that in Patriksson et al. (2015b) does, however, not verify this observation; see Sect. 9.3.) • Instance T2 (Patriksson et al. 2015b) was constructed to illustrate the performance of algorithms when the spread of the components' lives is small (i.e., high values of the shape parameter). Then the SORP tends to be more deterministic (like the ORP; see Almgren et al. 2012) and optimization methods are expected to outperform heuristics. • Instance T3 contains more components and time steps than T1 and T2. The ratio between a component's replacement cost c n and its expected life α n Γ (1 + 1/β n ) is similar for all n. The start-up cost d is within the range of the component replacement costs, c n , n ∈ N , so that each component will cause a significant portion of the total replacement costs. 15 T3 is intended to force the optimization policy to create many optimality cuts; hence, our improved lower bound on the recourse function yields a more efficient implementation of the optimization policy. 16 • In instance T4 the costs and expected lives of components are related as in T3, albeit the "regularity" is avoided according to the following: After deciding upon T , N , and d, for each n ∈ N , c n and β n are sampled from a uniform distribution on For each test instance, 100 different scenarios were generated and used for performance testing. These scenarios are called simulation scenarios to be distinguished from the scenarios used within the optimization policy (see Sect. 9.2).

The policies
For testing and illustration we have implemented and used the following four policies.  16 For a discretized current problem with N components in state ξ ξ ξ ∈ B N , at most |X ξ ξ ξ | = 2 N −1 T ξ ξ ξ optimality cuts can be generated, 1 T ξ ξ ξ being the number of failed components. Hence, for small values of N (such as in T1 and T2), differences in numbers of generated optimality cuts may seem insignificant.
Algorithm 2: Simulated annealing input : A set N of components with life distributions F n : R + → [0, 1] and initial guessesσ n for soft life thresholds, n ∈ N . A planning horizon S > 0, a number 0 < S, and set S of simulation scenarios. An upper limit r max on # iterations with no improvement. output: Soft life thresholds σ n , n ∈ N C best ← ∞; r ← 0; while r < r max do C ← average maintenance cost over all scenarios in S, given soft life thresholdsσ n ; if C < C best then C best ← C; σ n ←σ n , n ∈ N ; r ← 0; end σ n ← min max { N(σ n , σ n /(r + 1)), }, S + , n ∈ N , where N(·, ·) is a random number from the normal distribution; r ← r + 1; end return σ n , n ∈ N ; components which will fail within the next time step are also considered to be failed. Hence, several components may "fail" simultaneously.
• The optimization policy (the two-stage decomposition method; see Sect. 8). It occurs in an old implementation, based on the lower bound on the recourse function from Patriksson et al. (2015b), and a new implementation, based on the (tighter) lower bound on the recourse function given in Proposition 1(ii). Hence, the new implementation is expected to run faster than the old one. The optimization policy requires an arbitrary feasible starting solution. Both implementations are initialized at the solution provided by the run-tofailure policy, which is reasonable when the planning horizon is near. In general, however, the optimization policy may be improved through better starting solutions (obtained by some fast policy). • The expected value (EV) optimization policy, which is used by Almgren et al. (2012) to solve the SORP. It approximates the current problem by-for each componentassigning to the current individual its expected life, conditioned on its age, and by assigning to all future individuals the component's expected life. 17 • An age-based policy: at failure of any component, maintenance is performed and each component whose age surpasses its assigned soft life threshold is also replaced. This policy is akin to that of Crocker and Kumar (2000), in which maintenance is enforced whenever the age of any component reaches its assigned hard life threshold. Hard life thresholds are not needed in our set-up, in which failures are assumed either to be detectable only very shortly before occurring, or to cause no extra costs. The soft life thresholds are variables in an optimization problem to minimize the expected maintenance cost over the planning period; in our numerical experiments (see Sect. 9.3) it was solved using simulated annealing (Algorithm 2), with the components' expected lives used as initial guesses for the soft life thresholds. To prevent from getting stuck in poor local optima, the best solution from 100 separate runs of Algorithm 2 was used. Table 2 presents, for the instances T1-T4, lower bounds on the optimal maintenance costs and the maintenance costs resulting from the four policies described in Sect. 9.2. It also Table 2 Lower bounds on the expected value of costs (c x , x = lower bound) and average costs (100 simulation scenarios/instance) resulting from the four policies (c x , x = policy). The gaps are computed as (c policy /c lower bound − 1) · 100%. In the optimization policy, we used q n = 2, n ∈ N , for T1-T3, and q = (4, 2, 2, 3, 2, 6, 2) for T4   Instance  T1  T2  T3 T4 presents the relative differences in maintenance costs (so-called gaps) between each policy and the lower bound. For each instance, there is a significant gap between the costs resulting from the optimization policy and the respective lower bound, but as concluded in Sect. 5, the bound is not expected to be tight, so it is not likely that any policy will produce solutions with costs close to the lower bound. In our experiments, the age-based policy and the optimization policy resulted in approximately the same average cost for all four instances. A different age-based policy was employed by Patriksson et al. (2015b), and that age-based policy performed significantly worse for some instances, in particular instance T1.

The number of generated optimality cuts
While applying the old and the new implementation of the optimization policy on 100 simulation scenarios of each test instance, a vast number of discretized current problems was solved and for each of which a number of optimality cuts were generated; the corresponding frequencies are shown in Fig. 1a-d. Although the same simulation scenarios were used for both implementations of the optimization policy, they occasionally result in different maintenance decisions. One reason is that optimality cuts are typically generated in different orders, resulting in different scenarios being used for ORPIL, which leads to different approximations of the recourse functions. Another reason is that-in either implementation-the first optimal maintenance decision found is selected; hence, in case of multiple optima, the order in which the optimality cuts are generated has an impact on the optimal solution selected. Once the selected maintenance decisions within a scenario differ between the implementations, the subsequent current problems will also differ. In Fig. 1 only current problems encountered by both implementations are counted: a higher frequency of the scenarios are solved by fewer optimality cuts with the new optimization policy than with the old one, e.g. in instance T1 the new policy will most frequently require 2 optimality cuts (and slightly above 2 cuts on average), whereas the old policy will most frequently require 4 optimality cuts (and above 2.5 cuts on average). Each generated optimality cut results in one instance of ORPIL to be solved; with ORPIL being the computationally heaviest part of the algorithm, the computation time is approximately proportional to the number of generated optimality cuts.
The histograms in Fig. 1 represent tilted distributions of the numbers of optimality cuts generated: Current problems (occuring during the solution of the simulation scenarios) being excluded (due to the implementations proposing different maintenance decisions at earlier stages) are more prevalent when close to than far from the planning horizon. Hence, very few optimality cuts are generated when near to the horizon-often then an obvious decision is to replace only failed components (i.e., the starting solution in Algorithm 1).

Conclusions and future research
This article deals with a stochastic opportunistic replacement problem for a system of components with economic dependencies. In particular, we attempt to solve the so-called current problem, i.e., decide "which components to replace now", at a time when at least one component must be replaced. A simple but non-trivial lower bound on the expected maintenance cost is presented; it is based on the cost that would result from a perfect coordination of replacements (using a minimal number of expected replacement occasions) with no loss of component life (due to early replacements). The lower bound cannot be expected to be tight (i.e., equal the optimal value), except in some extreme cases. It provides, however, a restriction on the gain reachable by any optimization policy.
We improve-in terms of computing efficiency-a previously developed decomposition method for a two-stage approximation of the current problem: the improved lower bound on the recourse function value results in fewer optimality cuts being generated. Also the two-stage approximation is generalized, allowing for different numbers of non-identical life individuals for different components.
Experiments showing that the optimization policy usually outperforms simpler policies was presented by Patriksson et al. (2015b). In contrast, our experiments indicate that a fairly simple age-based policy performs on par with the optimization policy.