A compositional modeling framework for the optimal energy management of a district network
Daniele Ioli Alessandro Falsone Alessandro Vittorio Papadopoulos Maria Prandini
cs.SY






DEIB]Daniele Ioli
DEIB]Alessandro Falsone
MDH]Alessandro Vittorio PapadopoulosCA
DEIB]Maria Prandini

[CA]Corresponding author, Tel. +46 (0)21-1073 23
[DEIB]Politecnico di Milano, Piazza Leonardo da Vinci, 32, 20133 Milano, Italy
[MDH]Mälardalen University, Högskoleplan 1, 72123, Västerås, Sweden
[myfootnote]This work is partly supported by the European Commission under the UnCoVerCPS project, grant number 643921, and was performed when the third author was a post-doctoral researcher at Politecnico di Milano. 



This paper proposes a compositional modeling framework for the optimal energy management of a district network. The focus is on cooling of buildings, which can possibly share resources to the purpose of reducing maintenance costs and using devices at their maximal efficiency. Components of the network are described in terms of energy fluxes and combined via energy balance equations. Disturbances are accounted for as well through their contribution in terms of energy. Different district configurations can be built, and the dimension and complexity of the resulting model will depend on the number and type of components and on the adopted disturbance description. Control inputs are available to efficiently operate and coordinate the district components, thus enabling energy management strategies to minimize the electrical energy costs or track some consumption profile agreed with the main grid operator. 
 


	Smart grid modeling Compositional systems Energy management Building thermal regulation









§ INTRODUCTION


Building energy management, and temperature regulation in particular, has recently attracted the attention of various researchers (see, e.g., <cit.>).
Indeed, energy consumption in buildings represents approximately 40% of the worldwide energy demand, and more than half of this amount is spent for Heating, Ventilation and Air Conditioning (HVAC) systems <cit.>.
Energy management can be performed at the level of a single building, e.g., using the storage to shift in time the thermal energy request to time slots where the electricity costs are lower. As buildings started sharing equipments at the benefit of shared operating costs, increased flexibility, and overall performance improvement, energy management needs to be performed at the district network level, which calls for appropriate modeling and control strategies. Constructing models of interconnected systems is generally demanding, and here we propose a modular framework that simplifies this task and is also suitable for the application of different control design approaches.

The proposed modeling approach is oriented to energy management and compositional in that components are described in terms of thermal/electrical energy fluxes  and interact by exchanging energy, which makes easy to compose a district network configuration via energy balance equations.
Our modeling framework is built with a control-oriented perspective. It includes disturbances like, e.g., solar radiation, outside temperature, occupancy, and wind power production, as well as control inputs like, e.g., buildings temperature set-points, charge/discharge of storages, activation/deactivation of devices, that can be appropriately set so as to optimize performance at the district level.

Complexity and size of the model associated with a district configuration depend on number of components and type of description adopted per component. The model can be either deterministic or stochastic depending on the disturbance characterization as a deterministic or stochastic process, respectively. It can range from a low dimensional deterministic system  with continuous input and state that is convex in the control input, to a large dimensional  Stochastic Hybrid System (SHS) <cit.> with discrete and continuous input and state.

Given a certain configuration, one can then formulate energy management problems like the minimization of the cost of the electrical energy requested to the main grid or the tracking of some given electrical energy exchange profile that was agreed with the main grid operator according to a demand-response strategy. The district network in the latter case can be viewed as a user that actively participates to the electrical energy demand/generation balance of the overall grid, and, hence, to its stabilization.

Further contributions in the literature adopt a similar perspective. In <cit.>, the focus is on simulation so that the model dependence on the control input is not a concern. In <cit.>, the aim is the design of an energy management strategy via a simulation-based approach. The modeling effort is limited in this case, and the idea is to take an accurate model in the literature and run simulations to the purpose of policy design, with no concern of making explicit the dependence on the input and formally proving optimality.
The approach in <cit.> is the closest to our approach, in that it addresses energy management problems for a  microgrid that is built based on models of single components, combined via energy balance equations. Models are however simplified, in particular that of the building. Also, occupancy is not accounted for explicitly. A specific strategy for energy management is considered, whereas our framework is more comprehensive since it allows for the design of different strategies (certainty equivalence based, robust, stochastic) for the minimization of suitably defined (nominal, min-max, average) cost in presence of (nominal,  robust, probabilistic) constraints on comfort and actuation. Depending on the  network communication and computation capabilities, and on privacy issues like in the case of buildings not willing to disclose their consumption profile, a centralized, decentralized, or distributed optimization scheme can be conceived and implemented. Overall, our work is more general and it actually subsumes the approach in <cit.>.

It is worth noticing that other modeling frameworks have been developed in the literature <cit.>. However, the obtained models are typically more complex since they are based on partial differential equations, and require numerical optimization tools for solving the resulting nonlinear optimization problems <cit.>. 


This paper is based on our earlier work in <cit.>, which is extended in several directions. We provide a more detailed description of the district components, including a validation with respect to other commercial simulation tools of the building thermal model according to a norm defined by the American Society for Heating Refrigerating and Air-conditioning Engineers (ASHRAE). We show how to compose a network configuration and formulate an energy management problem as an optimization program. In particular, we show a simulation study of the results achieved in the case where nominal disturbances are present and computations are performed by a central unit. Finally, we suggest a multirate approach as a viable solution for allowing real-time computation of the control input, while retaining model accuracy.


The reminder of the paper is structured as follows. Section <ref> presents the models of the district network components, and Section <ref> shows how they can be connected to set up a network configuration, while defining objective and constraints of the optimal energy management problem. Section <ref> describes some configurations, providing examples of the kind of results one can achieve through the presented framework. Section <ref> shows how to deal with computational complexity, discussing a multirate approach, while Section <ref> concludes the paper. <ref> describes the procedure adopted for validating the model of the building.

 

§ DISTRICT NETWORK COMPONENTS



We consider a district network connected to the main grid that will provide the electrical energy needed to compensate for possible imbalance between demand and generation within the district. We model the evolution of the network over a finite time horizon [t_i,t_f], which is divided into M time slots of duration . The contribution in terms of energy requested/provided by the different components per time slot along the discretized control horizon is provided.
Components can consume (e.g., buildings), produce (e.g., renewable power generators), store (e.g. thermal storages and batteries), or convert energy (e.g., the chiller plants), and are combined via energy balance equations so as to build the overall model of the district. Each component may be affected by some inputs which can be either disturbances or control inputs. In the case when control inputs are available, a suitable strategy can be conceived to set them so as to efficiently manage the system along the time horizon [t_i,t_f].

In the rest of this section, we provide a model for the following components: building,  chiller, storage, combined heat and power unit, and wind turbine. Models are either derived from first principles or taken from the literature.  In the latter case, appropriate references are provided. Tables <ref>–<ref> summarize the main characteristics of the first 4 components. The last component provides an input to the network in terms of wind energy. Similarly to the wind energy contribution, one could consider the solar energy contribution provided by photovoltaic panel installations. Models partly derived from first principle and partly taken from the literature could be used to this purpose. This is not treated here, but the interested reader can refer to, e.g.,  <cit.>. Further components could also be added to the district network. The key idea when introducing our compositional framework is that if a component can be modeled in terms of energy, possibly depending on some control input and/or disturbance signal, then, it can be easily included in the network. When the dependence of the energy on the control input is convex, piecewise linear, or linear with additional binary variables, the problem of designing an energy management strategy can be reduced to a mixed integer linear or a convex optimization program for which efficient solvers exist. 





 §.§ Building

We consider a building as composed of n_z zones, where each zone is characterized by its own (average) temperature T_z,j, j=1,…,n_z.
The zones temperatures can be collected in a vector T_z = [T_z,1⋯ T_z,n_z] and we next determine the amount of cooling energy E_c needed for making them track a given profile. We say that the building is controllable if a control layer is present to this purpose. Suitable constraints will be imposed on the assigned profile to make the resulting tracking problem feasible while guaranteeing comfort conditions at the same time.

The cooling energy E_c,j requested by zone j can be derived based on the thermal energy balance within the zone, accounting for both thermal effects related to its structure and thermal phenomena related to occupancy, equipment, lights, etc, and solar radiation through windows.
More precisely, we have

    E_c,j = E_w,j + E_z,j+ E_p,j + E_int,j,

where E_w,j  is the amount of energy exchanged between  zone j and its adjacent walls, E_z,j is the contribution of the thermal inertia of zone j, and E_p,j and E_int,j is the heat produced by people and other heat sources within zone j, respectively.

The thermal model of the building is derived from first principles, following <cit.>.




  §.§.§ Walls contribution

For modeling the walls contribution we use a one-dimensional finite volumes model. Each wall is divided into vertical layers (`slices') that may differ in width and material composition. The area of each slice coincides with the wall area and each slice is assumed to have uniform density and temperature. The one-dimensional discretization is sensible since the heat flow is perpendicular to the crossed surface. Each internal slice exchanges heat only with nearby slices through conduction, whilst boundary slices are exposed towards either a zone or the outside of the building and exchange heat also via convection and thermal radiation. External surfaces are assumed to be gray and opaque, with equal absorbance and emissivity and with zero transmittance. Absorbance and emissivity are wavelength-dependent quantities, and here we shall consider two different values for shortwave and longwave radiation.

The heat transfer balance equation for the iþ slice of the wþ wall is given by:

    Ṫ_w,i =   1/C_w,i[ (k_w,i^i-1 + h_w,i^i-1)T_w,i-1 + (k_w,i^i+1 + h_w,i^i+1)T_w,i+1
       -(k_w,i^i-1 + h_w,i^i-1 + k_w,i^i+1 + h_w,i^i+1)T_w,i + Q_g,w,i + R_w,i],

where T_w,i denotes the temperature of the wall slice, C_w,i being its thermal capacity per unit area, and k_w,i^j and h_w,i^j, with j = i± 1, representing respectively the conductive and convective heat transfer coefficients between the i^th and the j^th slice of the same wall w. Q_g,w,i is the thermal power generation inside slice i and R_w,i represents radiative heat exchanges and is defined as

    R_w,i =
    		0,		    1 < i < m 
    α_w^S Q^S + α_w^L Q^L - ε_w,i Q_r(T_w,i),	   slice i facing outside
    ∑_w'=1,…,n_w
     j∈{1,M} F_(w,i)→(w',j)( ε_w',j Q_r(T_w',j) - ε_w,i Q_r(T_w,i) )	    slice i facing inside

where Q^S and Q^L denote the incoming shortwave and longwave radiation power per unit area, respectively, and α_w^S and α_w^L are the corresponding absorbance rates for wall w. Q_r(T_w,i) is the emitted radiation as a function of the slice temperature, ε_w,i < 1 being the emissivity and F_(w,i)→(w',j) the view factor that takes into account the fraction of radiation leaving slice i of wall w and reaching slice j of wall w'. Finally, n_w denotes the total number of walls.

Equation (<ref>) holds for every slice in every wall w.
If the wall is composed of m slices, we have m equations like (<ref>) with i=1,2,…,m.
When the superscript in the right-hand side of equation (<ref>) takes value 0 or m+1, reference is made to either a zone of the building (internal surface of the wall) or the outside of the building (external surface of the wall). Note that k_w,1^0 = k_w,m^m+1 = 0 as there is no thermal conduction on walls boundary surfaces, h_w,i^i-1 = 0 for i>1, h_w,i^i+1 = 0 for i<m, and ε_w,i = 0 for 1<i<m, since there is no thermal convection nor radiation between slices. As for the slice in contact with the ground, we assume that the energy exchange occurs via thermal conduction only (no convection nor radiation is considered), where the ground is considered as a thermal reservoir, and as such its temperature is constant.
Since we assume that each wall is a gray body, the power Q_r(T_w,i) radiated from each slice is governed by Q_r(T_w,i) = σ T_w,i^4, where σ is the Stefan-Boltzmann constant. This expression is approximately linear around the slice mean operating temperature T_w,i so that it can be replaced by

    Q_r(T_w,i) = 4 σT_w,i^3 T_w,i - 3 σT_w,i^4.

Then, the evolution of the temperatures T_w = [T_w,1⋯ T_w,m] of the m slices composing wall w can be described in matrix form by

    Ṫ_w = A_wT_w + B_wT_z + W_wd,

where we recall that T_z is the vector containing the temperatures of the n_z zones.
Vector d = [T_out  T_gnd  Q^S   Q^L   1] is the disturbance input and collects the outdoor temperature T_out, the ground temperature T_gnd, and the incoming shortwave Q^S and longwave Q^L radiations. The constant 1 in d is introduced to account for the constant term in (<ref>). Finally, A_w, B_w and W_w are suitably defined matrices that are easily derived based on the scalar equation (<ref>), whose coefficients depend on the wall characteristics.


Equation (<ref>) refers to a single wall. If there are n_w walls in the building, then, we can collect all walls temperatures in vector T = [T_1⋯T_n_w], and write the following equation for the evolution in time of T:

    Ṫ = AT + BT_z + Wd,

where A is a block-diagonal matrix with A_w as wþ block, B = [   B_1     ⋯ B_n_w ] and
W = [ W_1⋯W_n_w ].

If we consider zone j and one of its adjacent wall w, then the thermal power transferred from wall w to zone j is given by

    Q_w → j = S_w h_w,b^b' (T_w,b - T_z,j),

where S_w is the wall surface and the pair (b,b') can be either  (1,0) or (m,m+1) according to the notation introduced for (<ref>).  The total amount of thermal power transferred from the building walls to zone j can be expressed as Q_b,j = ∑_w∈𝒲_j Q_w→ j, where 𝒲_j is the set of walls w adjacent to zone j. Defining Q = [Q_b,1⋯ Q_b,n_z], we obtain

    Q = CT + DT_z,

where C and D are suitably defined matrices derived based on equation (<ref>). From (<ref>) and (<ref>), we finally get

    Ṫ = AT + BT_z + Wd
    Q = CT + DT_z



The obtained model, though linear, can be quite large. However,  its order can be greatly reduced by applying the model reduction algorithm
based on Hankel Single Value Decomposition (HSVD), as suggested in <cit.>.
 



The zone temperature profile to track T_z is taken as a linear function of time within each time slot of length , defined by the values u(k) = T_z(k) at the time steps k=0,1,…,M.
By approximating the input d as a piecewise linear function of time as well, with values ω(k) = d(k) at k=0,1,…,M, an exact discrete time version of the linear model (<ref>) can be derived (see <ref>).
The evolution of y(k) = Q(k) over the finite time horizon can then be computed as

    y = [y(0) ⋯ y(M)]= FT(0) + Gu + Hω

where we set u = [u(0) ⋯ u(M)] and ω = [ω(0) ⋯ω(M)],
and F, G and H are suitably defined matrices.


The thermal energy E_w(k) = [E_w,1(k) ⋯ E_w,n_z(k)] transferred from the walls to all zones can be computed by integrating Q(t) on each time
slot, which leads to the following approximate expression:

    E_w(k) = /2 ( y(k-1) + y(k) ), k=1,…,M.

Finally, from (<ref>) and (<ref>) we can derive the enlarged energy vector E_w = [E_w(1) ⋯ E_w(M)]:

    E_w = F̃T(0) + G̃u + H̃ω,

where F̃, G̃, and H̃ are obtained from matrices F, G, and H in (<ref>) via (<ref>).







  §.§.§ Zones energy contribution


In order to decrease the temperature of zone j in the time frame from (k-1) to k, we need to draw energy from the zone itself. This energy contribution can be expressed as

    E_z,j(k) = - C_z,j (T_z,j(k) - T_z,j((k-1))),

where C_z,j is the heat capacity of the jþ zone.
If we account for all n_z zones, and all M time frames within the finite horizon [t_i,t_f], equation (<ref>) can be written in the following matrix form

    E_z = Zu,

where we set E_z = [E_z(1) ⋯ E_z(M)] with E_z(k) = [E_z,1(k) ⋯ E_z,n_z(k)],
and Z is a suitably defined matrix.





  §.§.§ People energy contribution

Occupancy implies heat production, which in crowded places can be actually significant <cit.>.
According to an empirical model in <cit.>, the heat rate Q_p,j produced by the n_p,j occupants of a zone j at temperature T_z,j is given by

    Q_p,j = n_p,j ( p_2 T_z,j^2 + p_1 T_z,j + p_0 ),

where p_2 = -0.22 W/K^2, p_1 = 125.12 W/K and p_0 = -1.7685 · 10^4 W. Expression (<ref>) is almost linear in a sensible operating temperature range and can thus be accurately approximated by linearization around some comfort temperature T_z,j:

    Q_p,j   = n_p,j( (2 p_2 T_z,j + p_1)(T_z,j-T_z,j) + p_2 T_z,j^2 + p_1 T_z,j + p_0 ) 
       = n_p,j( p̃_1 T_z,j + p̃_0 ).

Recall now that the zone temperature profile T_z,j to track is assumed to be linear in time. If we approximate the occupancy n_p,j as a linear function of time within each time slot as well, as suggested in <cit.>, then equation (<ref>) can be analytically integrated from (k-1) to k to obtain the energy transferred to zone j in the kþ time slot:
	
    E_p,j(k) =     q_2,k(n_p,j)T_z,j(k)  + q_1,k(n_p,j)T_z,j((k-1)) + q_0,k(n_p,j)

where we set

    q_2,k(n_p,j) = p̃_1 6( 2 n_p,j(k) + n_p,j((k-1))  ) 
       q_1,k(n_p,j) = p̃_1 6( n_p,j(k) + 2 n_p,j((k-1))  ) 
       q_0,k(n_p,j) = p̃_0 2( n_p,j(k) + n_p,j((k-1))  )


The total amount of energy transferred to all zones in each time slot can be packed in a vector E_p(k) = [E_p,1(k) ⋯ E_p,n_z(k)] and then, defining E_p = [E_p(1) ⋯ E_p(M)] and n_p = [n_p,1(0)  n_p,1() ⋯ n_p,1(M) ⋯ n_p,n_z(0)  n_p,n_z() ⋯ n_p,n_z(M)], one can write
that

    E_p = N(n_p)u + e(n_p),

where N(n_p) and e(n_p) depend on the coefficients (<ref>).

Note that occupancy profiles can be either obtained from data or derived from a stochastic model, like, e.g., the one in <cit.> which is based on Poisson arrival/departure processes <cit.>.

Further energy contributions of the building occupants, in terms for instance of blinds movement and setpoint override, are not modeled here. Recent works on human-building interaction discuss the impact of human intervention on energy management strategies. The interested reader is referred to <cit.>, where a possible strategy to limit human intervention is proposed, and to <cit.>, where a model predictive control solution is suggested for timely adjusting the control action to unpredicted human disturbances.





  §.§.§ Other internal energy contributions

There are many other types of heat sources that may affect the internal energy of a building, e.g., lighting, daylight radiation through windows, electrical equipment, etc. The overall heat flow rate produced within zone j can be expressed as the sum of three contributions, namely

    Q_int,j = α_j Q^S  + λ_j+ κ_j I_ℝ^+(n_p,j),

where α_j is a coefficient that takes into account the mean absorbance coefficient of zone j, the transmittance coefficients of the windows and their areas, sun view and shading factors, and radiation incidence angle. I_ℝ^+(·) denote the indicator function on the positive real values. The thermal energy contribution to zone j due to internal lighting and electrical equipment is composed of two contribution: a constant term λ_j, and an additional therm κ_j that represents the change in internal lighting and electrical equipment when people are present. Note that Q_int,j does not depend on Q^L because windows are usually shielded against longwave radiation.
The energy E_int,j(k) during the k^th slot is given by:

    E_int,j(k) =    2[Q^S(k) + Q^S((k-1))] + λ_j
        +2κ_j [ I_ℝ^+(n_p,j(k)) + I_ℝ^+(n_p,j((k-1))) ]

and is obtained by (<ref>), where the first (linear) and second (constant) terms have been analytically integrated, whereas the third term has been treated separately, due to the presence of the indicator function. In the cases when occupancy drops to zero or becomes nonzero in a time slot, the  energy contribution is set to a half of the contribution in the case when occupancy is nonzero at the beginning and at the end of the time slot.
We can collect the thermal energy of the zones in a single vector E_int(k) = [E_int,1(k) ⋯ E_int,n_z(k)], and  then define E_int = [E_int(1) ⋯ E_int(M)], which is finally given by:

    E_int = M ω + L(n_p).





  §.§.§ Overall building cooling energy request


Now we can finally compute the cooling energy demand of all zones in the building for tracking the piecewise linear zone temperature profiles T_z specified via the input u at the discrete time instants k=0,1,…,M during the time horizon [t_i,t_f].
Specifically, from (<ref>) it follows that E_c=[E_c(1) ⋯ E_c(M)] with E_c(k)=[E_c,1(k) ⋯ E_c,n_z(k)] is the sum of four contributions:

    E_c = E_w + E_z + E_p + E_int,

where E_w is given in (<ref>), E_z in (<ref>), E_p in (<ref>), and E_int in (<ref>). This leads to the following expression for the cooling energy demand:

    E_c   = F̃T(0) + (G̃ + Z + N(n_p)) u + (H̃+M)ω  + e(n_p)+ L(n_p)
       = A_cT(0) + B_c(n_p) u + W_cω + b(n_p)

where A_c, and W_c are constant matrices, whereas B_c(n_p) and b(n_p) depend on the occupancy.
Note that the input u defining the zone temperature profiles enters affinely the system dynamics if the occupancy n_p were fixed.






  §.§.§ Building block: interfaces and related constraints


The thermal model of the building can be considered as a block with the following input/output interfaces:
the control input vector u specifying the piecewise linear  zone temperature profiles T_z at the discrete time instants k=0,1,…,M, and
disturbance input vectors n_p and ω  representing the occupancy and the collection of outdoor temperature T_out and incoming shortwave Q^S and longwave Q^L radiations, respectively; and the output vector E_c of the cooling energy demand requested by the zones in the building to track T_z.


Notice that the cooling energy demand cannot be negative. Furthermore, a profile where the zone temperature is required to decrease with a steep slope cannot be tracked. This can be formulated as a constraint  on the maximum amount of energy E_c,j^max that can be requested by a zone j per time slot (from which the upper bounding vector E_c^max of the same size of E_c can be derived), and, possibly, a maximum amount E_c,b^max that can be requested by the building during the whole time horizon. This maps into the following actuation constraints:

    0≤E_c≤E_c^max,   1E_c≤ E_c,b^max,

where 1 denotes a column vector with all elements equal to 1 so that 1E_c is the total cooling energy requested by the building. Note that when a vector is compared with a scalar like in (<ref>), it means that each component of the vector is compared with that same scalar.

Table <ref> summarizes the relevant quantities related to the building model. The type attribute is introduced to denote possible different models that can be used, which eventually has some impact on the energy management problem formulation. Type A is the controllable building model where the zone temperature profiles can be optimized via the control input u, whereas Type B is the uncontrollable building model where the zone temperature profiles cannot be chosen but are already specified via some given u vector. In a network configuration, it is possible to include both controllable and uncontrollable buildings. Comfort and cooling energy bounds can then be enforced only in the case of Type A model, which contributes  to the network description with equations and inequalities that are linear in the control input.







 


 §.§ Chiller plant

A chiller plant is an electrical devices that reduces the temperature of a liquid, typically water, via vapor compression or absorption cycle. In this way, it
converts the electric power provided by the electrical grid into cooling power, which is then conveyed to either some cooling load or some thermal storage via the chilled water circuit.

Chillers can be modeled through the equation

    E_,ℓ=a_1T_oT_cw+a_2(T_o-T_cw)+a_4T_oE_,c/T_cw-a_3/E_,c-E_,c,  0≤ E_,c≤ E_,c^max,

where E_,ℓ is the electrical energy absorbed by the chiller in order to provide the cooling energy E_,c  in a time slot of duration , and E_,c^max is the corresponding maximum cooling energy production.
Note that E_,ℓ depends also on the outdoor temperature T_o and the temperature of the cooling water T_cw. The latter is typically regulated by low level controllers so that it is maintained almost constant at some prescribed optimal operational value, which also facilitates the stratification in the thermal storage.
The chiller description (<ref>) is derived from the original Ng-Gordon model <cit.> which is based on entropy and energy balance equations and accounts also for heat losses and  pump contribution to the electric energy consumption (E_,ℓ> 0 when E_,c=0)..
Coefficients a_1, a_2, a_3, a_4 characterize the chiller performance. Depending on their values, we can have different efficiency curves as given by the Coefficient Of Performance (COP), which is the ratio between the produced cooling energy and the corresponding electrical energy consumption:

    COP = E_ch,cE_ch,ℓ.

Figure <ref> shows an example of curves of the COP for three chiller units of different size, with their respective approximations presented in the following sections.





We next introduce simpler approximations of relation (<ref>), which preserve convexity in the control input E_,c. 



  §.§.§ Chiller model approximations

A convex biquadratic approximation

    E_,ℓ=c_1(T_o)E_,c^4+c_2(T_o)E_,c^2+c_3(T_o),  0≤ E_,c≤ E_,c^max,

of the nonlinear Ng-Gordon model (<ref>) can be derived by using weighted least square to best fit the most relevant points, i.e, those that correspond to zero energy request and to the maximum COP values.

Another possible convex approximation of (<ref>) is via a PieceWise Affine (PWA) function given by the following convex envelope of a finite number of affine terms

    E_,ℓ=max{m_c(T_o)E_,c+q_c(T_o)},  0≤ E_,c≤ E_,c^max,

where  the coefficients of the affine terms are collected in the two vectors m_c(T_o) and q_c(T_o), and the max operator is applied component-wise.
Note that, if E_,ℓ in expression (<ref>) is to be minimized, then (<ref>) can be easily translated as a set of linear constraints with an epigraphic reformulation.

The quality of the biquadratic and PWA approximations is compared in Figure <ref>.







  §.§.§ On-off switching

As shown in Figure <ref>, the chiller absorbs some amount of electrical energy even when no cooling energy is produced.
In order to have the possibility of switching the chiller on and off, one can introduce the binary variable δ_(k), k=0,…,M, that represents the on (δ_(k)=1) and off (δ_(k)=0) logical status of the chiller at time k, k=0,…,M. The cooling energy request E_,c(k) and on-off command δ_(k) are related via the logical condition

    δ_(k) = 1   ⇔   E_,c(k) > 0.

Let E_,c^max be the maximum value for E_,c and ε a small quantity, typically set equal to the machine precision.
Using the Conjunctive Normal Form in <cit.>, (<ref>) can be expressed as a mixed integer linear condition:

    εδ_(k) ≤ E_,c(k) ≤ E_,c^maxδ_(k),

which leads to δ_(k) = 0 ⇔ E_,c(k) =0 and δ_(k) = 1 ⇔ E_,c(k) ∈ [ε,E_,c^max], that are practically equivalent to (<ref>).
Depending on the adopted approximation, we can rewrite the model of the chiller including the on-off condition as

    E_,ℓ(k)=
     (c_1(T_o(k))E_,c(k)^4+c_2(T_o(k))E_,c(k)^2+c_3(T_o(k)))δ_(k)
    max{m_c(T_o(k))E_,c(k)+q_c(T_o(k))}δ_(k)

with 0≤ E_,c(k)≤ E_,c^max.
The PWA formulation is particularly convenient since the product between a (piecewise) affine function Mx+q and a discrete variable δ can be reduced to a mixed integer linear condition <cit.>, by introducing the auxiliary variable z=δ (M x+Q) subject to 0≤ z ≤min{M x+Q+(1-δ) ,  δ}, where  is an upper bound on Mx+q.




  §.§.§ Chiller block: interfaces and related constraints


The chiller block can be described with a static map between the cooling energy  E_,c = [E_,c(0) ⋯ E_,c(M)] that it produces and the corresponding absorbed electrical energy E_,ℓ = [E_,ℓ(0) ⋯ E_,ℓ(M)].

The cooling energy that the chiller can provide is subject to the following bound:

    0≤E_,c≤ E_,c^max,

which maps into a bound on the absorbed electrical energy

    0≤E_,ℓ≤ E_,ℓ^max.


When the on-off command δ_ = [δ_(0) ⋯δ_(M)] is introduced as an additional control input, the following further constraint enters the chiller model:

    εδ_≤E_,c≤ E_,c^maxδ_.



Table <ref> summarizes the relevant quantities of the chiller model, with Type A, B, C, and D representing possible modeling variants. The max operator is applied element-wise, and the symbol * is the element-wise multiplication.




 


 §.§ Storage

Thermal Energy Storages (TESs) are becoming widely used in medium size grids. TESs represent the most effective way, or even sometimes the only way, to take advantage of renewable energy sources. This is indeed the case for thermal solar energy and geothermal energy systems. In a smart grid context, they can be used as energy buffers for unbinding energy production from energy consumption. More specifically, in a district cooling scenario, a TES for cooling energy can shift the request of cooling energy production to off-peak hours of electrical energy consumption, make chillers operate in high-efficiency conditions, and smooth peaks of electrical energy request with benefits both for power production and distribution network systems, see e.g. <cit.>.

There are many different technical solutions to store thermal energy, the most widely used are fluid tanks and Phase Changing Materials (PCMs) storages.
We focus next on fluid tanks modeling, and add a note on how the model can be extended to PCMs storages in Remark <ref>.
From an energy management perspective we will make use of a black box model, derived based on system identification techniques, that uses the energy exchange (added or removed) as input and the thermal energy stored as output. The simplest model is a first order AutoRegressive eXogenous  (ARX) system

    S(k+1) = aS(k) - s(k),

where S(k) is the amount of cooling energy stored and s(k) is the cooling energy exchanged (s(k)>0 if the storage is discharged, and  s(k)<0 if it is charged) in the kþ time slot,  while a ∈ (0,1) is a coefficient introduced to model energy losses.


By unrolling the thermal storage dynamics in (<ref>) we can express the cooling energy stored along the look-ahead discretized time horizon [t_i,t_f] in a compact form as

    S = Ξ_0 S(0) + Ξ_1 s,

where we set S = [S(1) ⋯ S(M)], s = [s(0) ⋯ s(M-1)], and Ξ_0 and Ξ_1 are suitably defined matrices.

A more sophisticated model can be obtained by introducing dissipation effects through the efficiency coefficients β_C ∈ [0,1] and β_D ∈ [0,1] for the charge/discharge dynamics as follows:

    S(k+1)= aS(k)-((1-β_C)δ_C+(1+β_B)δ_D)s(k),

where δ_C(k) ∈{0,1} and δ_D(k) ∈{0,1} indicate the mode in which the storage is operated: δ_C(k) = 1 and δ_D(k) = 0, the storage is charged (s(k)<0), δ_C(k) = 0 and δ_D(k) = 1 the storage is discharged (s(k)>0), and δ_C(k) = δ_D(k) = 0 the storage is not used.
Notice that δ_C(k) and δ_D(k) are mutually exclusive, which can be coded via the constraint

    δ_D(k) + δ_C(k) ≤ 1.

It is possible to set minimum and maximum thresholds for the energy exchange rate in both the charging and discharging phases by constraining s(k)  as follows:

    δ_D(k) s_D^min + δ_C(k) s_C^max≤ s(k) ≤δ_D(k) s_D^max + δ_C(k) s^min_C

with s_C^max < s_C^min≤ 0 and 0 ≤ s_D^min < s_D^max.
Note that if δ_C(k) = δ_D(k) = 0 (storage not in use), inequalities (<ref>) degenerate to the condition s(k)=0.

Model (<ref>) is bilinear in the control inputs since  δ_C(k) and δ_D(k) are multiplied by s(k). However, we can reduce it to the linear model

    S(k+1)= aS(k)-(1-β_C)s_C(k)-(1+β_D)s_D(k)

by replacing s(k) with the new control variables s_C(k) = δ_C(k)s(k) and s_D(k) = δ_D(k)s(k).
Accordingly, the constraint (<ref>) becomes

    δ_C(k)s_C^max≤ s_C(k) ≤δ_C(k) s_C^min
       δ_D(k)s_D^min≤ s_D(k) ≤δ_D(k) s_D^max.

The energy exchange s(k) can then be recovered from s_C(k) and s_D(k) as s(k) = s_C(k)+s_D(k).

Model (<ref>) subject to constraints (<ref>) and (<ref>) is equivalent to model (<ref>)
subject to constraints (<ref>), (<ref>) and (<ref>). This latter model has the advantage of being linear so that it can be expressed in compact form along the look-ahead discretized time horizon [t_i,t_f] as follows:

    S = Ξ_0 S(0)+ Ξ_C s_C + Ξ_D s_D
       s = s_D+ s_C,

where  s_C = [s_C(0) ⋯ s_C(M-1)], s_D = [s_D(0) ⋯ s_D(M-1)], and Ξ_C and Ξ_D are suitably defined matrices.
Note that those elements of vectors s_C and s_D that correspond to a zero charge and discharge command in vectors σ_C = [σ_C(0) ⋯σ_C(M-1)] and σ_D = [σ_D(0) ⋯σ_D(M-1)] are set to zero (see (<ref>) and (<ref>).
Given that the charge and discharge commands are mutually exclusive, we have that  s_Cs_D=0.




The described thermal storage system is active in that it can be directly operated by charge/discharge commands.
Passive thermal storages are instead physical elements, like the walls of a building, that can accumulate and release thermal energy but are not directly charged or discharged. Even though it is more difficult in principle to take advantage of passive thermal storages since there is no direct way to control them, in Section <ref> we shall show how an optimal energy management strategy can exploit them.




Note that batteries for electrical energy storage can in principle be modeled in the same way <cit.>. However, charging/discharging efficiencies depend on the battery State Of Charge (SOC) and energy losses can be related to the exchanged energy (exchange efficiency), so that a more complex model has to be specifically introduced. Also, additional constraints as for example the minimum and maximum charging time should be added to obtain a feasible operation of the battery.
PCMs thermal storages can be modeled as electric batteries with the fraction of liquid in the storage playing the role of the SOC in determining the model coefficients.





  §.§.§ Storage block: interfaces and related constraints


The proposed model for the thermal storage has as control input the energy exchange s, eventually decomposed into the charge and discharge inputs s_C  and s_D activated by the mutually exclusive commands σ_C  and σ_D. The stored energy S is the output of the model in  both cases. Since the storage capacity is limited and the stored energy is a positive quantity, the following constraints apply

    0≤S≤S^max.

In addition, the amount of energy that can be exchanged per time unit is limited, and it cannot exceed certain thresholds, i.e., the bounds

    s^min≤s≤s^max

apply to the energy exchange s, or bounds (<ref>) and (<ref>) apply to the charge and discharge inputs s_C  and s_D.

Table <ref> summarizes the relevant quantities related to the storage model. 
The type attribute denotes possible different models for the storage.



 


 §.§ Combined Heat and Power unit: Microturbine

A Combined Heat and Power (CHP) unit is a device that jointly produces electricity and heat power while consuming primal energy (fossil fuels or hydrogen) with the purpose of reducing the amount of energy wasted in the environment. In most cases one of these two products is a byproduct. For example, modern power plants recover waste heat and deliver it for district heating purposes.
CHPs with large capacity are becoming widely used and highly performing. At the same time a number of micro-CHP solutions are being developed, the most promising ones being microturbines and fuel cells that convert gas or hydrogen into heat and electricity. Combined Cooling, Heat and Power (CCHP) devices are also available that convert part of the produced heat into cooling energy.

We consider a microturbine modeled through two static characteristics describing the electrical power production and the heat production, both as a function of the fuel volumetric flow rate. Figure <ref> represents the characteristics of the C30 microturbine produced by Capstone company <cit.>.
We can see that both curves are almost linear. The electrical energy E_,ℓ(k) and the heat E_,h(k) produced by this microturbine during the kþ time slot can then be expressed as affine functions of the fuel volumetric flow rate u_(k), that is supposed to be constant in each time slot, i.e.,

    E_,ℓ(k) = m_ℓ u_(k) + q_ℓ,
        E_,h(k) = m_h u_(k) + q_h,

where m_ℓ, q_ℓ, m_h, and q_h are positive coefficients.



If we include the possibility of switching on or off the microturbine, we need to introduce the binary variable δ_(k), k=0,…,M, and modify the model as follows:

    E_,ℓ(k) = δ_(k)(m_ℓ u_(k) + q_ℓ),
        E_,h(k) = δ_(k)(m_h u_(k) + q_h).


Note that we do not model the microturbine transient from on to off, as instead suggested in <cit.>. Yet, the static model that we adopt is accurate given that  a sensible choice of  when addressing energy management is typically larger than the time scale of the microturbine dynamics.




  §.§.§ CHP block: interfaces and related constraints


The CHP block represents a microturbine and is characterized by two control inputs that can be set in each time slot: the fuel volumetric flow rate u_ and the on-off status of the microturbine δ_. It provides as outputs the electricity E_,ℓ and the heat E_,h produced per time slot.

 Since the microturbine specifications require a minimum fuel volumetric flow rate u_^min for the unit to be operative, we need to include the following logical condition:

    u_(k) ≤ u_^min  ⇔  δ_(k)=0,

which can be rewritten as:

    δ_(k)(u_^min +ε)≤ u_(k)≤δ_(k)u_^max+(1-δ_(k)) u_^min,

where u_^max is the maximum flow rate and ε>0 is set equal to the machine precision.
The product between the (piecewise) affine function and a discrete variable δ is a nonlinear mixed integer expression that can be reduced to a mixed integer linear condition <cit.>.

The constraints related to the CHP are of three types:


  * Fuel inlet bounds of u_ = [u_(0) ⋯ u_(M)]:

    u_^min≤u_≤ u_^max


  * Maximum heat and electrical energy that can be produced by the microturbine:

    0 ≤E_,h≤ E_,h^max,
       0 ≤E_,ℓ≤ E_,ℓ^max

with E_,h = [E_,h(0) ⋯ E_,h(M)], and E_,ℓ = [E_,ℓ(0)⋯ E_,ℓ(M)].

  * Logical on-off bounds (<ref>) expressed over the finite horizon k=0,…,M with δ_ = [δ_(0) ⋯δ_(M)].



Table <ref> summarizes the main characteristics of the CHP model.
Type A and B are the possible variants of the CHP model.




 


 §.§ Wind turbine


A wind turbine is used to convert the wind kinetic energy into electrical energy.
Four different operational modes are typically defined for a controlled wind turbine (see  Figure <ref>): Mode 1, when the wind speed value is within the range from zero up to the cut-in wind speed v_in and there is no power produced by the wind turbine, which is turned off; Mode 2, below the rated power P_n, thus called below-rated, where the power captured from the wind is maximized; Mode 3, above the rated wind speed, thus called above-rated, where the wind turbine is saturated to the rated power P_n, and as the wind speed increases above the nominal turbine speed v_n, the blade pitch angle is adjusted so that local angles of attack acting on local airfoil sections become smaller, and hence the loads become relatively smaller and the power keeps constant; Mode 4, when the wind speed is above the cut-out wind speed v_out, and the wind turbine is shut down, due to load and fatigue issues.

A turbine that is optimally sized for the site where it is installed is operating most of the time at the transition point between Mode 2 and Mode 3, also called at-rated <cit.>. The power generated by the wind turbine P_wt can then be calculated as follows:

    P_wt =
    
    0,    v_wind≤ v_in or   v_wind≥ v_out
    
    P_m(v_wind),    v_in≤ v_wind≤ v_n
    
    P_n,    v_n≤ v_wind≤ v_out

where P_m(v_wind) is the maximum power that can be extracted from the wind kinetic energy when the wind speed is v_wind, while P_n is the rated power.

Notice that the wind speed v_wind is acting as a disturbance on the turbine. Therefore, the power produced by the wind turbine as output given the disturbance input v_wind is a disturbance as well.
To the purpose of the energy management of the district network, we consider the static model in Figure <ref> (solid line) for the power produced by the wind turbine as a function of the wind speed.  As for the wind speed prediction, both physical and statistical models, e.g., based on Markov chain, have been considered in the literature <cit.>. Combining (<ref>) with wind speed prediction models one can determine the energy contribution of the wind turbine by computing the average power produced within a time slot, and then multiplying it by the time slot duration . 
Note that the static modeling of the wind turbine is appropriate if the time slot duration  is sufficiently large compared to the involved inertia. In our set-up of a district network,  small scale wind turbines for roof installation could be included, compatibly with  of the order of minutes.
 


 

§ DISTRICT NETWORK COMPOSITIONAL MODELING AND OPTIMAL ENERGY MANAGEMENT



In this section, we show how the components previously introduced can be interconnected in order to define a certain district network configuration. 
We consider a network of buildings located in a neighborhood and do not model the distribution network. 
Since the input/output interfaces of each component have been described in terms of thermal or electrical energy received or produced, energy balance equations and energy conversion functions can be adopted to combine the network components. For instance, the sum of the cooling energy requests of the buildings in the network should be equal to the sum of the cooling energy provided by chillers and taken from/stored in the thermal storages; each chiller receives as input a cooling energy request and provides as output the corresponding electrical energy consumption; the sum of the electrical energy consumption should be equal to the electrical energy produced by the local power generators, i.e. the CHP units and the wind turbine, taken from/stored in the batteries, and provided by the main grid.
Depending on the number of components and the adopted model for each component, the overall model of the district network has a different size and complexity, the most general one being hybrid due to the presence of both continuous and discrete variables, and stochastic due to the disturbances (e.g., occupancy, outside temperature, solar radiation, wind velocity) acting on the system, <cit.>.


Figure <ref> shows a possible district network configuration and the energy fluxes between its components and the main grid. The district network may be composed of multiple  buildings that share common resources such as cooling and heat storages, chillers, CHP units, batteries and renewable energy generators. The three nodes appearing in the figure do not correspond to any physical component but are introduced to point out that fluxes associated with the same kind of energy (electrical, heat, and cooling energy) add up to zero.
Some energy contributions can be controlled (e.g., those related to storage units), some others can be controlled only indirectly (e.g., electrical energy requested by the chiller), or cannot be controlled (e.g., renewable energy production). This is pointed out using different arrowheads in Figure <ref>.
As for buildings, some of them are controlled in that their cooling energy request can be modulated to some extent via the zone temperature setpoint. If the zone temperature setpoints are fixed and given by some comfort profiles, then the building is uncontrollable.

We assume that the district network is connected to the main grid, which supplies the electrical energy needed to maintain the balance between electrical energy demand and generation within the district network.

The district network is “smart” if it is possible to appropriately set the controllable variables so as to optimize its behavior.
A typical goal is to minimize the overall cost while guaranteeing the satisfaction of the energy needs of the users in the district.
Costs are mainly due to the electrical energy requested to the main grid and additional costs related to device operation such as startup and fuel costs. The overall cost is then given by:

    J = C_ℓ + C_ + C_ + C_f,

where the first term is the electrical energy cost C_ℓ = ∑_k=1^M C_ℓ(k); C_ = ∑_k=1^M C_(k) is the cost for the chillers startup; C_ = ∑_k=1^M C_(k) and C_f=∑_k=1^M C_f(k) are the costs for the CHPs startup and fuel consumption.


It is worth noticing that startup costs also serve the purpose of favoring solutions that avoid continuous and unrealistic switching of devices.
Note also that additional logical conditions are needed to account for them. For example, a chiller startup cost can be modeled as C_^max{δ_(k) - δ_(k-1), 0 }, where C_^ is the actual startup cost which is accounted for at k only if the chiller was off at k-1 and is switched on at k.
Similarly, for the CHP, its startup cost at k is given by C_^max{δ_(k) - δ_(k-1), 0 }.
The fuel costs of a CHP are proportional to the amount of fuel consumption during the kþ time slot, i.e., ψ_f  δ_(k) u_(k), where ψ_f is the unitary fuel cost.

As for the electrical energy cost, the cost per time slot C_ℓ(k) is typically given by a PWA function of the electrical energy exchange E_L(k) with the main grid, i.e.,

    C_ℓ(k)=max{c_1,ℓ(k)E_L(k)+c_0,ℓ(k)},

where the coefficients of the affine terms are collected in vectors c_1,ℓ(k) and c_0,ℓ(k), and the max operator is applied componentwise. This expression allows us to adopt different values for revenues (E_L(k)<0) and actual costs (E_L(k)>0), and to account for penalties when the electrical energy consumption/production E_L(k) exceeds certain thresholds. 
Note that, if C_ℓ is to be minimized, an epigraphic reformulation can be adopted to rewrite (<ref>) in terms of a set of linear inequalities.

To describe E_L for an arbitrary configuration, we adopt in this section the following short-hand notations. Components correspond to energy contributions and are defined through letters (building , chiller , storage , CHP microturbine ) with a superscript that denotes the model type (symbols are given in Tables <ref>–<ref>) and the kind of energy (electrical ℓ, cooling c, and heating h) provided as output. This is important, e.g., to distinguish between a thermal storage (^c) and an electric battery (^ℓ), and also in the case when a component allows for multiple kinds of energy as output. For instance, ^B,h stands for the heating energy produced by a CHP described by a linear on-off model. The subscript possibly denotes the energy request received as input, as in the case of a chiller that has to provide the net cooling energy requested by buildings after deduction (addition) of that provided (requested) by the thermal storage units.

We can for example derive the expression of E_L  for the configuration in Figure <ref>:

    E_L= ^A,ℓ_←{^B,c + ^A,c + ^c}+^B,ℓ+^ℓ.


If we then plug (<ref>) into equation (<ref>) and (<ref>), we get the expression for the cost function J to be minimized.

Note that J may be uncertain if there are disturbance inputs acting on the system. In such a case, one can either neglect uncertainty and refer to some nominal profile for the disturbance inputs or account for uncertainty and formulate a worst case or an average cost criterion based on J.
Furthermore, when we compose a district network model plugging together all the elements, we also get a number of constraints associated with them.
Constraints express both technical limits (i.e., maximum cooling energy that a chiller can provide) and performance requirements (i.e., comfort temperature range). Additional constraints can be added if needed (e.g., the maximum amount of electrical energy that the main grid can provide). Yet, constraints might be uncertain due to the presence of disturbances, and, hence they might be enforced only for the nominal profile, thus neglecting uncertainty, or as robust or probabilistic constraints.

Different approaches can then be adopted to address the energy management of the district network, depending also on the choice of the cost criterion (nominal/worst-case/average) and the constraints (nominal/robust/probabilistic).
Uncertainty on the parameters values could also be explicitly accounted for in the design. For instance, one could assume that parameters take equally likely values in some range and impose that performance is optimized over almost all instances except for a small set.

Furthermore, different architectures (centralized, decentralized or distributed) can be conceived and implemented for the resulting optimization problem solution, depending on the actual communication and computation capabilities available in the network, and on possible privacy of information issues like in the case of a building that is not willing to share its own consumption profile, while still aiming at cooperating for reducing the overall district cost.

The formulation of the optimal energy management problem involves defining the following quantities:


  * Global parameters, i.e., sampling time , and number of M of time slots of the look-ahead time-horizon. 


  * Optimization variables, i.e., the decision variables to be set by the optimization problem. Notice that energy balances must always hold, and this may decrease the actual degrees of freedom of the system. For example, in the case of a controllable building with a chiller plant, the cooling energy request to the chiller cannot be set freely, since it has to match the cooling energy needed for the building to track the temperature setpoint that becomes effectively the only decision variable.


  * Cost function, i.e., the quantity that has to be minimized, e.g., the electric energy costs or the deviation of the energy consumption from some nominal profile agreed with the main grid operator.


  * Constraints, i.e., the feasibility conditions that limit the solution space of the optimization problem. Notice that constraints can be classified in three categories:


  * Single component constraints, which are enforced at the level of each component separately and are related due to its dynamics and capabilities.
    For example, the energy accumulated in a storage is jointly dictated by the storage capacity and dynamics of the storage system.


  * Interconnection constraints, which relate variables of different components and originate from their cooperative interaction in the district. For example, the temperature setpoint in a controllable building cannot result in a cooling energy request that is larger than the energy that the chiller can produce and the energy that can be taken from the storage.


  * Control constraints, which are enforced to achieve some desired property of the energy management strategy. These are, for instance, the comfort constraints imposed on the temperature in a building or the constraints enforced at the end of the control time-horizon on the energy in the storage to avoid its depletion and allow for repetitive use of the control strategy in  a periodic fashion.










 

§ SOME NUMERICAL EXAMPLES



In this section, we present some examples where a district network configuration is considered and a related energy management problem is defined and solved.
All examples refer to a centralized architecture, with known profiles for the disturbances. 
We consider a one-day time horizon since this is a commonly used time horizon for building energy management, especially temperature control. 
We enforce a periodic solution to cope with the myopic attitude of the finite horizon strategy, which would empty the battery/storage and drive the zone temperatures to the limit of their admissible range at the end of the time horizon in order to save money, without caring of the next day.

Examples were chosen to be simple but realistic enough to highlight the capabilities of the proposed framework.
Many more examples could be presented with reference to different set-ups in terms of either district network configuration or energy management problem formulation. 
Distributed energy management strategies could be adopted for easing computations and preserving privacy of information, as suggested in <cit.>.  The stochastic nature of disturbances could be accounted for via a randomized approach as in <cit.>, which however refers to a single building configuration. Stochastic periodic control solutions, <cit.>, could be implemented as well. 




 §.§ Cooling of a controlled building with a chiller plant



Inspired by the numerical example presented in <cit.>, we start by considering the simple district network configuration in Figure <ref>, which consists of a controlled building and a chiller unit.



The controlled building is a medium-sized three story office building of the following dimensions: 20m long, 20m wide, and 10m tall. Each facade of the building is half glazed and the roof is flat. The biquadratic approximation (<ref>) is used for the chiller model.


Disturbances are treated as deterministic signals. Figure <ref> shows the profiles adopted for the occupancy and internal energy contributions, solar radiation and outside temperature.


In Section <ref>, we consider a single-zone setup for the building, where the three floors are treated as a unique thermal zone, with the same temperature setpoint. In Section <ref>, we move to a multi-zone setup, where each floor is a thermal zone with its temperature setpoint.
In both cases we neglect energy exchanges through thermal radiation among internal walls and we consider the ground floor thermally isolated from the ground.



  §.§.§ Single zone building



The purpose of this example is twofold:


  * Showing the role of the building structure as a passive thermal storage, that can accumulate and release thermal energy;

  * Compare the energy management strategies obtained with two different control objectives.



The problem is then formulated as follows:


  * Global parameters, the sampling time is set to = 10 minutes, and the time horizon is set to 1 day, i.e., M=144.


  * Optimization variables, the only optimization variable is the temperature set-point of the single zone T_z as defined via the control input u over the considered finite horizon.


  * Cost function, we here consider the two different cost functions:


  * cooling energy provided by the chiller:

    J_1 = ∑_k=0^M E_,c(k).


  * electricity consumption:

    J_2 = ∑_k=0^M E_,ℓ(k),


where the electricity consumption is related to the cooling energy of the chiller static characteristic that maps one into the other according to a specific COP (see Section <ref>).


  * Constraints, the following constraints are included in the optimization problem:


  * Single component constraints, the constraints of the single components are given by

    E_c≥ 0
        0≤E_,c≤ E_,c^max.




  * Interconnection constraints: the chiller satisfies the cooling load demand, i.e.,

    E_,c(k) = E_c(k),    k=0,…,M.



  * Control constraints: zone temperature should be within some comfort range, and a periodic solution is enforced by setting the same value for the zone temperature setpoint at the beginning and end of the time horizon:

    u^min≤ u≤ u^max
        u(M) =  u(0) =T_z(0)
       T(0) = T(M)


We consider an ideal setting where both T_z(0) and T(0) can be set so as to obtain periodic solution.


The resulting optimization problem is a convex constrained program that can be solved, for example, with CVX[<http://cvxr.com/>] with the SDPT3 solver.



Figure <ref> shows the resulting optimal temperature profiles T_z for the two cases.
Both solutions stay within the prescribed comfort temperature bounds. Notice that the discrepancy between the two curves is at most of about 1.6^∘C. Despite this distance being small, from Figure <ref> one can notice a clear difference in the required cooling energy for the two cases.

In the case of minimization of the electricity consumption (J_2), a “precooling” phase occurs from time 18:00 to time 10:00 of the next day (if we think about the solution applied over multiple days), which leads to a larger cooling energy request.

Intuitively, the second policy stores some cooling energy in the building structure, ahead of time, thus smoothing the cooling energy request in the central part of the day, when occupancy is larger, to get the chiller operating with high efficiency. The “building thermal mass” is  exploited as a passive thermal storage to add further flexibility to the system <cit.>.

On the other hand, the first policy exploits the fact that at night the temperature is lower, comfort constraints are satisfied (indeed they are set to be larger because there are no occupants in the office building), and the chiller does not need to provide any cooling energy to the load. Figure <ref> shows the electric energy consumption in the two cases, highlighting that the chiller is working at its minimum for most of the time in the cooling energy minimization policy (J_1). The integral of the curves in Figure <ref> is the electricity consumption and is larger for the cooling energy minimization policy. Indeed, Figure <ref> shows that the second policy makes the chiller operate close to the maximum COP value of the chiller, thus requiring much less electrical energy.



In summary, depending on the cost function adopted in the energy management strategy design, one can have significantly different behaviors of the same district network configuration, with a different performance, even with a limited difference in the temperature setpoint profiles.







  §.§.§ Multi-zone building
 

In this example, the controlled building has three thermal zones, one per floor, with three zone temperature setpoints, and our goal is to investigate the impact of a time-varying electricity price over the day time.

The problem is formulated as follows:


  * Global parameters, the sampling time is set to = 10 minutes, and the time horizon is set to 1 day, i.e., M=144.


  * Optimization variables, the optimization variables are the temperature setpoints of the three zones T_z = [T_z,1 T_z,2 T_z,3], as defined via the control input u over the considered finite horizon.


  * Cost function, we minimize the cost of the electrical energy needed for cooling the building

    J = ∑_k=0^M p(k) E_ℓ(k),

where p(k) is the time-varying unitary cost (see Figure <ref>).


  * Constraints, the following constraints are included in the optimization problem:


  * Single component constraints, the constraints of the single components are given by

    E_c≥ 0
        0≤E_,c≤ E_,c^max.



  * Interconnection constraints, the chiller satisfies the cooling load demand, i.e.,

    E_,c(k) = E_c(k),    k=0,…,M.



  * Control constraints,
zone temperature should be within some comfort range, and a periodic solution is enforced by setting the same value for the zone temperature setpoint at the beginning and end of the time horizon:

    u^min≤ u≤ u^max
        u(M) =  u(0) =T_z(0)






Figure <ref> shows the optimal temperature setpoints for the three zones (top graph), the stack of the cooling energies associated with each zone and how they compose the cooling energy request E_,c  to the chiller (middle graph), and the unitary price for the electrical energy. The obtained solution shows that all the zones are precooled in the first hours of the day. However, all the cooling energy provided by the chiller is conveyed to Zone 2. This can be explained by observing that Zone 2 is the central story, and it exchanges thermal energy with the other zones through its ceiling and floor, thus cooling them. In addition, considering jointly the cooling energy request profile with the price of energy, it is possible to see that whenever the price of energy is decreasing the energy request increases, and vice versa. Finally, notice that at the end of the day, Zone 2 is being cooled again, to bring the temperature back to its value at the beginning of the day, with a larger amount of cooling energy request when the price gets lower.



The example highlights the advantages of defining multiple thermal zones in a building. Figure <ref> shows the total cost obtained for the single-zone and multi-zone cases, and the number of iterations needed to find the optimal solution in the two cases. The total cost decreases by 8% with the multi-zone configuration, but the computational complexity increases as witnessed by the higher number of iterations of the solver and the total CPU time. This analysis can be performed at design time to decide the number of zones to be included, based on the available computational resources and the economical advantage.




 §.§ Microgrid with uncontrolled building




In this example, we consider a microgrid configuration as represented in Figure <ref>, which includes a set of uncontrolled buildings, a cooling and a heat storage, three chillers of different sizes (small, medium, and large, indexed by 1, 2 and 3, respectively, and whose COP curves are presented in Figure <ref>), a microturbine, and an electrical storage (a battery).



Buildings are modeled via their requests of cooling energy E_c (for air temperature control) and heating energy E_h (for getting warm water). The chillers and microturbine can be switched on and off by setting the associated logical variables  δ_,i, i =1,2,3 and δ_ to 1 and 0.
We adopt the PWA approximation (<ref>) for the three chillers, with 10 knots, obtained from their Ng-Gordon models with T_o = 22^∘C, T_cw = 10^∘C, and the following parameters:

    Chiller 1 (small):     a_1=0.0056,    a_2=10.11,   a_3=7.00,   a_4=0.9327, 
       Chiller 2 (medium):     a_1=0.0109,    a_2=20.22,   a_3=3.80,   a_4=0.9327, 
       Chiller 3 (large):     a_1=0.0230,    a_2=40.44,   a_3=1.98,   a_4=0.9327.

The measure units of the coefficients are [a_1] = W/K, [a_2] = W, [a_3] = K/W, [a_4] = 1.

The problem is then formulated as follows:


  * Global parameters, the sampling time is = 1 hour, and the time horizon is 1 day, hence the number of time lots is M=24.



  * Optimization variables, the optimization variables are the chillers and the microturbine on/off status δ_,i, i= 1,2,3, and δ_, the cooling energy requested to the chiller E_,c, the CHP fuel inlet u_, and the energy exchanged with the storage units s_i, i ∈{c,h,ℓ}.


  * Cost function, composed of different components:


  * The cost for energy trading:

    C_ℓ(k) = p(k) (E_ℓ(k) - E_,ℓ(k) -s_ℓ(k))

where p(k) is the unitary cost for trading electrical energy (the same used for the previous experiment), and E_ℓ(k) is the electrical energy demand;


  * The fuel cost for the microturbine given by:

    C_f(k) = ψ_f(k)  δ_(k) u_(k)

where ψ_f is the unitary cost for the fuel, here considered constant and unitary: ψ(k) = 1.


  * The startup costs of the microturbine:

    C_(k) = C_^max{δ_(k) - δ_(k-1), 0 }

where C_^ =1.


  * The startup costs of the chillers:

    C_(k) = ∑_i=1^3 C_,i^max{δ_,i(k) - δ_,i(k-1), 0 }

where C_,1^ =0.05, C_,2^ =0.1, C_,3^ =0.2.

Therefore, the overall cost function becomes:

    J = ∑_k=0^M C_ℓ(k) + C_f(k) + C_(k) + C_(k)



  * Definition of constraints, the following constraints are included in the optimization problem:


  * Single component constraints given by:

    δ_(k)(2 +ε)≤ u_(k) ≤ 10δ_(k)+ 2(1-δ_(k)) ,
       0 ≤ E_,h(k) ≤ 2160MJ,    0 ≤ E_,ℓ(k) ≤ 1080 MJ,
       -500MJ≤ s_h(k) ≤ 500 MJ,   
    0≤ S_h(k) ≤ 1500 MJ,
       -360MJ≤ s_c(k) ≤ 360 MJ,    0≤ S_c(k) ≤ 1800 MJ,
       -250MJ≤ s_ℓ(k) ≤ 250 MJ,   0≤ S_ℓ(k) ≤ 1500 MJ,
       εδ_,i(k) ≤ E_,i(k) ≤δ_,i(k) · 180 MJ,   ∀ i={1,2,3}.



  * Interconnection constraints, given by thermal energy balance equations:

    ∑_i=1^3 E_,c,i(k) + s_c(k) = E_c(k),
       E_,h(k) + s_h(k) = E_h(k)



  * Control constraints, a periodic solution is enforced by setting the same value for the on/off status of the devices and for the stored energies: 
    δ_(0) = δ_(M)  δ_,i(0) = δ_,i(M),    i = 1,2,3
        S_h(0)= S_h(M)   S_c(0)= S_c(M)    S_ℓ(1)= S_ℓ(M)




The resulting optimization problem is a Mixed Integer Linear Programming (MILP) problem, and can be solved using YALMIP to formulate the optimization problem and CPLEX as a solver.




Figure <ref> shows the cooling energy demand E_c and the heating energy demand E_h of the uncontrolled buildings over the considered time horizon (top left graph), the storage of cooling, heat and electrical energy (top right graph), with the respective control inputs (middle right graph), the fuel inlet to the microturbine (middle left graph), the energy demand for the three chillers (bottom left graph), and finally, the bottom right graph shows the logical on/off status of the chillers and of the microturbine.

Notice that the small and medium size chillers (associated with δ_,1 and δ_,2) are never switched off, while the large chiller is used to charge the cooling storage before the peak of cooling energy demand, and right after that, when the storage got empty.
The microturbine is switched on only during the peaks of heat energy demand (see the profile of δ_).


Figure <ref> shows the performance of the three chillers in terms of COP, as a function of time. It is possible to see that the small and medium chiller are always operating close to their maximum COP, while the large chiller either operates at its maximum COP or it is in the off mode, where its consumption is zero.

The cooling energy storage is charged at the beginning of the day when the cooling load is low, and then discharged when the load increases. It is finally recharged to meet the periodicity imposed via the control constraints.
On the other hand, the heating storage is used to meet the heating energy demand, and the microturbine is activated only when the heating storage is empty. Also in this case, part of the energy produced when the CHP was active, is stored to meet the energy demand of the next day.

 

§ MULTIRATE CONTROL



Increasing the number buildings and/or thermal zones per building necessarily leads to a larger computational effort for solving the energy management control problem since the number of optimization variables increases as well.
This may become an issue when a receding horizon strategy is adopted and optimization is performed online at every control instant. Indeed, real-time constraints can hamper the applicability of the approach.

A possible way to avoid this issue is to use larger values of the sampling time  for the discretization of the model, thus yielding a lower number of optimization variables for the same time horizon. As a side (positive) effect, discretizing with a larger sampling time make less stringent
the real time constraints since it gives more time to perform the computations and apply the solution.


Unfortunately, using a larger sampling time degrades the model accuracy, thus eventually deteriorating the control performance.
This issue can be tackled by taking a multirate control approach, where model and controller operate with different sampling periods. Specifically, if we let  be the sampling period of the model and introduce the rate M_R∈ℕ, then, in multirate control, the control action is only set every M_R time slots of length , or, equivalently, = M_R is the sampling period for the controller.
This choice allows for an accurate representation of the model dynamics, while still decreasing the number of optimization variables, and, as a consequence, the computational complexity, by a factor M_R. 
Clearly, the reduction of the number of optimization variables has some impact on the achievable performance in terms of cost and also reactiveness to possible disturbances with fast dynamics. The choice of the rate M_R must compromise between computational effort reduction and performance degradation, compatibly with the available resources.




 §.§ Example

Let us focus on the example presented in Section <ref>, with cost function given by the electrical energy

    J = ∑_k=0^M E_ℓ(k).

We sample the model with = 10 minutes, and we study the effects of employing different rates M_R for applying the control input, namely M_R = 1,6,12,24,36,48, corresponding to = 1/6,1,2,4,6,8 hours, thus progressively reducing the number of optimization variables.

Figure <ref> shows the optimal temperature profiles for the different rates. Notice that the curves associated with = 10 minutes (M_R = 1) and with = 1 hour (M_R = 6) are practically indistinguishable, but in the latter case we have reduced the number of optimization variables by a factor 6. 
The reduction of the optimization variables causes an increase of the costs, as shown in Figure <ref>. However, this increase is negligible up to = 2 hours, while the computation effort is almost constant for values of  larger than or equal to 1 hour, if evaluated in terms of total CPU time[The total CPU time presented in Figure <ref> is computed as the average total CPU time over 100 experiments for each considered , for the sole solver.].



We can also analyze the chiller performance presented in Figure <ref>.

The higher the rate M_R, the lower is the flexibility of the control input to do a fine adjustment of the temperature setpoint and compensate the disturbance variability. This is why the chillers are not constantly operating at a high efficiency levels when M_R is larger. This results in a less performing chiller, so one should look for a tradeoff among computational effort and cost efficiency.

Finally, we can conclude that the adoption of multirate control solutions is problem specific, depending on the available computational power, and on the complexity of the optimization problem to be solved. 
 

§ CONCLUSION



In this paper we presented a  modeling framework for the optimal operation of a district network, with reference in particular to the cooling of multiple buildings that are sharing resources like chillers or storages. Various components have been introduced and modeled in terms of energy fluxes so as to ease their composition via energy balance equations. A control-oriented perspective is adopted in that control and disturbance inputs are explicitly accounted for, in terms of their energy contribution.
We also described how to formulate an optimal energy management problem as a constrained optimization program where control inputs are the optimization variables  and need to be set so as to minimize some energy-related function (e.g., electric energy cost, deviation from some nominal profile of electric energy consumption), while satisfying comfort and actuation constraints.
Finally, a multirate approach was proposed to reduce the number of optimization variables while preserving the model accuracy. This has potential for real time applicability of the method when implemented according to the receding horizon strategy of model predictive control.
This will allow to compensate for unpredictable human-building interactions as discussed in <cit.>. 

Some numerical examples were also presented to show the versatility of the proposed framework.
Currently, we are addressing optimal energy management of a district network in presence of stochastic disturbances, the key challenge being how to account for them when embedded in a distributed setting with limited communications capabilities. The approach in <cit.> could be useful to this purpose.

 


§ MODEL VALIDATION




Reliability of the model is crucial when adopting model-based control design strategies. At the same time, if a model is accurate but very complex, then, design can become impractical.

As for what concerns the network district modeling for energy management purposes, the more difficult component to model is the building, since various factors need to be accounted for, including size and structure of the building, walls composition, presence of electrical devices, occupancy, and environmental conditions, like outdoor temperature and solar radiation. Also, model complexity grows as the size of the system increases.

Models and modeling frameworks for buildings have been proposed in the literature <cit.>. Most of them include a detailed characterization of the fluid dynamics phenomena, e.g., the evolution of the temperature and humidity of the thermal zones, and they typically require specialized Computational Fluid Dynamics (CFD) tools for simulation. Even though these approaches provide very accurate simulation results, they are difficult to use for control design purposes, due to their complexity. In this paper we adopted a control-oriented perspective and presented a simple model of the building where thermal zone temperatures act as control inputs and enters linearly the system dynamics.

Validation of a model of the building dynamics against experimental data is quite challenging, also because setting up a measurement facility for a building can be complex and expensive. In order to validate the presented model, we hence resort to the methodology introduced by the American Society for Heating Refrigerating and Air-conditioning Engineers (ASHRAE), and, more specifically, the validation method defined in the ANSI-ASHRAE 140 standard.
This standard specifies test cases and procedures for the evaluation of the technical capabilities and range of applicability of computer programs that compute the thermal performance of buildings and their HVAC systems. The current set of tests included in the standard consists of (i) comparative tests that focus on building thermal envelope and fabric loads, and mechanical equipment performance, and (ii) analytical verification tests that focus on mechanical equipment performance. Different building energy simulation programs, with different levels of modeling complexity, can be tested. For all tests included in the specifications, results provided by other certified simulation tools are presented, and they represent the baseline for validating new modeling and simulation software. A detailed description of the simulation tools included in the specification can be found in <cit.>.

We here provide the results of some of the main tests defined in the ANSI-ASHRAE 140 standard, and compare them with the baseline provided in the standard. Let us first introduce the test case and then describe the validation procedure.

We consider a building located at an altitude of 1649m above the sea level, and weather data series resuming the weather conditions for a whole year are available and provided by the standard. The data set contains: external dry bulb temperature, wind speed, wind direction, direct and diffuse solar radiation.


The building has a 48m^2 floor area, a single story with rectangular-prism geometry, and two south-facing windows, 6m^2 each (see Figure <ref>). Two set-up are considered, which differs in materials composition and walls thickness: lightweight (case 600 in the standard) and heavyweight (case 900 in the standard). The standard specifies the composition in terms of thickness, density, thermal conductivity and specific heat capacity of all layers of each wall, for both the lightweight and heavyweight cases. These values are listed in Tables <ref> and <ref>, respectively. According to the specification, density and specific heat of the underfloor insulation have the minimum admissible value allowed in the building model that is tested, and in any case a value not smaller  than zero. Also, the contribution of the internal loads and people, within the thermal zone is constant over the year and equal to Q_int + Q_p = 200W.





The standard also provides the values for the internal and external solar absorption and infrared emission coefficients
 α_i^S = α_i^L = 0.6 and ε_i = 0.9, i=1,…,m, and for the interior and exterior combined radiative and convective heat transfer coefficients, from which the radiative and convective coefficients can be recovered. Finally, the standard contains also the windows properties, the values of incidence angle-dependent optical properties, and the interior solar distribution. The reader is referred to the ANSI-ASHRAE 140 standard  for a complete list of building properties.

We focus on two procedures for validation described in the standard. The first one is denoted as Free Float (FF) in that the heating and cooling equipment is switched off and the zone temperature evolves freely subject to internal/external disturbances. The purpose of this test is to validate the physical model without the effect of any control action, and such validation is performed comparing some statistics (maximum, minimum, and annual average) of the zone temperature T_z over a year against other simulation tools. The second procedure prescribes to simulate the building together with the heating/cooling system by applying a simple control strategy: the controller has to maintain the air temperature inside the building between 20^∘C and 27^∘C. Specifically, the control strategy is:


  * Heat = ON if temperature <20^∘C; otherwise, Heat = OFF.

  * Cool = ON if temperature >27^∘C; otherwise, Cool = OFF.

The ANSI-ASHRAE 140 standard specifies that the air conditioning system produces only pure heating load and sensible cooling load outputs. That is, all equipment is 100% efficient with no duct losses and no capacity limitations. In this controlled case, the validation is performed comparing the hourly-integrated peak of the cooling and heating power provided to the building.
The thermostat was implemented as two saturated PI controllers with antiwindup, where the control variable is the amount of cooling and heating power to be injected in the system, equivalently to the implementation adopted in <cit.>. 

In the following we will denote as 600FF and 900FF the case when the free float validation procedure is applied to the lightweight and heavyweight buildings, and as 600 and 900 the case when the control is applied.




 §.§ Validation results


In this section we present the numerical results obtained in the 600FF and 900FF and 600 and 900 test cases.
For running the validation process, it is necessary to rewrite the model with the heat flow rate Q as control input, and the temperature of the zone T_z as the output of the system. To this aim, we consider a simulation model composed of a state vector including the temperature of the different slices of the walls as described in (<ref>), and the temperature of the zone T_z. The evolution of T_z is governed by the continuous-time version of (<ref>), made explicit with respect to Ṫ_z:

    Ṫ_z = -C_z^-1Q_z = -C_z^-1(Q - Q_w - Q_p - Q_int),

with Q_w, Q_p, and Q_int being the heat flow rates towards the zone of the walls, the occupancy, and of other internal equipment producing heat. Considering the expressions (<ref>), (<ref>), and (<ref>), one can write the expression of Q_z as a function of the states T and T_z, of the input Q and of the disturbances. This continuous time model is implemented in Modelica[<https://modelica.org/>], in order to carry out the validation process.

Table <ref> reports the obtained results in terms of maximum, minimum, and mean annual temperature. 
Our model (last column of the table) is compared  with the other ones provided in the standard under the free float validation procedure.


In the 600FF test case, the results obtained with the model considered herein are comparable with the ones obtained with the other building simulation models. As for 900FF, only the minimum temperature is comparable with the other results, while the maximum temperature is slightly higher than the values obtained with the other models, whereas the mean annual temperature is lower. Overall the obtained statistics produce reasonable results in the free float case, even though the model adopted for the presented framework is much simpler than the other simulation models.


Table <ref> summarizes the validation results when the presented control strategy is in place.


The hourly peak of cooling and heating power are comparable with those of the other tools in both the test cases.

Overall, the obtained results in this validation phase show that state-of-the-art simulation tools provide similar results to those obtained with our model, which has the key advantage of being simpler and hence more suitable for design purposes. 

 

§ DISCRETIZATION OF THE WALLS TEMPERATURE DYNAMICS



To solve the discrete-time optimal energy management problem, we need to consider a discretized version of (<ref>). Given the linearity of (<ref>), it holds that

    T((k+1))    = e^AT(k) +  ∫_k^(k+1) e^A((k+1)-τ)(BT_z(τ) + Wd(τ)) dτ.

If we assume that T_z (i.e., our control variables) and d are linearly varying within each time slot, then the integral in (<ref>) can be computed analytically. Formally, given

    T_z(τ) = T_z,k+1-T_z,k/(τ - k) + T_z,k,

where T_z,k = T_z(k), and

    d(τ) = d_k+1-d_k/(τ - k) + d_k,

where d_k = d(k), ∀τ:   k≤τ < (k+1) and k=1,…,M. If we set T_k = T(k) and Q_k = Q(k), k=1,…,M, then the dicretized system can be expressed as follows

    T_k+1   = Γ_xT_k + Γ_u,1T_z,k+1
    			+ (Γ_u,0-Γ_u,1)T_z,k + 
                + Γ_ω,1d_k+1 + (Γ_ω,0-Γ_ω,1)d_k 
    Q_k    =  CT_k   + DT_z,k

where

    Γ_x = e^A
       Γ_u,1 = 1(∫_0^ e^As(-s)ds) B
       Γ_u,0 = (∫_0^ e^As ds) B
       Γ_ω,1 = 1(∫_0^ e^As(-s)ds) W
       Γ_ω,0 = (∫_0^ e^As ds) W.



Applying the transformation ξ_k = T_k-Γ_u,1T_z,k-Γ_ω,1d_k we
obtain

    ξ_k+1   = Γ_xξ_k
    		+ ((Γ_x-I)Γ_u,1 + Γ_u,0)T_z,k +
                 + ((Γ_x-I)Γ_ω,1 +Γ_ω,0)d_k 
    Q_k    =  Cξ_k   + (CΓ_u,1 + D)T_z,k
    		+ CΓ_ω,1d_k

Dropping the bold notation for vectors and matrices, (<ref>) can be rewritten as the following discrete-time system

    {   x(k+1) = Ãx(k) + B̃u(k) + W̃ω(k) 
       y(k) = C̃x(k) + D̃u(k) + Ṽω(k) 
    .

where x(k) = ξ_k, u(k) = T_z,k, ω(k) = d_k, y(k) = Q_k, and the matrices are:

    [                 Ã = Γ_x, B̃ = (Γ_x-I)Γ_u,1+Γ_u,0,  W̃ = (Γ_x-I)Γ_ω,1+Γ_ω,0;                  C̃ = C,         D̃ = CΓ_u,1 + D,              Ṽ = CΓ_ω,1. ]


From (<ref>) one can derive the expression of x(k) and y(k) as a function of the initial state and the inputs up to  k:

    x(k) = Ã^k x(0) + ∑_h=0^k-1Ã^k-1-h( B̃u(h)+W̃ω(h) ) 
       y(k) = C̃Ã^k x(0) + ∑_h=0^k-1C̃Ã^k-1-h( B̃u(h)+W̃ω(h) ) + D̃u(k) + Ṽω(k).


By recalling that x(0)=ξ_0 = T_0-Γ_u,1T_z,0-Γ_ω,1d_0, we obtain the following expression of y(k) as a function of the original state variables and the inputs up to time k:

    y(k) =
        C̃Ã^kÃT_0+C̃Ã^k-1( (Γ_u,0-Γ_u,1)u(0) + (Γ_ω,0-Γ_ω,1)ω(0) ) 
       + ∑_h=1^k-1C̃Ã^k-1-h( B̃u(h)+W̃ω(h) ) + D̃u(k) + Ṽω(k).


 

elsarticle-num



