AMultiscale Differential Evolution Algorithm-Based Maintenance Plan Optimization for Building Energy Retrofitting

This paper addresses a maintenance plan optimization problem for building energy retrofitting with a novel differential evolutionbased algorithm. The problem manifests substantial complexity due to the parallel optimization of the maintenance time scale, maintenance instants, and maintenance rates, with multiple objectives and time scales. The problem is formulated based on an impulsive and switched phenomenon, with objectives of maximizing the energy performances and economic indicators for the retrofitting project over a finite horizon. A multiscale differential evolution algorithm is accordingly proposed to be the numerical solver. The proposed algorithm is able to optimize the decision variables simultaneously, especially the maintenance time scale that can influence the number of other decision variables. Simulation results illustrate the effectiveness of the proposed approach.


Introduction
While buildings result in a significant amount of energy consumptions [1], the building energy efficiency now becomes a topic that receives massive attention.The building energy retrofitting is considered to be a major method to surpass the energy efficiency in existing buildings.One major research focus for the building retrofitting is the retrofitting planning [2][3][4].In addition, studies also reveal that the maintenance of retrofitted items makes significant contributions to the final performance of a building energy retrofitting project, due to the negative impact from item malfunctions and deterioration [5][6][7].
The building energy retrofitting refers to activities in existing buildings that apply energy conservation interventions to improve efficiency or conserve energy or water, or manage demand [8].A common category of interventions is the replacement of less efficient equipment, such as lighting devices, air conditioners, and geysers [9,10].In a retrofitting project, a large amount of retrofitted items are involved, as well as substantial magnitude of capital investment.During operation, malfunctions and deteriorations are inevitable to the retrofitted items.Such a deterioration results in loss of energy performance according to the measurement and verification principle [8].As the maintenance can reverse the deterioration, it plays an important role in the sustainability of building energy efficiency.Therefore, the maintenance plan optimization is incorporated into the retrofitting planning.
The maintenance plan optimization problem introduced in [6] involves selecting the maintenance rate for a number of retrofitted items over fixed time scheduling, such that as much energy as possible can be saved and satisfying economic performance can be achieved, based on limited amount of capital investment.The term "maintenance rate" describes such a scenario that only a proportion of the failed items are maintained or repaired.The joint effects of item failures and maintenance actions result in dynamics of the overall energy performances over time.Therefore, the maintenance plan optimization problem is modelled to be an optimal control problem.Given that the aggregated performances of the retrofitted items are the main concern, the number of the working items is selected to be the state variable, and the maintenance rate is selected to be the control variable.
Wang et al. [7] further take into account the performance deterioration and failures simultaneously in the preceding maintenance plan optimization problem.A multistate system approach is employed, where retrofitted items are considered to have multiple working states and a failure state [11,12].The transition between working states describes the performance degradations over time.
To reverse such a degradation, two types of maintenance actions are introduced, namely, the preventive maintenance (PM) and corrective maintenance (CM).The PM is defined as a maintenance action that restores the system to a specified better condition when the system is still functioning, and CM refers to any action that restores the failed system to its working state [11].Similarly, the multistate-based maintenance plan optimization problem is interpreted into an optimal control problem, where numbers of items subject to each working state are the state variables and maintenance rates for PM and CM are the control variables.The maintenance time scheduling is predecided; i.e., the time instants where maintenance actions take place are known a priori and fixed.This also suggests that the time scale of the maintenance plan is bounded and fixed.The maintenance plan optimization problem is thereby a finite dimensional problem that is solved by varying the PM and CM rates.
Obviously, the aforementioned investigation is limited by the constraint of fixed time scheduling.In fact, the time scheduling is a common decision variable in the reliability engineering area [13].It is straightforward to take into account time scheduling optimization in addition to the maintenance plan optimization for building energy retrofitting, for the sake of further improvement of the aggregate performances.However, such an idea can lead to two difficulties.Firstly, the preceding control system modelling has to be improved.Due to the joint effects of the maintenance rate, maintenance time scheduling, and maintenance time scale, it is necessary to collectively optimize the three categories of decision variables to achieve the optimum performance.A control system modelling that includes the involved decision variables is required.Secondly, taking into account maintenance scale develops the optimization problem to be a multiscale one; i.e., the candidate solutions can have different scales.Such a fact results in a difficulty to employ conventional numerical solvers.Although multiscale optimization has been investigated in literatures [14][15][16][17], the existing approaches often develop reduced models and require detailed understanding of specific problems.As a result, the investigated problem is a complex one, and there lacks a proper approach to interpret and solve such a maintenance plan optimization.
In this paper, the preceding control system modelling is extended, where an impulsive and switched approach is employed and solved.The maintenance actions are considered to be a sequence of impulsive effects at variable instants that switch the retrofitted items to different working states.The numbers of items under the different working states, i.e., the system states, jump with such impulsive effects, whose impacts are indicated by the maintenance rates of PM and CM.The time scheduling describes the sequence of the impulsive effects, where the number of time instants is also variable.In this way, a hybrid system with impulsive effects and switching phenomena [18] is employed to model the system dynamics, namely, an impulsive and switched maintenance planning.The retrofitting items are categorized into several groups; each consists of items that are considered to be homogeneous ones, i.e., with the same inherent energy and reliability performances, the same operating schedules, and similar operational environment.Thereafter, the performances of the totality of the item groups can be characterized by the population of each item group.Further discussions about the item grouping can be found in [19].Such a problem is translated into a parameter optimization problem.There are two objectives: maximizing the overall energy savings and internal rate of return (IRR) over the sustainability period, according to [7].A weighted sum of the two objectives is employed to allow the numerical solver.The evolutionary algorithms (EAs) were a good choice to many complicated engineering problems [20,21].However, the conventional approach is less efficient as a solver to the investigated problem.It has to apply a trial-and-error process to find out a proper scale, based on which the maintenance rates and instants are optimized.Such a process can repeat substantial iterations to figure out a satisfying solution.Given the difficulties, a multiscale differential evolution (MSDE) approach is hereby proposed as the numerical solver.The MSDE employs a multiple population mechanism that has been employed in other studies of EA algorithms for the sake of better performances [22].In MSDE, the major concern is the parallel searching of the maintenance scale, rates, and instants, namely, the optimal combination.It divides the candidate solutions into several small subgroups, each subject to a different maintenance scale.A candidate solution consists of several pairs of maintenance rates and instants.The number of pairs is decided by the maintenance scale.The competitions simultaneously take place among the candidate solutions within same group and among the subgroups that pertain to different maintenance scales.As a random heuristic search algorithm, MSDE can hopefully, although not guaranteed, discover a satisfactory solution with optimum scale after sufficient iterations.Simulation results illustrate the effectiveness of the proposed approach.
The remainder of the paper consists of four sections.Section 2 gives the maintenance plan optimization problem formulation on an impulsive and switched basis.Section 3 explains the MSDE algorithm and the corresponding numerical solver design.Section 4 introduces the case study, the simulation results, and analysis.Section 5 draws conclusion and discusses future researches.

Maintenance Plan Optimization Formulation
2.1.Problem Statement.The maintenance plan optimization problem is established based on the retrofitting problem.Given a retrofitting project, the performances over time are estimated on a life-cycle cost analysis basis; i.e., the varying performances, operational costs, and maintenance costs are all taken into account.More details about such a retrofitting planning problem can be referred to [23].Accordingly, maintenance planning aims at selecting proper maintenance actions to further optimize the performances of retrofitted items over the sustainability period.
The maintenance plan optimization problem is modelled on a discrete time system basis.Given a series of time instants t k = kS with k = 0,1,2, … , T , let t k denote the sampling instants where the system state is obtained, S the sampling interval.A finite time horizon 0, TS is identified, i.e., the sustainability period.As introduced above, the decision variables of the maintenance plan optimization are the combination of the maintenance time scale, maintenance time scheduling, and maintenance rates.Let t k p and t k c denote the preventive and corrective maintenance instants, where k p = p 1 , p 2 , … , p N p and t k c , k c = c 1 , c 2 , … , c N c .Accordingly, N p and N c denote the maintenance scale.The maintenance time schedule is decided via identifying the maintenance instants.Let denote the indices of preventive and corrective maintenance instants, respectively.The elements of Q p and Q c are selected from k = 0,1,2, … , T , implying that the maintenance instants are commensurate with sampling instants t k .
As aforementioned, the performances reveal a dynamic nature under the joint impacts of failures (malfunctioning) and maintenance.As a result, an idea of the control system approaches is employed for the maintenance plan optimization.Generally, we employ x t k to represent the variable population of the working items, and x t k the maintenance rates.The optimization is to find proper N p , N c , Q p , and Q c , and corresponding u t k , such that a set of optimal x t k can be achieved to maximize the performance indicators, e.g., the energy savings and IRR.The contents of x t k and u t k are explained based on the optimization problem formulation as the following.
2.2.System Dynamics Modelling.The mathematical description of the system dynamics is hereby explained.Consider the following discrete dynamical system over 0, TS : where f • is a continuously differentiable function.Let the system state x t k i jump at instants t k i , i = 1, 2, … , N , namely, the switching instants.Such jumps are represented by the following: where △x t k i denotes an impulsive effect at t k i that incurs a jump of the system state.t k i is employed to describe the moment right after the impulsive effect.Accordingly, x t k i denotes the system state prior to the impulsive jump and x t k i the system state after.The system dynamics of a maintenance planning problem is established as the following: Equation (3) describes such a homogeneous retrofitted item group l that contains N l items with M l working states; M l denotes the best working states, and 1 denotes the worst.The system state x l,i t k , i = 1, 2, … , M l denotes the numbers of working items that are corresponding with group l and state i. f l i,j x l,i t k , t k denotes the proportion of items that degrades from state i to a worse state j over interval t k , t k+1 .g l,i x l,i t k , t k denotes the malfunctioning items from state i, namely, the population decay model.For the convenience of discussion, (3) is abbreviated to The impulsive jumps are described as the following: Similarly, △x l,i t k denotes the impulsive effects and t k the moment post to jumps.u l,i x l,i t k p , t k p , i = 1, 2, … , M l − 1 denotes the preventive maintenance rate and u c l x l,0 t k c , t k c the corrective maintenance rate.x l,0 t k denotes the number of malfunctioning items.The control of ( 3) and ( 4) is to identify proper N p and N c , Q p and Q c , and u c l and u l,i , namely, a maintenance plan.Let u p l t k denote the overall rate of preventive maintenance and u c l t k the overall rate of corrective maintenance at t k : According to (3), (4), and ( 5), there is a dependence of the maintenance plan and system states x l,i t k ; i.e., (3), (4), and ( 5) describe a closed-loop system.
2.3.Performance Indicator Formulation.The estimation of the aggregate performances requires two parts of information, the population and individual item performances.The population is estimated by the preceding system dynamics modelling, and the item performances are characterized by a series of performance indicators that are explained as the following.Firstly, an assumption is made such that the retrofitted items can be categorized into several groups; items from the same group are considered as homogeneous ones with the same energy and reliability performances.Let N denote the number of groups of homogeneous items.There are three categories of performance indicators, including the average energy savings, cost savings, and maintenance costs per item.For an item group l, the performance characteristics are represented by the following: where a l t k denotes the average energy savings over sam- and the corresponding cost savings are obtained: and the maintenance cost at each time instant is obtained: where M c t k denotes a constant cost that is charged when maintenance action takes place at t k .Such a cost refers to the expenditure of activate maintenance actions.The total expenditure of the retrofitting project is represented by where h 0 denotes the initial investment to the retrofitting project.The profit of the project is then obtained by P = B| all − h| all .The net present value (NPV) of the project over 0, TS is formulated as following: i.e., the IRR is the discount rate that makes NPV = 0 over horizon 0, TS .The IRR cannot be represented analytically, but still a bounded value subject to initial investment h 0 [24].3), (4), and ( 5) and performance indicator formulations ( 6), ( 7), ( 8), ( 9), ( 10), (11), and ( 12), the maintenance plan optimization problem is established as the following.
subject to system dynamics (3), (4), and (5) and the following constraints, In the above equations, α is the minimal acceptable energy saving value, namely, the targeted energy saving.α is usually given by the energy performance contract of a retrofitting project.d| T R is the IRR.λ 1 and λ 2 denote the weighting factors.The objective function formulation indicates that the energy and economic performances are simultaneously taken into account, where the weighting factors adjust the importance of the two objectives.β is the maintenance budget limit over 0, T , and T ′ is the payback period limit.T p is the payback period, which is computed from the last t n that makes NPV < 0. Let M N denote the total number of working states for the N item groups, u p • ∈ ℛ M N −N ×S p and u c • ∈ ℛ N×S c ; i.e., the minimization problem ( 13) and ( 14) is a finite dimensional problem.Equations ( 13) and ( 14) are different with the conventional optimal control problem.Firstly, the control inputs to be solved are a set of parameters.Secondly, the objective function is a weighted sum of two different performance indicators, instead of the quadratic performance index.Thirdly, d R | T is nonanalytic; as a result, the gradient method becomes infeasible to ( 13) and (14).Furthermore, N p and N c can change upon searching the optimum, implying that multiple scales can exist in such one optimization problem.As a result, a multiscale differential evolution (MSDE) algorithm is hereby proposed as the numerical solver.

A Multiscale DE Solver
A general form of the multiscale problem is employed as the following: to represent the investigated optimal control problem.x and y constitute a candidate solution.x indicates an implicit vector with a variant scale N x y that is determined by another vector y .y is an explicit vector with a fixed scale D. To the minimization problem ( 13) and ( 14 ω is a large positive coefficient that adjusts the penalty from the violation of constraints.The penalty functions P n with n = 1,2,3 are given as the following: To solve such a multiscale minimization problem, we propose a multiscale-based approach, namely, multiplescale DE (MSDE).The MSDE is actually an improvement of the conventional DE algorithm.Therefore, the concepts of basic DE will be introduced firstly and thereafter the implementation of MSDE.
3.1.Basic Differential Evolution.The differential evolution (DE) is a popular category of the evolutionary algorithm that is suitable to solve the continuous optimization problem.The DE algorithm can hopefully, although not guaranteed, 5 Complexity discover a satisfactory solution after sufficient iterations for most types of optimization problems.In the DE algorithm, the term "individuals" are adopted to refer to the candidate solutions that represent the possible values of the decision variables.These individuals are moved around in the search-space via a series of mathematical operations including Mutation, Crossover, and Selection.
In order to implement the DE algorithm, a population of NP individuals is generated at the initial state to cover as more search space as possible.Therefore, a uniform randomization is usually adopted at the initial stage.As fixed-scale problems can be regarded as a special case of multiscale problems, the following minimization problem min f x , x = x 1 , x 2 , … , x N ∈ ℛ N is adopted for discussion.The search space is constrained in N j=1 L j , U j , L j ≤ U j , where L j is the lower bound of search space and U j the upper bound.
At each generation, an individual by the mutation operation.There are many possible mutation operators from literatures, such as the following: More mutation operators be found in [25].In the above equations, the indices r 1 and r 2 are distinct integers uniformly randomized from the set i = 1, 2, … , NP .x r1,G − x r2,G is a difference vector to mutate corresponding target vector x i,G .The parameter F is the mutation factor which usually ranges on the interval (0,1).x best,G is the best solution at the current generation G. v i,G is the obtained mutant vector of x i,G .
After the mutation phase, DE performs a binomial crossover operation to x i,G and v i,G to generate a trial vector u i,G = u 1,i,G , u 2,i,G , … , u D,i,G for each i = 1, 2, … , NP, according to the following: where j = 1, 2, … , D. j rand is a random integer that ranges between 1, D .j = j rand implies that an element in vector u j,i,G is chosen, i.e., the acceptance of the mutant vector v i,G .rand 0, 1 is a random real number that ranges between [0,1].For each pair of j and i, rand 0, 1 is randomized independently and compares with a constant crossover rate CR, CR ∈ 0, 1 .With such an operation, the trail vector u i,G is guaranteed to be different with the target vector x i,G on at least one dimension.Such operation can result in violation of the boundary limits.In that case, the following repair is applied: The offspring are selected from the target vectors and the trail vectors, which is called selection operation.DE employs a greedy "one-to-one" replacement scheme.For each pair of x i,G and u i,G , the fitness of solutions f • is compared to find out a better solution that is selected to enter the next generation:

Complexity
The pseudocode of MSDE is given in Algorithm 1, and a flowchart is given in Figure 1.The idea of MSDE is further explained as the following.By decomposing the solution population into several subpopulations (islands), the individuals can evolve separately.There are double levels of evolution in MSDE.On the first level, individuals from the same subpopulation live in a local environment and interact only with the neighbors.On the second level, the subpopulations are connected by comparing the fitness of its local best.
The better fitness of the local best reveals a trend to find a better solution subject to the corresponding scale.Therefore, such an explicit variable is kept after the competition, and the other explicit variables have to be modified.Furthermore, a growth of the subpopulation with best explicit variable is encouraged; as a result, the worst Definition: NP denotes the population size; n denotes the quantity of subpopulations; NS denotes the subpopulation size; m denotes the shuffling period; D denotes the dimension of the explicit variable; N denotes the dimension of the implicit variable; F denotes the mutation factor; CR denotes the crossover rate; S F,CR denotes the set of mutation factor F and CR for each subpopulation.1: BEGIN 2: G = 1; 3: Initialize n subpopulations (including n explicit variables and NP implicit vectors); 4: Evaluate the fitness for each individual; 5: Set m, S F,CR ; 6: while the stopping criterion is not satisfied do 7: while k = 1 to n do 8: Randomly choose (F k , CR k ) from S F,CR ; 9: Find the local best y k,G , x best,k,G ; 10: end while 20:  The major differences between the basic DE and MSDE are as the following.The MSDE has several subpopulations and accordingly multiple scales over its iterations, while basic DE has only one population and one fixed scale.During iterations of MSDE, the proper scale is identified simultaneously upon searching the optimum of other decision variables, while original DE lacks a mechanism to optimize the problem scale as well as the implicit vectors.Therefore, computation burden of the MSDE base solution is less than original DE to a multiscale problem.

Simulation Results and Analysis
In order to verify the effectiveness the proposed approach, a case study is established and a simulation is performed based on the case study data.The case study is a part of a practical office building retrofitting project.Two categories of retrofitted items are involved.The first category is the compact fluorescent lamp (CFL).A type of 15 W CFLs is applied to replace the old inefficient lights.The CFLs are considered to be unrepairable items, with only one working state.The population decay model of the CFLs is given by where x cf l denotes the population of the CFL group, and x 0 cf l denotes the initial group, x cf l t 0 = x 0 cf l .b cf l and c cf l are the parameters of the decay model.The parameters can be estimated by experimental data fitting.More details about decay   22) can be found from [26].The second category is the air conditioner.More efficient air conditioners are installed to replace the old ones.The air conditioners are considered to be repaired items, with three working states: good, average, and bad.Based on the multistate model ( 3), the working state transition as well as the population decay of the air conditioner group is estimated by f a,c i,i−1 x ac,i t k , t k = x ac,i t k 1 − e ς p ac,i , i = 3, 2, g ac,i x ac,i t k , t k = x ac,i t k 1 − e −ς c ac,i , i = 3, 2, 1,

23
where x ac,i denotes the number of air conditioners under working state i. i = 3 denotes the good state, i = 2 the average state, and i = 1 the bad state.x ac,3 t 0 = x 0 ac , where x 0 ac denotes the initial population of the air conditioner group.x ac,3 t 0 = x ac,1 t 0 = 0. f ac i,i−1 denotes the amount of air conditioners that transit from working state i to i − 1, and g ac,i denotes the number of items that become malfunctioning from state i.The parameters ζ p ac,i and ς c ac,i are small positive constants that indicate the interval of the state transition.Equation ( 23) is obtained from the constant failure rate model, a common model from the reliability engineering.Detailed discussions of such a model can be found from [27].Table 1 gives the values of b cf l , c cf l , ζ p ac,i , and ς c ac,i in the case study.Table 2 gives the average performance characteristics per item.There are five categories of performance characteristics.The unit price is the initial cost of applying such an intervention, i.e., installing an item in this case study.The unit energy savings are the estimated monthly savings against the preretrofit items.Such a saving is verified following the M&V principles.Accordingly, the monthly cost savings, i.e., the reduced operational costs from energy savings, are estimated.
The preventive costs and corrective costs refer to the expenditure of applying a preventive maintenance or corrective maintenance to a failed item.Such expenditure includes the cost of both appliances and manpower.Furthermore, the initial quantities of the respective item group are given as well.At the initial stage, all the air conditioners are considered to be under good working state; therefore, quantities under average (i = 2) and bad (i = 1) are 0. From the unit prices and quantities, it can be identified that the initial capital investment h 0 is $20,692.Some configurations of the simulation are given as the following.For this case study, the sustainability period is set to be 10 years.The sampling interval is one month; i.e., there are 120 time instants over the sustainability period.The baseline energy performance is 8,685,311.7 kWh, which covers the pre-retrofit consumptions of lights and air conditioners.The baseline performance is extracted from a three-year energy bill of the inspected building.A 12% targeted energy saving is to be achieved by the retrofitting; as a result, the minimal energy saving is 1,042,237.4kWh.Furthermore, for the economical performance calculation, the payback period time limit is set to be 24 months, which is a common acceptable payback period limit in practice.For the maintenance cost calculation, a constant cost at each maintenance 9 Complexity instant is set to be $200.Several maintenance budgets are provided to establish different maintenance scenarios: $20,000 indicates a scenario with limited budget, $40,000 indicates a relatively insufficient one, and $65,000 indicates a sufficient one.As the comparison, a set of fixed maintenance time schedules are employed as the following: Q p = 7, 13, 19, 25, … , 115 and Q c = 13, 25, 37, 49, … , 109 .The fixed time schedules are also applied to the three scenarios of different maintenance budgets.The weights of the objective function ( 13) are λ 1 = 0 5 and λ 2 = 0 5, which is an equal consideration of the two objectives.For the MSDE configuration, 30 subpopulations are initialized, each containing 60 individuals.The mutation factor F linearly decreases from 1.0 to 0.2 over the iterations.It is selected such that the individuals can better converge to an optimum at the end of iterations.The crossover rate CR is 0.7 that encourages the crossover of individuals.The shuffling period is set to be 100 iterations, and the stopping criterion is 1000 iterations.Some simulation results are illustrated as the following.Table 3 shows the simulation results and Table 4 the obtained maintenance time schedules.Figure 2 the maintenance rates over time, and Figure 3 depicts the timely energy savings of the optimal maintenance plan with three different budget limits.The numbers are obtained by an average of 5-run results.From the tables, the optimal solution outperforms the fixed schedules; especially when the budget is sufficient, such an increment becomes larger.From Table 4, the maintenance actions take place more frequently than the fixed schedule.This reveals that the time scheduling can 10 Complexity further improve the maintenance planning performances, as a proper scale can imply a better solution in many cases.
From the simulation, the most important difference between the MSDE and the comparative results are the maintenance scale.A proper selection of the maintenance scale can result in significant improvement in the performances.A preliminary investigation reveals that, with a fixed time maintenance time scale, the population size can be reduced to half to achieve a similar performance.However, identifying a proper scale takes quite a number of trail-and-error processes.The conventional method can spend a lot more time than the MSDE approach.In general, the effectiveness of employing MSDE to solve maintenance plan optimization problems with multiple categories of decision variables is verified by the simulation, and the MSDE is promising to solve impulsive and switched optimal control problems with formulation.

Conclusions
In this paper, a maintenance plan optimization problem for building energy retrofitting that simultaneously involves variable maintenance time scale, maintenance time scheduling, and maintenance rates, is investigated and addressed.An impulsive and switched optimal control formulation is employed to model the novel problem based on existing maintenance plan optimization studies.By adopting a weighted sum of two objectives, the aggregate energy savings and internal rate of return over a finite time horizon, an optimal control problem that aims at finding the optimal combination of the maintenance time scale, instants, and rates is formulated.Given that the maintenance time scale can vary upon searching optimum, the problem becomes a multiscale one.A multiscale differential evolution algorithm is proposed to find the numerical solution to the multiscale optimal control problem with acceptable computational burdens.We believe that this approach can be applied to a general category of minimization problems with variant scales and facilitates the solution of other impulsive and switched optimal control problems.The effectiveness of the proposed approach is verified by simulation results.The current stage work calls for further studies from two major aspects: firstly, the maintenance plan optimization problem with variable time scheduling needs more investigation, such as the controller design and discussion on the robustness of the impulsive and switched modelling; secondly, the MSDE requires detailed studies on the mutation mechanism and shuffling period, such that the performance of MSDE can be further improved.

2. 4 .
Optimization Problem.Let P = P 1 , P 2 , … , P N denote the initial populations of the respective item groups, u p • denote the collection of u p l k p , and u c • denote the collection of u c l kc , l = 1, 2, … , N. Taking advantage of the preceding impulsive and switched system modelling (

2. 4 . 1 .
Optimization Problem ℳP O.For a dynamical system consisting of N subsystems in (3), (4), and (5) with state variablesx t = x 1 t , … , x l t , … , x N tand initial state x 0 , find the proper maintenance scale N p and N c , maintenance time schedule Q p and Q c , PM rates u p ⋅ = u p 1 kp , … , u p l kp , … , u p N kp , and CM rates u c ⋅ = u c 1 kc , … , u c l kc , … , u c N kc , to minimize the following performance index ), x denotes the time schedule Q p , Q c and maintenance rate u p • , u c • , y denotes the maintenance scale N p , N c .The maximum possible scale of x can be a large number denoted as N N ≠ D .Apparently, x ∈ ℛ M N −N+1 ×N p × N+1 ×N c and y ∈ ℛ 2 .The constraints are realized by a stationary penalty function in the objective function:

7
Complexitysubpopulation has to adopt a same explicit variable of the best subpopulation.

Figure 2 :Figure 3 :
Figure 2: The maintenance rates over time with different maintenance budgets.
t k and b l t k the cost savings.a l,i t k and b l,i t k are corresponding with x l,i t k .C l t k denotes the maintenance costs per item, where C 1,i t k is corresponding with preventive maintenance rate u l,i t k and C c l t k the corrective maintenance rate u c l t k .The overall energy savings over the sustainability period 0, TS are then represented by Complexitywhere d is the discount rate, d > 0. ϕ t k = 1, 2, 3 … indicates that the sampling instant t k lies within a specific year after the initial stage.Let d R | T denote the IRR.It is solved by letting the following equation equals to 0: .2. Multiple-Scale Differential Evolution.In MSDE, several subpopulations constitute the whole population of candidate solutions.Each subpopulation corresponds with a specific explicit variable y k,G , and individuals from one subpopulation share the same explicit variable.While individuals, i.e., the implicit vectors from the same subpopulation, are moved around to search better performances over iterations, the explicit variable is also changed by comparing the performances between subpopulations.
(1) A population of NP individuals is initialized into n subpopulations.Explicit variables are firstly initialized as y k,G = k = 1, 2, … , n .For the k − th subpopulation, NS individuals are then initialized with a specific scale determined by y k,G iables of other subpopulation other than the best and worst ones.For subpopulations with changing scale of explicit variables, all implicit variables have to be randomly generated in the search space in accordance with the new scale.If stopping criteria are not satisfied, go to step 2.

Table 1 :
Parameters for population decay and state transition models.

Table 2 :
Characteristics of retrofitted items.

Table 3 :
Performances of maintenance plans with different budget limits.

Table 4 :
Optimal maintenance time schedules.