# A Multiscale Differential Evolution Algorithm-Based Maintenance Plan Optimization for Building Energy Retrofitting

Authors: Bo Wang; Zhou Wu; Lei Liu
Journal: Complexity
 (2018)

---

## Abstract

This paper addresses a maintenance plan optimization problem for building energy retrofitting with a novel differential evolution-based 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.


---

## Body

## 1. 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–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–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–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. Section2 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.


## 2. 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 instantstk=kS with k=0,1,2,…,T, let tk 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 tkp and tkc denote the preventive and corrective maintenance instants, where kp=p1,p2,…,pNp and tkc, kc=c1,c2,…,cNc. Accordingly, Np and Nc denote the maintenance scale. The maintenance time schedule is decided via identifying the maintenance instants. Let Qp=p1,p2,…,pNp and Qc=c1,c2,…,cNc denote the indices of preventive and corrective maintenance instants, respectively. The elements of Qp and Qc are selected from k=0,1,2,…,T, implying that the maintenance instants are commensurate with sampling instants tk.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 employxtk to represent the variable population of the working items, and xtk the maintenance rates. The optimization is to find proper Np, Nc, Qp, and Qc, and corresponding utk, such that a set of optimal xtk can be achieved to maximize the performance indicators, e.g., the energy savings and IRR. The contents of xtk and utk 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 over0,TS:
(1)xtk+1=fxtk,k∈0,T;k≠ki,i=1,2,…,N,where f· is a continuously differentiable function. Let the system state xtki jump at instants tki,i=1,2,…,N, namely, the switching instants. Such jumps are represented by the following:
(2)xtki→=xtki+△xtki,i=1,2,…,N,xtki+1=fxtki→,where △xtki denotes an impulsive effect at tki that incurs a jump of the system state. tki→ is employed to describe the moment right after the impulsive effect. Accordingly, xtki denotes the system state prior to the impulsive jump and xtki→ the system state after. The system dynamics of a maintenance planning problem is established as the following:
(3)xl,Mltk+1=xl,Mltk−fMl,Ml−1lxl,Mltk,tk−gl,Mlxl,Mltk,tk,tk≠tkc,kc=c1,c2,…,cNcxl,Ml−1tk+1=xl,Ml−1tk+fMl,Ml−1lxl,Mltk,tk−fMl−1,Ml−2lxl,Ml−1tk,tk−gl,Ml−1xl,Ml−1tk,tk,tk≠tkp,kp=p1,p2,…,pNp⋮xl,1tk+1=xl,1tk+f2,1lxl,2tk,tk−gl,1xl,1tk,tk,tk≠tkp=p1,p2,…,pNpxl,nt0=x0l,n,1≤n≤Ml.Equation (3) describes such a homogeneous retrofitted item group l that contains Nl items with Ml working states; Ml denotes the best working states, and 1 denotes the worst. The system state xl,itk, i=1,2,…,Ml denotes the numbers of working items that are corresponding with group l and state i. fi,jlxl,itk,tk denotes the proportion of items that degrades from state i to a worse state j over interval tk,tk+1. gl,ixl,itk,tk denotes the malfunctioning items from state i, namely, the population decay model. For the convenience of discussion, (3) is abbreviated to xltk+1=Fxltk, where xltk=xl,Mltk,xl,Ml−1tk,…,xl,1tk. The impulsive jumps are described as the following:
(4)△xl,Mltk=ulcxl,0tk,tk,tk=tkc,△xl,Mltk=∑i=1Ml−1ul,ixl,itk,tk,tk=tkp,△xl,Ml−1tk=−ul,Ml−1xl,Ml−1tk,tk,tk=tkp,⋮△xl,1tk=−ul,1xl,1tk,tk,tk=tkp,xl,itk→=xl,itk+△xl,itk,1≤i≤Ml,xltk+1=Fxltk→,xl,0tk=Nl−∑i=1Mlxl,itk,kc∈Qc,kp∈Qp.Similarly,△xl,itk denotes the impulsive effects and tk→ the moment post to jumps. ul,ixl,itkp,tkp, i=1,2,…,Ml−1 denotes the preventive maintenance rate and ulcxl,0tkc,tkc the corrective maintenance rate. xl,0tk denotes the number of malfunctioning items. The control of (3) and (4) is to identify proper Np and Nc, Qp and Qc, and ulc and ul,i, namely, a maintenance plan. Let ulptk denote the overall rate of preventive maintenance and ulctk the overall rate of corrective maintenance at tk:
(5)ulptk=ul,Ml−1xl,Ml−1tk,tk,…,ul,1xl,1tk,tk,tk=tkp,ulptk=0,tk≠tkp,ulctk=ulcxl,0tk,tk,tk=tkc,ulctk=0,tk≠tkc,According to (3), (4), and (5), there is a dependence of the maintenance plan and system states xl,itk; 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. LetN¯ 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:
(6)altk=al,1tk,al,2tk,…,al,Mltk,bltk=bl,1tk,bl,2tk,…,bl,Mltk,Cltk=C1,1tk,Cl,2wtk,…,Cl,Ml−1tk,Clctk,where altk denotes the average energy savings over sampling interval tk−1,tk and bltk the cost savings. al,itk and bl,itk are corresponding with xl,itk. Cltk denotes the maintenance costs per item, where C1,itk is corresponding with preventive maintenance rate ul,itk and Clctk the corrective maintenance rate ulctk. The overall energy savings over the sustainability period 0,TS are then represented by
(7)EStk=∑l=1N¯∑i=1Mlal,itkxl,itk,ES∣all=∑k=1TEStk,and the corresponding cost savings are obtained:
(8)Btk=∑l=1N¯∑i=1Mlbl,itkxl,itk,B∣all=∑k=1TBtk,and the maintenance cost at each time instant is obtained:
(9)htk=∑l=1N¯∑i=1Ml−1Cl,itkul,itk+∑l=1N¯Clctkulctk+Mctk,where Mctk denotes a constant cost that is charged when maintenance action takes place at tk. Such a cost refers to the expenditure of activate maintenance actions. The total expenditure of the retrofitting project is represented by
(10)h∣all=h0+∑k=1Thtk,where h0 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:
(11)NPV=∑k=1TBtk−htk1+dϕtk−1−h0,where d is the discount rate, d>0. ϕtk=1,2,3… indicates that the sampling instant tk lies within a specific year after the initial stage. Let dR∣T denote the IRR. It is solved by letting the following equation equals to 0:
(12)∑k=1TBtk−htk1+dR∣Tϕtk−1−h0=0,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 h0 [24].


### 2.4. Optimization Problem

LetP=P1,P2,…,PN¯ denote the initial populations of the respective item groups, up· denote the collection of ulpkp, and uc· denote the collection of ulckc, l=1,2,…,N¯. Taking advantage of the preceding impulsive and switched system modelling (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.


#### 2.4.1. Optimization ProblemℳPO

For a dynamical system consisting ofN¯ subsystems in (3), (4), and (5) with state variables xt=x1t,…,xlt,…,xN¯t and initial state x0, find the proper maintenance scale Np and Nc, maintenance time schedule Qp and Qc, PM rates up⋅=u1pkp,…,ulpkp,…,uN¯pkp, and CM rates uc⋅=u1ckc,…,ulckc,…,uN¯ckc, to minimize the following performance index
(13)JP,Qp,Qc,up⋅,uc⋅=−λ1ES∣allα−λ2dR∣T,subject to system dynamics (3), (4), and (5) and the following constraints,
(14)ES∣all≥α,∫t=0Thtdt≤β,Tp≤T′.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∣RT 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. Tp is the payback period, which is computed from the last tn that makes NPV<0. Let MN¯ denote the total number of working states for the N¯ item groups, up·∈ℛMN¯−N¯×Sp and uc·∈ℛN¯×Sc; 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, dR∣T is nonanalytic; as a result, the gradient method becomes infeasible to (13) and (14). Furthermore, Np and Nc 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.


## 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 instantstk=kS with k=0,1,2,…,T, let tk 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 tkp and tkc denote the preventive and corrective maintenance instants, where kp=p1,p2,…,pNp and tkc, kc=c1,c2,…,cNc. Accordingly, Np and Nc denote the maintenance scale. The maintenance time schedule is decided via identifying the maintenance instants. Let Qp=p1,p2,…,pNp and Qc=c1,c2,…,cNc denote the indices of preventive and corrective maintenance instants, respectively. The elements of Qp and Qc are selected from k=0,1,2,…,T, implying that the maintenance instants are commensurate with sampling instants tk.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 employxtk to represent the variable population of the working items, and xtk the maintenance rates. The optimization is to find proper Np, Nc, Qp, and Qc, and corresponding utk, such that a set of optimal xtk can be achieved to maximize the performance indicators, e.g., the energy savings and IRR. The contents of xtk and utk 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 over0,TS:
(1)xtk+1=fxtk,k∈0,T;k≠ki,i=1,2,…,N,where f· is a continuously differentiable function. Let the system state xtki jump at instants tki,i=1,2,…,N, namely, the switching instants. Such jumps are represented by the following:
(2)xtki→=xtki+△xtki,i=1,2,…,N,xtki+1=fxtki→,where △xtki denotes an impulsive effect at tki that incurs a jump of the system state. tki→ is employed to describe the moment right after the impulsive effect. Accordingly, xtki denotes the system state prior to the impulsive jump and xtki→ the system state after. The system dynamics of a maintenance planning problem is established as the following:
(3)xl,Mltk+1=xl,Mltk−fMl,Ml−1lxl,Mltk,tk−gl,Mlxl,Mltk,tk,tk≠tkc,kc=c1,c2,…,cNcxl,Ml−1tk+1=xl,Ml−1tk+fMl,Ml−1lxl,Mltk,tk−fMl−1,Ml−2lxl,Ml−1tk,tk−gl,Ml−1xl,Ml−1tk,tk,tk≠tkp,kp=p1,p2,…,pNp⋮xl,1tk+1=xl,1tk+f2,1lxl,2tk,tk−gl,1xl,1tk,tk,tk≠tkp=p1,p2,…,pNpxl,nt0=x0l,n,1≤n≤Ml.Equation (3) describes such a homogeneous retrofitted item group l that contains Nl items with Ml working states; Ml denotes the best working states, and 1 denotes the worst. The system state xl,itk, i=1,2,…,Ml denotes the numbers of working items that are corresponding with group l and state i. fi,jlxl,itk,tk denotes the proportion of items that degrades from state i to a worse state j over interval tk,tk+1. gl,ixl,itk,tk denotes the malfunctioning items from state i, namely, the population decay model. For the convenience of discussion, (3) is abbreviated to xltk+1=Fxltk, where xltk=xl,Mltk,xl,Ml−1tk,…,xl,1tk. The impulsive jumps are described as the following:
(4)△xl,Mltk=ulcxl,0tk,tk,tk=tkc,△xl,Mltk=∑i=1Ml−1ul,ixl,itk,tk,tk=tkp,△xl,Ml−1tk=−ul,Ml−1xl,Ml−1tk,tk,tk=tkp,⋮△xl,1tk=−ul,1xl,1tk,tk,tk=tkp,xl,itk→=xl,itk+△xl,itk,1≤i≤Ml,xltk+1=Fxltk→,xl,0tk=Nl−∑i=1Mlxl,itk,kc∈Qc,kp∈Qp.Similarly,△xl,itk denotes the impulsive effects and tk→ the moment post to jumps. ul,ixl,itkp,tkp, i=1,2,…,Ml−1 denotes the preventive maintenance rate and ulcxl,0tkc,tkc the corrective maintenance rate. xl,0tk denotes the number of malfunctioning items. The control of (3) and (4) is to identify proper Np and Nc, Qp and Qc, and ulc and ul,i, namely, a maintenance plan. Let ulptk denote the overall rate of preventive maintenance and ulctk the overall rate of corrective maintenance at tk:
(5)ulptk=ul,Ml−1xl,Ml−1tk,tk,…,ul,1xl,1tk,tk,tk=tkp,ulptk=0,tk≠tkp,ulctk=ulcxl,0tk,tk,tk=tkc,ulctk=0,tk≠tkc,According to (3), (4), and (5), there is a dependence of the maintenance plan and system states xl,itk; 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. LetN¯ 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:
(6)altk=al,1tk,al,2tk,…,al,Mltk,bltk=bl,1tk,bl,2tk,…,bl,Mltk,Cltk=C1,1tk,Cl,2wtk,…,Cl,Ml−1tk,Clctk,where altk denotes the average energy savings over sampling interval tk−1,tk and bltk the cost savings. al,itk and bl,itk are corresponding with xl,itk. Cltk denotes the maintenance costs per item, where C1,itk is corresponding with preventive maintenance rate ul,itk and Clctk the corrective maintenance rate ulctk. The overall energy savings over the sustainability period 0,TS are then represented by
(7)EStk=∑l=1N¯∑i=1Mlal,itkxl,itk,ES∣all=∑k=1TEStk,and the corresponding cost savings are obtained:
(8)Btk=∑l=1N¯∑i=1Mlbl,itkxl,itk,B∣all=∑k=1TBtk,and the maintenance cost at each time instant is obtained:
(9)htk=∑l=1N¯∑i=1Ml−1Cl,itkul,itk+∑l=1N¯Clctkulctk+Mctk,where Mctk denotes a constant cost that is charged when maintenance action takes place at tk. Such a cost refers to the expenditure of activate maintenance actions. The total expenditure of the retrofitting project is represented by
(10)h∣all=h0+∑k=1Thtk,where h0 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:
(11)NPV=∑k=1TBtk−htk1+dϕtk−1−h0,where d is the discount rate, d>0. ϕtk=1,2,3… indicates that the sampling instant tk lies within a specific year after the initial stage. Let dR∣T denote the IRR. It is solved by letting the following equation equals to 0:
(12)∑k=1TBtk−htk1+dR∣Tϕtk−1−h0=0,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 h0 [24].


## 2.4. Optimization Problem

LetP=P1,P2,…,PN¯ denote the initial populations of the respective item groups, up· denote the collection of ulpkp, and uc· denote the collection of ulckc, l=1,2,…,N¯. Taking advantage of the preceding impulsive and switched system modelling (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.


### 2.4.1. Optimization ProblemℳPO

For a dynamical system consisting ofN¯ subsystems in (3), (4), and (5) with state variables xt=x1t,…,xlt,…,xN¯t and initial state x0, find the proper maintenance scale Np and Nc, maintenance time schedule Qp and Qc, PM rates up⋅=u1pkp,…,ulpkp,…,uN¯pkp, and CM rates uc⋅=u1ckc,…,ulckc,…,uN¯ckc, to minimize the following performance index
(13)JP,Qp,Qc,up⋅,uc⋅=−λ1ES∣allα−λ2dR∣T,subject to system dynamics (3), (4), and (5) and the following constraints,
(14)ES∣all≥α,∫t=0Thtdt≤β,Tp≤T′.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∣RT 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. Tp is the payback period, which is computed from the last tn that makes NPV<0. Let MN¯ denote the total number of working states for the N¯ item groups, up·∈ℛMN¯−N¯×Sp and uc·∈ℛN¯×Sc; 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, dR∣T is nonanalytic; as a result, the gradient method becomes infeasible to (13) and (14). Furthermore, Np and Nc 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.


## 2.4.1. Optimization ProblemℳPO

For a dynamical system consisting ofN¯ subsystems in (3), (4), and (5) with state variables xt=x1t,…,xlt,…,xN¯t and initial state x0, find the proper maintenance scale Np and Nc, maintenance time schedule Qp and Qc, PM rates up⋅=u1pkp,…,ulpkp,…,uN¯pkp, and CM rates uc⋅=u1ckc,…,ulckc,…,uN¯ckc, to minimize the following performance index
(13)JP,Qp,Qc,up⋅,uc⋅=−λ1ES∣allα−λ2dR∣T,subject to system dynamics (3), (4), and (5) and the following constraints,
(14)ES∣all≥α,∫t=0Thtdt≤β,Tp≤T′.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∣RT 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. Tp is the payback period, which is computed from the last tn that makes NPV<0. Let MN¯ denote the total number of working states for the N¯ item groups, up·∈ℛMN¯−N¯×Sp and uc·∈ℛN¯×Sc; 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, dR∣T is nonanalytic; as a result, the gradient method becomes infeasible to (13) and (14). Furthermore, Np and Nc 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.


## 3. A Multiscale DE Solver

A general form of the multiscale problem is employed as the following:(15)minJx→,y→,y→=y1,y2,y3,…,yD∈ℛD,to represent the investigated optimal control problem. x→ and y→ constitute a candidate solution. x→ indicates an implicit vector with a variant scale Nx→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), x→ denotes the time schedule Qp, Qc and maintenance rate up·,uc·,y→ denotes the maintenance scale Np, Nc. The maximum possible scale of x→ can be a large number denoted as NN≠D. Apparently, x→∈ℛMN¯−N¯+1×Np×N¯+1×Nc and y→∈ℛ2. The constraints are realized by a stationary penalty function in the objective function:
(16)J=−λ1ES∣allα−λ2dR∣T+ω∑n=13max0,Pn.ω is a large positive coefficient that adjusts the penalty from the violation of constraints. The penalty functions Pn with n=1,2,3 are given as the following:
(17)Pn=α−ES∣all,n=1,∑k=1Thtk−β,n=2,−NPV∣0Tp,n=3.To solve such a multiscale minimization problem, we propose a multiscale-based approach, namely, multiple-scale 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, 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 includingMutation, 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 problemminfx→,x→=x1,x2,…,xN∈ℛN is adopted for discussion. The search space is constrained in ∏j=1NLj,Uj, Lj≤Uj, where Lj is the lower bound of search space and Uj the upper bound.At each generation, an individualx→i,G=x1,i,G,x2,i,G,…,xD,i,G (namely, a target vector) is mutated to create a mutant vector v→i,G=v1,i,G,v2,i,G,…,vD,i,G by the mutation operation. There are many possible mutation operators from literatures, such as the following:
(18)v→i,G=x→i,G+F⋅xbest,G−x→i,G+F⋅x→r1,G−x→r2,G.More mutation operators can be found in [25]. In the above equations, the indices r1 and r2 are distinct integers uniformly randomized from the set i=1,2,…,NP. x→r1,G−x→r2,G is a difference vector to mutate the corresponding target vector x→i,G. The parameter F is the mutation factor which usually ranges on the interval (0,1). xbest,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 tox→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:
(19)u→j,i,G=v→j,i,G,ifrand0,1≤CRorj=jrand,x→j,i,Gotherwise,where j=1,2,…,D. jrand is a random integer that ranges between 1,D. j=jrand implies that an element in vector u→j,i,G is chosen, i.e., the acceptance of the mutant vector v→i,G. rand0,1 is a random real number that ranges between [0,1]. For each pair of j and i, rand0,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:
(20)u→j,i,G=minUj,2Lj−u→j,i,G,ifu→j,i,G<Lj,maxLj,2Uj−u→j,i,G,ifu→j,i,G>Uj.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 ofx→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:
(21)x→i,G+1=u→i,G,iffu→i,G≤fx→i,G,x→i,G,otherwise.


### 3.2. Multiple-Scale Differential Evolution

In MSDE, several subpopulations constitute the whole population of candidate solutions. Each subpopulation corresponds with a specific explicit variabley→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. Thereby, the MSDE employs a parallel processing mechanism to improve the basic DE and is able to address multiscale problems.The procedure of MSDE can be illustrated in the following steps:(1)
A population ofNP 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
(2)
At each iteration, within subpopulationk, the basic DE is implemented to search a better implicit vector x→k,G subject to given explicit variable y→k,G, following (18), (19), (20), and (21). Let x→best,k,G denote the solution with best fitness value, namely, the local best
(3)
After everym iteration, namely, the shuffling period, comparing the local best from subpopulations. Such a number refers to the period that allows each subpopulation to evolve. The fitness of the local best is also identified to be the fitness of the corresponding explicit variable. Copy the best-performing explicit variable to the worst-performing subpopulation, and mutate explicit variables 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.The pseudocode of MSDE is given in Algorithm1, 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.Algorithm 1: Pseudocode of MSDE.
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; SF,CR denotes the set of mutation factor F and CR for each subpopulation.
1:BEGIN
2:G=1;
3: Initializen subpopulations (including n explicit variables and NP implicit vectors);
4: Evaluate the fitness for each individual;
5: Setm, SF,CR;
6:while the stopping criterion is not satisfied do
7:while k=1 to n do
8:    Randomly choose (Fk,CRk) from SF,CR;
9:    Find the local besty→k,G,x→best,k,G;
10:while i=1 to NS do
11:    Randomly choosex→r1,k,G and x→r2,k,G from the k-th subpopulation, x→r2,k,G≠x→r1,k,G≠x→i,k,G;
12:v→i,k,G=x→i,k,G+Fk⋅xbest,k,G−x→i,k,G+Fk·x→r1,k,G−x→r2,k,G;
13:while j=1 to D do
14:if rand0,1<CRk then
15:u→j,k,i,G=v→j,k,i,G;
16:else
17:u→j,k,i,G=x→j,k,i,G;
18:end if
19:end while
20:if fu→i,k,G,y→k,G≤fx→i,k,G,y→k,G then
21:x→i,k,G+1=u→i,k,G;
22:else
23:x→i,k,G+1=x→i,k,G;
24:end if
25:if fu→i,k,G,y→k,G≤fx→i,k,G,y→k,G then
26:x→i,k,G+1=u→i,k,G;
27:else
28:x→i,k,G+1=x→i,k,G;
29:end if
30:end while
31:end while
32:G=G+1;
33:if modG,m==0 then
34:  Find the best subpopulationkb,G and worst one kw,G;
35:y→kw,G,G=y→kb,G,G;
36:  Mutatey→k,G where k≠kb,G≠kw,G;
37:  Initialize subpopulationsk≠kb,G according to new y→k,G;
38:end if
39:end while
40: Find global best solutionybest, xbest;
41:Return xbest;
42:ENDFigure 1
The flowchart of MSDE.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 subpopulation has to adopt a same explicit variable of the best subpopulation.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.


## 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, 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 includingMutation, 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 problemminfx→,x→=x1,x2,…,xN∈ℛN is adopted for discussion. The search space is constrained in ∏j=1NLj,Uj, Lj≤Uj, where Lj is the lower bound of search space and Uj the upper bound.At each generation, an individualx→i,G=x1,i,G,x2,i,G,…,xD,i,G (namely, a target vector) is mutated to create a mutant vector v→i,G=v1,i,G,v2,i,G,…,vD,i,G by the mutation operation. There are many possible mutation operators from literatures, such as the following:
(18)v→i,G=x→i,G+F⋅xbest,G−x→i,G+F⋅x→r1,G−x→r2,G.More mutation operators can be found in [25]. In the above equations, the indices r1 and r2 are distinct integers uniformly randomized from the set i=1,2,…,NP. x→r1,G−x→r2,G is a difference vector to mutate the corresponding target vector x→i,G. The parameter F is the mutation factor which usually ranges on the interval (0,1). xbest,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 tox→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:
(19)u→j,i,G=v→j,i,G,ifrand0,1≤CRorj=jrand,x→j,i,Gotherwise,where j=1,2,…,D. jrand is a random integer that ranges between 1,D. j=jrand implies that an element in vector u→j,i,G is chosen, i.e., the acceptance of the mutant vector v→i,G. rand0,1 is a random real number that ranges between [0,1]. For each pair of j and i, rand0,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:
(20)u→j,i,G=minUj,2Lj−u→j,i,G,ifu→j,i,G<Lj,maxLj,2Uj−u→j,i,G,ifu→j,i,G>Uj.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 ofx→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:
(21)x→i,G+1=u→i,G,iffu→i,G≤fx→i,G,x→i,G,otherwise.


## 3.2. Multiple-Scale Differential Evolution

In MSDE, several subpopulations constitute the whole population of candidate solutions. Each subpopulation corresponds with a specific explicit variabley→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. Thereby, the MSDE employs a parallel processing mechanism to improve the basic DE and is able to address multiscale problems.The procedure of MSDE can be illustrated in the following steps:(1)
A population ofNP 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
(2)
At each iteration, within subpopulationk, the basic DE is implemented to search a better implicit vector x→k,G subject to given explicit variable y→k,G, following (18), (19), (20), and (21). Let x→best,k,G denote the solution with best fitness value, namely, the local best
(3)
After everym iteration, namely, the shuffling period, comparing the local best from subpopulations. Such a number refers to the period that allows each subpopulation to evolve. The fitness of the local best is also identified to be the fitness of the corresponding explicit variable. Copy the best-performing explicit variable to the worst-performing subpopulation, and mutate explicit variables 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.The pseudocode of MSDE is given in Algorithm1, 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.Algorithm 1: Pseudocode of MSDE.
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; SF,CR denotes the set of mutation factor F and CR for each subpopulation.
1:BEGIN
2:G=1;
3: Initializen subpopulations (including n explicit variables and NP implicit vectors);
4: Evaluate the fitness for each individual;
5: Setm, SF,CR;
6:while the stopping criterion is not satisfied do
7:while k=1 to n do
8:    Randomly choose (Fk,CRk) from SF,CR;
9:    Find the local besty→k,G,x→best,k,G;
10:while i=1 to NS do
11:    Randomly choosex→r1,k,G and x→r2,k,G from the k-th subpopulation, x→r2,k,G≠x→r1,k,G≠x→i,k,G;
12:v→i,k,G=x→i,k,G+Fk⋅xbest,k,G−x→i,k,G+Fk·x→r1,k,G−x→r2,k,G;
13:while j=1 to D do
14:if rand0,1<CRk then
15:u→j,k,i,G=v→j,k,i,G;
16:else
17:u→j,k,i,G=x→j,k,i,G;
18:end if
19:end while
20:if fu→i,k,G,y→k,G≤fx→i,k,G,y→k,G then
21:x→i,k,G+1=u→i,k,G;
22:else
23:x→i,k,G+1=x→i,k,G;
24:end if
25:if fu→i,k,G,y→k,G≤fx→i,k,G,y→k,G then
26:x→i,k,G+1=u→i,k,G;
27:else
28:x→i,k,G+1=x→i,k,G;
29:end if
30:end while
31:end while
32:G=G+1;
33:if modG,m==0 then
34:  Find the best subpopulationkb,G and worst one kw,G;
35:y→kw,G,G=y→kb,G,G;
36:  Mutatey→k,G where k≠kb,G≠kw,G;
37:  Initialize subpopulationsk≠kb,G according to new y→k,G;
38:end if
39:end while
40: Find global best solutionybest, xbest;
41:Return xbest;
42:ENDFigure 1
The flowchart of MSDE.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 subpopulation has to adopt a same explicit variable of the best subpopulation.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.


## 4. Simulation Results and Analysis

In order to verify the effectiveness of 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(22)xcfltk+1=bcflxcfltk−bcflccflxcfl2tkxcfl0,where xcfl denotes the population of the CFL group, and xcfl0 denotes the initial group, xcflt0=xcfl0. bcfl and ccfl are the parameters of the decay model. The parameters can be estimated by experimental data fitting. More details about decay model (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
(23)fi,i−1a,cxac,itk,tk=xac,itk1−eςac,ip,i=3,2,gac,ixac,itk,tk=xac,itk1−e−ςac,ic,i=3,2,1,where xac,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. xac,3t0=xac0, where xac0 denotes the initial population of the air conditioner group. xac,3t0=xac,1t0=0. fi,i−1ac denotes the amount of air conditioners that transit from working state i to i−1, and gac,i denotes the number of items that become malfunctioning from state i. The parameters ζac,ip and ςac,ic 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 bcfl, ccfl, ζac,ip, and ςac,ic in the case study.Table 1
Parameters for population decay and state transition models.
Retrofits
bcfl
ccfl
ζac,ip
ζac,ic
15 W retrofit CFL
0.0947
0.775
N/A
N/A
Air conditioneri=3
N/A
N/A
0.095
0.029
Air conditioneri=2
N/A
N/A
0.04
0.033
Air conditioneri=1
N/A
N/A
N/A
0.05Table2 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 pre-retrofit 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 h0 is $20,692.Table 2
Characteristics of retrofitted items.
Retrofits
Quantities
Unit price ($)
Unit energy saving (kWh)
Unit cost saving ($)
Preventive cost ($)
Corrective cost ($)
15 W retrofit CFL
338
14
105.6
11.9
N/A
14
Air conditioneri=3
42
380
4320
486.65
N/A
175
Air conditioneri=2
0
380
3542.3
397.95
52
N/A
Air conditioneri=1
0
380
2651.75
278.35
70
N/ASome 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.4 kWh. 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 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. Table3 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 further improve the maintenance planning performances, as a proper scale can imply a better solution in many cases.Table 3
Performances of maintenance plans with different budget limits.
Cases
Budget limit ($)
Energy savings (kWh)
Percentage saved over target
IRR
Payback period (months)
NPV ($)
Maintenance cost ($)
Total investment ($)
Fixed schedule
20,000
1,167,215
11.99%
72.15%
18.72
51427.78
19,986
40,678
Optimal schedule
20,000
1,190,464
14.22%
76.23%
16.92
52949.3
19,998
40,690
Fixed schedule
40,000
1,670,773
60.30%
79.20%
18.12
68393.86
39,987
60,679
Optimal schedule
40,000
1,702,041
63.31%
79.27%
17.64
70076.16
39,986
60,678
Fixed schedule
65,000
1,829,741
75.56%
78.81%
18.64
70608.34
52,382
73,074
Optimal schedule
65,000
2,001,191
92.01%
80.44%
17.28
77670.84
58,046
78,738Table 4
Optimal maintenance time schedules.
Time schedule
Budget limit ($)
Maintenance number
Maintenance instants
Qp
20,000
6
{14,26,38,50,62,74,86}
Qc
20,000
5
{14,26,38,50,62,74}
Qp
40,000
12
{8,18,28,…,98}
Qc
40,000
12
{8,18,28,…,98}
Qp
65,000
34
{5,8,11,…,104}
Qc
65,000
34
{8,11,14,…,107}Figure 2
The maintenance rates over time with different maintenance budgets.Figure 3
Energy saving trajectories of the optimal time schedule and maintenance plan with different maintenance budgets.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 similar formulation.


## 5. 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.



---

*Source: 3237471-2018-10-15.xml*