Distributed State Estimation Using a Modified Partitioned Moving Horizon Strategy for Power Systems

In this paper, a distributed state estimation method based on moving horizon estimation (MHE) is proposed for the large-scale power system state estimation. The proposed method partitions the power systems into several local areas with non-overlapping states. Unlike the centralized approach where all measurements are sent to a processing center, the proposed method distributes the state estimation task to the local processing centers where local measurements are collected. Inspired by the partitioned moving horizon estimation (PMHE) algorithm, each local area solves a smaller optimization problem to estimate its own local states by using local measurements and estimated results from its neighboring areas. In contrast with PMHE, the error from the process model is ignored in our method. The proposed modified PMHE (mPMHE) approach can also take constraints on states into account during the optimization process such that the influence of the outliers can be further mitigated. Simulation results on the IEEE 14-bus and 118-bus systems verify that our method achieves comparable state estimation accuracy but with a significant reduction in the overall computation load.


Introduction
Power system state estimation (PSSE) plays an indispensable part in the power industry [1]. One common centralized approach named the weighted least squares (WLS) has been widely used for PSSE, employing a nonlinear measurement model. In recent years, phasor measurement units (PMUs) have drawn much attention because they can provide voltage and current phasors and the measurement model becomes linear [2,3]. The computation complexity will become simpler.
Even though PMUs can result in measurements with higher accuracy, the PSSE still constitutes a major challenge due to the presence of bad data or outliers [3][4][5]. Such outliers that are far away from the expected measuring data create the potential risk of misleading the estimated result [6]. The WLS, primarily based on a single snap shot of measurements, is not a robust method and it may lead to a biased estimated result even when a single bad measurement occurs [3,7]. One common way to alleviate the influence of outliers is to use more snap shots of measurements. Moreover, in order to reduce the influence of outliers further, the largest normalized residuals (LNR) test [8] is usually used in WLS to deal with bad data. Many other centralized algorithms based on the WLS have been proposed in [9][10][11]. Due to the increasing number of PMUs installed in substations, it is reasonable to assume that the power systems will be only measured by PMUs in the near future [12]. The authors in [13] use an adaptive approach in updating the accuracies of the PMUs while employing local decision metrics and an Internet of Things (IoT) paradigm. It allows considering adaptive values of the accuracies of different measurement devices in the estimation process. However, the algorithm is basically solving a centralized state estimation problem for the distribution system. In [3,4], the least absolute value (LAV) has been proposed to handle the outliers where the measurements are collected from PMUs. However, the algorithms described so far are mostly based on the measurement model only. A robust estimator via the moving horizon strategy has been developed for PSSE in [14]. The re-weighted moving horizon estimation (MHE) aims to solve an optimization problem at each time instant by using a limited amount of the most recent information, and the objective function of re-weighted MHE includes the measurement model error, the process model error and the error in the state estimate at the beginning of the window. Moreover, the constraints on states have also been exploited. By having these constraints in the optimization process, it is more robust to the outliers. The price to pay is an increased complexity since the computational load will increase. A modified MHE (mMHE) has been developed in [15]. In this reference, the objective function of the mMHE consists only of the sensor model error and the error from the prior estimate. The state at the beginning of the window is obtained in the optimization, instead of all states in the window under the scheme of MHE. Therefore, the mMHE will have a higher computing speed than the MHE. The tradeoff is the estimated accuracy of mMHE is a little smaller than that of MHE, but their difference is insignificant.
Due to the fact that power systems become larger and increase in complexity, the centralized estimators that process the measurements from the whole grid may no longer be feasible [16,17]. Based on the rapid growth of usage of the wide-area monitoring systems for modern power grids, many approaches on distributed state estimation are proposed. A fully distributed state estimation approach is presented whereby each local area solves the system-wide states [18]. The authors in [9] propose a fully decentralized adaptive re-weighted state estimation algorithm for hybrid PSSE, where both measurements from the supervisory control and data acquisition (SCADA) system and PMUs are used. However, for the two algorithms mentioned above, a large amount of iterations are required when the power systems are large. A distributed state estimation for power systems with linear models has been proposed and the alternating direction method of multipliers (ADMM) is used to solve the optimization problem [19]. A distributed robust bilinear state estimation method is proposed to multi-area power systems with nonlinear measurements where interregional communication is required but a central coordinator is not necessary [20]. One drawback of the distributed algorithms described so far is that they are largely based on just the measurement model and fall into the WLS category. The authors in [21] propose a new distributed framework where the process model is used. The average consensus algorithm is applied such that the substations can maintain the global state through information exchange with neighbors. However, for all the approaches mentioned above, the constraints on states are not considered during the optimization process. This may lead to a suboptimal solution. A distributed moving horizon estimation (DMHE) for power systems has been proposed in [22]. It is suitable for advanced applications such as wide-area monitoring and control that require the system-wide states to be available to all the regional transmission organizations (RTOs) [23]. However, the computational load for each local area is still heavy. The authors in [24] propose a distributed state estimation method for linear systems, and it is known as the partitioned moving horizon estimation (PMHE). The PMHE is more reasonable and suitable for large-scale systems monitoring because each local area (or subsystem) solves for the local states via a smaller optimization problem and the computational load is smaller whereby measurements are only sent to the local estimator but not to the centralized estimator, so a large amount of communication burden would be saved. Meanwhile, information is exchanged among neighboring areas only so the communication load is also small.
In this paper, considering the high speed of mMHE, the accuracy of MHE, and the advantage of PMHE to implement the MHE in a distributed way, a distributed state estimation method named the modified PMHE (mPMHE) is proposed and implemented for PSSE. Under the scheme of mPMHE, the error from the process model is ignored in the objective function. The proposed mPMHE has the following advantages:

•
Instead of all the states in the window estimated by PMHE, the mPMHE only estimates the state vector at the beginning of the window so it is faster than PMHE. The estimated precision of mPMHE is slightly lower than that of PMHE, but their difference is insignificant. The mPMHE achieves comparable state estimation accuracy but with a significant reduction in the overall computation load.

•
It is a distributed algorithm and is suitable for large-scale PSSE. Each local area solves for its own local states by using the local measurements and the estimated results from neighboring areas, so the computation load is small. In addition, the communication load is also small because the information is exchanged among neighboring areas only.

•
Constraints are taken into account during the optimization process and it is robust to outliers. Hence, good estimated results could be obtained.
This paper is organized as follows. The centralized state estimation is briefed in Section 2 and the mPMHE approach is proposed in Section 3. The simulations on the IEEE 14-bus and 118-bus systems are shown in Section 4 and conclusions are made in Section 5 respectively.

Measurement Model and State Equation
The linear measurement model based on the PMUs [3] is given as where H is the measurement matrix and k is the time step. This paper uses rectangular coordinates. z ∈ R m is the measurement vector composed of the real and imaginary components of the bus voltage (or the line current) phasors. The state vector x ∈ R n includes the real and imaginary parts of the voltage phasors. v ∈ R m is assumed to be noise with covariance R. The standard deviation of i-th measurement noise is denoted as σ i . In this paper, the following assumption is held: Assumption 1. The power system is observable so the matrix G = H T H is full rank.
The following simplified process model is considered for the state estimation [25,26]: where A is assumed to be an identity matrix according to [26] and w k represents the zero-mean disturbance with covariance Q.

Weighted Least Squares (WLS)
The WLS is an iterative algorithm that is applied to measurements including the power injections and power flows are usually collected from the SCADA system. However, the iterations are not necessary when PMUs are used. A traditional power system may be considered as a quasi-static system [25] because load demands change slowly and hence the state changes slowly. The sampling time of PMU measurements is usually in the order of milliseconds while the estimates are usually updated once every few minutes if the measurements are collected from SCADA [27]. In order to alleviate the influence of bad measurements, a total horizon length of N + 1 PMU measurements are used. During this interval, it is assumed that the system state is constant. The state x t can be estimated by solving the following cost function: where e i,k is the i-th measurement residual at time step k, Equation (4) gives Differentiating the above cost function (3) with respect tox t , where Using (4), Equation (5) can be written as Under Assumption 1, the matrixH T WH is also full rank and is an invertible matrix. Next, set Ψ(E) = 0, then from (6) the estimation of x t is given bŷ Note that the WLS is not a robust estimator and the largest normalized residuals (LNR) method [8] is usually used in WLS to deal with bad data. The normalized residuals are calculated as follows: The normalized residuals e norm i,k are calculated according to the residual covariance matrixΩ and measurement residual e i,k . If the normalized residuals e norm i,k are larger than a pre-determined threshold, the largest one will correspond to the bad measurement Z bad i . Once the largest normalized residual is found, the corresponding measurement is updated: The states will then be recalculated based on the updated measurements Z new i . Several iterations may be needed in order to make sure that all normalized residuals are less than the pre-determined threshold, for example, 3.0 [3].

Moving Horizon Estimation (MHE)
The common WLS estimator for power system state estimation has been discussed in previous subsection. However, the constraints on the states are not considered and this may lead to suboptimal estimates. In this subsection, the MHE algorithm in [28] is applied to the PSSE problem.
where X is a set of constraints defined by linear inequalities. In the following, the notation t − N|t − N − 1 refers to the time step for prediction from step t − N − 1 to t − N. Denote · as the Euclidean norm of a vector and · 2 S as the square of the weighted Euclidean norm of a vector, The objective function of MHE in (8) at time step t is given by where the arrival cost Φ t−N is given as For the arrival cost (10), we calculate P t−N|t−N−1 from P t−N−1|t−N−2 using the equation derived in [29]: The objective function of MHE includes three error terms: (i) the error between the measurement and sensor model prediction by (2); (ii) the error between the estimated state and its process model prediction by (1); (iii) and the error between the initial statex t−N in the horizon and the a priori state estimatex t−N|t−N−1 .

Modified Moving Horizon Estimation (mMHE)
In this subsection, we briefly review the mMHE method presented in [15]. The mMHE is an approach that compromises between the computational complexity and the estimated accuracy.
The objective function for mMHE includes two terms, a prior term given in (10) and the measurement error term given in (3), Note that the MHE approach needs to solve the problem with a higher dimension (N + 1)n and also acquires all the states in the current window,x t−N , . . . ,x t . The objective function of WLS given in (3) only includes the first error term of MHE and its dimension is n. The mMHE solves a problem with the same dimension as (3) so the computational time is also faster than MHE. Even though the errors from the process model are not included in the objective function of mMHE, there still exists some uncertainty and Q should be combined in (11).
Driven by the increasing demand of wide-area system monitoring and the high speed of mMHE, a distributed approach named the modified partitioned MHE (mPMHE) is proposed in the next section.

Modified Partitioned Moving Horizon Estimation (mPMHE)
In the previous section, the WLS, MHE and mMHE are implemented under the centralized setup, where all PMU measurements are collected and then sent to a control center. However, it may not be feasible in practice when a power grid is large-scale. In this section, on the basis of PMHE proposed in [24], we develop a distributed state estimation approach named the modified PMHE (mPMHE) for large-scale power system monitoring, where each local area estimates its local states based on local measurements and information exchanges among the neighboring areas. Moreover, the mPMHE takes less time than PMHE. The tradeoff is a decrease in the estimated accuracy. The convergence of mPMHE follows that of PMHE presented in [24].

mPMHE Problem Formulation
Let models (1) and (2) be partitioned into control areas with non-overlapping states [24]: . The global vectors are We assume that the power system partitioning is based on the following assumption: In this paper, the proposed mPMHE-i at time t is defined as: subject tox The local cost function in (18) is given by whereH is quite sparse and the second term on the right-hand side of (20) depends only on the neighboring areas.

Update Matrices
Denote the local observability matrix as O T ] T and 0 as the matrix of zero elements. We follow the method given by [24] to update Π .

Simulation Results
In this section, simulations on the IEEE 14-bus and 118-bus systems using the mMHE and mPMHE algorithms will be presented. We also illustrate the effect of the number of PMUs installed in the power systems and consider two scenarios: one with redundant observations on selected buses and one with a minimum number of PMUs for full topological observation.

Redundant Observations
In this example, the IEEE 14-bus system is shown in Figure 1 where the PMUs are placed according to [26]. Fifty-eight measurements, z i , i = 1, . . . , 58, consisting of 12 voltages (i = 1, . . . , 12) and 46 currents (i = 13, . . . , 58) are taken at each time instant k. When the horizon length N + 1 increases, the computation load of MHE and PMHE will also increase due to the increasing dimension of the optimization problem arising from MHE and PMHE. However, the level of estimation accuracy will be higher when the horizon increases. It should be noted that the horizon length can be chosen according to the required estimation accuracy. In this paper, for simplicity, the measurements with horizon length 3, i.e., k = t − 2, t − 1, t are used to give one set of estimates. According to [1], the WLS estimator usually uses one snap shot of measurements (the horizon length is 1) to estimate, in this paper we also show the result for comparison and it is represented by "WLS(1)" in the following figures. The measurement matrix H is calculated according to the parameters in [30]. There are n = 28 states in the vector T of (1) in which V r i and V im i (j = 1, . . . , n 2 ) are the real and imaginary parts of the bus voltage phasors, respectively. A is simplified as an identity matrix according to [25,26]. The IEEE 14-bus system is partitioned into four non-overlapping areas where each local area is locally observable. Due to the partitioned type, information is only exchanged with its immediate neighbors and the communication scheme is given in Figure 2. The measurements allocated for each area are given in Table 1, where I r i,j and I im i,j define the real and imaginary components of current phasors from bus i to bus j respectively.
and denote the Average of Mean Square Error (AMSE) in local area i for k ∈ [t c , t] to evaluate the estimation accuracy: where t c is the converging time step. An n-dimensional column vector comprising of all ones (or zeros) is denoted as 1 n (or 0 n ). I n denotes an n × n identity matrix. The initialization parameters of mPMHE algorithm are listed as follows: • The initial state vectors x The initial covariance matrices: P 1|0 = 10 3 I 28 , Π The noise covariances: Q The horizon length: . . , 14. Two different cases which consider the measurements including Gaussian and non-Gaussian noises are presented as follows: Case 1: Measurements with Gaussian noise In this case, the measurement noises are assumed to be Gaussian. The standard deviation of the voltage phasor measurement is arbitrary chosen as σ i = 0.005, i = 1, . . . , 12 and that of the current phasor measurements is set as σ i = 0.01, i = 13, . . . , 58, according to [31].
The simulations are performed using MATLAB version R2012b on a Windows 10 computer configured with Intel R Core TM , CPU i7-4500U, 1.80 GHz and 8 GB RAM, where the quadratic program problems arising from mMHE and mPMHE are solved by the alternating direction method of multipliers (ADMM), following the method presented in [14]. The constraints on the state variables are taken into account both in the MHE, mMHE, the PMHE in [24] and the proposed mPMHE to handle outliers. The comparison of AMSE values obtained from different estimators are given under the "Gaussian" column in Table 2. Even though the WLS (with horizon length 1) spends the least time, 0.3 ms, its AMSE value however is the largest. That is because the WLS estimator does not need to perform any iteration to obtain the estimated results. However, the WLS is not a robust estimator. When the horizon length is set as 3, the WLS takes more time, 1.6 ms, to obtain the estimated results. The AMSE value of WLS with LNR is 1.6 × 10 −3 and it takes 7.7 ms. The LAV estimator is built up according to the method provided by [3] and is conducted based on the matlab subfunction provided in GUROBI example. The AMSE value of LAV is equal to that of WLS with LNR, but the LAV spends more time than WLS with LNR in this case. It is significant to note that the MHE gets the highest accuracy (i.e., having the smallest AMSE 1.2 × 10 −3 ), compared with that of the WLS and WLS with LNR. The AMSE of mMHE with constraints is 1.3 × 10 −3 , which is quite close to the result of MHE. However, the mMHE with constraints only takes 6.6 ms while the MHE with constraints takes 11.8 ms. The mMHE takes 44% less time than that of MHE. The traditional PMHE gets better results than the mPMHE, but the mPMHE takes 2.6 ms and it is faster than PMHE, 4.9 ms. The time taken by the has a decreased percentage of 47% compared with PMHE, and a reduction of 60% compared with the centralized mMHE. Even though the mPMHE sacrifices slight estimated accuracy, its computation reduction is significant compared with that of PMHE. Figure 3 shows the convergence of mPMHE with constraints and the converging time step is around 20 time steps. The mPMHE and PMHE finally converges to the centralized results respectively. The convergence rate of mPMHE is quite close to that of PMHE but the computation load reduction of mPMHE is significant. Moreover, the estimated results of some buses from time step 10 to 80 are highlighted in Figure 4 in order to see the difference among different estimators. In order to test the proposed estimator affected by the high-magnitude outliers, suppose the outliers occur in the real part of current measurement I r 43 at time steps 40 and 50, as shown in Figure 5a. From Figure 5b, we can see that the WLS is seriously affected by the outliers, while the MHE, mMHE, PMHE and mPMHE can deal with the bad data and can still get good estimated results.

Case 2: Measurements with non-Gaussian noise
In order to verify that the proposed estimator can also deal with more outliers, the measurement noises are assumed to be non-Gaussian. For the voltage measurements z i,k , i = 1, . . . , 12, the noise v i is associated with the probability density function as follows: where σ i = 0.005. The first term accounted for 97% in f i (v i ) and the second term having a larger standard deviation is assumed to be the outliers [32]. For the current measurements z i,k , i = 13, . . . , 58, the noise v i is associated with the probability density function as follows [33] (22) where σ i = 0.01. The 3% of uniform distribution in (22) is useful for modeling initial conditions, disturbances and measurement errors that are equally likely to occur anywhere within a given interval.
According to the results shown in Figure 6, both PMHE and mPMHE also converge to the centralized results. However, the MSE is a little larger and it is not as smooth as that shown in Figure 3. That is because many outliers have been included in the measurements. The convergence time seems slightly longer and it is around 30 time steps according to Figure 7. Referring to the AMSE and average time under the "non-Gaussian" column in Table 2, the AMSE values are larger than those under Gaussian assumption, due to the occurrence of outliers. The WLS is not a robust estimator so its AMSE is the highest, 3.7 × 10 −3 . The WLS with LNR can detect the bad data, updated the relevant measurements and then the weights of the outliers can be mitigated. More iteration steps are needed for the WLS with LNR so the time increases to 14.0 ms. The AMSE of WLS with LNR is 1.8 × 10 −3 . With the inclusion of the constraints, the MHE can handle the outliers and gets the minimum value, 1.7 × 10 −3 . The AMSE of the mMHE with constraints is 1.8 × 10 −3 , while it only takes 7.3 ms but the MHE takes 15.9 ms to get one set of estimated results. The decreasing percentage of mMHE with constraints is 54%. The AMSE obtained from the distributed approaches finally converge to the centralized results, i.e., 1.7 × 10 −3 for the PMHE with constraints and 1.8 × 10 −3 for the mPMHE with constraints. Their convergence rates are still quite similar. The estimated error of mPMHE is slightly larger than that of PMHE but their difference is insignificant. However, the average computation time of the distributed approach (mPMHE with constraints) is reduced further, 3.7 ms, and it has a reduction of 43% compared with the PMHE with constraints (6.5 ms). The mPMHE with constraints has much bigger reduction percentages, 74% and 69%, compared with the common methods, WLS with LNR (14.0 ms) and LAV (12.0 ms), respectively.

Full Observation with Minimum Number of PMUs
In this example, the IEEE 14-bus system has been installed with minimum number of PMUs. If the PMUs are only installed at Buses 2, 6, 7 and 9 according to [26], Assumption 2 is not guaranteed when the local areas is partitioned as the one shown in Figure 1. A new partitioned type is required, as shown in Figure 8a and the relevant communication scheme is shown in Figure 8b.
A total number of 38 measurements, consisting of eight voltages and 30 currents are taken at each time and the measurements with horizon length 3 are used to give one set of estimates. The standard deviations of the voltage and current measurement noises follow those given in the previous subsection, and the MSE of different estimators under Gaussian and non-Gaussian assumptions are shown in Figure 9a,b respectively. It is significant that the converging time of mPMHE and PMHE with constraints is around 35 time steps. The AMSE and average time are given in Table 3. Due to the smaller number of measurements, the AMSE values are larger than those in Table 2, while the average time is less than those in Table 2. No matter under the Gaussian or non-Gaussian noise assumption, the AMSE values of mPMHE is equal to that of mMHE and a little larger than that of MHE, but the average computation time is much smaller than those for LAV, MHE and mMHE. For example, under the non-Gaussian noise assumption, the average computation time of mPMHE is 2.8 ms, and it has reduction percentages of 47% and 35%, compared with the centralized estimator (LAV, 5.3 ms) and the distributed state estimation method (PMHE, 4.3 ms).

Redundant Observations
In order to verify that the proposed algorithm is also effective in large-scale system, the IEEE 118-bus system is used and is partitioned into six local areas (subsystems) as shown in Figure 10. The relevant communication scheme is simplified as shown in Figure 11. The PMUs are placed according to [26] where 54 PMUs are used. A total number of 108 voltage measurements with the noise pdf in (21) with σ i = 0.005, i = 1, . . . , 108 and 366 current measurements with the noise pdf in (22) with σ i = 0.01, i = 109, . . . , 474 are considered. Area 1 has 46 measurements, Area 2 has 78 measurements, Area 3 has 76 measurements, both Area 4 and 5 have 110 measurements respectively, and Area 6 has 54 measurements. Every local area is locally observable. The initialization parameters of mPMHE algorithm in the IEEE 118-bus system are listed as follows: • The initial state vectors x The initial covariance matrices: Π According to the MSE result shown in Figure 12 and the estimated results shown in Figure 13, it is clear that the estimated results of PMHE and mPMHE converge to the centralized results obtained from MHE and mMHE. The convergence rate of mPMHE is still quite close to that of PMHE and the convergence cannot be expected in fewer than 400 time-steps. The AMSE and average computation time (per step) are given in Table 4. The mMHE still takes less time than that of MHE. The mPMHE with constraints can still get good estimated result and it takes the least time, 32 ms, compared with other robust estimators. The mPMHE reduces about 95% execution time compared with the centralized method (mMHE). Moreover, the mPMHE with constraints has a percentage decrease of 42% compared with that of PMHE (55 ms).

Full Observation with Minimum Number of PMUs
In this case, a minimum number of 32 PMUs are placed in the IEEE 118-bus system according to [26], and the whole system is partitioned into six local areas (subsystems) as shown in Figure 10. In order to make sure that every local area is locally observable, Bus 13 needs to be placed into Area 1 (Sub 1) and Buses 19 and 33 should be returned to Area 3 (Sub 3). The communication scheme is the same as that shown in Figure 11. A total number of 64 voltage measurements with the noise pdf in (21) with σ i = 0.005, i = 1, . . . , 64 and 250 current measurements with the noise pdf in (22) with σ i = 0.01, i = 65, . . . , 314 are considered.
According to the MSE result shown in Figure 14, it is clear that the estimated results of PMHE and mPMHE converge to the centralized results obtained from MHE and mMHE. The convergence rate of mPMHE is still quite close to that of PMHE and it is faster than that shown in Figure 12. This verifies that the structure of matrix H affects the convergence rate and the details of convergence property can be found in [24]. Compared with those under the redundant observations, the larger AMSE values and smaller average computation time (per step) will be obtained under the "Observation with minimum number of PMUs" column in Table 4. The mPMHE with constraints can still get good estimated result and it reduces about 62% execution time compared with the centralized method (LAV). Moreover, the mPMHE with constraints has a percentage decrease of 36% compared with that of PMHE (33 ms).
According to the simulation results on the IEEE 14-bus and 118-bus systems, given in Tables 2 and 4, it is clear that the average computation time of mPMHE increases, due to the dimension of local optimization problem in the IEEE 118-bus system which is larger than the IEEE 14-bus system. Therefore, the computation burden depends on the size of partitioned local areas and it would be scalable even when the size of a power system is larger than the IEEE 118-bus system. The computational load can be reduced if the larger power system is partitioned into more local areas (subsystems) and the dimension of each local optimization problem is smaller than that in the IEEE 118-bus system.  Remark: The moving horizon estimation (MHE), modified MHE (mMHE), partitioned moving horizon estimation (PMHE) and modified PMHE (mPMHE) take constraints into account during the optimization process.

Conclusions
Based on the wide-area monitoring systems, a fully distributed state estimation (mPMHE) based on the moving horizon estimation is proposed for power system state estimation, where each local area solves a smaller optimization problem to estimate its own local states by using the local measurements and the estimated results from its neighboring areas. The computation load is reduced compared with that of the centralized methods. In contrast with PMHE, the error from the process model is ignored in our proposed method. The estimated precision of mPMHE is slightly lower than that of PMHE but their difference is insignificant. The mPMHE achieves comparable state estimation accuracy but with a significant reduction in the computation load. Acknowledgments: This work was supported by the Singapore National Research Foundation (NRF) under its Campus for Research Excellence And Technological Enterprize (CREATE) programme, and Cambridge Centre for Advanced Research in Energy Efficiency in Singapore (CARES), C4T project.
Author Contributions: Tengpeng Chen made substantial contributions in proposing the algorithm, running simulations, analysing the data and manuscript preparation. Yi Shyh Eddy Foo and Xuebing Chen helped with the editing of the writing, and analyzed the simulation tests. K.V. Ling supervised this research.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: