A distributed network-based economic model predictive control for networked smart energy systems considering communication network constraints

Generally, the modern networked energy systems are constructed from some coupled sub-systems interconnected via their states and share their information through a communication network featuring time delays. This article suggests an efﬁcient multi-agent based cooperative-distributed networked economic model predictive control scheme for energy consumption optimisation and reducing the cost of consumed electricity in network-based smart energy systems considering communication network delays. In this regard, a buffer-based moving horizon strategy with estimator is suggested to estimate the own states of every sub-system and the coupled states, that is, the states in which their information exchange among sub-systems and their values are not accessible due to communication time delays. Moreover, the boundedness of the estimation error as well as the stability of the closed-loop system are established. The usability and proﬁciency of the suggested scheme are proved by applying the developed approach on two case studies.

as control objectives [13]. As far as we know, there is no work on the designing of smart energy systems network based controllers considering both economical object and communication network non-ideal behaviours which is the main motivation of the research in the current study.
This article suggests a new distributed network-based economic model predictive control, the so-called DNEMPC, scheme for optimisation of the energy consumption and reduction of the cost of consumed electricity in the network-based smart energy systems. The proposed strategy is based on moving horizon estimation (MHE) strategy which can handle the system constraints and works as an estimator and a predictor simultaneously. To this end, the states which are coupled between sub-systems and their values are unknown, due to time delays, are considered as unknown disturbances. Then each local MHE estimates these coupled states and also the local state of each sub-system by solving an optimisation problem. Moreover, the performance function of the control problem is considered as an economic cost function and the distributed controller is designed so that some economic performances of the process are optimised. The major contributions of the article are listed as follows: (i) Proposing a new distributed network-based economic model predictive scheme for the management of energy consumption and the reduction of the cost of consumed electricity in the network-based smart energy systems. (ii) Developing an estimation procedure to estimate the states coupled between sub-systems in which their values are not available due to delays and also the local state of each subsystem by solving an economic optimisation problem. (iii) Presenting the stability analysis of the proposed strategy. (iv) Reflecting the proposed method for designing a networkbased controller for smart energy systems via two case studies.
The remainder of this article is presented as follows. In Section 2, the transmission network between sub-systems is explained and the formulation of problem is provided. Furthermore, some primary assumptions are represented. The proposed MHE-based cooperative DNEMPC scheme is explained in Section 3. Section 4 is assigned to the convergence analysis of the proposed estimator as well as the stability analysis of the closed-loop system. The simulation studies to verify the obtained theoretical results are given in Section 5 and eventually, some conclusions are drawn in Section 6.
Notations. Throughout the current article, symbol ‖.‖ stands for the induced 2-norm and Euclidean norm for matrices and vectors, respectively. Superscripts "−1" and ''T'' indicate the inverse of square matrices and the transpose operation, respectively. ℝ n and ℝ n×m , respectively, stand for Euclidean space with n dimension and the collection of every n × m real matrices. Furthermore, M refers to the collection of 1, 2, … , M integers. Considering A i is a vector or matrix, the notation col  for a closed ball of radius r > 0 created at origin. In addition, the interior of set is represented by int( ). Furthermore, x + predicates as a successor state for state vector x. Moreover, if ⊆ ℝ n is a convex and compact set in which the origin will be in its interior, then is a ℂ-set. Considering ⊆ ℝ n , if Ax + w ⊆ for all w ∈ , then is called a robust positively invariant (RPI) for x + = Ax + w; w ∈ ⊆ ℝ n and all x ∈ .

NETWORK AND SYSTEM DESCRIPTION
This article considers a networked large-scale smart energy system containing M linear discrete-time sub-systems ( Figure 1). The sub-systems are interacted through their states. Sub-system i ∈ M can exchange the data with its neighbouring sub-systems N i where N i ⊆ M and N i ≠ ∅. The sub-systems are supposed to have the capability for exchanging the information through the communication channels featuring induced time-varying delays. For sub-system i, the time delay of the exchanged data from its neighbouring sub-system j is presented by Remark 1. Owing to [14,15], if the properties of the communication network are known, then the upper bound of delays can be approximated. Thus, all considered bounds of delays are supposed to be known and are constant integers. Moreover, the time-varying communication delays i j k could be replaced with the upper bounds of delays. This does not affect the computation of the controller law.
The current article considers, for delay compensation, a buffer-based strategy described in [7]. According to [7] all DNEMPC controllers accumulate their predictive input sequence with length Nc into unique package with a time stamp and transfer them to their neighbouring sub-systems via the communication channels. It should be noted that, in order to ensure that the time-labelled data is accurate and precise, it is necessary that the clocks of all sub-systems' controllers and sensors are synchronised. For this, some real-time clocks and appropriate synchronisation algorithm can be employed [16].
On the other hand, every local controller includes at most (M − 1) buffers where each buffer is assigned to a neighbouring sub-system and the related packet is saved till the entrance of the next buffer. Once a newer packet data is received, its time label is compared to the time label of the existing stored packet data in the buffer. If the newly received packet data is newer than the stored packet, the buffer will be updated, otherwise the newly received packet data will be ignored and the buffer contents are moved to the left, and afterwards a zero is placed to the right side of the buffer. The buffer output data will be employed for the prediction of interacting states which are not accessible at the current time step, owing to the communication channels time delays.
For instance, assume that at time step (k − 1) the obtained predictive control sequence packet of the jth local controller, that is, {u j (k − 1), u j (k), … , u j (k + N c − 2)}, is written on the relating buffer of the ith local controller. If at time step k, a noncredible packet with time delay i j k is received after the buffer is shifted to the left and zero is added to the right of the sequence. Now{u j (k), … , u j (k + N c − 2), 0}is arranged into the buffer. In the worst case, this stage is repeated d max -1 times until the next credible control packet is received. Thus, if d max ≤ N c , the existence of at least one predicted control data is guaranteed on the output side of the buffer.

Assumption 2.
For each sub-system, the control horizon Nc is equal to or greater thand , whered = d max − d min + 1.
Remark 2. It is worth mentioning that this article can be extended for the case of simultaneously occurring of communication delays and packet dropout. In some studies, packet dropouts are considered as the extended delays [6,7]. In other words, when the time delay in a packet is greater than a prespecified upper bound of delays (here, d max ), it is treated as a loss packet. In the case of packet dropout, an upper bound (T p ) on the number of successive time steps with packet dropouts is considered. In this model, in the worst case, a sub-system will receive a new valid packet from another sub-system within d max + T p time steps. This implies that we have a network with a maximum delay of d max + T p . Generally, the communication channels induced time delays that include the influence of packet dropouts, data received as disordering and network communication time delays.
The dynamic model for sub-system i is described as the linear discrete-time formulation as follows [9]: x and y i ∈ ℝ n i y are the control input, state vector and measured system output, respectively. Furthermore, i and i are ℂ-sets and matrices A ii , B i , A i j , D i and C i are known and constant with appropriate dimensions. Also, i and i are polytopic and polyhedral constraint sets, respectively, which the origin is considered an interior point in them. Moreover, v i k ∈ i ⊆ ℝ n i y and w i k ∈ i ⊆ ℝ n i x stand for the unknown output and state disturbances, respectively. Also, i is assumed as a set of all state variables where there exists a feasible control command existing in i . Besides, x p i j k shows the state trajectory of the j th sub-system which is predicted and computed in ith sub-system, when the induced communication time delay occurs.

Assumption 3.
For i ∈ M , the pairs (A ii , C i ) and (A ii , B i ) are supposed to be detectable, stabilisable, respectively. The whole regular system, ignoring the communicationinduced time delays and disturbances, can be written as i are input and state constraint sets, respectively. The A matrix is introduced as follows [17]:

THE SUGGESTED COOPERATIVE DNEMPC WITH INDUCED COMMUNICATION CHANNEL DELAYS
In this section, firstly an efficient distributed moving horizon estimator-predictor is proposed. Then using the estimated and predicted states, a cooperative DNEMPC is designed such that the states of the closed-loop system converge to a neighbourhood of the origin.
For reducing the estimation error, a pre-estimator is developed for every ith sub-system i as follows: wherex i k andŷ i k are the current estimates of x i k and current output of pre-estimator, respectively. Also, L i shows the gain matrix of pre-estimator which should be computed in such a are the initial states. Moreover,x i k denotes the current estimation of ith sub-system local state trajectory andx p i j k stands for the current predicted coupled state between ith and j th sub-systems.
Next for sub-system i, a local constrained MHE is designed via using the constrained optimisation problem as follows: subject tô where i refers to a non-negative weight for the ith moving horizon estimator and estimation horizon is shown by (Ne + 1) . The termx i l |k demonstrates thex i value at time step l , which is computed at time step k. By solving (4), the optimal solution is obtained and illustrated with ( The optimum sequence of the local states is calculated through (4b), where the current optimum state of ith sub-system is shown byx o k i . In (4a), the term ‖x i k−Ne −x i k−Ne ‖ 2 shows the arrival criterion and this article employs the latter estimate to update thex i k−Ne as the form In sub-system i, the current jth predicted state is computed employing (6): in which u b i j k shows the j th buffer output, j ∈ ℕ i installed in the sub-system at time step k and ℕ j ∖i presents the neighbouring sub-systems of the jth sub-system without considering the ith sub-system. It should be noted that, the u b i j k input trajectory is equivalent to u i k , when no communication channels time delay occurs. Moreover, when delays occur, y j k−1 is considered as 0 in (4e) and (6).
In the following, the NDMPC is developed for system (1). In this regard, for every ith sub-system, at the time step k, the cooperative control criterion function is considered as in which u i . is the control trajectory which should be designed. Moreover, the prices of electricity consumption considered in the optimisation problem as the cost coefficients Co i and Co i j . For all i ∈ I M , Q i > 0, R i > 0 and P i > 0 are weighting matrices. It is noteworthy that if (A ii , Q i ) are detectable, then Q i ≥ 0. Parameter Nc stands for both control and prediction horizons which meetd ≤ Nc. Also, termsx i (⋅) andx p i j k+l are estimated local state trajectory x i (⋅) and the predicted interacted state between the ith and jth sub-systems which are computed through the ith local MHE.
For ith sub-system, at time-step k and every iteration ρ, the optimal control command is obtained through solving (7), which is an optimal regulating control problem with constraints After iterations, at time step k, the obtained optimal solution of mentioned problem is illustrated by ( The final solution is resulted from a convex combination between the current and last optimal solutions of the economic model predictive control problem (7), that is, where a i is the ith subsystems' weighting factor which meets ∑ i∈I M a i = 1. Now, u i k is accumulated into one packet data and transmitted to every interconnected sub-system j ≠ i, through a communication network featuring time delays. Also, the first value of u i k is exerted to the ith sub-system.

THE ANALYSIS OF ESTIMATION ERROR
To derive an analytic expression for the estimation error, (4) is represented as a convex quadratic program (QP) as follows: in which constant matrices G i and i , with appropriate dimensions, present the constraints of (4b)-(4e). r i is a constant term, and z i = col (x i k−Ne , X j ) shows an uncertain vector of optimisation in which X j Δ = col ). The corresponding matrices H i and F i in (8) are where ) .
In (9b), to update matrix F i , we need the past inputs, current and past measured outputs, predicted state and a prior estimate of computed statex i k−Ne . For the ith MHE, the estimation error is introduced as where ). Note that a sequence of the estimation error is obtained using (4b) and (4e) in (10) where the current estimation error is denoted by e i k|k . The boundedness of e i k will be established with the following theorem.
Theorem 1. In sub-system (1) with pre-estimator (3), consider the constrained distributed MHE problem described in (4) and the constrained DNEMPC problem presented in (7). If Assumptions 1-3 hold, and if i ≥ 0 and A L i is Schur, then the current estimation error e i k will be bounded and there exists a ℂ-set i so that for all k ≥ 0 if e i k=0 ∈ i then e i k ∈ i .
The proof is given in the Appendix. For the ith sub-system, the prediction error is defined as Let e p k and x p k are vectors that involve whole prediction errors and predicted states, at time step k, respectively. The nest corollary is implied from Theorem 1, which is needed to guarantee the feasibility of all sub-systems' obtained input trajectory. Also, there are two sets p and p , both including the origin so that p ⊕ p ⊆ and for every x p k =0 ∈ p and e p k =0 ∈ p , x p k ∈ for all k > 0. Furthermore, there are two setŝand , both including the origin so that̂⊕ ⊆ and for everyx k =0 ∈̂and e k =0 ∈ , x k ∈ for all k > 0.
The convexity of i , Corollary 1 and the initialisation procedure employed in the suggested cooperative guarantee that if, for every sub-system, there is a feasible input trajectory, then there are feasible input trajectories for all sub-systems, at future times.

SIMULATION STUDIES
In this section to verify the usability and proficiency of the proposed DNEMPC strategy, it is employed for design of distributed controller for smart energy systems. For this, the following two case studies are considered.

Example 1
In the first simulation study, the proposed DNEMPC method is exerted to control the heat pump for heating the smart residential buildings with the aim of optimising energy consumption and reducing the cost of consumed electricity. In fact, the objective of the economic forecasting control presented in this case study is to control the heat pump in the heating apparatus of the networked distributed buildings with floor heating apparatus by considering the perturbations of radiation of solar and ambient temperature so that by transmitting the pump power consumption time, at times when electricity prices are lower, the demands of the system are met (reaching room temperature to the optimum temperature) as well as the cost of consuming the pump is also minimised.
In this regard, firstly, a dynamic model of two-building floorheating system joined to a heat pump placed in the ground source is presented. The current article adopts the developed method in [18] to construct the model of considered networked building climate control system. Commonly, for a building with N area zone, the model of the ith sub-system (i = 1, 2, … , N )can be represented as follows. where T is the state vector, where T r , T w and T f stand for building air temperature, water temperature in the heating pipes placed in the floor and floor temperature, respectively. Furthermore, u shows the control input which is the power employed via the compressor connected to the heat pump and W i = [T a s ] T stands for the disturbances which are the sun radiation s and ambient temperature T a . Furthermore, the indoor temperature is considered as the output variable, y i = T i r . Moreover, matrices in (11) are as follows: C p, f , where Q ra and Q fr denote the heat which is transferred from the buildings' air to the environment and from the floor to the buildings' air, respectively. Moreover, Q wf shows the heat exchanged between the floor and water which is flowing in heating pipes, under the floor. Finally, the term U A stands for a product of the surface area of the layer between two heatexchanging media and the heat conductivity. The nomenclatures of systems' parameters are provided in Table 1 [18].
Here we consider two building climate control system (two sub-systems) with the same dynamic which interacted through states. The states are considered as T i f and T i w for i = 1, 2. It is supposed that, for every sub-system, the measured output is the indoor temperature, that is, T i r . The inputs to sub-systems are the power employed via the compressor in the heat pump [19] Moreover, the data is exchanged between two sub-systems among the non-ideal communication channels featuring random delays varying between 1 and 4 s.
To show the potential of the proposed DNEMPC, we need to know the electricity prices for a specific time period. This work considers the day-ahead prices of electricity consumption using the Nordic power exchange market presented in [20], as depicted in Figure 2. Furthermore, the disturbances ϕ s and T α are considered as diurnal cycles featuring added noise [20], as illustrated in Figure 3.
First a local pre-estimator is designed for each sub-system as the form of (4) using (19). In the context of DMHE design, Ne = 20. For the DMPC control problem for both subsystems, the control and prediction horizons are equal, and Nc = 30.
The main objective is to minimise the total cost of electricity consumption in a predefined time period as well as holding the indoor temperature, T r , in a specified limit in the presence of communication network time delay. The data is transmit-ted among the non-ideal communication channels with random delays varying between 1 and 12 s. Figure 4 shows time-varying communication delay sequences induced by the communication network in a day. Moreover, to have a better perceptiveness, a conventional economic model predictive (CEMPC) approach, presented in [21], is adopted and applied for comparison.
The predicted indoor temperature and obtained optimal schedule of compressor, for a 5-day horizon and over the predefined time-varying constraints are provided in Figure 5. The constraints show that the temperature can be lower during the night-time compared to daytime. Furthermore, Figure 6 illustrates the actual electricity prices and the obtained optimal heat pump power input. Figures 5  and 6 verify that the consumption of power is shifted to periods when the electricity is cheap and show that the house floor thermal capacity has the ability to save enough energy so that, during the daytime, the heat pump can be turned off. This issue implies that the slow heat dynamics of the floor can be employed to move the consumption of energy to the periods where the price of electricity is low as well as the indoor temperature is kept in the desired value. Furthermore, Figure 6 reveals that the proposed DNEMPC method has better performance than conventional CEMPC method in shifting the consumption of power to periods when the electricity is cheaper.

Example 2
To verify the usability and proficiency of the proposed DNEMPC strategy, it is applied for the design of load frequency controller for a four-area interconnected power system over communication network. The main aim of the load frequency controller in the modern multi-area power system is to keep the uniform frequency at each area and to maintain the tie-line power exchanges in a pre-defined tolerance in the presence of load disturbances or sudden changes in load demands. In this regard, a four-area network-based power system model, as depicted in Figure 7, is given for the design of network-based  Four-area interconnected power system LFC system. It is assumed that the ith control area (i = 1, 2, 3, 4) is coupled with the jth control area i≠j, via a tie line. Furthermore, each control area can be expressed with the following equation [22]: where x i ∈ ℜ n i , u i ∈ ℜ m i and y i ∈ ℜ l i stand for the state, control and output vectors of the ith control area, respectively. Also,x p i j ∈ ℜ n j is the state vector of the neighbouring control area. The output and state vectors for the ith area are presented as follows [22]: Moreover,A i , A i j and B i stand for real and constant matrices which are adopted from [22] and the values of entire the param- eters are presented in Table 2: For all areas, the nominal value of operating frequency, static load and the set point of load reference are considered same as presented in [22], where f = 60 Hz, ΔP L i = 0 and −0.3 ≤ ΔP re f i ≤ 0.3 (i = 1, 2, 3, 4). Also, the time-varying induced delays of communication channels are considered as 12 (k), 13 (k), 14 (k), 23 (k), 24 (k) and 34 (k).  To prove the superiority of the suggested method, a distributed networked model predictive control (DNMPC) and CEMPC control approach are adopted, respectively, from the [9,21] and applied for comparison. Furthermore, the efficiency of the developed network-based control technique is evaluated through the following scenarios: Plan 1. The data is transferred via the communication channels with random delays changing between 1 and 5 s (i.e. d max is chosen as 5 s).  Figures 8 and 9, for plans 1 and 2, respectively. Figure 8 reveals that the proposed DNEMPC technique and the applied DNMPC strategy [9] have much better performance in the presence of communication-induced delays than the CEMPC method [21]. In addition, Figure 9 proves the superiority of the suggested DNEMPC technique scheme in comparison to other methods when a large delay is considered. Moreover, the Euclidean norm of the estimation error for the developed DNEMPC and plans 1 and 2 are provided in Figure  10. From this figure, it can be concluded that the states of the system are regulated satisfactorily.

CONCLUSION
The current article suggests an efficient cooperative-based distributed economic model predictive control for optimising the consumption of energy and reducing the cost of consumed electricity in networked smart energy systems. The proposed approach is based on MHE method which can handle system constraints, and works as an estimator and a predictor simultaneously. Then each local MHE estimates these coupled states, and also the local state of each sub-system, by solving an optimisation problem. To verify the efficiency of the designed control scheme, the distributed economic model predictive control method is applied to control the heat pump for heating the smart residential buildings with the aim of optimising energy consumption and reducing the cost of the consumed electricity. The obtained results show the profitability of the proposed method. The proposed strategy can be developed and employed easily to other wide-area systems.
.λ iA , G iA and ξ iA stand for the active Lagrange multipliers and the relating active inequality matrices of (7), respectively.
The fact A i L is Schur implies that A i e and M i A i e (M i ) −1 are Schur matrices and there is a ℂ-set i which is robust positively invariant for (22) [17]. It follows that M i A i e (M i ) −1 i ⊕̄i e ⊂ i and if e i k=0 ∈ i , then e i k ∈ i , ∀k ≥ 0.