Demand Response Analogues for Residential Loads in Natural Gas Networks

Demand response for electrical power networks is a mature field that has yielded numerous efficiency and resilience benefits like managing the peaks and valleys of electricity usage, reduction in peak electricity usage, to name a few. However, only recently has the study of the counterpart to demand response for electric power in natural gas started to receive similar levels of attention. Natural gas systems are increasingly operating at or near capacity, which challenges these systems to meet all the needs for gas, especially during severe winter weather. However, unlike demand response programs in electrical networks, demand response in gas networks cannot shift peak usage. Here, we develop analogues to demand response that can help improve the resilience of natural gas systems by reducing peak consumption and thereby limiting potential disruptions such events can cause. This paper develops a mathematical formulation to support a residential-level demand response analogue for natural gas based on current and anticipated smart thermostat technologies. The mathematical formulation takes the form of an optimal control problem (OCP) that leverages physics constraints to model temperature changes, balance equitable service, and optimize the gas consumption for collections of houses. On test problems, the formulation demonstrates significant benefits, including the ability to cut peak demand by 15% while still ensuring equitable service to customers.


I. INTRODUCTION
One of the major challenges faced by modern energy systems is the problem of managing peak resource usage. During peak usage periods, energy systems operate at or near their capacity, which limits flexibility, increases green house gas emissions [1], [2], reduces reliability [3], decreases resilience, and, ultimately, drives costly infrastructure investments. To address this challenge, energy system operators are increasingly adopting and developing control solutions to help them better manage peak usage. One of the solutions adopted by the electric power industry is demand response (DR), which is used to shift energy usage away from peak consumption time periods. However, DR for other energy industries, such as natural gas, is relatively nascent. DR in the electric power industry generally relies on two technologies: load control devices and market strategies. Market strategies are used to price energy and provide incentives to consumers to use energy at times desirable to the system operator [4]. Load control devices (energy storage, power-heating switch devices, etc.), are used to directly flatten energy demand by moving consumption from times with high utilization to times with low utilization. One example of load control in electric power DR programs is home cooling systems, where thermal inertia is used to pre-cool buildings prior to peak demand hours.
Within the electric power industry, there are numerous studies that have demonstrated the benefits of DR programs. For example, in 2009, the United States Federal Energy Reg-ulatory Commission (FERC) announced a national DR action plan for maximizing the adoption of DR. The plan included a national communication program, the development of analytical tools, the outlining of regulations, and examples of model contracts [5], [6]. Subsequent reports estimated that adoption of DR has saved approximately 2,133 MW in wholesale markets (2018) and approximately 30,000 MW in retail markets (2017). While natural gas is an important energy resource [7], it has received little DR attention, as compared to electric power, and much of the existing work has focused on how electric DR can support natural gas by, for example, indirectly reducing natural gas used to produce electricity. Unfortunately, winter months see a considerable increase in non-electric natural gas usage (e.g., home heating) for which electric power DR offers no relief. Moreover, extreme weather events, like the 2021 cold weather snap in Texas, stress the natural gas infrastructure in ways that an electric power DR implementation cannot support. Hence, there is an important opportunity for gas specific DR programs, e.g., a DR program focused on residential home heating.
A DR program focused on residential home heating differs from similar programs for electric power DR for several reasons. In an electric power DR program, the operator introduces variable pricing to encourage greater use during offpeak hours; the consumer weighs the options and determines if it is viable to shift usage to periods with lower pricing, and usually, energy usage is shifted away from peak time periods. However, when considering a natural gas DR program for residential home heating, it is apparent that the aims and ideas of an electric power DR do not directly apply. For instance, variable pricing of gas may exasperate social inequities for a lifeline use of natural gas when it is most needed (e.g., cold temperatures).
This observation reveals the question at the heart of this paper-if demand response, as usually understood, does not translate to gas DR, is there an analogue that can be introduced? 1 While an operator cannot shift peak demand, it is still useful to reduce total and peak consumption, either at the level of an individual household, or for an ensemble of houses, especially in extreme winter weather conditions such as a polar vortex. For consumers, it requires agreement to allow control of heating equipment to be guided by the forecast of the ambient temperature and current internal temperature, along with a willingness to give up some amount of thermal comfort (in the form of deviation from set-point temperatures) for monetary incentives. A reduction in consumption could be the difference between total collapse of the distribution network (as happened in Texas recently [8]) and keeping that network functioning through challenging conditions with an equitable distribution of gas supply.
There are many challenges associated with natural gas DR [9] and this paper focuses on addressing the challenge associated with the formulation of mathematical models based on 1 In the subsequent discussion, the terms natural gas DR or gas DR are used with the tacit understanding that the terms do not correspond in a natural way to their counterpart in electric power DR. physics and thermodynamics for controlling natural gas consumption. To the best of our knowledge, this is the first effort to formulate and demonstrate the workings of such a model to support natural gas DR programs. The model assumes current and anticipated capabilities of smart thermostats and is focused on smaller consumers, who are a significant source of DR capability in electric power systems, but are yet to be studied in natural gas systems. Since home heating systems are the major source of natural gas consumption in residential households, the model developed in this paper optimizes the gas consumption (home heating) of a collection of houses to reduce peak gas demand during normal winter operation, as well as during extreme temperature events, such as those induced by severe weather events.
While a DR program has various aspects (incentive programs, objectives, thermostat controller options, etc.), all DR strategies are capable of being formulated as an optimal control problem (OCP) from the perspective of mathematical modeling. The OCP we develop is a mixed-integer optimization problem coupled with ordinary differential equations (ODEs) that model the dynamics of house temperature (heat gain from gas-fired furnaces and heat loss to the surroundings due to cooler external temperatures). The details of the OCPs are driven by two factors: (i) the implementation of the control policy on smart thermostats (centralized or decentralized) and (ii) the objective of the OCP (e.g., maintaining household temperatures as close to desired set points as possible while keeping peak consumption below a specified threshold, or reducing the overall gas consumption while keeping temperature deviations to a pre-specified limit). Overall, the goal of all the OCPs presented in this article is to assess (qualitatively and quantitatively) the merits of a DR policy through the viewpoint of both the natural gas system operator and residential customers. The models also serve as an important first step towards identifying requirements for DR technologies and market structures. We expect the methods developed in this paper to be used by stakeholders as a baseline for identifying potential benefits and determining the underlying requirements for a successful implementation of DR in a natural gas setting. In summary, while the crux of this article is the mathematical formulation of OCPs that can model gas DR for residential use, the key contributions of this paper include: • Formulation of an Ordinary Differential Equation (ODE) model for a home heating system inspired by dynamic models for Thermostatically Controlled Loads (TCLs) and leverage these to formulate OCPs that optimize natural gas load scheduling at a household level. • Illustration of the use of a receding-horizon based algorithm to solve the resulting OCPs when the discretized problem becomes large. • Presentation of case studies which highlight the aggregation benefits and the degree of thermostat flexibility required to achieve a certain level of peak reduction.
The remainder of this paper is organized as follows: Sec-tion II provides a literature review and relevant background for DR. Section III discusses a high level description of the problem. In Section IV, we present a baseline model to reflect current practice and formulate the physics-based, decentralized and centralized natural gas DR models and describe the solution method in Section V. Finally, Section VI presents the case studies and simulation results and Section VII summarizes the conclusions and offers future perspectives.

II. BACKGROUND AND LITERATURE REVIEW
Optimization is the core computational technology used for investigating and controlling DR. Within the literature, both deterministic and stochastic optimization methods have been developed. Deterministic models assume that all information is known by the users and the operators [10], [11]. Stochastic models introduce uncertainty associated with supply, demand, price, customer behavior, environmental variables, etc.
[12]- [19]. The objective functions consider economic factors like cost savings [20], [21], operation profits [18], [22]- [24], engineering, reliability [25], risk minimization [18], greenhouse gas emission reductions [19], [26], social welfare [17], [22] and user comfort [27]. Other lines of work include optimal control of ensembles of Thermostatically Controlled Loads (TCLs) [28]- [30] which are very closely related to the models and DR formulations that are presented in this paper. Ensemble-based optimal control problems reduce electricity consumption from TCLs (typically air-conditioners) during the summer months by utilizing thermal inertia. In this line of work [28], a simple Newtonian thermal model (an RC circuit model) for heat exchange between the surroundings and a house is utilized. We use similar thermal models in this article to develop a DR program in a natural-gas setting to cater to the heating needs of households in winter months. There are two fundamental differences between conventional (electric power based) TCL and the natural gas furnace based control considered here. First, TCLs operate in steady-state, while gas furnace operation is inherently transient. Second, TCLs have almost instantaneous response, while the gas furnaces operate on much longer time-scales, thus ruling out methods based on thermal-inertia that were successful for TCLs. While DR response has been a part of the electric power industry's technology adoption conversations for at least a decade, it is only recently that such technologies have been discussed in the context of other energy industries, such as natural gas. This interest is largely driven by two recent emerging drivers. First, historically low prices have greatly expanded natural gas's utilization and placed enormous strains on existing capacity [31]. At the same time, economic, political, and social factors have limited the ability of the industry to add capacity [32]. As a result, DR is being proposed as a technology for improving the utilization of existing capacity [33]. Second, extreme events like polar vortices and the recent 2021 cold snap in Texas have created sudden spikes in the natural gas consumption for home and commercial heating that contributed to loss of load to millions of homes [8]. Such events can create supply and capacity shortages that interrupt gas deliveries. During recent events, states like Minnesota have employed an extreme form of DR that asked customers to decrease thermostat set-point temperatures to 60 F to counter such demand spikes [34]. Despite the potential for DR to support the natural gas industry in both these contexts, work and approaches for modeling and controlling DR in natural gas systems has been limited (unlike electric power) but is expected to increase dramatically in the next several years [9].
There is also some research on natural gas as a supplemental energy resource for electric power DR. [23] developed a DR market model where natural gas based generation provides a supportive role for electric power DR. [21] developed a model of the combined electric and natural gas system in Great Britain and used it to study the influence of DR in combined system modeling. [19] used electric power DR technologies, like pumped storage and incentives, to improve the operational performance of coupled natural gas and electric systems. More recently, [35] developed a pricing strategy to maximize the social benefits of DR in coupled gas-grid systems. Finally, [36] described a price-based DR for optimally scheduling deferable gas loads (power plants) to maximize profits and [37] built on this approach to develop a combined pricing mechanism to schedule production and consumption of natural gas systems to maximize social welfare. In all these papers, the focus was largely on how DR in power systems could impact the operations of joint gaselectric systems with limited attention devoted to how to use gas assets in DR. Figure 1 shows a conceptual diagram of consumer level DR programs. While this is not the only way DR programs are structured, it is a common method and this structure is used to inform the OCP and solution methods developed in this paper. In this figure, the operator of the DR program seeks to reduce energy consumption during peak periods. Consumers (households) are given incentives to participate in the DR program to help meet the goals of the system's operator. Load aggregators (blue boxes in Fig. 1) are sometimes used to recruit participants, serve as an interface between the individual consumers and markets, and are used to coordinate consumer responses to DR requests.

III. DEMAND RESPONSE SETUP AND ASSUMPTIONS
A key component of the gas DR of this paper is ensuring relatively equal distribution of adverse conditions to participants. In this context, an adverse condition is thermal discomfort (modeled as average deviation from temperature set points or the maximum deviation from temperature set points) which has detrimental effects on human health and performance [38]- [40]. In the OCP, the requirement of equitable distribution for customer satisfaction is enforced by including a social welfare objective function which indirectly minimizes the thermal discomfort. Another component of DR is participant incentives as participants need incentives Here, an aggregator may refer to a utility company or other third party in a particular geographic location.
(typically monetary) to give up temperature comfort. This paper uses the OCP to identify key parameters (maximum temperature deviation and number of participants) to inform the choice of incentive levels but does not directly address the challenge of implementing incentives. The final component of DR is constraints that model the governing equations describing the dynamics and control of the system. Whereas electric power based DR is well developed, natural gas DR is relatively new and these governing equations have yet to be fully established. In this article, we utilize thermodynamic models of heat flow developed for electric power TCLs [28], [30] and contribute a natural gas DR formulation for these dynamic models.
To construct the OCP, we make the following assumptions about DR. These assumptions allow us to simplify the mathematical model, and investigate the worst-case guarantee that can be provided for the benefits a natural gas DR program can support. The assumptions include: 1) The natural gas utility operator aggregates the natural gas demand from all households in a locality. From here on, we use the term natural gas utility operator and aggregator interchangeably.
2) The only source of natural gas consumption for the households is home heating. This is justified because home heating systems are the major source of flexible natural gas consumption in most households during the winter months. 3) Each participating household in the DR program has a smart thermostat. In this article, we consider two common thermostats ("Type 1" and "Type 2") and their use in an OCP. The "Type 1" thermostat is capable of (a) obtaining a temperature forecast for the future, (b) allowing the customer to establish a temperature set-point for the thermostat, (c) allowing the customer to choose a parameter that adjusts the priority of the thermostat between maintaining a set-point temperature to saving on gas consumption, and (d) has minimal computing capability to solve optimization problems. The use of a "Type 1" thermostat permits a decentralized control policy as no information is shared between households or between the household and the aggregator.
The "Type 2" thermostat is capable of (a) sending the current temperature of the house and the set-point temperature of the thermostat to the aggregator and (b) receiving on-off control policy from the aggregator that it implements. These controllers support centralized control policies. Unlike the "Type 1" thermostat, the "Type 2" thermostat does not perform any computation and all of the computation is performed by the aggregator. In both cases, the underlying OCP is largely the same and relies on the same algorithm developed in this paper. 4) Effects of the distribution pipeline network used to deliver gas to the households are saved for future work. For this viability study, we argue that if one does not quantitatively or qualitatively find any benefit to using the OCP model for DR under these assumptions, then any DR formulation with less strict assumptions will only do worse and be unviable. Furthermore, the effectiveness of utilizing the proposed DR formulations are quantitatively corroborated using a combination of (i) the comfort for each participant examined in terms of statistics of temperature deviation from the set-point temperature and (ii) the overall reduction in the amount of gas consumed by the participant. On the other hand, benefits to the aggregator are quantified using the overall reduction in the gas consumption achieved when more participants sign-up for the DR programs.

IV. DR OPTIMAL CONTROL PROBLEMS
Before presenting Optimal Control Problem (OCP) formulations for the natural gas DR problem, we introduce some notation. A brief outline of the notation is provided in the nomenclature section below while a detailed description of the entries is provided in context as and when they are used to formulate the OCPs.

Sets and parameters:
H -set of houses indexed by h Θ(t) -forecast of ambient temperature α h , β h -thermal coefficients of house h Q h -burn rate of the gas furnace in the house h θ 0 h -initial temperature of house h θ h -set-point temperature of house h γ -load curtailment factor T -optimization time horizon λ h -thermostat parameter that controls its mode of operation D -max. total gas consumed by all the households in H V h -volume of house h κ h -heat transfer coefficient of the walls of the house h A h -surface area of the walls of the house exposed to the outside temperature E g -specific heating value of natural gas C a -isochoric specific heat of air ρ a -density of air State variables: θ h (t) -temperature of house h Control variables: q h (t) -on-off variable for the natural gas supply of house h Other variables: ∆ h -maximum deviation of the house temperature from the set-point temperature of the house G h -total gas consumption for the house from [0, T ] We next present the thermodynamics of heating inside a single house. For ease of readability, we present the explanation of the notation in detail.

B. THERMODYNAMICS OF HEATING
In a locality where a DR program is implemented, we denote by H the set of houses, and index each house by h. We let Θ(t) be the given time-varying forecast of the ambient/outside temperature. We let θ 0 h denote the initial temperature inside house h. Each house, h, has a gas furnace that burns fuel at the rate Q h (kg s −1 ). For each house h ∈ H, we let {α h , β h } denote its thermal coefficients. Also associated with each house h are the following state and control variables: θ h (t) -current temperature of the house (K) and q h (t) ∈ {0, 1} -a binary control variable that controls supply of natural gas. For each house, h ∈ H, we denote the differential algebraic equations (DAEs) that govern the thermodynamics of heating by the following equations: We also remark that Eq. (1a) has one state variable and one control variable and is similar to an RC electric circuit. The stability analysis of such RC equations is well studied in the literature and it is known that these systems are stable and controllable [41] and have good robustness properties [42]. Utilizing the above dynamic equations, we now present three OCP models that are reflective of current practice and thermostats available. The first OCP provides a baseline that models how thermostats are operated today. The second provides the OCP for Type 1 thermostats (inherently decentralized) and the third provides the OCP for Type 2 thermostats (inherently centralized). Structurally, the OCPs are nearly the same, yielding a common algorithm discussed later.

C. BASELINE OCP MODEL
In the baseline model, each house is equipped with a normal thermostat that implements an instantaneous control policy. At any time instant t, the thermostat decides to turn-on or turn-off the gas supply if the current temperature of the house is less or greater than (respectively) the set-point temperature of the thermostat. This is reflective of the thermostats that are currently deployed in the majority of households in the United States. For this baseline model, the control policy is defined by the following equation In Eq. (2), H(·) is the Heaviside step function andθ h is the set-point temperature of the thermostat. The model is a decentralized model because implementing the control policy in each household h ∈ H requires only local information from h. The baseline model is integrated over a finite time horizon [0, T ] to obtain a solution θ * h (t) and q * h (t) for every h ∈ H. The mass flow rate of natural gas consumed by all houses in the set H is computed from the solution as: In Eq. (3), d b models the total mass flow rate of gas consumed (in kg s −1 ) by all houses in the set H. We refer to this value as the baseline gas consumption for a set of households in H.

D. THERMOSTAT TYPE 1: DECENTRALIZED LOOK-AHEAD OCP
The decentralized look-ahead model uses "Type 1" smart thermostats. The output of this model is a decentralized onoff control policy for each house h. The control policy for each household is dependent only on current temperature and the parameters of the thermostat in that house. This model is a look-ahead model because, unlike the baseline model, this model accounts for temperature forecasts and computes a control policy accordingly. We next introduce additional notations to formulate this decentralized look-ahead OCP.

VOLUME 4, 2016
The Type 1 thermostat has the capability to obtain a temperature forecast, Θ(t), up to a set time T in the future (referred to as the optimization time horizon). The first thermostat parameter is the desired temperature setpoint,θ h , that the user wants to maintain inside the house h ∈ H. The second parameter, λ h ∈ [0, 1], represents the mode in which the thermostat operates. When the value of λ h = 1, the thermostat keeps the maximum deviation of the house temperature from its set point to a minimum for the entire optimization time horizon. In the other words, the DR OCP computes a control policy that minimizes the maximum deviation over the entire time horizon, defined by When λ h = 1, the thermostat computes a control policy that minimizes the total natural gas consumption over the entire optimization time horizon (see Eq. (4c)). For all other values of λ h ∈ (0, 1), the objective is to minimize a convex combination of ∆ h and G h with λ h as the multiplier and models a trade-off between these outcomes. However, this trade-off cannot be quantified a priori, in the sense that it is not possible to choose a value of λ h that will guarantee a certain outcome for either ∆ h or G h unless the pareto front is solved for. The OCP of the decentralized look-ahead model that the smart thermostat solves for each house h ∈ H is formulated as follows: Here, the optimal control policy for each house is obtained by solving the OCP in (4) for each house. In the above formulation, Eq. (4b) represents the dynamics of heating for the house h, Eq. (4c) computes the total gas consumption for house h over the entire optimization time-horizon, and Eq. (4d) provides the initial conditions for the dynamics in Eq. (4b). Finally, Eq. (4e) is the binary restrictions on the control variable q h (t) and Eq. (4f) provides the bounds on the state variables, θ h (t).
The decentralized control policy that is implemented by the thermostat allows the customer to save on natural gas consumption (incentive) without compromising on the temperature deviation from the set-point temperature. From an aggregator stand-point, it provides a reduction in peak (DR) without compromising on the comfort level of the participants. In that sense, this decentralized look-ahead model is a win-win for both the participants of the DR program and the aggregator. This formulation can be thought of as indirect DR, since conventional DR usually responds to prices or overall demand. Here we have an indirect price response, in the form of reduced overall consumption. The decentralized formulation explores the limits of what can be done in the context of DR using today's smart thermostats with access to temperature and forecasts.

E. THERMOSTAT TYPE 2: CENTRALIZED LOOK-AHEAD MODELS
We next discuss the OCP for "Type 2" thermostats. Unlike a Type 1 thermostat, the Type 2 thermostat has no computational capability. The thermostat sends the temperature set-point and the current temperature of the house to the aggregator and receives and implements an on-off control policy that is computed by the aggregator. The aggregator has a forecast of the ambient temperature for the duration of the optimization time horizon and uses it to compute the control policy based on look-ahead information. The model is centralized because the aggregator utilizes (i) set-point temperatures of each household, (ii) initial temperatures of each household, (iii) parameters of each house, (iv) ambient temperature forecasts, and (iv) additional constraints on consumption to solve a single OCP that computes a control policy for each household. The computed control policy is transmitted to the thermostat in each household for implementation.
Similar to the decentralized look-ahead models, there is an optimization time horizon is [0, T ]. We let, ∆ h , denote the maximum deviation between the set-point temperature and the actual temperature of house h over the entire optimization time horizon. To model peak shaving, we let D denote the maximum total amount of gas consumed by all the households in H and let γ ∈ (0, 1] denote the percentage of gas reduction that the aggregator achieves with the centralized DR program. With these notations, we present the following look-ahead model with two different objective functions (i) mean deviation objective and (ii) max deviation objective: In both the objective functions, the value of ∆ h is interpreted as discomfort of participants. The objective function in Eq. (5a) and Eq. (5b) minimizes the average and maximum discomfort, respectively, over all the households participating in the DR program. The centralized look-ahead OCP model with objective functions in Eq. (5) is subject to the following constraints: The constraints in Eq. (6e) enable the operator to enforce a peak mass flow rate of γ · D. Here, D is determined by the system operator by estimating of how much total gas is consumed by all the households during winter peaks. γ is a percentage factor that the operator can vary to perform reliability analysis in the case of extreme events. We refer to γ as the curtailment factor that indicates the percentage by which the operator wishes to curtail the maximum mass flow rate in the natural gas system. The centralized model offers an advantage over the decentralized one in that once γ is chosen, the operator already knows what to expect regarding peak consumption, and if the deviation performance corresponding to a chosen value of γ is unsatisfactory, then it is clear that a lower value of γ needs to be considered to achieve better deviation performance.

V. ALGORITHMS
This section presents a receding-horizon algorithm to solve the OCPs described in section IV. To solve the OCPs, we first discretize the state and control variables over time. It is useful to separate the temporal resolution of the control variables from the state variables to better model the capabilities of the controller and limit chattering [43]. Before we present the receding-horizon algorithm, we first present a discretized version of the decentralized and centralized look-ahead models.
In the algorithm, the optimization time horizon [0, T ] is discretized into n equispaced time instants K = {0 = t 0 , t 1 = ∆t, . . . , t n = n∆t = T } with a time step of ∆t = T /n, for every house h ∈ H. For OCP's based on Type 1 thermostats, the decentralized look-ahead model is a Mixed-Integer Linear Program (MILP) of the form: The MILPs for the Type 2 OCP are solved as a single optimization problem for all houses by the aggregator. The constraints of these MILPs after discretization, are as follows:

A. A RECEDING-HORIZON (RH) ALGORITHM
One issue that hinders solving the OCP MILPs to global optimality is the size of the MILPs. For the Type 1 thermostat there is one MILP for each house, h ∈ H, that is solved separately. The size of the MILP (Eqs. (7) and (8)

VI. CASE STUDIES
In this section, we present two case studies. The first case study analyzes the efficacy of the models for a typical winter day in the Chicago area and the second case study analyzes the impact of the 2019 Polar vortex in the Chicago area from temperature data. For both the case studies, the set of houses and the parameters for each house remain unchanged. The Julia Programming language [44] was used for all the implementations and Gurobi [45] was used to solve the MILPs that were obtained by discretizing the OCPs. All the computational experiments were run on an Intel Broadwell E5-2695 processor with a base clock rate of 2.10 GHz and a RAM of 128 GB.

A. HOUSE DATA
For all the case studies, the total number of houses in the set H is set to 140. 14 categories of houses were chosen ranging from single bedroom apartments to single family homes. The thermal coefficients for each house was calculated according to the formulae derived in Appendix A. The VOLUME 4, 2016 Solve MILP for [0, FIGURE 2: Schematic of the RH algorithm. IC denotes initial conditions. furnace burn rate, Q h , for each house was set in the range of [3×10 −5 , 12×10 −5 ] kg s −1 depending on the square footage of the corresponding house. The value of the parameters specific to each house category (volume, square footage, initial temperature, set-point temperature, etc.) are made available at https://github.com/kaarthiksundar/ng-dr-data. The other parameter values used for computing the thermal coefficients for each house are shown in Table 1.

Parameter
Value (SI Units)

B. AMBIENT TEMPERATURE DATA
For a typical winter day case study, historical ambient temperature data for the Chicago area on January 2, 2020 was chosen. For the polar vortex case, the temperature data from January 26, 2019 noon to January 27, 2019 noon was chosen. The temperature data for both cases were obtained from https://www.noaa.gov/weather. The plot of the temperature profiles for both the cases are shown in Figures 3a and 3b, respectively. Throughout the rest of the article, we refer to the two case studies as "typical-day" and "polar-vortex".

C. RESULTS FROM THE TYPE 1 THERMOSTAT MODEL
Here we present the results obtained for the the Type 1 thermostat model in Eq. (4) for both case studies and compare it with the results from the baseline model in Eq.
(1) and Eq. (2). For the decentralized look-ahead model, the value of λ h for all the houses is varied in the set {0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1} and the look-ahead time period, denoted by T rh (in hours) is chosen from the set {1, 3}. The look-ahead OCP for each house is solved separately but in parallel. The state and control discretization for the MILP reformulations for the look-ahead OCPs were set to 1 and 3 minutes, respectively. Table 2 shows the average computation time, in seconds, to solve the MILPs to optimality. The average is computed over all the houses since each house has its own MILP. As T rh increases, it is evident that the size of the MILP increases and run time increases.    Tables 3 and 4 show the statistics of temperature deviation and the total gas consumption obtained using the Type-1 thermostat for T rh ∈ {1, 3} hours. The deviation performance is a measure of discomfort caused to the households, i.e., the greater the deviation, the greater the discomfort. The level of discomfort for λ h = 0.65 provides ample evidence that experiments for λ h < 0.65 are unnecessary. In all the tables, the first row is computed by solving the baseline model. The results indicate that a lower value of the lookahead time period results in more short-sightedness of the solutions. Reducing the value of T rh degrades the deviation performance needed to achieve lower gas consumption. This trend is observed in both the Tables 3 and 4. For the "typical-day" case study results in Tables 3a and 3b, we observe substantial reduction in gas consumption    Table 3a, when the value of λ h = 0.75, the gas savings per house is 0.07 kg. When the number of houses aggregated is around 10000, this leads to a total gas consumption reduction of ≈ 700 kg, which is a substantial amount of reduction in consumption even on a typical winter day. Furthermore, it is important to note that this reduction in gas consumption does not come at the cost of compromising on the deviation performance (in comparison to the baseline model). In summary, utilizing smart thermostats with lookahead features reduces a substantial amount of consumption without compromising deviation performance for a typical winter day.
In the "polar-vortex" case study, we observe that though reduction in gas consumption is not as substantial as in the "typical-day" case study, there exist regimes (λ h values) for which the decentralized look-ahead models provide some reduction in gas consumption along with good deviation performance. Also, unlike the "typical-day" case, where consumption reduction is achieved for all λ h values that provide deviation performance comparable to the baseline model, the values of λ h closer to 1 cause an increase in the gas consumption. This is again expected due to the severity of the polar vortex. Nevertheless, if the operator can ensure all the smart thermostats for the decentralized model have appropriate values of λ h , it can still lead to some reduction in gas consumption.
Finally, the Tables 5a and 5b present the amount of reduction in gas consumption that can be achieved for both the "polar-vortex" and the "typical-day" case studies, respectively, for certain regimes (λ h values) when compared against the baseline model which is reflective of current practice. The value of T rh was set to 3 hours for the results in both the tables. The results show that even in an extreme event scenario like the polar vortex, there exists certain regimes of operation of the thermostats that provide an overall reduction in natural gas consumption while providing a better temperature deviation performance than what can be achieved by current practice (baseline model). In a typical winter day, this reduction is amplified even further with a wider range of λ h values providing the same type of benefit to the customers. From an aggregator stand-point, installation of smart Type 1 thermostats on a greater number of houses only increases the aggregation benefits as shown in the final column of the Tables in 5.

D. RESULTS FOR THE TYPE-2 THERMOSTAT MODELS
Given the Type 1 thermostat model results, we focus on Type 2 thermostat results for the polar-vortex case study. For these results, the following parameters are used: (i) D (peak mass flow rate) is set to the maximum mass flow rate of the baseline model and (ii) γ (curtailment factor) is chosen from the set {0.85, 0.90, 0.95}. Similar to the previous results, all the experiments for the model are run for T rh ∈ {1, 3} hours. In all the results, "min." is used to denote the Type-2 thermostat model with the mean deviation objective and "min. max." is used to denote the max. deviation objective. Finally, a time limit of 1.5 hours was placed on every computational experiment. Table 6 shows the average computation time and the average optimality gap for all the runs of the Type 2 thermostat model with the two different objectives. In Table 6, the averages are computed over each run of the centralized MILP in Sec. V-A. One immediate observation is the high computation time and optimality gap values for the "min." model when compared against the "min. max." model. This trend can be explained by degeneracy in the solution space of the OCPs, i.e., existence of multiple solutions with the same mean deviation when minimizing mean deviation. This degeneracy is removed when minimizing the maximum deviation which reduces the computation time. Also, as T rh increases, the number of variables and constraints increases which leads to increased computation times. Finally, as the value of the curtailment factor decreases, the control policy bounds natural gas consumption. In an extreme event setting, this can make finding optimal solutions more difficult. Physically, this is equivalent to constraining the capacity of the system. Consistent with results in Table 6, the box plots (Figure 4) also show that the "min." model produces a deviation performance much worse than the "min. max." model within the computational time limit imposed. The spread of the data in Figure 4 reveals two salient trends: (i) When using the "min." objective function, a 1 hour look-ahead leads to too big a spread in temperature deviation data. However, by using a 3 hour look-ahead, one can achieve peak reduction of 5%, 10% or even 15% (corresponding to curtailment factor values of 0.95, 0.90, and 0.85, respectively) with a deviation performance at least as good as the baseline model. (ii) When using the "min-max" objective function, a 1 hour look-   3: Statistics of temperature deviation and amount of gas consumed for the entire day for the "typical-day" case study. Note that reduction in gas consumption does not come at the cost of compromising on the deviation performance (in comparison to the baseline model). Thus, utilizing smart thermostats with look-ahead features can lead to substantial amount of reduction in consumption without compromising the deviation performance for a typical winter day.
ahead leads to unacceptably large temperature deviations even though the spread is limited. However, with a lookahead of 3 hours, even 15% peak reduction is feasible with temperature deviation values and spread that are much lesser than the baseline. This leads to a conclusion that if the technology and infrastructure is indeed set up for implementing the centralized models, it can provide the operator a more efficient way of providing equitable services to the customers without incurring a substantial increase in the natural gas consumption. Going further, the results from the centralized model allow one to conclude that equitable service can be provided while simultaneously reducing the peak natural gas consumption during an extreme event scenario.

VII. FUTURE RESEARCH DIRECTIONS AND CONCLUSION
This paper develops an OCP for natural gas DR based on physics and thermodynamics, and demonstrates significant benefits of DR in the field of natural gas. In order to demonstrate the feasibility of a DR implementation for natural gas focused on residential loads (home heating), the residential loads are treated as the sole consumers of natural gas. The natural gas pipeline network is ignored, and two kinds of smart thermostats are modelled. Our formulation of the baseline OCP is reflective of the current state of practice, and our formulation has been compared to the current state for two different instances -a typical winter day where the maximum and minimum temperature attained are 15 degrees apart, and an extreme weather event represented by a polar vortex where the temperature range is almost three times as large. It has been shown that our decentralized formulation can lead to appreciable savings in gas consumption for a typical day, and do slightly better than current practice for days representative of winter storms. However, our centralized formulation can save up to 15% gas consumption even on such days. For both the centralized and decentralized look-ahead OCP models of DR, we have shown in the case studies that substantial reduction in gas savings for end customers can be achieved without compromising on the comfort level, while at the same time reducing peak consumption for the aggregator. The only cost to the customers is the thermostat installation in the case of the decentralized OCP and the communication network to communicate temperature inside the house to a central operator in the case of a centralized OCP. Hence, if this cost is offset using incentives, it will attract customers for this DR program. Thus, the aggregator could use the formulations to determine how much benefit they would get, and then use this information to create an incentive that would get people to participate. Currently, household customers have fixed rate plans or variable rate plans (whose prices change on a monthly basis) so the only way to lower the cost is to reduce consumption. Hence, our DR formulations are reflective of today's technologies and market structures, and more sophisticated price based DR  We observe that though the reduction in gas consumption is not as substantial as in the "typical-day" case study, there exist regimes (λ h values) for which the decentralized look-ahead models provide some reduction in gas consumption along with good deviation performance.  will require the introduction of real-time gas rates.
To make the study more realistic in the future, current research can be extended in the following directions: (i) inclusion of network effects due to the pipeline network and   (b) Centralized look-ahead model with min. max. objective function FIGURE 4: Box plot of the simulated temperature deviation (with respect to the set point temperature) for all the houses, i.e., ∆ h values for all the houses in the set H. The spread of the data reveals two salient trends: (i) When using the "min." objective function, simulation with a 1 hour look-ahead leads to temperatures with too big a spread in temperature deviation. However, by using a 3 hour look-ahead, one can achieve peak reduction of 5%, 10% or even 15% (corresponding to curtailment factor values of 0.95, 0.90, 0.85 respectively) with a deviation performance at least as good as the baseline model. (ii) When using the "min-max" objective function, simulations with a 1 hour look-ahead lead to unacceptably large temperature deviations even though the spread is limited. However, with a look-ahead of 3 hours, even 15% peak reduction is feasible with temperature deviation values and spread that are much lesser than the baseline.
the study of its impact on the proposed DR models, (ii) formulation and solution of DR models that include uncertainty in the temperature forecast, (iii) data-driven estimation of the thermal coefficients of the house using measurements and the known uncertainty in the estimates of these coefficients, and finally (iv) quantify the link between tine-based rate/economic incentives and direct load control. In conclusion, our article represents the start of an interesting and important problem of designing DR for a natural-gas only system. This work is timely given the realities of climate change; climate change is causing an increased number of severe winter weather events and this in turn is causing a stress on the natural gas pipeline networks, impeding the ability of these systems to provide equitable services to all its customers. The lack of resiliency of the natural-gas grid results in inconvenience, disruption and even loss of lives. The proposed DR models are a simple technique to increase the resiliency of the natural gas systems, leading to more efficient operation in severe events as well as normal winter weather. Nevertheless, large scale implementation of these models requires major technological and infrastructure investment as well as policy changes that need to be initiated by the federal government. .

APPENDIX A COMPUTATION OF THERMAL COEFFICIENTS FOR A HOUSE
We consider a house h ∈ H and let θ 0 h denote the initial temperature inside house h. We let Θ(t) denote the timevarying forecast of the ambient (outside) temperature. Each house has a gas furnace that burns fuel at the rate Q h (kg s −1 ). Each house h is also defined by V h -the volume of the house (m 3 ), κ h -the heat transfer coefficient of the walls of the house based on insulation (W m −2 K −1 ), A h -the surface area of the walls of the house exposed to the outside temperature (m 2 ), andθ h -the constant temperature set-point of the house (set by the occupants) (K). Finally, we let E g (J kg −1 ) denote the specific heating value of natural gas. When the furnace is on, the heat flux produced by the furnace is given by Q h · E g (J s −1 ). Furthermore, the following statements are assumed about the physical system: the air inside the house undergoes an isochoric (constant volume) process, and hence the relevant specific heat of the gas C a (J kg −1 K −1 ) is calculated at constant volume. We assume that C a is constant over the temperature range of interest. We also neglect the thermal compressibility (density dependence on temperature) of air and assume a constant density ρ a (kg m −3 ) over the temperature range of interest. The only source of heat supply is the furnace, and heat is constantly lost through the walls of the house at a rate proportional to the surface area of the walls and the difference in temperature between the house and the external environment.
Given these notations and assumptions, the dynamics of heat flow within the house h is governed by: Comparing Eq. (9) and (1a), we can obtain expressions for the thermal coefficients α h and β h as: In all the computational experiments, the thermal coefficients for each house is computed using Eq. (10). KAARTHIK SUNDAR received the Ph.D. degree in mechanical engineering from Texas A&M University, College Station, TX, USA, in 2016. He is currently a Research Scientist in the Information Systems and Modeling Division at Los Alamos National Laboratory, Los Alamos, NM, USA. His research interests include problems pertaining to vehicle routing, path planning, and control for unmanned/autonomous systems; nonlinear optimal control, estimation, and large-scale optimization problems in power and gas networks; combinatorial optimization; and global optimization for mixed-integer nonlinear programs.

SHRIRAM SRINIVASAN earned his MS in
Mathematics and PhD in Mechanical Engineering, all from Texas A&M University, College Station. His research interests are in computational mechanics and reduced-order models of structured systems such as fracture networks and gas pipeline grids. He is a staff scientist in the Applied Mathematics and Plasma Physics Group of the Theoretical Division at Los Alamos National Laboratory, New Mexico. RUSSELL BENT received his Ph.D. degree in computer science from Brown University in 2005. Since then he has been a scientist at Los Alamos National Laboratory (LANL), NM, USA. He is currently in the Applied Mathematics and Plasma Physics Group (T-5), where he leads LANL's interorganizational Advanced Network Sciences Initiative (ANSI), https://lanl-ansi.github.io/. ANSI is an interdisciplinary initiative that enables fundamental and applied research to address long-term challenges in critical infrastructure design, operation, and security. The primary philosophy of ANSI is that combining insights from Theoretical Physics, Applied Mathematics, Computer Science, and Engineering can result in novel computational methods that address a variety of emerging challenges in infrastructure networks. Dr. Bent is the principal or co-principal investigator for DOE projects in critical infrastructure systems research and development that focus on improving robustness of infrastructure systems to extreme events, increasing resilience of distribution networks, modeling interdependencies between systems, managing disasters that impact critical infrastructure, modeling smart grid technologies, and developing methods for mixed-integer, non-linear optimization. Such projects include the Grid Modernization Laboratory Consortium (GMLC) Project Extreme Event Modeling, the GMLC and DOE Smart Grid Project LPNORM: A LANL, PNNL, and NRECA Optimal Resiliency Model, and the DOE Advanced Grid Modeling (AGM) Project Joint Power System and Natural Gas Pipeline Optimal Expansion Planning. He is also the lead developer for the software Alpine.jl, A Global Solver for Nonconvex MINLPs (https://github.com/ lanl-ansi/Alpine.jl) and the software GasModels.jl, a toolbox for modeling natural gas systems (https://github.com/lanl-ansi/GasModels.jl) He is the author of one book, Online Stochastic Combinatorial Optimization, and over 80 peer reviewed journal and conference publications. Dr. Bent is also an associate editor for the INFORMS Journal of Computing.