Non-Linear Threshold Algorithm for the Redundancy Optimization of Multi-State Systems

To improve system performance, redundancy is widely used in different kinds of industrial applications such as power systems, aerospace, electronic, telecommunications and manufacturing systems. Designing high performant systems which meet customer requirements with a minimum cost is a challenging task in these industries. This paper develops an efficient approach for the redundancy optimization problem of series-parallel structures modeled as multi-state systems. To reach the target system availability, redundancies are used for components among a list of products available in the market. Each component is characterized by its own availability, cost and performance. The goal is to minimize the total cost under a system availability constraint. Discrete levels of performance are considered for the system and its components. The extreme values of such performance levels correspond to perfect functioning and complete failure. A piecewise cumulative load curve represents consumer demand. System availability corresponds to the aptitude to fulfill this demand. The multi-state system availability evaluation uses the universal moment generating function technique. The proposed optimization algorithm is based on the non-linear threshold accepting metaheuristic, while using a self-adjusting penalty guided strategy. The obtained results demonstrate the approach efficiency for solving the redundancy optimization problem of multi-state systems. Its effectiveness is also tested using the classical redundancy optimization problem of binary-state systems. The algorithm is evaluated by comparison to the best known methods. For multi-state systems, it is compared to genetic algorithm and tabu search. For binary-state systems, it is compared to genetic algorithm, tabu search, ant colony optimization and harmony search. The obtained results demonstrate that the proposed approach outperforms these state-of-the-art benchmark methods in finding, for all considered instances, a high-quality solution in a minimum computational time.


Introduction
To improve system reliability, redundancy is widely used in different kinds of industrial applications including aerospace, power systems, electronic, telecommunications and manufacturing systems. The redundancy optimization problem (ROP) of a series-parallel system considers many sub-systems in series and for each sub-system, different component versions are connected in parallel. Thus, two design objectives are generally considered. The first objective is to find the number of redundant components and their location to ensure the required performance subject to resource constraints. The second objective is to minimize the cost under performance constraints. The current reliability optimization literature distinguishes between the ROP of multistate systems (MSS), and the ROP of binary-state systems (BSS).
The ROP of MSS is more recent. Multi-state reliability modeling considers that a system and its components may have more than two levels of performance varying from perfect functioning to complete failure. A review of MSS literature can be found, for example, in (Levitn, 2005;Levitin and Lisnianski, 2001;Lisnianski and Levitin, 2003;Ushakov and Levitin;2002). The ROP for series-parallel MSS was initially presented in (Ushakov, 1987), where the universal generating function method was used for reliability evaluation (Ushakov, 1986). Many papers dealt with the development of heuristic methods based on metaheuristics for the ROP of series-parallel MSS (Kuo et al., 2001). These methods include genetic algorithms (Levitin et al., 1997;Levitin et al., 1998;Lisnianski et al., 1996), ant colony (Nourelfath et al., 2003), harmony search in (Nahas and Thien-My, 2010), a heuristic algorithm in (Ramirez and Coit, 2004), a coupling of tabu search and genetic algorithms in (Ouzineb et al., 2010), and artificial bee colony algorithm (Yeh and Hsieh, 2011). The present paper contributes to the literature on solution methods for the ROP of MSS. An efficient approach is developed, and compared to the state-of-the-art benchmark methods in terms of solution quality, computation time and implementation effort.
In the ROP of BSS (Chern, 1992), the system and its components are assumed to be either in a good state or in a failed state (Misra and Ljubojevic, 1973). Numerous optimization approaches and formulations have been proposed to solve the ROP under the binary-state assumption (Kuo and Prasad, 2000). These techniques can be classified as exact optimization approaches and non-exact or heuristic approaches. Exact techniques include dynamic programming (Bellman and Dreyfus, 1958;Fyffe et al., 1968;Nakagawa and Miyazaki, 1981;Yalaoui et al., 2005), mixed-integer and nonlinear programming (Tillman et al., 1977), and integer programming (Bulfin and Liu, 1985;Gen et al., 1993;Ghare and Taylor, 1969;Misra and Sharma, 1991). More recent exact methods are (Billionnet, 2008;Caserta and VoB, 2015a;Caserta and VoB, 2016;Onishi et al., 2007). Nonexact methods based on metaheuristics include genetic algorithms (Coit and Smith, 1996b;Painton and Campbell, 1995;Yokota et al., 1995;Yokota et al., 1996), tabu search algorithm (Kulturel-Konak et al., 2003), artificial neural networks , simulated annealing (Kim et al., 2006), ant colony optimization (Liang and Smith, 2004;Nahas et al., 2007;Nahas et al., 2008;Zhao et al., 2007), and harmony search (Nahas and Thien-My, 2010). More recent heuristics are (Chang and Kuo, 2018;Liu, 2016). In (Agarwal and Gupta, 2006;Agarwal et al., 2010), the authors focused on solving highly constrained redundancy optimization problems in binary complex systems. As the ROP of binary-state series-parallel systems has been recently solved to optimality (see for example (Caserta and VoB, 2015b), it is solved in this paper just to illustrate the efficiency of our approach in comparison to other existing meta-heuristics approaches. The objective of proposed approach is rather to solve the more recent ROP of multistate systems, which has never been solved to optimality.
All the above-mentioned approaches were (independently) proposed to solve either the ROP of binary-state systems, or the ROP of multi-state systems. When the proposed approaches are based on the same metaheuristic (for BSS and MSS), the operators and the algorithm steps are either different, or inefficient if used for both cases due to different modeling assumptions. Consequently, a method that is initially designed to solve the ROP of BSS cannot be used for MSS, and viceversa. To the best of our knowledge, there is no existing method designed to solve both the ROP of series-parallel BSS and MSS. The present contribution develops a unified and efficient heuristic based on the non-linear threshold accepting (NLTA) metaheuristic , which consists in a simple deterministic iterative local search method. Recently, this algorithm was applied to efficiently solve the redundancy allocation problem considering multiple redundancy strategies (Nahas et al., 2018). The NLTA algorithm represents a variant of the conventional (linear) threshold accepting algorithm. Our approach also uses a self-adjusting penalty guided strategy. The resulting method is efficient and very simple to implement. Numerical results for various test problems from previous reported research are presented, both for multi-state and binary cases. The performance of our approach is verified by comparing its results with those of the stateof-the-art benchmark methods, using all the existing benchmark problems.
The manuscript is structured as follows. Section 2 develops the method used for the availability evaluation of multi-state systems. Section 3 presents a description, and a mathematical formulation for the studied redundancy optimization problem. Section 4 presents the non-linear threshold accepting algorithm and the self-adjusting penalty guided strategy. Section 5 details the application of this approach to solve the formulated problem. Section 6 presents the numerical results. Section 7 concludes the paper.

Availability Evaluation for Multi-State Systems
To optimize the multi-state system redundancy, it is mandatory to develop an efficient evaluation procedure to estimate the availability of each series-parallel configuration. We consider a system with states corresponding to different levels performance. For repairable MSS, availability is assessed by the probability to satisfy the demand (Levitn, 2005;Levitin et al., 1998;Ushakov, 1987;Xue and Yang, 1995): where, ( ) is the system performance at time , and is the demand.
The present paper uses universal moment generating function (UMGF) to calculate the MSS availability.
Consider an example of a power generation system. Reliability is considered as a measure of the aptitude of the system to meet the customer demand ( ) by supplying an adequate amount of electrical energy ( ). For repairable MSS, a steady-state availability A can be defined as ( ≤ ). The following equations present the distribution of states probabilities, and the steady-state availability: It is convenient to distinguish, for a period of operation, different intervals with specific durations ( 1 , 2 , . . . , ). To each interval, corresponds a demand level ( 1 , 2 , . . . , ), which needs to be fulfilled. In this context, the MSS availability is:

Definition and Properties
Let be a discrete random variable having accessible states. Each state is characterized by a performance level and its corresponding probability . The UMGF of is a polynomial function: The function ( ) can be used to find the probabilistic characteristics of . In particular, the availability is: such as the operative Ψ is distributive and defined by: As an illustrative example, let consider single components such as each component has a nominal functioning state with performance , and a total failure state with performance 0. The availability of each component is = ( = ), and ( = 0) = 1 − . The UMGF of such binary state component is: In the next subsection, we define some basic composition operators to evaluate the MSS availability of a series-parallel system. This is done by using simple algebraic operations on the individual UMGF of components depending on their interactions (Ushakov, 1986).

Composition Operators
Let consider a subsystem containing components connected in parallel. In this case, the performance of the system is the sum of the performances of its components. The following operator defines the UMGF of parallel components: The Γ operator is a product of the individual UMGF of components. We have: For binary state components whose individual UMGF are defined in equation (9), we have: (1 − + ).
In series system, the component having the least performance constitutes the bottleneck and defines the system performance. For example, the UMGF of two series components is evaluated using the following operator: It is possible to express the UMGF of the series-parallel configuration by recursively applying the operators Γ and on individual UMGF components (Levitn, 2005). The UMGF approach is convenient for numerical implementation and in high dimension combinatorial optimization contexts. It is advantageously used in this paper to solve the ROP of MSS. : total cost of series-parallel MSS (ℎ, , , )

The ROP of Multi-State Systems Notation
: availability of series-parallel MSS Each subsystem is composed of many identical components in parallel. Several components versions are available in the market. Each version ℎ has different characteristics: nominal capacity (ℎ ), availability (ℎ ) and cost (ℎ ). The subsystem is defined by the version number ℎ (1 ≤ ℎ ≤ ) and the number of parallel components (1 ≤ ≤ ). The total cost is = ∑ =1 × (ℎ ). We consider homogeneous structures, which means that only one component type is selected to provide redundancy. The decision variables are defined by the vector (h,r) specifying the MSS structure, with = (ℎ 1 , ⋯ , ℎ ) and = ( 1 , ⋯ , ).
The goal of MSS availability optimization is to minimize the system cost with respect to a required availability level 0 (Lisnianski et al., 1996): Subject to (ℎ, , , ) ≥ 0 In the classical binary-state case, each component version has a different reliability, cost and weight. The goal of the considered optimization problem is to maximize the system reliability with respect to weight and cost limitations.

The Non-Linear Threshold Accepting Algorithm
Notation : parameter which decreases during the NLTA process

The Conventional Threshold Acceptance Algorithm
The threshold accepting (TA) algorithm is a simple local search algorithm developed in (Dueck and Scheuer, 1990). Like any local search algorithm, TA algorithm starts by a feasible solution ( 0 ), explores the neighborhood of the current solution and moves to a neighbor solution ( ) if it improves the objective function or if the deterioration (i.e. ( ) − ( 0 )) is less than a threshold parameter . The parameter is decreasing by a constant value Δ at each iteration. The pseudo-code of the method, for a minimization model, is given below: Step 1. Construct an initial solution 0 Step 2. Initialize 0 and Step 3. Set = 0 Step 4. While < − 1 -Create randomly a neighbor solution ′ and set Δ = ( ′) − ( ) -If Δ < then set +1 = ′ else +1 = Step 5. Set +1 = − Δ and = + 1 Step 6. End While

Extension
A new acceptance rule is used to extend the conventional threshold acceptance algorithm. To introduce this rule, we first define the following acceptance function: where, 0 is a positive parameter, and is a parameter which decreases during the execution of the algorithm.
The resulting extension is similar to simulated annealing and threshold accepting algorithm because during the search process it may accept low-quality solutions according to a different acceptance criterion. Initially, parameter is set to adequately high value. Since, the acceptance function ( ) (equation (12)) is inversely proportional to , any decrease in will increase the value of acceptance function, thus making it hard for a bad solution to be accepted. For any new solution 0 generated, it may be accepted for one of two reasons. The new solution is improvement over the current one, i.e. ( ) < ( 0 ), or it is accepted based on the comparison between ( ) and 0 , where 0 is the ratio between the objective values of the current solution and the new one ( 0 = ( )/ ( 0)) . Whenever an iteration changes, or a new solution is accepted, it reduces the value of by a constant value Δ . The parameters 0 and Δ are used to regulate the decay rate and thereby the search process. Manipulating the values of these parameters controls the convergence of ( ). Having an acceptance function that decreases in a nonlinear way may allow for a better explorative ability of the algorithm. One can remark that Equation (12) resembles the transfer function of a low-pass filter (RC-filter) used in electronics to reduce the amplitude of signals with higher frequencies. The following steps present the proposed algorithm for a minimization problem: Step 1. Construct a solution 0 Step 2. Initialize 0 , , Δ and (maximum number of iterations) Note that the objective function of a given solution is penalized in our approach. In what follows, we demonstrate that this simple method, coupled with a penalty guided strategy, efficiently solves hard instances of the ROP of MSS and BSS. Before presenting the computational results, let us explain first the rationale behind using this non-linear function.

Why the Proposed NLTA Approach is Successful?
While inheriting the characteristics of TA and simulated annealing, our algorithm use a non-linear acceptance function. To answer the above question, we need then to better understand the role of the non-linear decay rate on the behavior of the algorithm.

High Jump Competition Metaphor
As a metaphor, let consider the process encountered in high jump competition. The high jump is a track and field athletics event in which competitors must jump over a horizontal bar placed at measured heights. In competition the winner is the person who cleared the highest height. An important remark that has inspired our work is that the bar is raised in a non-linear way. For example, it can be raised by 5 cm at the beginning, then by 3 cm after, etc. Let consider a set of competitors that are chosen from a large number of persons. A competitor is chosen if he performs the best jump or if his jump is within an acceptable range. This range is larger at the beginning but lowered in a non-linear way as indicated above. In our metaphor, the persons participating to the high jump competition represent the search space. Choosing between the competitors corresponds to a number of iterations. In this metaphor, it is more probable to select a good competitor if the bar is raised in a non-linear way. Also, for more adaptability to competitors performance, it should be better to evaluate the decay depending on previous performance of the competitors.

Non-linearity and Adaptability
The metaphor above highlights that non-linearity and adaptability are desirable in the search process of good solutions (competitors). The non-linear decay rate adapts to the quality of the best solution more easily than a linear one. In fact, the NLTA adjusts the decay rate at every iteration. On the one hand, substantial improvements are observed, the decrease of the decay rate is also high. On the other hand, slower improvements in the best solution imply slower reactions of the decay rate. Unlike the classical (linear) TA, the NLTA should never become greedy, since the decay rate is adapted to the quality of the best obtained solution. These characteristics and the penalization strategy have allowed for an efficient solution of the ROP of BSS and MSS as described in the next section.

Application to the ROP of Multi-State Systems
To apply the proposed NLTA approach, we need to define the neighborhood space and the way the objective function is penalized.

Neighboring Solutions
Since we deal with homogeneous series-parallel MSS, each subsystem contains a number of identical components. The following procedure is used to define the neighborhood solutions:

Penalizing the Objective Function
For this problem, infeasible solutions are also accepted to force the penalized search towards the feasible solutions area. That is, the proposed penalized objective function is: The value of is updated as follows: • If all solutions are infeasible, set to × 1 .
• If all previous solutions are feasible, set to × 2 where, 1 and 2 are positive numbers.
For the ROP of BSS, the move applied to the current solution consists in either adding or subtracting one component. A subsystem is randomly selected. Then, the version of the added of subtracted component is randomly selected from the available versions. Also, infeasible solutions are accepted to force the penalized search towards the feasible solutions area.
The penalty functions are introduced to encourage the algorithms to explore both feasible and infeasible regions. Searching near the boundaries have proven beneficial. Gendreau et al. (1994) introduced the concept of adding infeasible solutions in the initial population. However, to penalize these solutions, they used self-adjusting penalties. Weights increase if, during the last few iterations, only infeasible solutions are found. They decrease if all recent solutions are feasible. In (Coit and Smith, 1996a;Coit and Smith, 1996b), the authors proposed GA to solve the RAP and found that penalizing the infeasible solutions benefits in improving final feasible solutions. Similar techniques are applied in (Kulturel-Konak et al., 2003;Ouzineb et al., 2008). Using the same concept, the penalty weights have been adjusted in our proposed approach.

Computational Results
The algorithms were implemented in MATLAB on a computer with a processing speed of 2 Ghz.

ROP of Multi-State Systems 6.1.1 Test Problems
Four benchmarks are used. The first three examples were introduced in (Levitin et al., 1998), (Lisnianski et al., 1996) and (Levitin et al., 1997). The fourth example was introduced in (Ouzineb et al., 2008). Table 9 and 10 present the data of Example 1. Three different values of the availability index 0 have been used to generate three variations. Table 11 and 12 present the data of Example 2. The reliability index 0 has been set to eight different values to create eight variations of the problem. Considering possible discounts, the cost of the component is defined as a function of the number of components purchased: where, * = is the cost per unit, and 1 , 2 , 1 , 2 are shown in Table 13 for each component type Table 14 and 15 present the data of Example 3. The availability index A0 has been set to three different values to create three variations of the problem.

Parameter Settings
For this problem, the user-specified parameters Δ , 0 and have been tuned by following the same tuning procedure used for binary-state systems. The found values to be most appropriate were: 1/ 0 = 0.0085, = 50 and Δ = 0.0001. In addition, preliminary experiments showed that the best results are obtained when the initial value of is set to 10 with = 5, = 5 , 1 = 0.99 and 2 = 1.01.

Comparison Results
The algorithm is evaluated by comparison to the existing benchmark methods: genetic algorithms (GA) (Levitin et al., 1998;Levitin et al., 1997;Lisnianski et al., 1996) and tabu search (Ouzineb et al., 2008). Table 1 Table 2 shows the total cost of the best solution for each of the 17 instances. The average (Av.) and the standard deviations (St.D) of the best solutions found in 10 trials are also presented in Table 2. This table shows that:  In 6 of the 17 test cases, the NLTA outperforms the genetic algorithm, while it provides the same results in the other cases.  In all cases tested, the results of the NLTA are in par with the best known results. The above results clearly highlight that our approach outperforms the existing methods. The reasons of this performance are discussed in the next subsection.

ROP of Binary-State Systems 6.2.1 Test Problems
The results of the algorithm are compared to the results reported in the literature for the 33 problem instances introduced by (Fyffe et al., 1968;Nakagawa and Miyazaki, 1981). The system consists of 14 subsystems. In Appendix A, Table 17 provides all the input data for the examples studied in this paper, with , and are (respectively) the reliability, the cost and the weight of component for subsystem .

Parameter Settings
In meta-heuristic approaches, like simulated annealing and threshold accepting algorithm, to achieve high performance, it is essential to adjust a number of input parameters. In the case of NTAA, there are three user-specified input parameters: Δ , 0 and . Parameter values most favorable to the considered optimization model need to be explored. Therefore, similar to common practices used in the literature, 10 simulations were performed for each case, where value of one input parameter was explored while keeping others constant and some statistical information about the average evolution of the key parameters was collected. Multiple tests were carried out, setting with the available budget = 130 and the allowed weight constraint = 191, in order to identify the suitable values for the parameters. Maximum number of iterations (i.e. 2 × 10 6 ) was used as stopping criteria. By performing these experiments, the most suitable values for the parameters were: 1/ 0 = 0.04, = 30 and Δ = 1.4 × 10 −5 .

Comparison Results
Benchmark methods include genetic algorithm (GA) (Coit and Smith, 1996b), tabu search (TS) (Kulturel-Konak et al., 2003), the ant colony optimization (ACO) (Liang and Smith, 2004) and the Harmony search algorithm (HSA) (Nahas and Thien-My, 2010). The obtained results for 10 trials are shown in Table 3. This table presents the system reliability of the NLTA solution for each instance. The maximum ( ) for the best solutions found in 10 trials and the average ( ), the standard deviations (St.D) of the 10 solutions and the average computation times. From this table, we can see that the standard deviation is very low and the average computation time does not exceed 87 seconds for any instance. Figure 1 shows the convergence of the best solution of the NLTA for = 130 and = 191. From this figure, we can see that not more than 15 × 10 5 iterations are necessary for the algorithm to converge. Indeed, this number is much lower than the size of the search space (larger than 7.6 × 10 33 for this instance example). Table 4 presents the comparison results of the NLTA and the best solutions found in previous papers: GA (Coit and Smith, 1996b), TS (Kulturel-Konak et al., 2003), ACO (Liang and Smith, 2004) and HSA (Nahas and Thien-My, 2010). To illustrate the efficiency of the proposed approach (NLTA), the following errors are defined: where, ( ) represents the best result obtained by NLTA and (ALG) represents the best result obtained by another approach, i.e. GA, TS, ACO or HSA.  Table 4 and these figures clearly highlight that NLTA outperforms the other approaches.

Discussion
In general, it is not always easy to justify why a given metaheuristic provides good results to solve a given problem. It is however widely recognized that a potentially good approach should seeks to balance between diversification and intensification strategies. Diversification explores the search space by examining unvisited regions, while intensification exploits a local region to look more closely around good solutions identified during the past search. Our method tends to find the right balance between diversification and intensification.
For the redundancy optimization problems, the best known solution methods are based on metaheuristic approaches such as GA, TS, ACO and HSA. The algorithm proposed in this paper is based on a simple iterative local search method. The impressive results found by this algorithm give rise to an important question: Why such simple method performs very well for the ROP? Our method has two important characteristics: the use of a self-adjusting penalty strategy, and the nonlinearity of the acceptance rule. To analyze the effect of these characteristics, additional experiments were conducted to quantify their contributions to the algorithm efficiency.   Tables 5 and 6 present the results obtained without the penalty strategy for the ROP of binary and multi-state systems, respectively. Figures 3 and 4 present the errors between the NLTA versions with and without penalties and show that the used penalty strategy plays an important role in solving the ROP.
Furthermore, to analyze the importance of the non-linearity, we tested the conventional threshold accepting algorithm (Dueck and Scheuer, 1990), with and without penalty strategies. Tables 7 and  8 present the comparison results, while Figures 5 and 6 sketch the errors and highlight two important points. First, the linear threshold algorithm without penalty strategy provides relatively bad solutions. Second, the introduction of the penalty strategy into the linear threshold algorithm improves the solution quality, but the overall efficiency remains much less competitive in comparison to existing methods. This means that the efficiency of our method comes not only from the use of the self-adjusting penalty strategies, but from the non-linear form of the acceptance rules as well. Under such a non-linear control of the search process speed, the intensification is also encouraged, and the solution acceptance depends on the value of the fraction between the objective value of the current solution and the value of its neighbor.

Conclusion
In this paper, we developed an evaluation/optimization approach to solve the redundancy optimization problem of series-parallel multi-state systems. This approach provides a powerful tool to reliability designers, when different varieties of components can be considered to achieve highavailability for complex multi-state systems. Extensive numerical results for the test problems from previous research highlight the advantages of the proposed approach in terms of solution quality, computation time and implementation effort. Using existing benchmark problems, it was found that this approach out-performs the state-of-the-art benchmark methods dealing with the ROP of multistate systems. Not only the proposed approach is shown to be efficient, but it is also simple to implement. Two important features characterise the paper contribution and justify the efficiency of the proposed approach. First, the availability evaluation method uses the universal moment generating function technique, which is convenient for numerical implementation in high dimension combinatorial optimization. Second, the optimization algorithm is adaptable and uses a non-linear acceptance function. Such non-linearity and adaptability are desirable in the search process of good solutions. Furthermore, by using a self-adjusting penalty guided strategy, this algorithm is encouraged to explore both feasible and infeasible regions.
To achieve high availability levels for series-parallel multi-state systems, this paper considered that each subsystem was composed of identical components (i.e. homogeneous systems). It was also assumed that the characteristics of components are known and deterministic. While the extension to heterogeneous systems is straightforward, taking into account parameters uncertainties may require further methodological development. We are currently working on these issues, and on the extension of our algorithm as general-purpose method to deal with others combinatorial optimization problems encountered in the context of multi-state system reliability. The data for multi-state systems are presented in Tables 9-16 and those of binary-state systems  examples are in Table 17.

Confilct of Interest
The authors confirm that there is no conflict of interest to declare for this publication.