An effective approach to reducing strategy space for maintenance optimisation of multistate series–parallel systems

Reliability Engineering and System Safety 138 (2015) 40–53

Contents lists available at ScienceDirect

Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress

An effective approach to reducing strategy space for maintenance optimisation of multistate series–parallel systems Yifan Zhou a,n, Tian Ran Lin d, Yong Sun c, Yangqing Bian a, Lin Ma b a

School of Mechanical Engineering, Southeast University, Nanjing 211189, PR China School of Chemistry, Physics and Mechanical Engineering, Queensland University of Technology, 2 George St, Brisbane 4001, QLD, Australia c CSIRO, QCAT, 1 Technology Court, Pullenvale 4069, QLD, Australia d School of Mechanical Engineering, Qingdao Technological University, Qingdao 266520, PR China b

art ic l e i nf o

a b s t r a c t

Article history: Received 19 December 2013 Received in revised form 13 January 2015 Accepted 15 January 2015 Available online 23 January 2015

Maintenance optimisation of series–parallel systems is a research topic of practical significance. Nevertheless, a cost-effective maintenance strategy is difficult to obtain due to the large strategy space for maintenance optimisation of such systems. The heuristic algorithm is often employed to deal with this problem. However, the solution obtained by the heuristic algorithm is not always the global optimum and the algorithm itself can be very time consuming. An alternative method based on linear programming is thus developed in this paper to overcome such difficulties by reducing strategy space of maintenance optimisation. A theoretical proof is provided in the paper to verify that the proposed method is at least as effective as the existing methods for strategy space reduction. Numerical examples for maintenance optimisation of series–parallel systems having multistate components and considering both economic dependence among components and multiple-level imperfect maintenance are also presented. The simulation results confirm that the proposed method is more effective than the existing methods in removing inappropriate maintenance strategies of multistate series–parallel systems. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Multistate system Maintenance optimisation Strategy space reduction Linear programming Heuristic algorithm

1. Introduction Maintenance optimisation of series–parallel systems is an important research topic and has attracted strong research interest from many pioneering researchers in the last two decades [1–5]. Many practical industrial systems can be regarded as series– parallel systems. A typical example is a production line in a manufacturing factory. A production line can have multiple production phases. Each phase can have several production units organised in parallel to enhance the performance of the system [6,7]. These production units can have multiple performance levels according to their actual degradation states. Therefore, performance of a series–parallel system is a function of the degradation state of its components. Moreover, there are various types of dependence, such as economic dependence, structure dependence, and stochastic dependence, among the components [8]. A system-level maintenance strategy, which is a combination of the component-level maintenance strategies, is thus needed to reduce the component degradation and to enhance the production rate of a system. However, the system-level strategy space of a typical

n

Corresponding author. E-mail address: [email protected] (Y. Zhou).

http://dx.doi.org/10.1016/j.ress.2015.01.018 0951-8320/& 2015 Elsevier Ltd. All rights reserved.

series–parallel system is often large, which poses a problem for system maintenance optimisation. For example, there are more than six million possible combinations of component-level maintenance strategies for the case study described in Ref. [9], even though the system contains only 14 components and has less than five degradation states for each component. The heuristic algorithm is often employed to find a costeffective system-level maintenance strategy when the strategy space is large. Huang and Wang [10] used the genetic algorithm to optimise the imperfect maintenance strategy of a multistate system. In their paper, the component degradation process was described by a non-homogeneous continuous-time Markov model. The genetic algorithm was also employed by Liu and Huang [11] to solve the selective maintenance problem of a series–parallel system under imperfect maintenance. Their work [11] was further extended by Pandey et al. [12] who employed a hybrid model to describe the effects of imperfect maintenance. In a later publication, Pandey et al. [13] developed a more complex maintenance optimisation problem by considering multistate components in the system. The differential evolution algorithm was employed in both Refs. [12,13] to solve the maintenance optimisation problem. The genetic algorithm was also employed by Vu et al. [14] in their study of maintenance planning and dynamic grouping problem of complex systems with either positive or negative economic

Y. Zhou et al. / Reliability Engineering and System Safety 138 (2015) 40–53

Nomenclature UGF ACO N Nn NPRO NPRO,n NPRO;n Nn(i, j0 ) Kn θ θn θPR,n θOP,n,w

θn Θn Θn (λn)i,j (μn)i,j (cn)i,j CST,n γn,i rp R RP C M;n

universal generating function ant colony optimisation number of subsystems in the series–parallel system number of components in subsystem n number of possible production rate levels of the series–parallel system number of possible production rate levels of subsystem n number of possible production rate levels of the system excluding subsystem n number of components in state i, when subsystem n is in state j0 . number of states of components in subsystem n the parameter vector of the maintenance strategy for the series–parallel system the parameter vector of the maintenance strategy for subsystem n the preventive replacement threshold of a component in subsystem n the threshold of the imperfect opportunistic maintenance that can improve the condition of a component in subsystem n from state i to state (i  w) the parameter vector of the maintenance strategy for the system excluding subsystem n the maintenance strategy space of subsystem n the maintenance strategy space of the system excluding subsystem n the transition rate of a component in subsystem n from state i to state j the repair rate of a maintenance activity that restores a component in subsystem n from state i to state j the cost of a maintenance activity that restores a component in subsystem n from state i to state j the setup cost before a series of maintenance activities for subsystem n the production rate of a component in subsystem n when the state of the component is i reward per unit production rate of the series– parallel system the average revenue per unit time of the series– parallel system the average production revenue per unit time of the series–parallel system the average maintenance cost per unit time of subsystem n

dependence. Doostparast et al. [15] used the simulated annealing algorithm to identify the optimal maintenance plan to minimise the total maintenance cost with respect to a desired level of system reliability. Rami et al. [16] obtained the optimal sequence of maintenance actions of a series–parallel system using the harmony search optimisation method. The heuristic algorithm has also been employed in several recent papers to deal with more complex maintenance optimisation problems. For example, several recent papers used the heuristic algorithm to optimise the system structure and the maintenance strategy of a multi-state system simultaneously [17–21]. Other works formulated the maintenance optimisation of a multistate system as a multiobjective optimisation problem [22,23]. Additional to the maintenance cost, other factors such as the average product number awaiting processing, system reliability and availability have also

41

Cn(i, θn) the cost of the maintenance activity of a component in subsystem n when the state of the component is i and the strategy adopted is θn
S  λn ðθn Þ i0 ;j0 the transition rate of subsystem n from state i0 to state j0 under maintenance strategy θn ðπn ðθn ÞÞi0 the steady-state probability that subsystem n is in state i0 under maintenance strategy θn U S ðz Þ the UGF representing the production rate distribution of the series–parallel system un ðzÞ the UGF representing the production rate distribution of subsystem n gv level v of the production rate of the system pv the probability that the system production rate is gv g n;v level v of the production rate of subsystem n pn;v the probability that the production rate of subsystem n is g n;v g n;v level v of the production rate of the system excluding subsystem n pn;v the probability that the production rate of the system excluding subsystem n is g n;v pnn;v the optimal value of pn;v , which maximises the average system revenue difference under two strategies of subsystem n maxÞ pðn;v the upper bound of the production rate distribution of subsystem n in stochastic ordering minÞ pðn;v the lower bound of the production rate distribution of subsystem n in stochastic ordering maxÞ pðn;v the upper bound of the production rate distribution of the system excluding subsystem n in stochastic ordering minÞ pðn;v the lower bound of the production rate distribution of the system excluding subsystem n in stochastic ordering
 UBR θn ; θ0n the upper bound of the average revenue difference when the maintenance strategy of subsystem n changes from θ0n to θn The designations of the subscripts used in the terminologies described in the Nomenclature and in the subsequent text are also defined as follows: i and j i0 and j0 n and m l v and u

the the the the the

indices of a component state indices of a subsystem state indices of a subsystem index of a component indices of production rate levels

been included as objective functions in the above-mentioned papers. In contrast, a number of recent publications focused on improving the accuracy and efficiency of the heuristic algorithm during maintenance optimisation. For example, Moghaddam and Usher [24] used both the genetic algorithm and the simulated annealing algorithm for maintenance optimisation of a multi-component system. The effectiveness and efficiency of the two algorithms were compared and discussed in their paper. Bris et al. [4] employed the genetic algorithm to optimise the periodic preventive maintenance plan of a series–parallel system. Samrout et al. [5] utilised the ant colony optimisation (ACO) algorithm to further improve the optimisation result obtained in Ref. [4]. Recently, Lin and Wang proposed a hybrid genetic algorithm [25] and an improved particle swarm optimisation [26] to tackle the optimisation problem

42

Y. Zhou et al. / Reliability Engineering and System Safety 138 (2015) 40–53

described in Ref. [4]. The results presented in Refs. [25,26] showed that both the two methods can render a more cost-effective maintenance strategy than those of the previous published methods. Although the heuristic algorithm is an effective tool in finding an appropriate solution in a large strategy space, it could be computationally expensive, and the maintenance strategy derived by the algorithm could be a local optimal solution [27]. Besides the heuristic algorithm, alternative approaches have also been developed to deal with the large optimisation space problem by reducing the maintenance strategy space of a system. A system-level maintenance strategy is the combination of subsystem-level maintenance strategies. Therefore, the systemlevel strategy space can be reduced substantially if inappropriate subsystem-level strategies are excluded during the optimisation process. Only a few studies have been devoted to strategy space reduction of series–parallel systems thus far. Galante and Passannanti [28] developed an algorithm to identify and remove the dominated subsystem-level strategies incurring a high maintenance cost and low subsystem-level reliability. Cao et al. [27] proposed a similar method focusing on the redundancy allocation problem of a series–parallel system. A major limitation of these works is that only systems with binary-state components were considered. Furthermore, system reliability was the only performance indicator during optimisation. On the contrary, the component degradation of a series–parallel system is often a gradual process with multiple states. Moreover, the performance of a multistate system can be described more accurately by the distribution of production rates instead of a single scalar value (i.e., reliability). Zhou et al. [29] proposed a strategy space reduction method based on stochastic ordering to address this issue. In their work, the subsystem-level strategies causing a high maintenance cost and a low production rate in stochastic ordering were excluded during the optimisation process. Although the method presented in Ref. [29] can significantly reduce the strategy space, the criterion used to eliminate inappropriate subsystem-level strategies considers only the current subsystem. Nonetheless, influence of subsystem-level maintenance strategies on the objective function of maintenance optimisation also depends on the maintenance strategies of other subsystems. An improved criterion is then developed in this paper to remove the inappropriate subsystem-level strategies more effectively by considering the inter-connection between subsystems. The method developed in this paper uses the upper bound of expected system revenue difference under two different subsystem-level strategies as the criterion to reduce the strategy space. The revenue of a system is defined as the difference between the production revenue and the maintenance cost of the system. When the maintenance strategy of a subsystem varies, the system revenue will change. If the maximum possible increment of the average system revenue is below zero, the new maintenance strategy of the subsystem will be eliminated from the strategy space during maintenance optimisation. A method for reducing maintenance strategy space based on bounding values has been proposed by the same authors in Ref. [30]. However, the bounding values estimated in Ref. [30] were based on the property of stochastic ordering only, which are too conservative and hinder the effectiveness of strategy space reduction. To overcome this problem, this paper converts the calculation of the bounding values of the average revenue difference to a linear programming problem that can be solved efficiently. A theoretical proof is also given to validate that the bounding values calculated using linear programming are at least as accurate as those derived in Ref. [30]. The numerical study presented in the last part of the paper also confirms that the new method is much more effective than the existing methods, particularly when the number of subsystems is large.

The rest of the paper is organised as follows: Section 2 introduces the formulation of the maintenance optimisation problem and assumptions made in this paper. Section 3 develops a method to calculate the objective function of the maintenance optimisation problem. Section 4 describes the linear programming based strategy reduction method. A numerical example is presented in Section 5 to validate the method and its effectiveness. The main findings are summarised in Section 6.

2. Problem formulation and assumptions Existing research for imperfect maintenance of series–parallel systems often assumed that the components in a system have only two states (i.e., working and failed) [12,31,32]. Nevertheless, the degradation of a component is usually a gradual process and can have multiple states. A series–parallel system having multistate components under imperfect maintenance is considered in this study. The economic dependence among components within a subsystem is also included in the modelling. 2.1. Problem formulation The system considered in this study consists of N subsystems inter-connected in series. Subsystem n contains Nn identical parallel components. It is assumed that the degradation process of these components follows a continuous-time discrete-state Markov chain. Each component in subsystem n has Kn states, where state 1 represents the perfect condition, and state Kn is the failure state. The transition duration of the component in subsystem n from state i to state j follows an exponential distribution with a rate parameter ðλn Þi;j , where λn is the matrix of transition intensities of components in subsystem n. The production rate of a component in subsystem n is γ n;i when the state of the component is i. For a series–parallel system, the system production rate is the minimum production rate of all the subsystems, i.e., G ¼ min ðGn Þ, where Gn and G are respectively the production n ¼ 1;…;N

rates of subsystem n and the whole system. On the other hand, the production rate of a subsystem is the sum of the production rates P n of its components, i.e., Gn ¼ N l ¼ 1 Gn;l , where Gn;l is the production rate of component l in subsystem n. It is further defined that the reward per unit production rate of the series–parallel system is rp. A multi-threshold model similar to Ref. [33] is adopted in this study to describe the maintenance strategy. A component in subsystem n is preventively maintained to a faultless condition when its state reaches the threshold θPR,n. When the state of a component in subsystem n is θOP,n,w þ 1 4iZθOP,n,w and at least one other component in the same subsystem is under maintenance, the former will be imperfectly repaired to a state (i  w). For a failed component, corrective maintenance is a compulsory action to restore it to a brand new state. Consequently, the maintenance activity applied to a component depends on the component state and on whether other components in the same subsystem are under maintenance. A more cost-effective strategy could be obtained if the maintenance activity of a component also considers the condition of other components. However, this additional consideration will make it difficult to model the maintenance strategy and the system degradation process, which is thus not considered in this study. The durations of different maintenance activities are all assumed to follow the exponential distribution, and the repair rates of the maintenance activities for components in subsystem n are modelled uniformly by a matrix μn. The element (μn)i,j in matrix μn denotes the repair rate of a maintenance activity that restores a component from state i to state j. The costs of

Y. Zhou et al. / Reliability Engineering and System Safety 138 (2015) 40–53

maintenance activities are modelled by a matrix cn . The element ðcn Þi;j in the cost matrix denotes the cost of a maintenance activity that restores a component from state i to state j. The economic dependence is introduced here by a setup cost that is entailed to start the maintenance of a subsystem. The setup cost before a series of maintenance activities for subsystem n is CST,n. The multi-threshold maintenance strategy of the components in subsystem n can be  described by a vector θn ¼  θPR;n θOP;n;K n  2 θOP;n;K n  3 ⋯ θOP;n;1 , where K n Z θPR;n Z θOP;n;K n  2 Z θOP;n;K n  3 Z ⋯ Z θOP;n;1 Z2 and ( θOP;n;w þ 1 4 w þ 1 θOP;n;w Z w þ 1 : ð1Þ θOP;n;w ¼ θOP;n;w þ 1 θOP;n;w þ 1 r w þ 1 when θOP;n;w ¼ θOP;n;w þ 1 , the maintenance activity that restores the component in subsystem n by w states will not be adopted. Therefore, by setting the values of the thresholds, the multithreshold model can be used to describe different types of maintenance strategies, e.g., run to failure or preventive replacement with imperfect opportunistic maintenance. The system-level strategy θ¼ [θ1 θ2 … θN] is the combination of subsystem-level strategies. The objective function of this study is the expected revenue per unit time of the series–parallel system, which is described as: XN R ¼ RP  C ; ð2Þ n ¼ 1 M;n where RP is the expected production revenue per unit time of the system, and C M;n is the expected maintenance cost per unit time of subsystem n. 2.2. Assumptions The major assumptions made in this work are summarised below:

 Components in a subsystem are identical. This assumption is

 





  

reasonable because practical subsystems often use the same components to facilitate the training of maintenance technicians and to reduce the cost of spare parts. If a component is under maintenance, its degradation process stops and its production rate is zero. Other components, on the other hand, continue to work at their original production rates. If the degradation of a component crosses the opportunistic maintenance threshold during the maintenance of other components in the same subsystem, the component will be placed into maintenance immediately without additional setup cost. Once the preventive replacement starts in a subsystem, all the components in the subsystem whose degradation state crosses the opportunistic maintenance threshold will be maintained at the same time. The economic dependence among different subsystems is not considered, and each subsystem is maintained independently. This assumption is consistent with many practical industry applications in which subsystems are located in different places and maintained by different personnel. The states of the components are always known; thus, inspection scheduling is not considered. Resources to carry out maintenance (e.g., spare parts and maintenance personnel) are always available. The structure of the system is fixed, and the redundancy optimisation is not investigated.

3. Calculation of the objective function The objective function used in this study is the average revenue per unit time, which is defined as the difference between the

43

average production revenue and the average maintenance cost. A steady-state analysis of the Markov chain for the productivity analysis of a repairable system was developed by Gupta and Bhattacharya [34]. Some recent research further considers practical factors such as influence of buffer [35]. For multistate systems, the universal generating function (UGF) is often used to calculate the system performance by combining the steady-state analysis results of components [36]. The UGF is also employed in this paper to evaluate the system-level production rate. The average production revenue and the average maintenance cost are both assessed based on the steady-state distributions of subsystems. Subsequently, the steady-state analysis of subsystems has to be undertaken first to enable the evaluation of the average maintenance cost and the average production revenue. 3.1. Steady-state analysis of subsystems The degradation process of components in a subsystem depends on each other when opportunistic maintenance is adopted. Therefore, the degradation process of the entire subsystem needs to be modelled. Based on our assumptions, the components in a subsystem are identical. Thus, the sorted combi  nation of component states (i.e., X n;1 ðt Þ X n;2 ðt Þ ⋯X n;Mn ðt Þ , where X n;l ðt Þ Z X n;l þ 1 ðt Þ) can be used to represent the health state of the subsystem. Due to opportunistic maintenance, the change of the health state of a subsystem does not follow the Markovian property. For example, considering the case that the state of the most degraded component in subsystem n satisfies θOP;n;1 r X n;1 ðt Þ o θPR;n , if the other components in subsystem n whose health states are above θPR;n have just been replaced, subsystem n is under opportunistic maintenance. Conversely, subsystem n is not under opportunistic maintenance when preventive maintenance is not triggered. To retain the Markovian property, the state of subsystem n at time t, X n ðt Þ, is divided into three mutually exclusive sets according to the current component state and maintenance activity: 8 Y > < n;1 X n ðt Þ A Y n;2 > :Y n;3

X n;1 ðt Þ o θPR;n and no component is under maintenance; θOP;n;1 r X n;1 ðt Þ o θPR;n and some components are under maintenance; : X n;1 ðt Þ Z θPR;n :

ð3Þ The size of the three sets Y n;1 , Y n;2 , and Y n;3 are M n;1 , M n;2 , and M n;3 , respectively. Therefore, the number of possible states of subsystem n is M n;1 þ M n;2 þ M n;3 , i.e., X n ðt Þ ¼ 1; 2; …; M n;1 þ M n;2 þ M n;3 . To analyse the steady-state distribution of subsystem n under strategy θn, the transition intensity matrix of the subsystem, λSn ðθn Þ A R ðMn;1 þ Mn;2 þ Mn;3 ÞðMn;1 þ Mn;2 þ Mn;3 Þ , should be established first. During each state transition of a subsystem, only one component in the subsystem changes its state. Therefore, the i0 th row of λSn ðθn Þ is calculated using the process given in Table 1. Once the transition intensity matrix λSn ðθn Þ is obtained, the steady-state probability vector of subsystem n under strategy θn , πn ðθn Þ A R ðMn;1 þ Mn;2 þ Mn;3 Þ1 , can be derived by solving the balance equations of the Markov process. Here, ðπn ðθn ÞÞi0 is the steady-state probability that subsystem n is in state i0 . 3.2. Calculation of the average maintenance cost The average maintenance cost per unit time for each subsystem can be calculated separately due to the assumption that the subsystems are maintained independently. The maintenance cost of a subsystem is the sum of the setup cost and the cost incurred by various maintenance activities such as corrective replacement, preventive maintenance, and opportunistic imperfect maintenance.

44

Y. Zhou et al. / Reliability Engineering and System Safety 138 (2015) 40–53

Table 1 The process of calculating the i0 th row of the transition matrix. Assume the current state of subsystem n is i0 , Xn ¼ i0 , and let
S  λn ðθn Þ i0 ;j0 ¼ 0; j0 ¼ 1; 2; …; M 1n þ M 2n þ M 3n For each component l in subsystem n Step 1: Assume that the current state of component l is i, the next possible state of component l under strategy θn can be: 8 0 > < j A ði; K n  i A Y n;1 or i o θOP;n;1 j¼1 > : j ¼ iw

iZ θPR;n i0 A Y n;2 [ Y n;3 and θOP;n;w r io θOP;n;w þ 1

The first equation denotes that component l is not under maintenance; the second equation indicates that component l is preventively replaced; the last equation is the situation that component l is under imperfect opportunistic maintenance For each possible new state j of component l Step 2: Identify the corresponding state of subsystem n,X 0n ¼ j0 , after component l is changed to the new state j. Then, set 8
S
  < λn ðθn Þ i0 ;j0 þ μn i;j i0 A Y n;2 [ Y n;3 and iZ θOP;n;1
S   , λn ðθn Þ i0 ;j0 ¼
S : λn ðθn Þ i0 ;j0 þ ðλn Þi;j else
 where, ðλn Þi;j is the transition intensity of component l from state i to j, and μn i;j is the repair rate of component l from state i to state j End End

The setup cost of subsystem n occurs only when the state of the subsystem changes from a state in set Y n;1 to a state in set Y n;3 . Consequently, the average setup cost per unit time of subsystem n under strategy θn can be calculated by 0 1 X X
 S @ðπn ðθn ÞÞi0  ð4Þ C ST;n ðθn Þ ¼ C ST;n  λn ðθn Þ i0 ;j0 A: i0 A Y n;1

j0 A Y n;3

The maintenance activity may occur under the following two situations: (1) subsystem state changes from a state in set Y n;1 to another in set Y n;3 ; (2) a component crosses the threshold of a certain maintenance activity when the state of the subsystem is in set Y n;2 or Y n;3 . Thus, the average cost incurred by various maintenance activities can be calculated by: C R;n ðθn Þ ¼

X i0 A Y n;1

þ

0 @ðπn ðθn ÞÞi0 

X i0 A Y n;2 [ Y n;3

X j0 A Y n;3

0

KX n

@ðπn ðθn ÞÞi0  1




λSn ðθn Þ

KX n

 0 0

i ;j

X






!1 
 N n i; j0  C n ði; θn Þ A

i¼1




  I N n i; j0 ¼ N n i; i0 þ 1

i ¼ 1 j0 A Y n;2 [ Y n;3


  λSn i0 ;j0  C n ði; θn ÞA:

ð5Þ

Here, I ð U Þ is the indicator function given by  0 A is false I ðAÞ ¼ ; ð6Þ 1 A is true
 and the function N n i; j0 denotes the number of components in state i, when subsystem n is in state j0 . The notation C n ði; θn Þ is the cost of the maintenance activity of a component in subsystem n when the state of the component is i and the strategy adopted is θn , which is given by 8 0 i o θn;OP;1 > < i Z θn;PR C n ði; θn Þ ¼ ðcn Þi;1 : ð7Þ > : ðc Þ θn;OP;w r i oθn;OP;w þ 1 n i;i  w Using Eqs. (4) and (5), the average maintenance cost of subsystem n can be calculated as: C M;n ðθn Þ ¼ C ST;n ðθn Þ þ C R;n ðθn Þ:

ð8Þ

3.3. Calculation of the average production revenue The system-level production rate under a maintenance strategy θ is a discrete random variable whose distribution can be

represented using the UGF as U S ðz; θÞ ¼

N PRO X

pv ðθÞ  zgv ;

ð9Þ

v¼1

where gv is the vth level of the system production rate, and g 1 o ⋯ o g NPRO ; NPRO is the number of levels of the system production rate. The variable pv ðθÞ is the probability that the system production rate is gv when the system-level maintenance strategy is θ. The UGF is the z-transform of the moment generating function, which is widely used in the performance evaluation of complex systems [36]. The notation z is the independent variable in the z domain of the z-transform. The UGF of the system is the combination of the UGFs of its subsystems. The UGF of subsystem n is given by XNPRO;n un ðz; θn Þ ¼ p ðθ Þ  zgn;v : ð10Þ v ¼ 1 n;v n Here, the variable gn,v is the vth level of the production rate of subsystem n and g n;1 o …o g n;NPRO;n ; NPRO,n is the number of possible production rate of subsystem n. The variable pn;v ðθn Þ denotes the probability that the production rate of subsystem n is gn,v when the maintenance strategy of the subsystem is θn. The probability, pn;v ðθn Þ, can be calculated based on the steady-state distribution of subsystem n under strategy θn. The UGFs of the subsystems are then combined to form the UGF of the system by applying the composition operator  as:   U s ðz; θÞ ¼  u1 ðz; θ1 Þ; …; uN ðz; θN Þ (N ) NX PRO;1 PRO;N X p1;v ðθ1 Þ  zg1;v ; …; pN;v ðθN Þ  zgN;v ¼  v¼1

¼

NX PRO;1 v1 ¼ 1

⋯

NX PRO;N vN ¼ 1

v¼1






∏ pn;vn ðθn Þ  zϕ g1;v1 ;…;gN;vN : N

ð11Þ

n¼1

This polynomial formulation considers all the possible mutually exclusive combinations of subsystem-level production rates. The combination function ϕð U Þ computes the system production rate using the production rates of its subsystems according to the structure of the system. Because the subsystems are interconnected in series, the function ϕð U Þ can be written as:
   φ g 1;v1 ; …; g N;vN ¼ min g 1;v1 ; …; g N;vN : ð12Þ The average production revenue in Eq. (2) can now be calculated using the distribution of the system production rate as: XNPRO RP ðθÞ ¼ r p  p ðθÞ  g v : ð13Þ v¼1 v

Y. Zhou et al. / Reliability Engineering and System Safety 138 (2015) 40–53

can be estimated by the following three steps:

4. Strategy space reduction of multistate systems This section presents a new method using linear programming for the strategy space reduction of multistate systems. A theoretical proof is also provided in the section to illustrate that the new method is at least as effective as the existing methods in strategy space reduction.

Step 1: Calculate the upper and lower bounds of the cumulative maxÞ minÞ probability of subsystem-level production rates, F ðm;v and F ðm;v
m ¼ 1; 2; …; n 1; n þ1; …; N; v ¼ 1; 2; …; N PRO;m using maxÞ F ðm;v ¼ minθm A Θm F m;v ðθm Þ

The method developed in this study for strategy space reduction is based on the following consideration - a change of the maintenance strategy for a subsystem can cause variation in the average revenue of the system. When the upper bound of this variation is below zero, the new maintenance strategy for the subsystem should be removed from the search space of maintenance optimisation. The upper bound of the average revenue difference when the maintenance strategy of subsystem n changes from θ0n to θn is described by
 

 UBR θn ; θ0n Z R θn ; θn  R θ0n ; θn for 8 θn A Θn ; ð14Þ

0  where the notations R θn ; θn and R θn ; θn represent the average revenue per unit time when the maintenance strategy of subsystem n is θn and θ0n , respectively. Their difference,


 R θn ; θn  R θ0n ; θn , is never larger than UBR θn ; θ0n for 8 θn A Θn . Here, Θn is the maintenance strategy space of the system excluding subsystem n. It is obvious that strategy θn is not a part of the optimal system-level maintenance strategy and should be
 removed when UBR θn ; θ0n r 0. Since subsystem n connects to the rest of the system in series, the average revenue difference can be written as:


 R θn ; θn  R θ0n ; θn ¼ C M;n θ0n C M;n ðθn Þ þ r p ! N PRO;n PRO;n X
 NX


   :  pn;v θn  pn;u ðθn Þ  pn;u θ0n  min g n;v ; g n;u

and

maxÞ maxÞ maxÞ ¼ F ðm;v  F ðm;v pðm;v 1

and

v ¼ 1; …; NPRO;n :

ð17Þ

In addition, the upper and lower bounds of the production rate ðmaxÞ minÞ distributions in stochastic ordering (i.e., pn;v and pðn;v , v¼ 1 ; 2; …; NPRO;n ) which satisfy Xu Xu
 pðmaxÞ r p θn ð18Þ v ¼ 1 n;v v ¼ 1 n;v and Xu v¼1

ðminÞ pn;v Z

Xu v¼1


 pn;v θn u ¼ 1; …; N PRO;n

for

8 θn A Θn ð19Þ

ð22Þ

and minÞ minÞ minÞ ¼ F ðm;v F ðm;v pðm;v  1:

ð23Þ
 and v ¼ 1; 2; …; N PRO;n according to Step 3: Calculate the following lemma (see Appendix for the proof of the lemma): ðmaxÞ pn;v

minÞ pðn;v

maxÞ and Lemma 1. The bound distributions of the production rates pðn;v  ðminÞ
pn;v v ¼ 1; 2; …; N PRO;n that satisfy the inequalities Xu Xu
 pðmaxÞ r p θn ð24Þ v ¼ 1 n;v v ¼ 1 n;v

and Xu

ðminÞ pn;v Z

v¼1

Xu v¼1


 pn;v θn ; u ¼ 1; 2; …; N PRO;n

for

8 θn A Θn : ð25Þ

can be derived by combining the bound distributions of the subsystem-level production rates as maxÞ pðn;v ¼

u¼1


 The variable pn;v θn is the probability that the system excluding subsystem n has a production rate g n;v under maintenance strategy θn . Eq.
(15) can be regarded as a linear function of the variables pn;v θn , v ¼ 1; …; N PRO;n , when θn and θ0n are fixed. Therefore, the maximum value of Eq. (15) for the two given strategies θn and θ0n can be obtained using linear programming. Several constraints are
 applied to the variablespn;v θn , v ¼ 1; …; N PRO;n during linear programming. It can be derived according to the properties of discrete random variables that XNPRO;n

 ð16Þ pn;v θn ¼ 1 v¼1

ð21Þ 




Pv

where the variable F m;v ðθm Þ ¼ u ¼ 1 pm;u ðθm Þ v ¼ 1; 2; …; N PRO;m is the cumulative probability of the production rate of subsystem m under strategy θm.
 maxÞ minÞ Step 2: Convert F ðm;v and F ðm;v v ¼ 1; 2; …; N PRO;m into the ðmaxÞ ðminÞ bounds of probability distribution pm;v and pm;v using

NX PRO;1

⋯

v1 ¼ 1

N PRO;n Xþ 1 X 1 NPRO;n

⋯

vn  1 ¼ 1 vn þ 1 ¼ 1

NX PRO;N vN ¼ 1

N

∏

m ¼ 1; m a n

pðmaxÞ m;vm

  I min g 1;v1 ; …; g n  1;vn  1 ; g n þ 1;vn þ 1 ; …; g N;vN ¼ g n;v

ð15Þ


 0 r pn;v θn r1

ð20Þ

minÞ ¼ maxθm A Θm F m;v ðθm Þ; F ðm;v

4.1. Strategy space reduction using linear programming

v¼1

45

!

and minÞ pðn;v ¼

NX PRO;1 v1 ¼ 1

⋯

N PRO;n Xþ 1 X 1 NPRO;n vn  1 ¼ 1 vn þ 1 ¼ 1

⋯

NX PRO;N vN ¼ 1

N

∏

m ¼ 1; m a n

pðminÞ m;vm

!   I min g 1;v1 ; …; g n  1;vn  1 ; g n þ 1;vn þ 1 ; …; g N;vN ¼ g n;v : Therefore, the linear programming problem has N PRO;n variables and 4N PRO;n þ 1 constraints. Here, Eq. (16) can be regarded as one single equality constraint; Eq. (17) can be represented by 2NPRO;n inequality constraints, while Eqs. (18) and (19) each has NPRO;n inequality constraints. 4.2. A comparison study The effectiveness of the new method on strategy space reduction of series–parallel systems is evaluated theoretically in this section. 4.2.1. The new method versus the method in Ref. [29] The new method calculates the average system-level revenue difference caused by the change of subsystem-level strategies, while the method in Ref. [29] aims to identify the dominated subsystem-level strategies incurring a high maintenance cost and

46

Y. Zhou et al. / Reliability Engineering and System Safety 138 (2015) 40–53

a low production rate in stochastic ordering. In another word, the new method considers the whole system, while the method in Ref. [29] focuses only on the current subsystem. The new method therefore, will be more effective than that in Ref. [29]. The following lemma is provided to support this claim.

The following lemma is developed to prove that the new method can remove at least the same number of strategies as the method in Ref. [30].
 Lemma 3. The upper bound of revenue difference UBR θn ; θ0n r 0, if

0 0 0 UBR θn ; θn ¼ C M;n θn  C M;n ðθn Þ

Lemma 2. The upper bound of
the averageu revenue difference u
P P UBR θn ; θ0n r 0, if C M;n ðθn Þ Z C M;n θ0n and pn;v ðθn Þ Z pn;
0 v¼1 v¼1 v θn , 8 u ¼ 1; 2; …; N PRO;n .

þ rp











UBR θn ; θ0n ¼ C M;n θ0n  C M;n ðθn Þ þ r p

N PRO;n

X

pnn;v 

v¼1




 pn;u θ0n





 min g n;v ; g n;u





NX PRO;n





u¼1

u¼1

X 

 ðmaxÞ   I pn;u ðθn Þ Z pn;u θ0n pn;v min g n;v ; g n;u

u¼1




¼ C M;n θ0n



v¼1

0

N PRO;n

 C M;n ðθn Þ þ r p  E@

0

v¼1

N PRO;n

 rp  E@

X

X

v¼1

X   
  pnn;v  min g n;v ; g n;u ¼ C M;n θ0n C M;n ðθn Þ þ r p

n

pn;v

v



NX PRO;n u¼1





  pn;u ðθn Þ pn;u θ0n  E min Gn θnn ; g n;u ;


 where the random variable Gn θnn is the production rate of the system excluding subsystem n, which follows the probability mass

  
 function Pr Gn θnn ¼ g n;v ¼ pnn;v
. The  term, min Gn θn ; g n;u g, is a non-decreasing function of Gn θn . Using the properties of stochastic ordering, the following inequalities can be obtained as:

1   min g n;v ; Gn ðθn Þ A

1 
 pnn;v min g n;v ; Gn θ0n A;

N PRO;n  X

    ðminÞ pn;v  min g n;v ; g n;u E min Gn θnn ; g n;u Z

where the random variable Gn ðθn Þ represents the production rate
 of subsystem n under strategy θn, and pn;v ðθn Þ ¼ Pr Gn ðθn Þ ¼ g n;v . PNPRO;n n Because the term, v ¼ 1 pn;v min g n;v ; Gn ðθn Þ is a
non-decreasing  Pu Pu 0 function of Gn ðθn Þ and v ¼ 1 pn;v ðθn Þ Z v ¼ 1 pn;v θn , 8 u ¼ 1; 2; …; N PRO;n , the following inequality can be derived using the properties of stochastic ordering as: 0 1 0 1 NPRO;n N PRO;n X X   
0  n n @ A @ rE E pn;v min g n;v ; Gn ðθn Þ pn;v min g n;v ; Gn θn A: v¼1

Proof:

u¼1

N PRO;n
 X   pn;u θ0n pnn;v min g n;v ; g n;u v¼1

v¼1



 ðminÞ   i þ I pn;u ðθn Þ o pn;u θ0n pn;v min g n;v ; g n;u r0:

NPRO;n

u¼1

NX PRO;n


 pn;u ðθn Þ  pn;u θ0n

NX PRO;n



 UBR θn ; θ0n ¼ C M;n θ0n  C M;n ðθn Þ þr p  pn;u ðθn Þ  pn;u θ0n

pn;u ðθn Þ

N PRO;n NX PRO;n X
   ¼ C M;n θ0n  C M;n ðθn Þ þ r p pn;u ðθn Þ  pnn;v min g n;v ; g n;u

 rp




N PRO;n

Proof. The linear programming algorithm in this study
proposed  can obtain the specified values of pn;v θn , v ¼ 1; …; N PRO;n that
 maximise the system revenue difference Rðθn Þ  R θ0n . These specified values are denoted by pnn;v , v ¼ 1; …; N PRO;n . The upper bound of the average revenue difference of a system can then be calculated by:

NX PRO;n h

v¼1


0 Because
 C M;n ðθn Þ Z C M;n θn , it can be concluded that 0 UBR θn ; θn r 0. Lemma 2 shows that the dominated strategies identified by the method in Ref. [29] can also be excluded by the new method. Therefore, the current method is at least as effective as the method in Ref. [29] for maintenance strategy reduction. 4.2.2. The new method versus the method in Ref. [30] The new method and the one presented in Ref. [30] both remove inappropriate strategies according to the bounds of the average revenue variation. Though in Ref. [30], the upper bound of the average system revenue difference when subsystem n is under strategies θn and θ0n is calculated by NX PRO;n



 ½ pn;u ðθn Þ  pn;u θθn UB0R θn ; θ0n ¼ C M;n θ0n C M;n ðθn Þ þ r p  u¼1

X 

   ðmaxÞ I pn;u ðθn Þ Z pn;u θ0n  pn;v  min g n;v ; g n;u  N PRO;n v¼1



   ðminÞ  min g n;v ; g n;u Þ: þ I pn;u ðθn Þ o pn;u θ0n  pn;v ð26Þ

v

and N PRO;n  X

    ðmaxÞ pn;v  min g n;v ; g n;u : E min Gn θnn ; g n;u r v

Therefore,

 UBR θn ; θ0n ¼ C M;n θ0n  C M;n ðθn Þ þr p 

NX PRO;n u¼1







  pn;u ðθn Þ  pn;u θ0n  E min Gn θnn ; g n;u

NX PRO;n


 ½ pn;u ðθn Þ pn;u θ0n r C M;n θ0n  C M;n ðθn Þ þr p  u¼1

X 

   ðmaxÞ  I pn;u ðθn Þ Z pn;u θ0n  pn;v  min g n;v ; g n;u N PRO;n v¼1



   ðminÞ þ I pn;u ðθn Þ o pn;u θ0n  pn;v  min g n;v ; g n;u  r0:

5. Numerical study The performance of the new method is evaluated in this section by comparing the strategy reduction outcomes with those obtained using other strategy reduction methods and the ACO algorithm. 5.1. Description of the studied system The system used in this numerical simulation is shown in Fig. 1, which has six subsystems inter-connected in series. The state transition intensities of components in different subsystems are listed in Table 2. The repair rates and costs for different types of maintenance activities are shown in Tables 3 and 4, respectively. The

Y. Zhou et al. / Reliability Engineering and System Safety 138 (2015) 40–53

47

Fig. 1. The system structure. Table 2 The transition intensities of components in different subsystems. n

(λn)1,2

(λn)1,3

(λn)1,4

(λn)1,5

(λn)1,6

(λn)2,3

(λn)2,4

(λn)2,5

(λn)2,6

(λn)3,4

(λn)3,5

(λn)3,6

(λn)4,5

(λn)4,6

(λn)5,6

1 2 3 4 5 6

0.2 0.3 0.25 0.4 0.3 0.2

0.2 0.24 0.2 0.3 0.2 0.15

0.12 0.12 0.1 0.2 0.1 0.1

0.2 0.1 – 0.1 – 0.04

– – – 0.05 – 0.02

0.4 0.32 0.4 0.35 0.4 0.3

0.34 0.3 0.2 0.3 0.2 0.25

0.3 0.15 – 0.2 – 0.1

– – – 0.1 – 0.05

0.42 0.4 0.3 0.4 0.5 0.4

0.3 0.3 – 0.35 – 0.3

– – – 0.2 – 0.15

0.45 0.6 – 0.5 – 0.4

– – – 0.3 – 0.3

– – – 0.6 – 0.5

Table 3 The repair rates of components in different subsystems. n

(μn)2,1

(μn)3,2

(μn)3,1

(μn)4,3

(μn)4,2

(μn)4,1

(μn)5,4

(μn)5,3

(μn)5,2

(μn)5,1

(μn)6,1

1 2 3 4 5 6

4.5 5 4 5.5 5 4.9

3.5 4.5 3.8 4.5 4.8 4.5

3 4 3 4 3.8 3.5

3 4 – 4.3 – 4

2.5 3.4 – 3.5 – 3.2

2 3 2 3 2 3

– – – 4 – 3.5

– – – 3.5 – 3

– – – 3 – 2.5

1 2 – 2.5 – 2

– – – 1 – 1

Table 4 The maintenance costs of components in different subsystems. n

(cn)2,1

(cn)3,2

(cn)3,1

(cn)4,3

(cn)4,2

(cn)4,1

(cn)5,4

(cn)5,3

(cn)5,2

(cn)5,1

(cn)6,1

1 2 3 4 5 6

0.15 0.25 0.15 0.3 0.2 0.4

0.075 0.125 0.075 0.15 0.2 0.15

0.15 0.25 0.15 0.3 0.1 0.4

0.075 0.125 – 0.15 – 0.15

0.1 0.2 – 0.2 – 0.2

0.15 0.25 0.5 0.3 0.6 0.4

– – – 0.15 – 0.12

– – – 0.22 – 0.2

– – – 0.25 – 0.25

0.6 1 – 0.3 – 0.4

– – – 1.5 – 1.2

Table 5 The production rates of components in different subsystems.

5.2. Performance evaluation of the strategy reduction method using linear programming

n

γn,1

γn,2

γn,3

γn,4

γn,5

γn,6

1 2 3 4 5 6

0.4 0.5 0.3 0.8 0.4 0.5

0.3 0.35 0.2 0.6 0.3 0.45

0.2 0.2 0.1 0.4 0.1 0.35

0.1 0.1 0 0.3 0 0.2

0 0 – 0.15 – 0.1

– – – 0 – 0

units used for the time and cost are day and 103 $. The production rates of components in different subsystems are listed in Table 5, and the reward per unit system-level production rate is assumed to be rp ¼10  103 $. According to the maintenance strategy structure considered in this study, the size of the system-level strategy space is 22  22  8  64  8  64¼126,877,696, which is too large to be processed by the enumeration method.

In this section, the performance of the new method is evaluated by comparing the results with those obtained using other optimisation algorithms. A sensitivity analysis is also carried out to investigate the relationship between the effectiveness of the new algorithm and the values of system parameters.

5.2.1. Maintenance optimisation using the ACO algorithm For the use of the ACO algorithm, the system-level maintenance strategy θ is divided into six steps, i.e., θn, n ¼1,…,6. The ants then randomly choose a subsystem-level maintenance strategy at each step. The probability that an ant selects the sth strategy for subsystem n is:

 α τn;s þ β ηn;s
; pðn; sÞ ¼ P

 t α τn;t þ β ηn;t

ð27Þ

48

Y. Zhou et al. / Reliability Engineering and System Safety 138 (2015) 40–53

strategy is θn1 ¼ [4,4,3,2], θn2 ¼[3,3,3,3], θn3 ¼[3,3,2], θn4 ¼[3,3,3,3,3], θn5 ¼[3,3,3], and θn6 ¼[4,4,4,3,3]. The optimisation process is executed using MATLAB 7.14 on a desktop computer with an Intel i7 3770 CPU and an eight Gigabytes of RAM. The elapsed time of the optimisation process in this simulation is 201 s.

where τn;s quantifies the pheromone trail at the sth strategy for subsystem n, and ηn;s is the visibility of that strategy. The parameters, α and β, represent the importance of the pheromone trail and the visibility. The visibility of a strategy is defined as the difference between the average subsystem-level production revenue and maintenance cost: ηn;s ¼ r p 

NX PRO;n

  pn;v θðnsÞ  g n;v  C M;n θðnsÞ ;

5.2.2. Maintenance optimisation using the method based on linear programming The new method of strategy space reduction using linear programming is also deployed for the maintenance optimisation of the same system. The results are listed in the first row of Table 6. At first, the new method is applied directly to all the subsystems where six compact subsystem-level strategy spaces are obtained as Θn , n ¼1,2,…,6. Thus, the size of the strategy space is reduced significantly from 126,877,696 to ∏6n ¼ 1 jΘn j¼ 4. After that, maintenance strategies of two adjacent subsystems are combined and the new method of the three new strategy sets

is applied; the sizes

are Θ1;2 ¼ 2, Θ3;4 ¼2, and Θ5;6 ¼1. Then, the strategy space Θ1;2 and Θ3;4 are combined again to become Θ1;2;3;4 that contains only two possible strategies. Finally, Θ1;2;3;4 is further combined with Θ5;6 to form Θ1;2;3;4;5;6 , so that only two candidate strategies are left. During each step, the dominated subsystem-level strategies are removed using the method in Ref. [29] before the new method based on linear programming is applied. Furthermore, the new method is performed twice in each step of strategy reduction to improve the estimation of the bound of average revenue difference, as inappropriate strategies are removed. It is noted that the optimal strategy obtained by the new method is the same as that identified by the ACO algorithm. Nevertheless, the total computation time for the strategy reduction and optimal strategy searching is just over 1 s, which is much faster than the ACO algorithm.

ð28Þ

v¼1

θðnsÞ

where denotes the sth optional maintenance strategy for subsystem n. Once an ant has chosen a subsystem-level strategy for each step, a system-level maintenance strategy is generated. After all the ants have identified their system-level strategy, the pheromone trail intensities are updated as: XNa τ0n;s ¼ ð1  ρÞ  τn;s þ Δτnn;sa ; ð29Þ n ¼1 a

where Na is the number of ants, and ρ(0 oρ r1) is the evaporation rate of the pheromone. The notation Δτnn;sa denotes the pheromone that is deposited on the sth strategy of subsystem n by ant na, which is calculated by: Δτnn;sa



8 < Q  R
θðna Þ   min
R
θðtÞ  t ¼ 1;…;N ¼ a : 0

if edge ðn; sÞ is visited by ant na

;

otherwise

ð30Þ ðna Þ

is the system-level strategy where Q is a constant, and θ selected by ant na. Based on the trial and error test results, the parameters of the ACO algorithm are selected as: Q¼ 100,000, ρ¼ 0.2, α¼0.9, β¼ 0.1, and Na ¼22 þ22 þ8 þ64 þ 8þ 64 ¼188. At the end of this process, most artificial ants concentrate on the maintenance strategies that generate high revenue. The
maximum average revenue obtained by the ACO algorithm is R θn ¼ 3133.7 $, while the corresponding

5.2.3. Maintenance optimisation using other strategy space reduction methods The strategy reduction methods developed in Refs. [29,30] are also applied to this maintenance optimisation problem for comparison. It is shown in Table 6 that the method of Ref. [29] can reduce the size of the strategy space to 380,916. The optimal solution is then identified using the enumeration method. The computation time used in this exercise is 373 s. Table 6 further shows that the method of Ref. [30] can reduce the number of strategies to 42 and the total computation time is 24 s. The numerical simulation result confirms that the new method using linear programming provides a better and more efficient approach to the maintenance optimisation of series–parallel systems.

Table 6 Results of strategy reduction using different methods.









∏6i ¼ 1 jΘi j Θ1;2 Θ3;4 Θ5;6 Θ1;2;3;4 Θ1;2;3;4;5;6 Elapsed time (s) The proposed 4 method Method in 2646,000 [29] Method in 1815,450 [30]

2

2

1

2

2

1

109

72

108

3527

380,916

373

84

34

82

7

42

24

Upper bound of average revenue difference

3

Upper bound obtained by linear programming Upper bound obtained by the method in [30]

2.5 2 1.5 1 0.5 0

20

40

60

80

100

120

140

160

180

200

The indices of strategy pairs

Fig. 2. Upper bound of average revenue difference identified by the new method and the method in Ref. [30].

Y. Zhou et al. / Reliability Engineering and System Safety 138 (2015) 40–53

Because the new method and the method in Ref. [30] both used the upper bound of the average revenue difference as the strategy reduction criterion, the upper bounds identified by these two algorithms are compared in this study. The strategy reduction process of subsystem 1 is discussed in the following. After removing the dominated strategies in stochastic ordering, the number of potential strategies for subsystem 1 is 15. The upper bounds of the average revenue difference between 15  14 ¼210 pairs of strategies are then calculated, which are plotted in Fig. 2. Fig. 2 confirms that linear programming can always produce a smaller bounding value than the method in Ref. [30], which is consistent with Lemma 3. All the 210 bounding values derived by the method in Ref. [30] are larger than zero, which indicates that no strategy can be further excluded. On the other hand, 56 bounding values calculated by using linear programming are less than zero, and 11 maintenance strategies are further excluded.

Table 7 The strategy reduction results under different values of reward per system-level production unit, rp. rp 0.1 1 10 50 100 1,000 10,000

jΘ1 j

j Θ2 j

j Θ3 j

j Θ4 j

jΘ5 j

j Θ6 j



Θ1;2



Θ3;4



Θ5;6



Θ1;2;3;4;5;6

1 1 1 2 2 2 2

1 2 2 2 2 3 3

1 1 2 2 2 2 2

1 1 1 2 2 2 2

1 1 1 2 3 3 3

1 1 1 1 1 2 2

1 2 2 3 2 2 2

1 1 2 3 2 2 2

1 1 1 2 2 2 2

1 2 2 4 4 4 4

Table 8 The strategy reduction results when the production rates of subsystem four is changed. γ 04;i =γ 4;i

j Θ1 j

j Θ2 j

j Θ3 j

j Θ4 j

jΘ5 j

j Θ6 j



Θ1;2



Θ3;4



Θ5;6



Θ1;2;3;4;5;6

0.1 0.2 0.5 1 2 5 10 100

1 1 1 1 2 2 2 2

1 2 1 2 1 1 1 1

1 1 1 2 2 2 2 2

1 1 2 1 2 1 1 1

1 1 1 1 1 1 1 1

1 1 2 1 2 2 2 2

1 2 1 2 2 2 2 2

1 1 2 2 2 2 2 2

1 1 2 1 2 2 2 2

1 2 4 2 2 4 4 4

49

5.2.4. Sensitivity analysis The reward of unit system-level production rate, rp, is one of the critical parameters affecting the effectiveness of the new method developed in this study. The number of remaining strategies after reduction under different values of rp is listed in Table 7 for comparison. It is shown that only one strategy remains when rp is at the low end of the value range used in this simulation. This is because when the rp value is small, the factor dominating the maintenance optimisation is the maintenance cost. Thus, only the subsystem-level strategies incurring the lowest maintenance cost can survive after strategy reduction. On the other hand, when the value of rp is large, the maintenance cost can be ignored. Therefore, the objective function of maintenance optimisation becomes the average system-level production rate. In this case, a further increase of the rp value will not change the strategy reduction result. Table 7 shows that the strategy reduction method using linear programming is effective for all cases using different values of rp. Another important factor affecting the strategy reduction of series–parallel systems is the uneven production rates among subsystems. In this analysis, the production rate of subsystem 4 is changed from γ 4;i to γ 04;i . The strategy reduction results for different ratio of γ 04;i to γ 4;i are listed in Table 8. It is shown that when the ratio is small, say 0.1, only one candidate strategy remains. In this case, subsystem 4 is the bottleneck of the system, and the only candidate strategy for all the other subsystems is the one incurring the lowest maintenance cost. On the other hand, when the value of γ 04;i is large, changing the maintenance strategy of subsystem 4 does not affect the system-level production rate significantly. Therefore, the strategy incurring the lowest maintenance cost is selected for subsystem 4. In this case, a further increase of γ 04;i will not change the result of the strategy reduction. Table 8 shows that the effect of unbalanced production rates among subsystems on the strategy reduction is insignificant. The system strategy space increases rapidly when the number of subsystems grows. Consequently, the effectiveness of the new method on systems comprising different number of subsystems is compared with those of the ACO algorithm and the method in Ref. [30]. The number of subsystems is increased by multiplying the base system described in Section 5.1, which connects to each other in series. To ensure that the ACO algorithm can stop in a reasonable time duration, the maximum number of ants used in the simulation is set at 200, and the maximum number of iterations is set at 1000. The other parameters of the ACO

Table 9 Computation time (s) of optimisation using the new method and that of ACO algorithm for different number of subsystems. N

6 12 18 24 30 36

rp ¼ 10

rp ¼20

rp ¼40

rp ¼60

rp ¼ 80

rp ¼ 100

rp ¼ 120

LP

ACO

LP

ACO

LP

ACO

LP

ACO

LP

ACO

LP

ACO

LP

ACO

1 4 29 27 14 27

201 610 1174 1754 2368 3475

1 5 31 71 25 29

189 793 1280 2016 3764 4393

2 4 13 111 60 28

176 1058 2083 2667 4001 4352

2 5 13 103 49 32

158 767 2591 3100 3573 7328

2 4 14 133 50 32

155 685 1457 1718 3920 5549

2 5 11 112 54 36

158 898 1432 2573 5776 7167

2 5 8 227 51 33

192 714 1310 2427 5714 6762

LP denotes the new strategy reduction method using linear programming.

Table 10 Computation time (s) of optimisation using the method in Ref. [30] for different number of subsystems. N

rp ¼ 10

rp ¼20

rp ¼ 40

rp ¼ 60

rp ¼ 80

rp ¼ 100

rp ¼120

6 12

24 436,000

10 436,000

8 436,000

7 436,000

6 4 36,000

6 436,000

6 436,000

50

Y. Zhou et al. / Reliability Engineering and System Safety 138 (2015) 40–53

Table 11 Comparison of the optimal average revenue (103 $) obtained by the new method and the ACO algorithm for different number of subsystems. N

6 12 18 24 30 36

rp ¼10

rp ¼ 20

rp ¼ 40

rp ¼ 60

rp ¼80

rp ¼ 100

rp ¼ 120

LP

ACO

LP

ACO

LP

ACO

LP

ACO

LP

ACO

LP

ACO

LP

ACO

3.134  2.459  7.128  11.35  15.34  19.23

3.134  2.493  7.128  11.35  15.42  19.23

10.60 3.166  2.654  7.848  12.59  16.98

10.60 3.166  2.654  7.852  12.63  17.04

25.80 15.04 7.204 0.4819  5.520  11.07

25.80 15.04 7.204 0.2830  5.556  11.33

41.19 27.11 17.35 9.288 2.111  4.419

41.19 26.74 17.35 9.275 1.496  4.761

56.59 39.21 27.68 18.18 9.986 2.528

56.59 39.21 27.56 18.18 9.974 2.477

71.99 51.38 38.02 27.12 17.92 9.645

71.99 51.20 37.82 26.88 17.91 9.110

87.39 63.59 48.35 36.06 25.88 16.79

87.39 63.07 48.35 35.97 25.88 16.76

LP denotes the new strategy reduction method using linear programming. When the ACO fails to converge to the global optimal solution, the cell of the table is filled with light grey.

Fig. 3. The system structure of the case study in Ref. [9].

algorithm are the same as those used in Section 5.2.1. The computation times of the new algorithm and the ACO algorithm for different number of subsystems and rp values are listed in Table 9. It is shown that the new method is much more efficient than the ACO algorithm for all the cases. E.g., the new method can still identify the optimal maintenance strategy in less than 50 s when there are 36 subsystems, or by another word, 126,877,6966 E4.1717  1048 candidate strategies. In contrast, the ACO algorithm takes more than 1 h to solve the same problem. Another interesting observation is that the computation time of the new method does not always increase with the number of subsystems, as the computation time also depends on the number of strategies removed at each step. The method in Ref. [30], on the other hand, is sensitive to the size of the system. As shown in Table 10, the computation time of the method in Ref. [30] is reasonable when there is a small number of subsystems. However, when the number of subsystems increases to 12, it takes more than 10 h to obtain the optimal maintenance strategy. Table 11 compares the optimal average revenue values obtained by the new method and the ACO algorithm. It is shown that the ACO algorithm can only find the global optimal solution for a system with small number of subsystems. When the number of subsystems is large (i.e., more than 24), the ACO algorithm is most likely to converge to a local optimum. Although the solution of the ACO can be improved by adjusting its parameters, it can take longer time to converge. 5.3. Application of the strategy reduction methods to the case studies of the existing publications [9,29] Two additional case studies extracted from the existing papers are also simulated in this study to illustrate the applicability and effectiveness of the new method. 5.3.1. The case study from Ref. [9] The system investigated in Ref. [9] is also used in this numerical simulation. The structure of the system is shown in Fig. 3. It is

Table 12 Results of strategy reduction for the case study in Ref. [9] using different methods.







∏5i ¼ 1 jΘi j Θ1;2 Θ3;4 Θ1;2;3;4 Θ1;2;3;4;5 Elapsed time (s) Method in [29] Method in [30] The proposed method

3020,544 3020,544 16

94 93 2

95 28 2

3599 8 2

662,216 1,256 1

774 103 4

Table 13 Computation time (s) of optimisation using the current method and that of ACO algorithm for different number of subsystems of the case study in Ref. [9]. Number of subsystems

5 10 15 20 25 30

ct ¼1000

ct ¼4000

ct ¼ 7000

ct ¼ 10,000

LP

ACO

LP

ACO

LP

ACO

LP

ACO

4 9 27 73 167 230

118 553 873 1722 4728 4632

4 10 16 51 104 131

165 748 1209 2513 2552 3064

4 12 18 41 129 152

271 1001 1911 2737 4762 6908

4 12 19 44 116 165

219 784 3026 3411 2770 6920

LP denotes the new strategy reduction method using linear programming.

Table 14 Computation time (s) of optimisation using the method in Ref. [30] for different number of subsystems of the case study in Ref. [9]. Number of subsystems

ct ¼1,000

ct ¼4,000

ct ¼ 7,000

ct ¼ 10,000

5 10

103 436,000

116 436,000

112 436,000

114 4 36,000

noted that only the case where the refund is not involved is considered in this study. Based on the assumption made in Ref. [9], each component only has a preventive replace threshold,

Y. Zhou et al. / Reliability Engineering and System Safety 138 (2015) 40–53

51

Table 15 Comparison of the optimal average revenue (103 $) obtained by the new method and the ACO algorithm for different number of subsystems of the case study in Ref. [9]. ct ¼1000

Number of subsystems

5 10 15 20 25 30

ct ¼4000

ct ¼ 7000

ct ¼ 10,000

LP

ACO

LP

ACO

LP

ACO

LP

ACO

53.296  11.064  59.787  104.02  146.32  186.60

53.074  11.406  60.467  105.19  146.54  187.19

353.72 197.81 101.71 27.529  34.342  88.615

353.72 197.25 101.03 25.964  35.381  90.909

656.37 416.70 275.15 170.48 85.463 12.242

656.27 416.60 274.23 170.26 83.398 9.6936

959.38 636.86 450.81 317.34 209.24 120.40

959.38 636.72 450.58 314.72 208.31 117.20

LP denotes the new strategy reduction method using linear programming. When the ACO fails to converge to the global optimal solution, the cell of the table is filled with light grey.

Table 16 Results of strategy reduction for the case study in Ref. [29] using the method in Ref. [30]. nso

nrb

∏4i ¼ 1 jΘi j



Θ1;2



Θ1;2;3



Θ1;2;3  jΘ4 j

Elapsed time (s)

1 1 1

1 2 3

31,680 7,560 3,744

13 7 7

5 3 3

25 15 15

1.4 1.3 1.4

Table 17 Results of strategy reduction for the case study in Ref. [29] using the new method based on linear programming. nso

nlp

∏4i ¼ 1 jΘi j



Θ1;2



Θ1;2;3



Θ1;2;3  jΘ4 j

Elapsed time (s)

1 1 1

1 2 3

576 16 16

2 2 2

2 2 2

8 8 8

2.9 2.6 2.6

5.3.2. The case study from Ref. [29] To further illustrate the effectiveness of the new method, the case study from Ref. [29] is also simulated. The number of candidate strategies was effectively reduced from 60,840,000 to 45,009 in Ref. [29], and was further reduced to 15 in Ref. [30] (see Table 16). In Table 16, nso and nrb denote the implementation numbers of the methods in Refs. [29,30], respectively. Table 17 shows the result of the same case study using the new method, where the notation nlp denotes the implementation number of the new method in each step. For all the three situations considered, the new method reduces the number of strategies from 60,840,000 to 8, which is a better outcome than that of Ref. [30]. Although the reduction method using linear programming takes slightly longer time to identify the upper bound of the average revenue difference than the method of Ref. [30], the overall computation time is similar to that of Ref. [30]. This is because most subsystem-level strategies are excluded by the new method, and linear programming is only invoked for a small number of times.

6. Conclusions and the number of possible system-level strategies is 27  4  27  9  256 ¼6718,464. The simulation results of the three strategy reduction methods are listed in Table 12 for comparison. It is shown that by using the method in Ref. [29], the number of possible system-level maintenance strategies can be reduced from 6718,464 to 662,216, and the computation time is more than 770 s. The method in Ref. [30] is then applied to the same system. It takes 103 s to identify 1256 candidates and find the optimal strategy within these candidates. On the contrary, the new method takes less than 5 s to find the optimal system-level strategy. Sensitivity analysis about the relationship between the number of subsystems and the efficiency of optimisation algorithms is also performed. The investigated system is expanded by connecting different number of original systems in Fig. 3 in series. The ACO algorithm, the method in Ref. [30], and the new method are all adopted. The computation time and the optimisation results are respectively listed in Tables 13–15. In these tables, ct is the reward for unit production rate of the system. Table 13 shows that the new method can obtain the optimal strategy much more efficiently than the ACO algorithm. As shown in Table 14, although the method in Ref. [30] is more efficient than the ACO algorithm when the number of subsystems is five, it cannot obtain the optimal strategy within 10 h when the number of subsystems is more than 10. Table 15 reveals that the ACO algorithm can hardly find the global optimal solution when there are more than five subsystems.

This paper has presented a new method to reduce strategy space using linear programming for efficient maintenance optimisation of series–parallel systems with multistate components. A theoretical proof has been provided in the paper to illustrate that the method is at least as effective as the existing methods in strategy space reduction of series–parallel systems. The method has also been successfully applied to the maintenance optimisation problem of a series–parallel system considering both economic dependence among components and imperfect maintenance activities with different effects. The numerical examples presented in the study confirmed that the method can obtain the optimal maintenance strategy more efficiently than the existing strategy reduction methods and the ACO algorithm. A sensitivity analysis of the critical system parameters showed that the algorithm was equally effective when a large range of system parameters were used in the simulation. It is worthy to point out that the method developed in this paper can also be used to efficiently process a system with a large number of subsystems; this is in direct comparison with the long computation time taken when other strategy reduction methods are used.

Acknowledgements The research work is supported by the National Natural Science Foundation of China (Grant no. 71201025) and the Research Fund

52

Y. Zhou et al. / Reliability Engineering and System Safety 138 (2015) 40–53

for the Doctoral Program of Higher Education of China (Grant no. 20110092120007).

Appendix The proof of Lemma 1: Proof. Only Eq. (24) is proved as follows, because Eq. (25) can be proved in the same way. It can be derived that u X v¼1

ðmaxÞ pn;v ¼

u NX PRO;1 X v ¼ 1 v1 ¼ 1

⋯

N PRO;n X 1 NPRO;n Xþ 1

⋯

vn  1 ¼ 1 vn þ 1 ¼ 1

NX PRO;N vN ¼ 1

 ðmaxÞ ðmaxÞ ðmaxÞ pðmaxÞ 1;v1 …pn  1;vn  1 pn þ 1;vn þ 1 ⋯pN;vN I   min g 1;v1 ; …; g n  1;vn  1 ; g n þ 1;vn þ 1 ; …; g N;vN ¼ g n;v

¼

NX PRO;1

N PRO;n Xþ 1 X 1 NPRO;n

⋯

v1 ¼ 1

vn  1 ¼ 1 vn þ 1 ¼ 1

NX PRO;N

⋯

vN ¼ 1

 ðmaxÞ ðmaxÞ ðmaxÞ pðmaxÞ 1;v1 …pn  1;vn  1 pn þ 1;vn þ 1 ⋯pN;vN I   min g 1;v1 ; …; g n  1;vn  1 ; g n þ 1;vn þ 1 ; …; g N;vN r g n;u     Þ ¼ E I min G1 θð1maxÞ ; …; Gn  1 θðnmax 1 ;   Þ ðmax Þ Gn þ 1 θðnmax rg n;u þ 1 ; …; GN θN and u X




  pn;v θn ¼ Pr Gn θn r g n;v

v¼1

¼

NX PRO;1 v1 ¼ 1




⋯

N PRO;n Xþ 1 X 1 NPRO;n vn  1 ¼ 1 vn þ 1 ¼ 1

⋯

NX PRO;N vN ¼ 1

p1;v1 ðθ1 Þ…pn  1;vn  1 ðθn  1 Þ

pn þ 1;vn þ 1 ðθn þ 1 Þ⋯pN;vN ðθN Þ


I min g 1;v1 ; …; g n  1;vn  1 ; g n  þ 1; vn þ 1 ; …; g N;vN r g n;u ÞÞÞ ¼ EðI ðminðG1 ðθ1 Þ; …; Gn  1 ðθn  1 Þ; Gn þ 1 ðθn þ 1 Þ; …; GN ðθN ÞÞ r g n;u :
 Here, the random variableGm θðmmaxÞ which denotes the production rate of subsystem m follows bound distribution pðmaxÞ m;vm

the   ðmaxÞ 8 vm ¼1,…,NPRO,m, i.e., pðmaxÞ ¼ Pr G θ ¼ g , whereas the m m m;vm m;vm random variable Gm ðθm Þ is the production rate of subsystem m
 under strategy θm , i.e., pm;vm ðθm Þ ¼ Pr Gm ðθm Þ ¼ g m;vm . Because
 I minðG1 ðθ1 Þ; …; Gn  1 ðθn  1 Þ; Gn þ 1 ðθn þ 1 Þ; …; GN ðθN ÞÞ r g n;u is a non-increasing function of G1 ðθ1 Þ; …; Gn  1 ðθn  1 Þ; Gn þ 1 ðθn þ 1 Þ; …;
 Pu Pu ðmaxÞ 8 u ¼ 1; 2; …; N PRO;n , GN ðθN Þ and v ¼ 1 pn;v θn Z v ¼ 1 pn;v it can be obtained using the property of stochastic ordering      Þ ðmax Þ that E I min G1 θð1maxÞ ; …; Gn  1 θðnmax  1 ; Gn þ 1 θn þ 1 ; …;  GN θðNmaxÞ Þ r g n;u ÞÞ rEðI ðminðG1 ðθ1 Þ; …; Gn  1 ðθn  1 Þ; Gn þ 1 ðθn þ 1 Þ; …; GN ðθN ÞÞ r g n;u ÞÞ:
 P P ðmaxÞ Consequently, uv ¼ 1 pn;v r uu pn;v θn .

Appendix A. Supporting information Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.ress.2015.01.018.

References [1] Levitin G, Lisnianski A. Importance and sensitivity analysis of multi-state systems using the universal generating function method. Reliab Eng Syst Saf 1999;65:271–82. [2] Levitin G, Lisnianski A. Joint redundancy and maintenance optimization for multistate series–parallel systems. Reliab Eng Syst Saf 1999;64:33–42. [3] Nourelfath M, Ait-Kadi D. Optimization of series–parallel multi-state systems under maintenance policies. Reliab Eng Syst Saf 2007;92:1620–6. [4] Bris R, Châtelet E, Yalaoui F. New method to minimize the preventive maintenance cost of series–parallel systems. Reliab Eng Syst Saf 2003;82:247–55. [5] Samrout M, Yalaoui F, Châtelet E, Chebbo N. New methods to minimize the preventive maintenance cost of series–parallel systems using ant colony optimization. Reliab Eng Syst Saf 2005;89:346–54. [6] Aziz A, Jarrahi F, Abdul-Kader W. Modeling and performance evaluation of a series–parallel flow line system with finite buffers. INFOR 2010;48:103–20. [7] Massim Y, Yalaoui F, Chatelet E, Yalaoui A, Zeblah A. Efficient immune algorithm for optimal allocations in series–parallel continuous manufacturing systems. J Intell Manuf 2010:1–17. [8] Nicolai RP, Dekker R. Optimal maintenance of multi-component systems: a review. Complex system maintenance handbook. London: Springer; 2008. p. 263–86. [9] Liu Y, Huang H Optimization of multi-state elements replacement policy for multi-state systems. In: Reliability and maintainability symposium (RAMS) proceedings—annual; 2010. p. 1-7. [10] Huang C-H, Wang C-H. Optimization of preventive maintenance for a multistate degraded system by monitoring component performance. J Intell Manuf 2014:1–20. [11] Liu Y, Hong-Zhong H. Optimal selective maintenance strategy for multi-state systems under imperfect maintenance. IEEE Trans Reliab 2010;59:356–67. [12] Pandey M, Zuo MJ, Moghaddass R, Tiwari MK. Selective maintenance for binary systems under imperfect repair. Reliab Eng Syst Saf 2013;113:42–51. [13] Pandey M, Zuo MJ, Moghaddass R. Selective maintenance modeling for a multistate system with multistate components under imperfect maintenance. IIE Trans 2013;45:1221–34. [14] Vu HC, Do P, Barros A, Bérenguer C Maintenance planning and dynamic grouping for multi-component systems with positive and negative economic dependencies. IMA J Manage Math 2014. [15] Doostparast M, Kolahan F, Doostparast M. A reliability-based approach to optimize preventive maintenance scheduling for coherent systems. Reliab Eng Syst Saf 2014;126:98–106. [16] Rami A, Hamdaoui H, Sayah H, Zeblah A. Efficient harmony search optimization for preventive-maintenance-planning for nuclear power systems. Int J Simul Multisci Des Optim 2014;5:A17. [17] Peng R, Xie M, Ng SH, Levitin G. Element maintenance and allocation for linear consecutively connected systems. IIE Trans 2012;44:964–73. [18] Xiao H, Peng R. Optimal allocation and maintenance of multi-state elements in series–parallel systems with common bus performance sharing. Comput Ind Eng 2014;72:143–51. [19] Nourelfath M, Châtelet E, Nahas N. Joint redundancy and imperfect preventive maintenance optimization for series–parallel multi-state degraded systems. Reliab Eng Syst Saf 2012;103:51–60. [20] Moghaddass R, Zuo MJ, Pandey M. Optimal design and maintenance of a repairable multi-state system with standby components. J Stat Plan Infer 2012;142:2409–20. [21] Liu Y, Huang H, Wang Z, Li Y, Yang Y. A joint redundancy and imperfect maintenance strategy optimization for multi-state systems. IEEE Trans Reliab 2013;62:368–78. [22] Sheikhalishahi M, Ebrahimipour V, Farahani MH. An integrated GA-DEA algorithm for determining the most effective maintenance policy for a kout-of-n problem. J Intell Manuf 2014;25:1455–62. [23] Wang C-H, Tsai S-W. Optimizing bi-objective imperfect preventive maintenance model for series–parallel system using established hybrid genetic algorithm. J Intell Manuf 2014;25:603–16. [24] Moghaddam KS, Usher JS. A new multi-objective optimization model for preventive maintenance and replacement scheduling of multi-component systems. Eng Optimiz 2010;43:701–19. [25] Lin T-W, Wang C-H. A hybrid genetic algorithm to minimize the periodic preventive maintenance cost in a series–parallel system. J Intell Manuf 2012;23:1225–36. [26] Wang C-H, Lin T-W. Improved particle swarm optimization to minimize periodic preventive maintenance cost for series–parallel systems. Expert Syst Appl 2011;38:8963–9. [27] Cao D, Murat A, Chinnam RB. Efficient exact optimization of multi-objective redundancy allocation problems in series–parallel systems. Reliab Eng Syst Saf 2013;111:154–63. [28] Galante G, Passannanti G. An exact algorithm for preventive maintenance planning of series–parallel systems. Reliab Eng Syst Saf 2009;94:1517–25. [29] Zhou Y, Zhang Z, Lin TR, Ma L. Maintenance optimisation of a multi-state series–parallel system considering economic dependence and statedependent inspection intervals. Reliab Eng Syst Saf 2013;111:248–59. [30] Zhou Y, Zhang Z, Bian Y, Gao H. An efficient maintenance optimisation method for series–parallel systems using stochastic ordering and revenue difference boundaries. In: 2013 International conference on quality, reliability, risk, maintenance, and safety engineering (QR2MSE). Emeishan, Sichuan, China: IEEE; 2013. p. 658–63.

Y. Zhou et al. / Reliability Engineering and System Safety 138 (2015) 40–53

[31] Levitin G, Lisnianski A. Optimization of imperfect preventive maintenance for multi-state systems. Reliab Eng Syst Saf 2000;67:193–203. [32] Nahas N, Khatab A, Ait-Kadi D, Nourelfath M. Extended great deluge algorithm for the imperfect preventive maintenance optimization of multi-state systems. Reliab Eng Syst Saf 2008;93:1658–72. [33] Castanier B, Grall A, Bérenguer C. A condition-based maintenance policy with non-periodic inspections for a two-unit series system. Reliab Eng Syst Saf 2005;87:109–20.

53

[34] Gupta S, Bhattacharya J. Discrete Markov chains—an analytical tool for productivity analysis of surface mining systems. Int J Surf Min Reclam Environ 1999;13:111–6. [35] Pascual R, Godoy D, Louit DM. Throughput centered prioritization of machines in transfer lines. Reliab Eng Syst Saf 2011;96:1396–401. [36] Levitin G. Universal generating function in reliability analysis and optimization. London: Spring-Verlag London Limited; 2005.