A Discussion of Mixed Integer Linear Programming Models of Thermostatic Loads in Demand Response

In smart grids, it is expected that electricity retailers will offer time-differentiated tariffs with significant price variations in short periods. Consumers are then encouraged to engage in demand response strategies by making the most of their flexibility in the operation of end-use appliances to minimize the electricity bill without compromising on comfort. Air conditioning systems are particularly suited to be controlled, by profiting from the thermal inertia of building units. This paper presents novel mixed integer linear programming formulations to optimize the operation of a thermostatically-controlled air conditioning system, thoroughly discussing their main features and associated computational difficulties resulting from their combinatorial nature.


Introduction
Buildings in the commercial and residential sectors account for about 40% of the total energy consumption in European Union countries [1,2]. Heating, Ventilation and Air Conditioning (HVAC) systems represent one of the most significant loads contributing to electrical energy consumption. HVAC systems are particularly suited to be controlled by making the most of: (i) the thermal inertia of spaces to be conditioned (leading to some time dissociation between the energy consumption and the energy service provided), and (ii) the flexibility of consumers' preferences regarding the provision of energy services in face of dynamic time-of-use tariffs (displaying some willingness to bear the indoor temperature somewhat farther from the reference settings for a limited period of time). Since time-differentiated tariff schemes are expected to become relevant commercial offers in electricity retail markets in smart grids, demand response actions should be developed to achieve optimal operation of HVAC systems. Optimal HVAC control can be beneficial for consumers (by reducing the energy bill without jeopardizing comfort), retailers (by managing buying and selling prices) and grid operators (by contributing to release the strain in distribution networks and enhancing the utilization of renewable sources). Consumer engagement in demand response programs can be made operational by means of automated energy management systems endowed with optimization algorithms and the capability to control appliances. Therefore, adequate and tractable models for optimal thermostat programming should be developed as operational tools (to be embedded with sensors and control) for demand response programs.
The potential of thermostatic loads for demand response actions has been exploited in the literature, due to their significant consumption and the possibility of being controlled. For this purpose, mathematical models have been proposed, in particular mixed integer linear programming (MILP) formulations, i.e. with integer and continuous variables. The authors in [3] consider heat pumps for heating residential buildings with a floor heating system using a linear state space model in an Economic Model Predictive Control framework. Thermal models are developed in [4] for a smart house for determining the value of thermal inertia in demand response dynamic pricing using a MILP formulation. The authors in [5] present a user-centric demand response control for scheduling the electric space heating load under price and load uncertainty to minimize a weighted-sum of the expected payment, loss of comfort, and financial risk of a customer while considering the end-user's preferences. The household thermal behavior is modeled by means of a two-capacity building model. An HVAC system is considered in [6] that is controlled (in combination with other type of loads, PV generation and storage) by a home energy management system, thus enabling residential consumers to participate in demand response programs. A price-based demand response strategy for multi-zone office buildings to optimize the energy cost of HVAC units and the thermal discomfort of occupants is formulated in [7] as a MILP model. The authors in [8] develop an approach based on a partial-differential equation model of thermal diffusion to determine the thermostat settings to minimize the electricity bill for a consumer with energy time-of-use and power prices, in which the optimal thermostat programming for HVAC is formulated as a constrained dynamic optimization problem. The authors in [9] consider the optimization of the investment and operation planning of a decentralized energy system, which is subject to different sources of uncertainties, encompassing photovoltaic generators and load flexibility using heat pumps in combination with thermal energy storage units for space heating and domestic hot water, which is tackled by two-stage stochastic programming.
This paper presents novel comprehensive MILP models in which the physical characteristics of a thermostatically-controlled air conditioning system are considered. These models have the common aim of determining the optimal thermostat operation using different forms of thermostat control. It will be discussed that the combinatorial nature of these models makes it difficult to solve them by a state-of-the-art exact solver. This fact impairs embedding realistic mathematical models in automated building energy management systems when a fine-grain time discretization of the planning period is considered.
The study is made in the perspective of the building owner or tenant, i.e. the one who aims to minimize the energy bill. These models can also be useful to optimize demand response strategies that can be of interest for the grid operator (possibly through the mediation of an aggregator of end-users' flexibility that uses this asset to participate in ancillary service provision markets). Thus, this paper adopts an Operational Research methodology focusing on the development of mathematical programming models.
The paper begins by developing a simple thermodynamic model to derive a timediscretized equation expressing the instantaneous indoor temperature (at instant t i ) as a function of the indoor temperature, the external temperature and the air conditioning system operation at the preceding instant (t i−1 ). MILP models are then presented accounting for the hysteresis behavior of the thermostat between minimum and maximum temperatures defining the dead-band around a reference temperature (set-point). It is shown that, if just modelling this behavior, the model is solved in an acceptable computation time. However, if a minimum ON or OFF period is considered (also to avoid excessive switching), offering increased operation flexibility, the computation effort becomes impracticable.
In Sect. 2 the building thermal model is developed. MILP models of a thermostatically-controlled air conditioning system with different features are described in Sect. 3. Section 4 reports extensive computational experiments of the different MILP models. Conclusions are drawn in Sect. 5.

Development of the Building Thermal Model
A space, a building unit, is considered for being heated/cooled by means of an air conditioning (AC) system. Assuming the building unit as a (homogeneous) control volume, whose instantaneous indoor temperature is uniform and denoted by θ in (t), the overall thermal energy balance on a heat rate basis (at some generic instant t) is [10]: for AC in heating mode, and (2) for AC in cooling mode. In these expressions: q ext (t) is the instantaneous external load [kW]; q g (t) is the instantaneous internal load (rate of heat generation, by occupants, lighting, appliances) [kW]; C = ρc p V is the overall thermal capacity [kJ/ • C], in which ρ and c p are, respectively, the mass density [kg/m 3 ] and the specific heat at constant pressure [kJ/(kg. • C)] of indoor air, if no other thermal inertia of the building is considered, or some weighted values for indoors, and V is the volume of the building unit [m 3 ]; q AC (t) is the instantaneous air conditioning heating/cooling load [kW].
Considering a finite-difference approach over a time-step Δt to integrate equation (1) in a fully explicit time discretization -all variables at the "previous" instant t i−1 , except for temperature θ in (t i ) at the "present" instant of Δt -it can be assumed that the rate of the energy storage, i.e. the right-hand side in (1) is given by U is the (weighted-average) overall heat transfer coefficient of the building unit envelope [kW/(m 2 · • C)] and A is the surface area of the envelope [m 2 ]. Therefore, U.A represents the overall thermal conductance of the unit envelope [kW/ • C].
Introducing the expressions above in (1), equations (3), (4), (5) are obtained: In a dynamic simulation of buildings, the internal thermal load q g [kW] is usually defined according to the operation profile (e.g., daily/weekly profile of occupation, of lighting, etc.), i.e. all internal loads (thermal power) associated with the type of utilization of the building. Without loss of generality for the application of the model, and for the sake of simplification in this context, the following assumptions are considered: 1. the internal heat loads are neglected: where COP is the AC coefficient of performance, P demand and P AC nom are the power demand and the AC nominal power, respectively, and s t i−1 is the AC control variable  (5) can be rewritten as: Therefore, since COP , P AC nom , Δ t and C are always positive quantities, the sign in equation (10) is positive as this equation has been derived from (1) for heating mode. Otherwise, γ = − Δ t C COP (negative) for cooling mode. The coefficients α, β and γ derive from the building characteristics (area, envelope, etc.) and the COP of the AC. Considering a building unit with a 225 m 2 floor area and 3 m height, and the properties of the air inside (density of 1.19 kg/m 3 and specific heat of 1.007 kJ/(kg. • C)), the air thermal capacity is C ≈ 810 kJ/ • C. Taking U -values of 0.35 and 0.30 W/(m 2 . • C) for the building facades and roof, respectively, a value of U.A ≈ 0.129 kW/ • C is obtained. Assuming a COP = 2.5, typical values of parameters α, β and γ are displayed in Table 1 for different discretization intervals of the planning period.

Mathematical Models of a Thermostatically-Controlled AC System
This section presents different MILP models aimed at determining the optimal operation of the AC within a planning period to minimize energy costs, considering distinct forms of control. Although other types of loads could be considered (such as shiftable and interruptible appliances), the goal herein is to focus on the mathematical modelling of the AC operation. Shiftable loads include dishwashers, laundry machines and clothes driers, i.e. appliances for which the operation cycle cannot be interrupted once initiated. The supply of interruptible loads, which include electric water heaters and the battery of electric vehicles, can be interrupted within preferred time slots provided a certain amount of energy is supplied. In general, shiftable and interruptible loads require MILP models, which can be solved very fast by state-of-the-art solvers (e.g. Cplex), although requiring many binary variables. Examples of such models are the ones presented in [11] and [12]. However, the MILP modelling of the thermostat behavior is more complex and imposes a higher computational effort to obtain the optimal solution. In particular, MILP models are presented accounting for the hysteresis behavior of the thermostat between minimum and maximum temperatures defining the dead-band around a reference temperature (set-point) established by the user according to his comfort preferences.
The planning period consists of T time intervals t = 1, . . . , T of a given duration; this discretization can be, for instance, 15-, 5-or 1-min. The length of the time interval is denoted by Δt, in hours. For instance, if a 1-min discretization is used then Δt = 1/60 h. Figure 1 displays the behavior of the thermostat of an AC in heating mode. The hysteresis operation of the thermostat prevents excessive switching when the indoor temperature is around the set-point, which may be established as the midpoint within the dead-band defined by the comfort range of indoor temperature [θ min ,θ max ] specified by the user. These comfort bounds may depend on t, i.e. may be different during the planning period ([θ min t , θ max t ], t = 1, . . . , T ); for the sake of clarity, we will consider the comfort bounds constant throughout the planning period in the formulations presented below. Hysteresis determines the AC control variable: s t = 1, i.e. the AC continues ON as the indoor temperature θ in t increases above θ min until it reaches θ max (a in Fig. 1); s t = 0, i.e. the AC continues OFF as the indoor temperature decreases below θ max until it reaches θ min (b in Fig. 1).

Modelling the Thermostat Behavior
Model 1A forces s t = 1 when the indoor temperature (θ in t ): -is below the lower bound of the comfort band (θ in t < θ min ) OR Fig. 1 Behavior of an AC thermostat for heating mode -is within the comfort band (θ min ≤ θ in t ≤ θ max ) AND in the previous instant the AC was also ON (s t −1 = 1) Model 1A (heating mode): s.t.: The binary variables s t control the thermostat switching: The auxiliary binary variables y t and w t establish the consistency conditions associated with thermostat switching to the ON position (s t = 1): M is a positive large number, which is used in (13) and (14) to enforce the logical conditions. M must be larger than any temperature value; for instance, M = 100 (or any higher value) is suitable as temperatures are in • C. The constraints of Model 1A ensure that: 1. if θ in t < θ min , then y t = 1 and w t = 1 by (13) and (14); hence, by y t +w t −s t ≤ 1 (16), s t = 1; 2. if θ min ≤ θ in t < θ max and s t −1 = 1, then y t = 1 and w t = 1 by (13) and (15); hence, by y t + w t − s t ≤ 1 (16), s t = 1; 3. in other cases, s t is free to be 0 or 1. The objective function pushes these variables to 0, in order to turn the AC to the OFF state and minimize cost.
Additional input information to be specified include: -The temperatures defining the thermostat dead-band, i.e. the minimum temperature θ min ( • C) below which the thermostat should be ON (s t = 1) and the maximum temperature θ max ( • C) above which the thermostat should be OFF (s t = 0). -The initial indoor temperature θ in 0 ( • C) and the initial status (0, 1) of the thermostat s 0 . This model assumes that the "natural state" of the control variable s t is 0 (since it has a positive coefficient in the objective function to be minimized). The objective function forces the binary variables to 0 whenever the constraints do not impose these variables to be 1. However, in this model it may happen that, when indoor temperature drops below the upper bound temperature (θ max ) after having been above it (and hence s t = 0), the control variable switches from 0 to 1 within the thermostat comfort band (note that there are no constraints forcing s t = 0 because this would be the normal state of the variable). This may be beneficial for the objective function by avoiding switching on at some later time intervals when the electricity cost is higher. Therefore, Model 1A does not replicate exactly the thermostat hysteresis behavior. This model guarantees the maintenance of the ON state within the comfort dead-band (a in Fig. 1) but not the maintenance of the OFF state (b in Fig. 1). This freedom to switch ON to benefit the cost objective function is the reason why this model takes so much computation time. A model forcing the control variables to 0 or 1 according to the thermal balance equation and thermostat hysteresis, i.e. a rule-based system modelling the thermostat switching, is almost instantaneously solved. This is implemented in Model 1B. In addition to Model 1A, Model 1B (also for heating mode) explicitly forces s t = 0 when the indoor temperature (θ in t ): -is above the upper bound of the comfort band (θ in t > θ max ) OR -is within the comfort band (θ min ≤ θ in t ≤ θ max ) AND in the previous instant the AC was also OFF (s t −1 = 0) (b in Fig. 1 The binary variables s t keep the same meaning as in Model 1A, i.e. they control the thermostat switching. The auxiliary binary variables y t and z t establish the consistency conditions associated with thermostat switching to the ON (s t = 1) and OFF (s t = 0) positions: Constraints (21) impose that s t = 1 if θ in t < θ min . Constraints (22), (23) together with (25), (26) impose that s t keeps the value of s t −1 within the comfort band. Constraints (24) impose that s t = 0 if θ in t > θ max . However, this model is just the mathematical formulation of logical implications, which establish the values that the control variable should have according to the indoor temperature θ in t and thus no optimization is really at stake. That is, s t variables are rigidly determined by the conditions that establish θ in t and the thermostat operation (Fig. 1). In this case, there is a single feasible solution (the one that complies with the thermostat hysteresis switching rules), so obtaining the solution to this model is almost instantaneous.

Modelling Minimum ON/OFF Duration Periods
A new model (Model 2) has been developed to offer the possibility of controlling the AC to make the most of time-differentiated tariffs, i.e. being ON during lower electricity price periods when it is not strictly necessary to heat the building space so that thermal comfort is achieved even limiting the time ON when electricity prices are higher. In face of this freedom given to the model, it is necessary to inhibit excessive switching, which is accomplished by imposing a new set of constraints. In the comfort band, the AC control variables s t are determined to minimize costs in face of dynamic tariffs, but once there is a switch from 1 to 0 or from 0 to 1, then the new state ON/OFF should be kept for at least d time intervals (these minimum ON/OFF periods could be different for each state). These new constraints establish: This model is also quite demanding from the computation time point of view due to offering the possibility of switching from 1 to 0 or from 0 to 1 in the comfort band. Model 2 (heating mode): s.t.:

Modelling an Inverter AC
Inverter technology AC appliances have growing acceptance because of increased efficiency, extended life and elimination of abrupt load and temperature variations. In this type of appliances, a variable-frequency drive adjusts the compressor and the cooling/heating output. The previous models can be extended to accommodate inverter units, for instance considering they can be operated at different levels with respect to the nominal power. This can be modelled by introducing additional binary decision variables: For instance, for R = 5 with the power levels 20-40-60-80-100% of the nominal power, the AC operation is modelled in Model 3. In this model, θ min and θ max should be interpreted as absolute bounds the user is willing to endure. Note that the sets of constraints in the previous models could also have been used here (for simplification purposes they were omitted). p AC t is the power required to operate the AC at time t of the planning period, i.e. the AC is either OFF (p AC t = 0) or supplied at one of the specified power levels.

Model 3 (heating mode):
be combined with the thermostat modelling of Models 1 and 2 and/or with the cost vs. comfort trade-off.
In the next section, illustrative results of the exact solution (whenever possible with a given computation budget) of these models using a MIP solver are presented drawing attention to the computational effort.

Data
A planning period of 24 h with a time discretization of 1 min is adopted, i.e. 1440 time steps, which imposes a heavier computational burden. The thermostat deadband is defined by θ min = 20 • C and θ max = 24 • C. The initial indoor temperature θ in 0 = 18 • C in all models, except in model 3 in which it is θ in 0 = 20 • C. The initial thermostat status s 0 = 0. The external temperatures θ ext t ( • C) for a winter day (from a meteorological station in Coimbra, central Portugal, on January 1 st , 2012) are available at http://www.uc.pt/en/org/inescc/publications/files/temp_out_ winter.xlsx.
The nominal power of the AC, P AC nom = 1.5 kW. The time-differentiated energy prices c t (e/kWh), considering 10 sub-periods, are given in Table 2. Experiments have been made on an Intel(R) Core(TM) i5, 3.33 GHz computer with GAMS ver. 24.0.2 and Cplex ver. 12.5.0.0.

Results
A computational budget of 7200 s or relative MIP gap of 0.5% was established. The behavior of thermostat switching as well as the evolution of indoor temperature are displayed in Figs. 2 and 3 for Model 1A and Model 1B, respectively. The cost of the corresponding optimal solutions for Model 1A and Model 1B, as well as the dimension of the models and the Cplex time to obtain these solutions are given in Table 3 Table 4. Considering that the AC can operate at power levels 20-40-60-80-100% of P AC nom in Model 3, the cost of the optimal solution as well as the dimension of the model and the Cplex time to obtain these solutions are given in Table 5. The corresponding solution of the AC operation and indoor temperature variation is displayed in Fig. 6.  The frequency of switching in Model 1A is higher than in Model 1B due to the freedom given by that model to minimize cost. In Model 1A, the AC is ON for a longer time in pricing sub-period P 3 to account for the increase in price in subsequent sub-periods (Fig. 2). Note that the solution presented in Fig. 2 and Table 3 for Model 1A still presents a MIP gap of 43.47% after 2 h of computation, due to the combinatorial difficulties of this model. Model 1B, as expected, is solved to optimality instantly.
The temperature peaks in Figs. 4, 5 and 6 are due to the capability of the models to accommodate for comfort by making the most of lower price sub-periods. The MIP gaps in Model 2 (Table 4) are still significant after 2 h of computation.
The optimal solution to Model 3 is better than the solutions obtained for the other models due to the possibility of the inverter AC to work at different power levels.

Conclusion
This paper presented a set of mixed integer linear programming models in order to optimize the operation of a thermostatically-controlled air conditioning system (AC), aiming at minimizing the energy costs. These models capture different physical control characteristics: -Model 1A models the hysteresis operation of the thermostat (in heating mode), preventing from excessive switching when the indoor temperature is around a setpoint. This model gives, however, some freedom to thermostat switching when the indoor temperature is within the dead-band, enabling to change from OFF to ON if this benefits the cost objective function. -Model 1B removes the flexibility given by the previous model, forcing the AC to keep its state when the indoor temperature is within the dead-band. This model implements, in a MILP formulation, a rule-based system of logical implications.
The feasible region is reduced to one solution and the model can be solved very quickly, unlike Model 1A which is computationally very demanding. -Model 2 gives some flexibility while inhibiting excessive switching, by considering that, when the indoor temperature is within the dead-band, the thermostat must keep its state (ON or OFF) for at least a predefined number of time intervals. This model is also quite computationally demanding. -Model 3 introduces a new feature that can be used to extend the previous models: it models inverter units, considering that the AC can operate at different levels of nominal power.
Modelling cost reduction at the expense of accepting worsening the indoor temperature by some degrees was also introduced in this paper. It has been shown that the combinatorial nature of these models imposes a severe computation burden, which in some cases impairs obtaining the optimal solutions in a practical, acceptable computational time by a state-of-the-art solver. Therefore, this work lays the foundations for understanding those computational difficulties and gives insights for the development of other approaches, namely dedicated meta-heuristics customized for the features of the models, in order to offer good solutions with a suitable computational effort.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made. The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.