Building control system with heat disturbance estimation and prediction

An environmental control system for a building including heating, ventilation, or air conditioning (HVAC) equipment that operates to affect a zone of the building and a controller including a processing circuit. The processing circuit is configured to estimate a thermal resistance between air of the zone and of an external space using values of a temperature of the zone air, a temperature of the external space air, and a heat transfer rate of the HVAC equipment, each value corresponding to a different time step within a time period. The processing circuit is configured to use the thermal resistance, time step specific values of the temperatures, and time step specific values of the heat transfer rate to estimate corresponding values of a heat disturbance. The processing circuit is configured to operate the HVAC equipment using a model-based control technique based on the heat disturbance values.

BACKGROUND

The present disclosure relates generally to control systems for buildings. The present disclosure relates more particularly to system identification for controlling building equipment.

System identification refers to the determination of a model that describes a system. For example, system identification may be used to identify a system describing environmental conditions. Because the physical phenomena that govern such systems are often complex, nonlinear, and poorly understood, system identification requires the determination of model parameters based on measured and recorded data from the real system in order to generate an accurate predictive model. However, various issues such as rounding of sensor measurements to a nearest whole number, high-frequency inputs, etc. can affect accuracy of the predictive model if not properly accounted for.

SUMMARY

One implementation of the present disclosure is an environmental control system for a building, according to some embodiments. The system includes heating, ventilation, or air conditioning (HVAC) equipment that operates to affect a variable state or condition of a zone of the building, according to some embodiments. The system includes a controller including a processing circuit, according to some embodiments. The processing circuit is configured to estimate a thermal resistance between air of the zone and air of an external space outside the zone using values of a temperature of the air of the zone, a temperature of the air of the external space, and a heat transfer rate of the HVAC equipment that provides heating or cooling to the zone, according to some embodiments. Each of the values correspond to a different time step within a time period, according to some embodiments. The processing circuit is configured to, for multiple time steps within the time period, use the thermal resistance between the air of the zone and the air of the external space, a time step specific value of the temperature of the air of the zone, a time step specific value of the temperature of the air of the external space, and a time step specific value of the heat transfer rate of the HVAC equipment to estimate a corresponding time step specific value of a heat disturbance, according to some embodiments. The processing circuit is configured to operate the HVAC equipment over the time period using a model-based control technique based on the time step specific values of the heat disturbance, according to some embodiments.

In some embodiments, the processing circuit is configured to perform a system identification process based on the time step specific values of the heat disturbance to identify a predictive model. The predictive model is used in the model-based control technique, according to some embodiments.

In some embodiments, the system identification process includes identifying a scaling parameter that scales the time step specific values of the heat disturbance.

In some embodiments, the processing circuit is configured to identify a Kalman gain and a stochastic model of the heat disturbance based on the time step specific values of the heat disturbance or predicted heat disturbance values based on the time step specific values of the heat disturbance. The Kalman gain and the stochastic model are used in the model-based control technique, according to some embodiments.

In some embodiments, identifying the Kalman gain and the stochastic model includes performing a multi-step system identification process.

In some embodiments, the processing circuit is configured to filter the time step specific values of the heat disturbance through at least one of an anti-spike filter or a smoothing filter.

In some embodiments, the processing circuit is configured to fit the time step specific values of the heat disturbance to at least one of a Gaussian function, a sinusoid function, or a user-defined function.

Another implementation of the present disclosure is a method for operating heating, ventilation, or air conditioning (HVAC) equipment of a building, according to some embodiments. The method includes estimating a thermal resistance between air of a zone of the building and air of an external space outside the zone using values of a temperature of the air of the zone, a temperature of the air of the external space, and a heat transfer rate of the HVAC equipment that provides heating or cooling to the zone, according to some embodiments. Each of the values correspond to a different time step within a time period, according to some embodiments. The method includes, for multiple time steps within the time period, using the thermal resistance between the air of the zone and the air of the external space, a time step specific value of the temperature of the air of the zone, a time step specific value of the temperature of the air of the external space, and a time step specific value of the heat transfer rate of the HVAC equipment to estimate a corresponding time step specific value of a heat disturbance, according to some embodiments. The method includes operating the HVAC equipment over the time period using a model-based control technique based on the time step specific values of the heat disturbance, according to some embodiments.

In some embodiments, the method includes performing a system identification process based on the time step specific values of the heat disturbance to identify a predictive model. The predictive model is used in the model-based control technique, according to some embodiments.

In some embodiments, the system identification process includes identifying a scaling parameter that scales the time step specific values of the heat disturbance.

In some embodiments, the method includes identifying a Kalman gain and a stochastic model of the heat disturbance based on the time step specific values of the heat disturbance or predicted heat disturbance values based on the time step specific values of the heat disturbance. The Kalman gain and the stochastic model used in the model-based control technique, according to some embodiments.

In some embodiments, identifying the Kalman gain and the stochastic model includes performing a multi-step system identification process.

In some embodiments, the method includes filtering the time step specific values of the heat disturbance through at least one of an anti-spike filter or a smoothing filter.

In some embodiments, the method includes fitting the time step specific values of the heat disturbance to at least one of a Gaussian function, a sinusoid function, or a user-defined function.

Another implementation of the present disclosure is a controller for operating heating, ventilation, or air conditioning (HVAC) equipment of a building, according to embodiments. The controller includes one or more processors, according to some embodiments. The controller includes one or more non-transitory computer-readable media storing instructions that, when executed by the one or more processors, cause the one or more processors to perform operations, according to some embodiments. The operations include estimating a thermal resistance between air of a zone of the building and air of an external space outside the zone using values of a temperature of the air of the zone, a temperature of the air of the external space, and a heat transfer rate of the HVAC equipment that provides heating or cooling to the zone, according to some embodiments. Each of the values correspond to a different time step within a time period, according to some embodiments. The operations include, for multiple time steps within the time period, using the thermal resistance between the air of the zone and the air of the external space, a time step specific value of the temperature of the air of the zone, a time step specific value of the temperature of the air of the external space, and a time step specific value of the heat transfer rate of the HVAC equipment to estimate a corresponding time step specific value of a heat disturbance, according to some embodiments. The operations include operating the HVAC equipment over the time period using a model-based control technique based on the time step specific values of the heat disturbance, according to some embodiments.

In some embodiments, the operations include performing a system identification process based on the time step specific values of the heat disturbance to identify a predictive model. The predictive model is used in the model-based control technique, according to some embodiments.

In some embodiments, the system identification process includes identifying a scaling parameter that scales the time step specific values of the heat disturbance.

In some embodiments, the operations include performing a multi-step system identification process based on the time step specific values of the heat disturbance or predicted heat disturbance values based on the time step specific values of the heat disturbance to identify a Kalman gain and a stochastic model of the heat disturbance. The Kalman gain and the stochastic model are used in the model-based control technique, according to some embodiments.

In some embodiments, the operations include filtering the time step specific values of the heat disturbance through at least one of an anti-spike filter or a smoothing filter.

In some embodiments, the operations include fitting the time step specific values of the heat disturbance to at least one of a Gaussian function, a sinusoid function, or a user-defined function.

DETAILED DESCRIPTION

Overview

Referring generally to the FIGURES, systems and methods for using system identification to estimate heat disturbances in building systems with complex sensor measurements and/or complex heating, ventilation, or air conditioning (HVAC) dynamics are shown and described. The systems and method described herein provide improved system models and therefore improved estimations of heat disturbance to more accurately heat/cool a building while optimizing costs related to said heating/cooling.

As described in the present disclosure, heat disturbance (also referred to as ({dot over (Q)}other) can refer to internal heat load/gains due to solar radiation, occupancy, electrical equipment, etc. A deterministic piece of heat disturbance can describe a portion of a total heat disturbance that can be determined based on parameter values and initial conditions of a heat disturbance estimation. A stochastic piece of heat disturbance can describe some inherent randomness in the heat disturbance. Estimating the stochastic piece can be quite difficult, but may be critical for a control system to generate accurate decisions. In performing control processes such as model predictive control (MPC), accurate estimations of the deterministic piece and the stochastic piece of heat disturbance can be critical for generating control decisions that optimize (e.g., reduce) costs, maintain occupant comfort, etc.

In some embodiments, a building zone group thermal model captures crucial system dynamics for MPC and control. The building zone group thermal model can be described by the following differential equations:

In the above building zone group thermal model, {dot over (Q)}HVAC, Toa, and {dot over (Q)}othercan be inputs to the thermal model. {dot over (Q)}HVACcan be measured and controlled, Toacan be measured but cannot be controlled, and {dot over (Q)}othercan be neither measured nor controlled. The states of the above thermal model can be seen as Tia, a measured state, and Tm, a non-measured state. An output to the thermal model can be Tia. To obtain values of the resistances and capacitances (i.e., Cia, Cm, Rmi, and Roi) of the thermal model, a system identification process can be performed. To perform the system identification, data for the building zone can be collected and used to solve for the parameters that provide accurate predictions of the building thermal dynamics.

In some embodiments, the system identification process is performed by a two-step process where a first step puts the building thermal model into a state space form and identifies the resistances and capacitances (i.e., Cia, Cm, Rmi, and Roi) using {dot over (Q)}HVACand Toaas inputs. Based on the first step, a measured output Tia, can be obtained. In a second step, a disturbance model can be augmented to the system to estimate the heat disturbance values. In general, the augmented model can be given by the following state space representation:

[x.⁡(t)d.⁡(t)]=[AcBd0Ad]⁡[x⁡(t)d⁡(t)]+[BcBd⁢d]⁢u⁡(t)y⁡(t)=[Cc⁢⁢Cd]⁡[x⁡(t)d⁡(t)]+Dc⁢u⁡(t)
where the disturbance model is characterized by the Ad, Bd, Bdd, and Cdmatrices. Further, the system can be converted into a discrete time model and an observer gain can be identified resulting in the following the closed-loop state estimation system:

However, augmenting the disturbance model and/or estimating the disturbance values may fail if an inaccurate disturbance model is used, if sensors utilize quantized measurements (i.e., measurements are rounded to a nearest integer), or if high frequency inputs permeate the data. To mitigate these issues, various regressions based on steady state assumptions can be performed to calculate an initial profile of the heat disturbance. The initial profile can be used in conjunction with various measurements of environmental conditions to obtain a model that accurately predicts changes in the heat disturbance over time. Particularly, the model can be used to predict scaled values of the heat disturbance over time and can be used in control processes to operate building equipment to optimize (e.g., reduce) overall costs and maintain occupant comfort. These and other features of the systems and methods are described in detail below.

Building HVAC Systems

The BMS that serves building10includes a HVAC system100. HVAC system100can include a plurality of HVAC devices (e.g., heaters, chillers, air handling units, pumps, fans, thermal energy storage, etc.) configured to provide heating, cooling, ventilation, or other services for building10. For example, HVAC system100is shown to include a waterside system120and an airside system130. Waterside system120may provide a heated or chilled fluid to an air handling unit of airside system130. Airside system130may use the heated or chilled fluid to heat or cool an airflow provided to building10.

HVAC system100thereby provides heating and cooling to the building10. The building10also includes other sources of heat transfer that the indoor air temperature in the building10. The building mass (e.g., walls, floors, furniture) influences the indoor air temperature in building10by storing or transferring heat (e.g., if the indoor air temperature is less than the temperature of the building mass, heat transfers from the building mass to the indoor air). People, electronic devices, other appliances, etc. (“heat load”) also contribute heat to the building10through body heat, electrical resistance, etc. Additionally, the outside air temperature impacts the temperature in the building10by providing heat to or drawing heat from the building10.

HVAC System and Model

Referring now toFIG. 2, a block diagram of the HVAC system100with building10is shown, according to an exemplary embodiment. More particularly,FIG. 2illustrates the variety of heat transfers that affect the indoor air temperature Tiaof the indoor air201in zone200of building10. Zone200is a room, floor, area, etc. of building10. In general, the primary goal of the HVAC system100is to maintain the indoor air temperature Ea in the zone200at or around a desired temperature to facilitate the comfort of occupants of the zone200or to meet other needs of the zone200.

As shown inFIG. 2, the indoor air temperature Tiaof the zone200has a thermal capacitance Cia. The indoor air temperature Tiais affected by a variety of heat transfers {dot over (Q)} into the zone200, as described in detail below. It should be understood that although all heat transfers {dot over (Q)} are shown inFIG. 2as directed into the zone200, the value of one or more of the heat transfers {dot over (Q)} may be negative, such that heat flows out of the zone200.

The heat load202contributes other heat transfer {dot over (Q)}otherto the zone200. The heat load202includes the heat added to the zone by occupants (e.g., people, animals) that give off body heat in the zone200. The heat load202also includes computers, lighting, and other electronic devices in the zone200that generate heat through electrical resistance, as well as solar irradiance.

The building mass204contributes building mass heat transfer {dot over (Q)}mto the zone200. The building mass204includes the physical structures in the building, such as walls, floors, ceilings, furniture, etc., all of which can absorb or give off heat. The building mass204has a temperature Tmand a lumped mass thermal capacitance Cm. The resistance of the building mass204to exchange heat with the indoor air201(e.g., due to insulation, thickness/layers of materials, etc.) may be characterized as mass thermal resistance Rmi.

The outdoor air206contributes outside air heat transfer {dot over (Q)}oato the zone200. The outdoor air206is the air outside of the building10with outdoor air temperature Toa. The outdoor air temperature Toafluctuates with the weather and climate. Barriers between the outdoor air206and the indoor air201(e.g., walls, closed windows, insulation) create an outdoor-indoor thermal resistance Roito heat exchange between the outdoor air206and the indoor air201.

The HVAC system100also contributes heat to the zone200, denoted as {dot over (Q)}HVAC. The HVAC system100includes HVAC equipment210, controller212, an indoor air temperature sensor214and an outdoor air temperature sensor216. The HVAC equipment210may include the waterside system120and airside system130ofFIG. 1, or other suitable equipment for controllably supplying heating and/or cooling to the zone200. In general, HVAC equipment210is controlled by a controller212to provide heating (e.g., positive value of {dot over (Q)}HVAC) or cooling (e.g., a negative value of {dot over (Q)}HVAC) to the zone200.

The indoor air temperature sensor214is located in the zone200, measures the indoor air temperature Tia, and provides the measurement of Tiato the controller212. The outdoor air temperature sensor216is located outside of the building10, measures the outdoor air temperature Toa, and provides the measurement of Toato the controller212.

The controller212receives the temperature measurements Toaand Tia, generates a control signal for the HVAC equipment210, and transmits the control signal to the HVAC equipment210. The operation of the controller212is discussed in detail below. In general, the controller212considers the effects of the heat load202, building mass204, and outdoor air206on the indoor air201in controlling the HVAC equipment210to provide a suitable level of {dot over (Q)}HVAC. A model of this system for use by the controller212is described with reference toFIG. 3.

In the embodiments described herein, the control signal provide to the HVAC equipment210by the controller110indicates a temperature setpoint Tspfor the zone200. To determine the temperature setpoint Tsp, the controller212assumes that the relationship between the indoor air temperature Tiaand the temperature setpoint Tspfollows a proportional-integral control law with saturation, represented as:
{dot over (Q)}HVAC,j=Kp,jεsp+Ki,j∫E0tεsp(s)ds(Eq. A)
εsp=Tsp,j−Tia(Eq. B)
where j∈{clg, hlg} is the index that is used to denote either heating or cooling mode. Different parameters Kp,jand Ki,jare needed for the heating and cooling mode. Moreover, the heating and cooling load is constrained to the following set: {dot over (Q)}HVAC,j∈[0, {dot over (Q)}cig,max] for cooling mode (j=clg) and {dot over (Q)}HVAC,j∈[−{dot over (Q)}htg,max0] for heating mode (j=htg). As discussed in detail below with reference toFIG. 4, the controller212uses this model in generating a control signal for the HVAC equipment210.

Referring now toFIG. 3, a circuit-style diagram300corresponding to the zone200and the various heat transfers {dot over (Q)} ofFIG. 2is shown, according to an exemplary embodiment. In general, the diagram300models the zone200as a two thermal resistance, two thermal capacitance, control-oriented thermal mass system. This model can be characterized by the following system of linear differential equations, described with reference toFIG. 3below:

Indoor air node302corresponds to the indoor air temperature Tia. From indoor air node302, the model branches in several directions, including down to a ground304via a capacitor306with a capacitance Cia. The capacitor306models the ability of the indoor air to absorb or release heat and is associated with the rate of change of the indoor heat transfer {dot over (T)}ia. Accordingly, the capacitor306enters Eq. C on the left side of the equation as Cia{dot over (T)}ia.

From indoor air node302, the diagram300also branches left to building mass node310, which corresponds to the thermal mass temperature Tm. A resistor312with mass thermal resistance Rmiseparates the indoor air node302and the building mass node310, modeling the heat transfer {dot over (Q)}mfrom the building mass204to the indoor air201as

1Rm⁢i⁢(Tm-Ti⁢a).
This term is included on the right side of Eq. C above as contributing to the rate of change of the indoor air temperature {dot over (T)}ia.

The diagram300also branches up from indoor air node302to outdoor air node314. A resistor316with outdoor-indoor thermal resistance Roiseparates the indoor air node302and the outdoor air node314, modeling the flow heat from the outdoor air206to the indoor air201as

1Ro⁢i⁢(To⁢a-Ti⁢a).
This term is also included on the right side of Eq. C above as contributing to the rate of change of the indoor air temperature {dot over (T)}ia.

Also from indoor air node302, the diagram300branches right to two {dot over (Q)} sources, namely {dot over (Q)}HVACand {dot over (Q)}other. As mentioned above, {dot over (Q)}othercorresponds to heat load202and to a variety of sources of energy that contribute to the changes in the indoor air temperature Tia. {dot over (Q)}otheris not measured or controlled by the HVAC system100, yet contributes to the rate of change of the indoor air temperature {dot over (T)}ia. {dot over (Q)}HVACis generated and controlled by the HVAC system100to manage the indoor air temperature Tia. Accordingly, {dot over (Q)}HVACand {dot over (Q)}otherare included on the right side of Eq. C above.

The second differential equation (Eq. D) above focuses on the rate of change {dot over (T)}min the building mass temperature T. The capacity of the building mass to receive or give off heat is modelled by capacitor318. Capacitor318has lumped mass thermal capacitance Cmand is positioned between a ground304and the building mass node310and regulates the rate of change in the building mass temperature Tm. Accordingly, the capacitance Cmis included on left side of Eq. D. Also branching from the building mass node310is resistor312leading to indoor air node302. As mentioned above, this branch accounts for heat transfer {dot over (Q)}mbetween the building mass204and the indoor air201. Accordingly, the term

1Rm⁢i⁢(Ti⁢a-Tm)
is included on the right side of Eq. D.

As described in detail below, the model represented by diagram300is used by the controller212in generating a control signal for the HVAC equipment210. More particularly, the controller212uses a state-space representation of the model shown in diagram300. The state-space representation used by the controller212can be derived by incorporating Eq. A and B with Eq. C and D, and writing the resulting system of equations as a linear system of differential equations to get:

As described in detail below, the controller212uses a two-step process to parameterize the system. In the first step, the controller212identifies the system parameters θ={θ1, θ2, θ3, θ4, θ5, θ6} (i.e., the values of Cia, Cm, Rmi, Roi, Kp,j, Ki,j). The disturbance state d is then introduced into the model and an Kalman estimator gain is added, such that in the second step the controller212identifies the Kalman gain parameters K. In some embodiments, the temperature setpoint Tspis not used as a system input, rather, {dot over (Q)}HVACis used as a direct input for controller212in generating a control signal for the HVAC equipment210.

As used herein, the term ‘variable’ refers to an item/quantity capable of varying in value over time or with respect to change in some other variable. A “value” as used herein is an instance of that variable at a particular time. A value may be measured or predicted. For example, the temperature setpoint Tspis a variable that changes over time, while Tsp(3) is a value that denotes the setpoint at time step3(e.g., 68 degrees Fahrenheit). The term “predicted value” as used herein describes a quantity for a particular time step that may vary as a function of one or more parameters.

Controller for HVAC Equipment with System Identification

Referring now toFIG. 4, a detailed diagram of the controller212is shown, according to an exemplary embodiment. The controller212includes a processing circuit400and a communication interface402. The communication interface402is structured to facilitate the exchange of communications (e.g., data, control signals) between the processing circuit400and other components of HVAC system100. As shown inFIG. 4, the communication interface402facilitates communication between the processing circuit400and the outdoor air temperature sensor216and the indoor air temperature sensor214to all temperature measurements Toaand Tiato be received by the processing circuit400. The communication interface402also facilitates communication between the processing circuit400and the HVAC equipment210that allows a control signal (indicated as temperature setpoint Tsp) to be transmitted from the processing circuit400to the HVAC equipment210.

The processing circuit400is structured to carry out the functions of the controller described herein. The processing circuit400includes a processor404and a memory406. The processor404may be implemented as a general-purpose processor, an application-specific integrated circuit, one or more field programmable gate arrays, a digital signal processor, a group of processing components, or other suitable electronic processing components. The memory406, described in detail below, includes one or more memory devices (e.g., RAM, ROM, NVRAM, Flash Memory, hard disk storage) that store data and/or computer code for facilitating at least some of the processes described herein. For example, the memory406stores programming logic that, when executed by the processor404, controls the operation of the controller212. More particularly, the memory406includes a training data generator408, a training data database410, a model identifier412, a model predictive controller414, and an equipment controller416. The various generators, databases, identifiers, controllers, etc. of memory406may be implemented as any combination of hardware components and machine-readable media included with memory406.

The equipment controller416is configured to generate a temperature setpoint Tspthat serves as a control signal for the HVAC equipment210. The equipment controller receives inputs of the indoor air temperature Tiafrom the indoor air temperature sensor214via the communication interface402and {dot over (Q)}HVACfrom the model predictive controller414(during normal operation) and the training data generator408(during a training data generation phase described in detail below). The equipment controller uses Tiaand {dot over (Q)}HVACto generate Tspby solving Eq. A and Eq. B above for Tsp. The equipment controller416then provides the control signal Tspto the HVAC equipment210via the communication interface402.

The model predictive controller414determines {dot over (Q)}HVACbased on an identified model and the temperature measurements Tia, Toa, and provides {dot over (Q)}HVACto the equipment controller416. The model predictive controller414follows a model predictive control (MPC) approach. The MPC approach involves predicting future system states based on a model of the system, and using those predictions to determine the controllable input to the system (here, {dot over (Q)}HVAC) that bests achieves a control goal (e.g., to maintain the indoor air temperature near a desired temperature). A more accurate model allows the MPC to provide better control based on more accurate predictions. Because the physical phenomena that define the behavior of the system (i.e., of the indoor air201in the building10) are complex, nonlinear, and/or poorly understood, a perfect model derived from first-principles is generally unachievable or unworkable. Thus, the model predictive controller414uses a model identified through a system identification process facilitated by the training data generator408, the training data database410, and the model identifier412, described in detail below.

System identification, as facilitated by the training data generator408, the training data database410, and the model identifier412, is a process of constructing mathematical models of dynamic systems. System identification provides a suitable alternative to first-principles-derived model when first principles models are unavailable or too complex for on-line MPC computations. System identification captures the important and relevant system dynamics based on actual input/output data (training data) of the system, in particular by determining model parameters particular to a building or zone to tune the model to the behavior of the building/zone. As described in detail below, the training data generator408, the training data database410, and the model identifier412each contribute to system identification by the controller212.

The training data generator408is configured to generate training data by providing an excitation signal to the system. In some embodiments, the training data generator408provides various {dot over (Q)}HVACvalues to the equipment controller416for a number N of time steps k, and receives the measured output response of the indoor air temperature Tiaat each time step k from the air temperature sensor214. The various {dot over (Q)}HVACvalues may be chosen by the training data generator408to explore the system dynamics as much as possible (e.g., across a full range of possible {dot over (Q)}HVACvalues, different patterns of {dot over (Q)}HVACvalues, etc.). In some embodiments, the training data generator408provides various Tspvalues to the equipment controller416instead of the various values of {dot over (Q)}HVAC.

If the equipment controller416receives the various {dot over (Q)}HVACvalues, various control inputs Tspcan be generated in response. The temperature setpoint Tspfor each time step k is provided to the HVAC equipment210, which operates accordingly to heat or cool the zone200(i.e., to influence Tia). In some embodiments, the temperature setpoints Tspare used by the training data generator408to be included in the training data. The training data generator receives an updated measurement of the indoor air temperature Tiafor each time step k and may also receive the outdoor air temperature Toafor each time step k. The training data generator408thereby causes the states, inputs, and outputs of the system to vary across the time steps k and generates data corresponding to the inputs and outputs.

The inputs and outputs generated by the training data generator408are provided to the training data database410. More particularly, in the nomenclature of the model of Eq. E and Eq. F above, the training data generator408provides inputs Tspand Toaand outputs {dot over (Q)}HVACand Ea for each time step k to the training data database410.

The training data database410stores the inputs and outputs for each time step k provided by the training data generator408. Each input and output is tagged with a time step identifier, so that data for the same time step can be associated together. The training data database410thereby collects and stores input and output data for each time step k, k=0, . . . N, or, more specifically, Tsp(k), Toa(k), Tia(k), and {dot over (Q)}HVAC(k), for k, k=0, . . . , N. This data is grouped together in the training data database410in a set of training data ZN. In the notation of Eq. G and Eq. H, ZN=[y(1), u(1), y(2), u(2), . . . , y(N), u(N)].

In some embodiments, the training data is refined using a saturation detection and removal process. System and methods for saturation detection and removal suitable for use to refine the training data ZNare described in U.S. patent application Ser. No. 15/900,459, filed Feb. 20, 2018, incorporated by reference herein in its entirety. For example, as described in detail therein, the training data may be filtered by determining whether the operating capacity is in a non-transient region for a threshold amount of a time period upon determining that an error for the building zone exists for the time period, and in response to a determination that the operating capacity is in the non-transient region for at least the threshold amount of the time period, indicating the time period as a saturation period. Data from the saturation period can then be removed from the training data.

The model identifier412accesses the training data database410to retrieve the training data ZNand uses the training data ZNto identify a model of the system. The model identifier412includes a system parameter identifier418and a gain parameter identifier420. As shown in detail inFIG. 5and discussed in detail with reference thereto, the system parameter identifier418carries out a first step of system identification, namely identifying the model parameters, while the gain parameter identifier420carries out the second step, namely determining a Kalman gain estimator. The model parameters and the Kalman gain estimator are included in an identified model of the system, and that model is provided to the model predictive controller414. The model predictive controller can thus facilitate the control of the HVAC equipment210as described above.

Referring now toFIG. 5, a detailed view of the model identifier412is shown, according to an exemplary embodiment. As mentioned above, the model identifier412includes the system parameter identifier418and the gain parameter identifier420. The system parameter identifier418is structured to identify the matrices A, B, C, D of Eqs. G and H, i.e., the values of θ={θ1, θ2, θ3, θ4, θ5, θ6}. In the embodiment described herein, this corresponds to finding the values of Cia, Cm, Rmi, Roi, Kp,j, and Ki,j.

The model framework identifier422thereby determines that the system parameter identifier418has the goal of determining a parameter vector {circumflex over (θ)}Nfrom the set of θ∈⊂d, whereis the set of admissible model parameter values. The resulting possible models are given by the set: M={(θ), θ∈}. The goal of the system parameter identifier418is to select a parameter vector {circumflex over (θ)}Nfrom among possible values of 0 that best matches the model to the physical system (i.e., the vector θ is a list of variables and the vector {circumflex over (θ)}Nis a list of values), thereby defining matrices A, B, C, and D. The model framework identifier422also receives training data ZNand sorts the training data (i.e., Tsp(k), Toa(k), Tia(k), and {dot over (Q)}HVAC(k), for k, k=0, . . . , N) into the notation of Eq. G-H as input/output data ZN=[y(1), u(1), y(2), u(2), . . . , y(N), u(N)].

The prediction error function generator424receives the model framework M={(θ), θ∈} and the training data ZNfrom the model framework identifier422. The prediction error function generator424applies a prediction error method to determine the optimal parameter vector {circumflex over (θ)}N. In general, prediction error methods determine the optimal parameter vector {circumflex over (θ)}Nby minimizing some prediction performance function VN(θ,ZN) that is based in some way on the difference between predicted outputs and the observed/measured outputs included in the training data ZN. That is, the parameter estimation θNis determined as:
{circumflex over (θ)}N={circumflex over (θ)}N(ZN)=argVN(θ,ZN).

The prediction error function generator424use one or more of several possible prediction error approaches to generate a prediction performance function VN(θ,ZN). In the embodiment shown, the prediction error function generator applies a simulation approach. In the simulation approach, the prediction error function generator424uses the model(θ), the input trajectory [u(1),u(2), . . . ,u(N)], and an initial state x(0) to produce predicted outputs in terms of 0. That is, the prediction error function generator424predicts:
[ŷ(1|0,θ),ŷ(2|0,θ) . . .ŷ(k|0,θ) . . . ,ŷ(N|0,θ)],
where ŷ(k|0, θ) denotes the predicted output at time step k given the training data from time0and the model(θ). The prediction error function generator424then calculates a prediction error at each time step k is given by ε(k, θ):=y(k)−ŷ(k|0, θ). The prediction error function generator424then squares the two-norm of each prediction error ε(k θ0) and sums the results to determine the prediction performance function, which can be written as:
VN(θ,ZN)=Σk=1N∥y(k)−ŷ(k|0,θ)∥22(Eq. I).

In an alternative embodiment, the prediction error function generator424applies a one-step-ahead prediction error method to generate the prediction performance function VN(θ, ZN). In the one-step-ahead prediction error method, the prediction error function generator424uses past input-output data and the model(θ) the model to predict the output one step ahead in terms of θ. That is, in the one-step ahead prediction error method, the prediction error function generator424generates one-step ahead predictions ŷ(k|k−1, θ), which denotes the predicted output at time step k given the past input-output sequence Zk−1and using parameters θ. The one-step ahead prediction ŷ(k|k−1, θ) is then compared to the measured output y(k) by the prediction error function generator424to determine the prediction error at k, defined as ε(k,θ):=y(k)−ŷ(k|k−1, θ). The prediction error function generator424then squares the two-norm of the prediction errors for each k and sums the results, generating a prediction performance function that can be expressed in a condensed form as:

In other alternative embodiments, the prediction error function generator424uses a multi-step ahead prediction error approach to generate the prediction performance function. The multi-step ahead prediction error approach is described in detail below with reference to the gain parameter identifier420andFIGS. 7-8.

The prediction error function generator424then provides the performance function VN(θ,ZN) (i.e., from Eq. I or Eq. J in various embodiments) to the optimizer426.

The optimizer426receives the prediction error function generated by the prediction error function generator424and optimizes the prediction error function in θ to determine {circumflex over (θ)}N. More specifically, the optimizer426finds the minimum value of the prediction error function VN(θ,ZN) as θ is varied throughout the allowable values of θ∈. That is, the optimizer426determines {circumflex over (θ)}Nbased on:
{circumflex over (θ)}N={circumflex over (θ)}N(ZN)=argVN(θ,ZN).

The optimizer426then uses {circumflex over (θ)}Nto calculate the matrices A, B, C, and D. The system parameter identifier418then provides the identified matrices A, B, C, D to the gain parameter identifier420.

The gain parameter identifier420receives the model with the matrices A, B, C, D (i.e., the model parameters) from system parameter identifier418, as well as the training data ZNfrom the training data database410, and uses that information to identify the gain parameters. The gain parameter identifier420includes an estimator creator428, a prediction error function generator430, and an optimizer432.

The estimator creator428adds a disturbance model and introduces a Kalman estimator gain to account for thermal dynamics of the system, for example for the influence of {dot over (Q)}otheron the system. The estimator creator428generates an augmented model with disturbance state d, given by:

[x.⁡(t)d.⁡(t)]=[AcBd00]⁡[x⁡(t)d⁡(t)]+[Bc0]⁢u⁡(t);y⁡(t)=[Cc⁢⁢Cd]⁢{x⁡(t)d⁡(t)}+Dc⁢u⁡(t)
where the parameters Ac, Bc, Cc, and Dcare the matrices A, B, C, D received from the system parameter identifier418and the disturbance model is selected with

The estimator creator428then converts the model to a discrete time model, for example using 5-minute sampling periods, resulting in the matrices Adis, Bdis, Cdis, Ddisand the disturbance model discrete time matrix Bddis. The estimator creator428then adds a parameterized estimator gain, resulting in the following model:

The matrix K(ϕ) is the estimator gain parameterized with the parameter vector ϕ where:

In this notation, {circumflex over (x)}(t+1|t) is an estimate of the state at time t+1 obtained using the Kalman filter and made utilizing information at sampling time t. For example, with a sampling time of five minutes, {circumflex over (x)}(t+1|t) is an estimate of the state five minutes after the collection of the data that the estimate is based on. The goal of the gain parameter identifier is to identify parameters {circumflex over (ϕ)}N(i.e., a vector of for each of ϕ1. . . ϕ8) that make the model best match the physical system.

The estimator creator428then provides the discrete time model with estimator gain (i.e., Eqs. K-L) to the prediction error function generator430. The prediction error function generator receives the model from the estimator creator428as well as the training data ZNfrom the training data database410, and uses the model (with the estimator gain) and the training data ZNto generate a prediction performance function.

The prediction error function generator430follows a multi-step ahead prediction error method to generate a predication performance function VN(ϕ, ZN). The multi-step ahead prediction error method is illustrated inFIGS. 7-8and described in detail with reference thereto. As an overview, in the multi-step-ahead prediction error method, the prediction error function generator430uses past input-output data and the model(θ) the model to predict the output multiple step ahead in terms of ϕ. That is, in the multi-step ahead prediction error method, the prediction error function generator430generates multi-step ahead predictions ŷ(k+h|k−1, ϕ), which denotes the predicted output at time step k+h given the past input-output sequence Zk−1and using parameters ϕ. The index h corresponds the number of steps ahead the prediction is made, and for each time step k predictions are made for h=0, . . . , hmax(i.e., when h=2, the prediction is three steps ahead because h is indexed from zero).

Each multiple multi-step ahead prediction ŷ(k+h|k 1, ϕ) is then compared to the corresponding measured output y(k) by the prediction error function generator430to determine the prediction error at k, defined as ε(k, θ):=y(k)−ŷ(k+h|k−1, ϕ). The prediction error function generator430then squares the two-norm of the prediction errors for each k and sums the results, in some embodiments using an weighting function w(h). The prediction error function generator430thereby generates a prediction performance function that can be expressed in a condensed form as:

The multi-step ahead prediction error method is described in more detail below with reference toFIGS. 7-8. In alternative embodiments, the prediction error function generator430follows the simulation approach or the one-step ahead prediction error approach discussed above with reference to the prediction error function generator424.

The prediction error function generator430then provides the prediction performance function (i.e., Eq. M) to the optimizer432. The optimizer432receives the prediction error function VN(ϕ, ZN) generated by the prediction error function generator430and optimizes the prediction error function in ϕ to determine {circumflex over (ϕ)}N. More specifically, the optimizer426finds the minimum value of the prediction error function VN(ϕ, ZN) as ϕ is varied throughout the allowable values of ϕ. In some cases, all real values of ϕ are allowable. That is, the optimizer426determines {circumflex over (ϕ)}Nbased on:
{circumflex over (θ)}N={circumflex over (ϕ)}N(ZN)=arg minϕVN(ϕ,ZN).

The optimizer432then uses {circumflex over (ϕ)}Nto calculate the matrices Kx(ϕ) and Kd(ϕ), resulting in a fully identified model. The gain parameter identifier420provides the identified model to the model predictive controller414.

In some embodiments, the prediction error function generator430reconfigures the multi-step ahead prediction problem by defining augmented vectors that allow the multi-step ahead prediction performance function (Eq. M) to be recast in an identical structure to the single-step ahead prediction performance function (Eq. J). Existing software toolboxes and programs (e.g., Matlab system identification toolbox) configured to handle the single-step ahead prediction error approach can then be used to carry out the multi-step ahead prediction error approach. To reconfigure the problem for that purpose, the prediction error function generator430considers, the system model of the form:
x(k+1)=Ax(k)+Bu(k);
y(k)=Cx(k)+Du(k).
where the one-step prediction of {circumflex over (x)}(k+1|k) using a steady-state Kalman gain is:
{circumflex over (x)}(k+1|k)=A{circumflex over (x)}(k|k−1)+Bu(k)+K(y(k)−C{circumflex over (x)}(k|k−1)−Du(k));
ŷ(k|k−1)=C{circumflex over (x)}(k|k−1)+Du(k).

In the multi-step prediction Kalman gain system identification problem, the complete pattern of the algebraic manipulations is shown by the 4-step prediction. The prediction error function generator430considers a case with four input data points and four output data-points starting from time h=0 to time h=3, so that hmax=3. The one-step prediction (with the prediction error function generator430given x0) is given by the equation:
{circumflex over (x)}(1|0)=Ax0+Bu(0)+K(y(0)−Cx0−Du(0));
ŷ(0|0)=Cx0+Du(0).

The prediction of the second step is
{circumflex over (x)}(2|0)=A{circumflex over (x)}(1|0)+Bu(1)=A(Ax0+Bu(0)+K(y(0)−Cx0−Du(0)))+Bu(1);
ŷ(1|0)=C{circumflex over (x)}(1|0)+Du(1)=C(Ax0+Bu(0)+K(y(0)−Cx0−Du(0)))+Du(1).

The prediction of the third step is

The forth step prediction is

With these 4-step predictions, the pattern needed to cast the multi-step prediction problem as a 1-step prediction is revealed. By aggregating the matrices multiplying x0, y(0), u(0), u(1), u(2), and u(3), the pattern revealed is:
{circumflex over (x)}(1|0)=Ax0+Bu(0)+K(y(0)−Cx0−Du(0));
{circumflex over (x)}(2|0)=(A2−AKC)x0+(AB−AKD)u(0)+Bu(1)+AKy(0);
{circumflex over (x)}(3|0)=(A3−A2KC)x0+(A2B−A2KD)u(0)+ABu(1)+Bu(2)+A2Ky(0);
{circumflex over (x)}(4|0)=(A4−A3KC)x0+(A3B−A3KD)u(0)+A2Bu(1)ABu(2)+Bu(3)+A3Ky(0);
ŷ(0)=Cx0+Du(0);
ŷ(1|0)=(CA−CKC)x0+(CB−CKD)u(0)+Du(1)+CKy(0);
ŷ(2|0)=(CA2−CAKC)x0+(CAB−CAKD)u(0)+CBu(1)+Du(2)+CAKy(0);
ŷ(3|0)=(CA3−CA2KC)x0+(CA2B−CA2KD)u(0)+CABu(1)+CBu(2)+Du(3)+CA2Ky(0).

Based on that pattern, the prediction error function generator430defines the following vectors:

The new system that has the 4-step prediction casted into a one-step prediction which can be analyzed by the prediction error function generator430using an existing system identification software product as:

In order to have the general formulation at time k for predicting hmaxstep ahead in time, this four-step example can be extrapolated to define the general augmented input and output vectors as:

With these definition, the general formulation at time k for predicting hmaxsteps ahead in time is:
{circumflex over (x)}(k+1|k)=A{circumflex over (x)}(k|k−1)+[B0]ũ(k)+[K0 . . . 0]({tilde over (y)}(k)−{circumflex over ({tilde over (y)})}(k).

As described above, in the multi-step ahead prediction error method the prediction error function generator430generates a function of the form:

If w(h)≡1 for all h, and using the augmented input and output vectors defined above, the multi-step ahead prediction performance function can be reconfigured into the following one-step ahead prediction performance function by the prediction error function generator430:

The prediction error function generator430then uses this reconfigured format of the prediction performance function with existing software toolboxes suited for the one-step ahead prediction error approach. The prediction error function generator430may include machine-readable media storing computer code executable to apply such software.

System Identification Methods

Referring now toFIG. 6, a flowchart of a process600for system identification is shown, according to an exemplary embodiment. The process600can be carried out by the controller212ofFIGS. 2 and 4.

At step602, the controller212applies an excitation signal to the HVAC equipment210. For example, the training data generator408may vary the {dot over (Q)}HVACvalues supplied to the equipment controller416, causing an excitation signal to be generated in the temperature setpoint Tspinputs provided to the HVAC equipment210. In general, the excitation signal is designed to test the system in a way to provide robust data for use in system identification.

At step604, training data is collected and stored by the controller212. Training data includes measureable temperature readings, i.e., Toaand Tia, controller-determined values {dot over (Q)}HVACand Tspfor each of a plurality of time steps k, k=0, . . . , N. The training data therefore includes inputs u(k) and the outputs y(k) for the time period. The training data is received from temperature sensors214,216, training data generator408, and/or equipment controller416and stored in training data database410.

At step606, the controller212identifies the model parameters for the system. That is, as discussed in detail above, the controller212determines the matrices A(θ), B(θ), C(θ), and D(θ) that minimize a prediction performance function VN(ZN, θ) for the model:
{dot over (x)}+(t)=Ac(θ)x(t)+Bc(θ)u(t);  (Eq. G)
y(t)=Cc(θ)x(t)+Dc(θ)u(t);  (Eq. H).

In identifying the model parameters, a simulation approach or a one-step-ahead prediction error approach is followed in some embodiments. These approaches are described in detail above with reference to the prediction error function generator424ofFIG. 5. In other embodiments, the model parameters are determined at step606using a multi-step ahead prediction error method, described in detail with reference toFIGS. 7-8.

At step608, the controller212identifies the gain estimator parameters. That is, the controller212determines the matrices Kxand Kdof Eq. K above. In preferred embodiments, the controller212uses the multi-step ahead prediction error method to find the matrices Kxand Kd. The multi-step ahead prediction error method is described in detail below with reference toFIGS. 7-8. In alternative embodiments, a simulation approach or a one-step-ahead prediction error approach is followed to find the matrices Kxand Kd.

At step610, the identified model is validated by the controller212. The controller212uses the identified model to generate control signal inputs Tspfor the HVAC equipment210using model predictive control. The controller then monitors the temperature measurements Toaand Tiafrom temperature sensors214,216, the input Tsp, and the value {dot over (Q)}HVACto determine how well the model matches system behavior in normal operation. For example, the training data database410may collect and store an addition set of training data that can be used by the model identifier412to validate the model. If some discrepancy is determined, the identified model may be updated. The identified model can thereby by dynamically adjusted to adjust for changes in the physical system.

Referring now toFIGS. 7-8the multi-step ahead prediction error approach for use in system identification is illustrated, according to an exemplary embodiment. InFIG. 7, a flowchart of a process700for identifying system parameters using the multi-step ahead prediction error approach is shown, according to an exemplary embodiment.FIG. 8shows an example visualization useful in explaining process700. Process700can be carried out by the system parameter identifier418and/or the gain parameter identifier420ofFIG. 5. In the embodiment described herein, the process700is implemented with the gain parameter identifier420.

Process700begins at step702, where the gain parameter identifier420receives training data ZN=[y(1), u(1), y(2), u(2), . . . , y(N), u(N)] from the training data database410. The training data includes measured outputs y(k) (i.e., Tia(k) and {dot over (Q)}HVAC(k)) and inputs u(k) (i.e., Toa(k) and Tsp(k)) for each time step k, k=1, . . . , N. N is the number of samples in the training data. The gain parameter identifier420also receives the system model from the system parameter identifier418.

At step704, the prediction error function generator430uses the training data for a time step k to predict outputs ŷ for each subsequent time step up to k+hmax. The value hmaxcorresponds to the number of steps ahead the predictions are made, referred to herein as the prediction horizon. Because hmaxis indexed from zero in Eq. M above, the prediction horizon is one more than the value of hmax. For example, in the case shown inFIG. 8and described below, predictions are made three steps ahead, corresponding to hmax=2 in the notation of Eq. D and a prediction horizon of three. The prediction horizon may be any integer greater than one, for example four or eight. The prediction horizon can be tuned experimentally, to determine an ideal prediction horizon. For example, too long of a prediction horizon may lead to poor prediction while too short of a prediction horizon may suffer the same limitations as the one-step ahead prediction error method mentioned above. In some cases, a prediction horizon of eight is preferred.

More specifically, at each step704the predicted outputs [ŷ(k|k−1), ŷ(k+1|k−1), . . . ŷ(k+hmax|k−1)] are predicted based on the past training data (i.e., through step k−1), denoted as Zk−1, along with future inputs [u(k), u(k+1) . . . u(k+hmax)]. These predictions are made using the model M(ϕ), such that predicted outputs ŷ depend on ϕ.

To illustrate the predictions of step704,FIG. 8shows a simplified visualization in which y(k) and ŷ(k) are depicted as scalar values for the sake of simplified explanation. InFIG. 8, the graph800plots the values of y and ŷ over time t for five time steps past a starting time t=0. The solid circles802represent measured outputs y(t) from the training data. The unfilled boxes804represent predicted outputs ŷ(t|0), that is, the outputs predicted for each time step based on the input/output data available at time t=0 (e.g., y(0)). The dashed lines represent the propagation of the predictions; for example, graph800includes three unfilled boxes804connected by a dashed line to the solid circle802corresponding to y(0). This shows that the predictions ŷ(t|0), 1≤t≤3, represented by the unfilled boxes804were based on the measured value of y(0).

At step706, the prediction error function generator430compares the predicted outputs ŷ to the measured outputs y for each future step up to k+hmax(i.e., for all predicted outputs ŷ generated at step704). More specifically, an error term for each step may be defined as y(k+h)−ŷ(k+h|k−1,ϕ). Because y and ŷ are vectors, the two-norm of this error term may be taken and squared to facilitate comparison between prediction errors as scalars, such that the error term becomes ∥y(k+h)−ŷ(k+h|k−1, ϕ)∥23. This term appears in Eq. M above.

As shown inFIG. 8, step706can be understood as measuring the distance between, for example, each unfilled box804and the corresponding solid circle802(i.e., the unfilled box804and the solid circle802at the same time t). Thus, in the example ofFIG. 8, step706includes calculating three error terms.

At step708, the error terms are weighted based on a weighting function w(h). The weighting function w(h) allows the prediction errors to be given more or less weight depending on how many steps ahead the prediction is. The weighting function w(h) is preferably a monotonically decreasing function of h, so that farther-out-in-time predictions have less influence on the prediction error. In some embodiments, the weighting function w(h)=1. Step708thereby corresponds the w(h) term in Eq. M above.

The process700then returns to step704to repeat steps704-706for each value of k, k=1, As illustrated inFIG. 8, repeating step704corresponds to generating the predictions represented by the unfilled circles808and the unfilled triangles810. The unfilled circles808chart the predictions based on the output data available at time t=1, i.e., ŷ(t|1), for t=2, 3, 4. The unfilled triangles chart the predictions based on the output data available at time t=2, i.e., ŷ(t|2), for t=3, 4, 5. Process700therefore involves making multiple predictions for most time steps: for example,FIG. 8shows three separate predictions for time t=3.

At step706, the prediction error function generator430again compares the predicted outputs ŷ for the new value of k to the measured outputs y for each future step up to k+hmaxto define the error term ∥y(k+h)−ŷ(k+h|k−1, 0)∥22as included in Eq. M. At step708, the terms are again weighted by the weighting function w(h). The weighting function w(h) may be the same for each k.

In the notation of Eq. M, each iteration of steps704-708thus corresponds to steps necessary to generate the values used by the inner (right) summation indexed in h, while repetition of the steps704-708corresponds to the iteration through k represented in the outer (left) summation. At step710, then, these summations are executed. In other words, the system identification circuit108sums the weighted error terms generated by steps704-708to generate a prediction performance function as:

The prediction performance function is a function of the input data ZNand the parameter variable ϕ. Typically, the input data ZNis given (i.e., received by the model identifier412and used in the calculation of error terms as described above). Thus, the prediction performance function is primarily a function of ϕ.

At step712, the prediction performance function VN(ϕ, ZN) is minimized to find an optimal parameter vector {circumflex over (θ)}N=argVN(ϕ, ZN). Any minimization procedure may be followed. The result of step712is a vector {circumflex over (ϕ)}Nof identified model parameters that tune the model({circumflex over (ϕ)}N) to accurately predict system evolution multiple steps ahead. At step714, the model identifier412provides the identified system model (i.e.,({circumflex over (ϕ)}N)) to the model predictive controller414for use in generating control inputs for the HVAC equipment210.

According to various embodiments, process700is run once at set-up to establish the system model, run periodically to update the system model, or run repeatedly/continuously to dynamically update the system model in real time.

Experimental Results

To illustrate the advantages of the systems and methods described above, the following experimental results are included and shown inFIGS. 9-18. The HVAC system100and the building10were put through two experiments: a heating and a cooling experiment.

Heating Experiment

In the heating experiment, a simulated HVAC system100is in a heating mode to heat a simulated building10. Because the system is simulated the actual values of the system parameters and the unmeasured time-varying disturbances ({dot over (Q)}other) are known in the experiment for sake of comparison to the identified model.

To start, the controller212provides excitation signal to the HVAC equipment210. The excitation signal902illustrated on graph900inFIG. 9varies the temperature setpoint Tspusing a pseudorandom binary signal that varies between the maximum and minimum allowable temperatures in the comfort zone (Tmax=25° C., Tmin=23° C.).

While the excitation signal is applied to the HVAC equipment210, training data is collected and stored in the controller212for each time step k as described above. Each time step kin the heating experiment corresponds to five minutes of time (i.e., a data sample is recorded every five minutes). The training data is used by the system parameter identifier418to identify the model parameters as described above. In the heating experiment, the following results were found:

The first step of model parameterization, carried out by the system parameter identifier418, thereby determined the parameters to the correct order of magnitude, but some differences are present due to the time-varying disturbances (i.e., {dot over (Q)}other).

Next, the Kalman gain parameters are identified by the gain parameter identifier420. In the experiment, the gain parameters are identified using a one-step ahead prediction error method, a two-step ahead prediction error method, a five-step ahead prediction error method, a ten-step ahead prediction error method, and a fifty-step ahead prediction error method. As an example of the results,FIG. 10shows a graph1000of the actual indoor temperature and the predicted indoor temperature over time as generated by the five-step ahead prediction error method.FIG. 11shows a graph1100of the actual {dot over (Q)}HVACand the predicted {dot over (Q)}HVACover time as generated by the five-step ahead prediction error method. As shown inFIGS. 10 and 11, the predicted values of Tiaand {dot over (Q)}HVACconsistently track the actual values.

The different number of steps (i.e., hmaxvalues) were included to allow comparison between the parameters identified using different numbers of steps. The Kalman gains identified using the various numbers of steps are presented in the following table:

FIG. 12shows a graph1200of the estimated building mass temperature Tmover time for the one-step prediction error method, the two-step prediction error method, and the fifty-step prediction error method, as well as the actual Tmof the simulated building10. As the number of steps increase, the Tmestimates improve, following the actual Tmline on graph1200closer.

Cooling Experiment

In the cooling experiment, a simulated HVAC system100is in a cooling mode to cool a simulated building10. As above, because the system is simulated the actual values of the system parameters and the unmeasured time-varying disturbances ({dot over (Q)}other) are known in the experiment for sake of comparison to the identified model.

A similar procedure as the heating experiment is followed to generate models, with the Kalman gain generated using the multi-step prediction error method with a variety of number of steps (i.e., various prediction horizons hmax) (e.g., one step, two steps, eight steps, twelve steps, twenty steps).FIG. 13shows output predictions generated using the one-step prediction error method compared to actual outputs, with Tiagraphed over time on graph1300and {dot over (Q)}HVACgraphed over time on graph1350. Similarly,FIG. 14shows output predictions generated using the two-step prediction error method compared to actual outputs, with Tiagraphed over time on graph1400and {dot over (Q)}HVACgraphed over time on graph1450.

To compare the results of the various multi-step prediction error methods (i.e., various number of steps), several metrics are used. First, a weighted mean absolute prediction error (WMAPE) metric is an exponentially weighted average of the absolute prediction error at each time step and given by:

WMAPE⁡(k)=∑i=kk+Nh-1⁢e-i/Nh⁢y⁡(i)-y^⁡(i❘k)∑i=kk+Nh-1⁢e-i/Nh,k=0,1,2,…⁢.
where Nh∈>0is the prediction horizon, y(i) is the actual output at time step i and ŷ(i|k) is the predicted output with the identified model given a measurement at time step k and the input sequence u(k), u(k+1), . . . , u(i−1). In the WMAPE equation, y is used to refer to a scalar (i.e., one of the two outputs), and the WMAPE is computed separately for both outputs. The horizon used to calculate the WMAPE in the cooling experiment was twelve.FIG. 15shows a graph1500of the WMAPE for Tiafor the one-step ahead prediction error method and a 12-step ahead prediction error method for comparison.FIG. 15also shows a graph1550of the WMAPE for {dot over (Q)}HVACfor the one-step ahead prediction error method and a 12-step ahead prediction error method for comparison.

Another metric that can be used to evaluate the results of the cooling experiment is the root mean squared prediction error (RMSPE). RMSPE is calculated for a range of values of q from zero-step ahead prediction up to Nh-step ahead prediction. That is, given a set of measured output values {y(0), . . . , y(M)} for M∈≥0, the RMSPE is given by:

RMSPE⁡(q)=∑i=q+1M⁢(y⁡(i)-y^⁡(i|i-q))2M-q
for all q∈{0, . . . ,Nh−1}. The RMSPE helps identify the prediction error over the prediction horizon. In the example here, the RMSPE is calculated for 288 steps (i.e., Nh=288).FIG. 16shows a graph1600of the RMSPE for Ea for the one-step ahead prediction error method and a two-step ahead prediction error method for comparison.FIG. 16also shows a graph1650of the WMAPE for {dot over (Q)}HVACfor the one-step ahead prediction error method and the two-step ahead prediction error method for comparison.

A third way to compare across different numbers of steps is to visualize the power of prediction.FIG. 17andFIG. 18shows examples visualizations1700and1800of this third metric. To generate the visualizations1700,1800, ten lines of N-steps-ahead predictions are plotted using the Kalman gain generated by each multi-step ahead prediction method. That is, a first line starts x0(i.e., an initial state) and plots the N step ahead prediction, from {circumflex over (x)}(1|0) all the way to {circumflex over (x)}(N|0). The second line takes {circumflex over (x)}(1|0) and plots N steps ahead, and so on, until ten lines are plotted. The closer the lines are to being on top of each other, the better the output multi-step prediction. In the examples ofFIGS. 15 and 16, the lines are plotted for twelve steps ahead (N=12).

The visualization1700ofFIG. 17is thereby generated for each of the one-step ahead prediction error method, the eight-step ahead prediction error method, the twelve-step ahead prediction error method, the twenty-step ahead prediction error method for the output Tia. The visualization1800ofFIG. 18is generated in the same way for each of the one-step ahead prediction error method, the eight-step ahead prediction error method, the twelve-step ahead prediction error method, the twenty-step ahead prediction error method for the output {dot over (Q)}HVAC.

In both visualization1700and visualization1800, the eight-step prediction error method is shown to have the best results (i.e., the lines are plotted closest together), even though the lines were plotted twelve steps ahead. Thus, in some embodiments, an eight-step ahead prediction error method may be preferred (i.e., hmax=7). Because each time step is five minutes in the experiment, this implies that a prediction horizon of forty minutes in the Kalman gain identification is well suited for generating a model that predicts one hour (12 steps) into the future.

Heat Load Modeling

Overview

Referring generally toFIGS. 19-21, methods for estimating heat disturbance Qotherare shown, according to some embodiments. Heat disturbance refers to heat in a building (or any space) that originates from sources beyond measurement and/or control of an environmental control system of the building. For example, heat disturbance may result from sunlight, heat radiating from electrical equipment, body heat radiation, etc. Accurately estimating heat disturbance can increase accuracy of estimations made during a model predictive control process. Without estimations of heat disturbance, a significant source of heat in a building may go unaccounted for, thus reducing accuracy of model predictive control and increasing energy usage and/or occupant discomfort.

As explained in greater detail below, heat disturbance can be modeled as a summation of a deterministic heat disturbance prediction and a stochastic heat disturbance prediction. The deterministic heat disturbance can describe a piece of a total heat disturbance that can be determined based on parameter values and initial conditions of a heat disturbance estimation problem. In some embodiments, the deterministic heat disturbance is calculated using a process for estimating deterministic load as described in U.S. patent application Ser. No. 14/717,593 filed May 20, 2015, incorporated by reference herein in its entirety. However, the stochastic heat disturbance, a piece of the total heat disturbance that describes some inherent randomness of the heat disturbance, can be difficult to calculate.

In some embodiments, the deterministic heat disturbance model is obtained using a pattern recognition and linear regression strategy as described in U.S. patent application Ser. No. 14/717,593 filed May 20, 2015, incorporated by reference herein in its entirety. In some embodiments, the stochastic heat disturbance model is obtained through identification of an autoregressive (AR) model separate from a system state space model used in model predictive control (MPC). In some embodiments, the stochastic heat disturbance model is obtained through identification of a model that is part of an overall state space model used in estimation and prediction during an MPC process.

Processes for Heat Load Modeling

Referring now toFIG. 19, a process1900for estimating historical heat disturbance is shown, according to some embodiments. By estimating historical heat disturbance based on prior data, a model can be trained to predict current and future heat disturbances for use in performing model predictive control. The historical heat disturbance can be used to train a model for determining how heat disturbance affects a space (e.g., a zone in a building). In some embodiments, some and/or all steps of process1900are performed by controller212described with reference toFIG. 2.

Process1900is shown to include receiving training data for a building system (step1902), according to some embodiments. The training data for the building system may describe measurements of various environmental information taken at previous time steps and control inputs of the system, for example as described above with reference to training data generator408and training data database410. For example, the measurements may include information such as building temperature measurements, building humidity measurements, occupancy measurements, electrical usage measurements, time of day, day of week, external weather conditions, etc. As heat disturbance can originate from many sources, large amounts of training data may be necessary to properly estimate historical heat disturbance. In some embodiments, the training data received in step1902includes time-series data. Time-series data can further indicate how heat disturbance changes over time. For example, in the middle of a day when occupants are present in the building system, a higher heat disturbance due to thermal radiation of people may be present than in the middle of the night when fewer occupants are present. In some embodiments, the training data is refined using a saturation detection and removal process. System and methods for saturation detection and removal suitable for use to refine the training data are described in U.S. patent application Ser. No. 15/900,459, filed Feb. 20, 2018, incorporated by reference herein in its entirety. In some embodiments, step1902is performed by controller212.

Process1900is shown to include performing system identification to identify system matrices A, B, C, and D (step1904), according to some embodiments. Based on the system matrices A, B, C, and D, a system model determined via the system identification can be used when estimating heat disturbance. In some embodiments, the system identification performed in step1904to identify the system model is similar to and/or the same as the system identification process as described in U.S. patent application Ser. No. 16/240,028, filed Jan. 4, 2019, incorporated by reference herein in its entirety. By identifying the system matrices, a system model can be captured to be used in determining heat disturbance. In general, the system matrices describe dynamics of the system. In step1904, the system model identified can illustrate dynamics (e.g., thermal dynamics) of the building system based on the training data. In some embodiments, step1904is similar to and/or the same as process600described with reference toFIG. 6. In some embodiments, step1904is performed by controller212.

Process1900is shown to include augmenting the system matrices with a disturbance model having two or more dimensions (step1906), according to some embodiments. In general, an augmented system of the system identified in step1904is given by the following state space representation:

[x⁡(t).d⁡(t).]=[AcBd0Ad]⁡[x⁡(t)d⁡(t)]+[BcBd⁢d]⁢⁢u⁡(t);(Eq.⁢N)y⁡(t)=[Cc⁢⁢Cd]⁡[x⁡(t)d⁡(t)]+Dc⁢u⁡(t);(Eq.⁢O)
where Ad, Bd, Bdd, and Cdare matrices characterizing the disturbance model and the parameters Ac, Bc, Cc, and Dcare the matrices A, B, C, D identified in step1904. In some embodiments, the disturbance model is a parameterized model such that the parameterized model can generate accurate models of heat disturbance over repeated usages of the model given various parameters.

In some embodiments, the disturbance model is determined based on a known higher order disturbance model. The disturbance model is described in greater detail in Eq. T and Eq. U below with reference to step1908. If the disturbance model is determined based on a known higher order disturbance model, an augmented system of the system identified in step1904may have the following form:

where θ1through θ8are parameters that can be identified or set to a prospected value, and all other variables being the same as above. In general, the augmented system model can be used for determining historic heat disturbances. In some embodiments, step1906is performed by controller212.

Process1900is shown to include identifying and/or selecting parameters of the disturbance model and identifying a Kalman gain (step1908), according to some embodiments. In some embodiments, the parameters are determined using a known higher order disturbance model. Using the known higher order disturbance model, a rate of change in estimated historical heat disturbance values can be captured. Further, periodic heat disturbances can be represented using the known higher order disturbance model. For example, solar irradiance may follow a periodic schedule with a peak around noon each day. In some embodiments, the known higher order disturbance model is a second order disturbance model of an oscillator system having two states, d1and d2, where d1can estimate (or calculate) values of heat disturbance and d2can estimate (or calculate) a rate of change in values of the heat disturbance. In general, the oscillator system may have the following form:

[d·1d·2]=Ad⁡[d1d2]+Bd⁢d⁢u=[01-w2-2⁢γ]⁡[d1d2]+Bd⁢d⁡[Ts⁢pTo⁢a];(Eq.⁢T)
where w is a frequency tuning parameter, γ is a damping tuning parameter, Bddis a matrix mapping a forcing input, Tspis an indoor air temperature setpoint, and Toais an outdoor air temperature. For example, the values of w and y can be selected to provide a user-selected period or frequency for the oscillator system, for example a period of one day that reflects oscillations in solar irradiance as described above. As such, the tuning parameters can be set to

w=2⁢π2⁢4×3⁢6⁢0⁢0,γ=0,
and the matrix that maps the forcing input Bddcan be set w zeroes

(i.e.,Bd⁢d=[0000]).
In this example, a pure oscillator system with zero damping and a frequency corresponding to a one-day period is achieved.

In some embodiments, the disturbance model is determined via a system identification using the training data received in step1902. Similar to using the known higher order disturbance model, two disturbance states, d1and d2, represent values of heat disturbance and a rate of change in the values of the heat disturbance respectively. In general, the disturbance model determined via system identification has the following form:

Process1900is shown to include estimating historical heat disturbances using the disturbance model and the training data (step1910), according to some embodiments. To estimate the historical heat disturbance, the disturbance model including the augmented system matrices as well as the Kalman gain can use the training data as input. In general, accuracy of the estimated historical heat disturbance hinges on accuracy of the system identification in step1904and parameters of the disturbance model identified/selected in step1908. Accordingly, one result of process1900is a set of estimated historical heat disturbance values that includes a heat disturbance value for a plurality of time steps for which training data is available. In some embodiments, step1910is performed by controller212.

Referring now toFIG. 20, a process2000for predicting a forecasted heat disturbance {circumflex over (Q)}otherforecastand performing model predictive control based on the forecasted heat disturbance is shown, according to some embodiments. Process2000uses an autoregressive (AR) model to model a stochastic heat disturbance. In some embodiments, the AR model is used online (i.e., while the system is operating) to correct predictions of a deterministic heat disturbance by accounting for residuals (i.e., prediction errors), thereby predicting the stochastic heat disturbance of Qother. In some embodiments, some and/or all steps of process2000are performed by controller212.

Process2000is shown to include estimating historical heat disturbance (step2002), according to some embodiments. In some embodiments, step2002is accomplished by performing process1900described with reference toFIG. 19. Process2000may require the system matrices, the disturbance model, the Kalman gain, and the estimated historic heat disturbance as determined in process1900to estimate {circumflex over (Q)}otherforecast. In some embodiments, step2002is performed by controller212.

Process2000is shown to include receiving an augmented system model, a disturbance model, and estimated historical heat disturbances (step2004), according to some embodiments. In general, the disturbance model and the estimated historical heat disturbances are as determined in step2002. In some embodiments, step2004is performed by controller212.

Process2000is shown to include training a deterministic heat disturbance model for online predictions of the deterministic portion of heat disturbances (step2005), according to some embodiments. In some embodiments, the deterministic heat disturbance model uses a linear regression and data fitting and/or pattern recognition to estimate the deterministic heat disturbance value similar to and/or the same as a deterministic load as described in U.S. patent application Ser. No. 14/717,593 filed May 20, 2015, incorporated by reference herein in its entirety. The disturbance model can be trained using the training data of process1900and the estimated historical heat disturbance so that the disturbance model can properly estimate deterministic heat disturbance values based on environmental conditions, a time of day, a day of the week, etc. In some embodiments, step2005is performed by controller212.

Process2000is shown to include training an autoregressive (AR) model to predict residuals of the disturbance model (step2006), according to some embodiments. In some embodiments, the AR model is trained using the estimated historical heat disturbances gathered in step2002. In some embodiments, the AR model is a first order model that captures the stochastic part of the heat disturbance Qother. In some embodiments, the AR model is a higher order model. The AR model can be determined off-line and may not be part of a state space model used during a model predictive control process. In general, a first order AR model can have the following form:
ê(k+1)=a1e(k)  (Eq.V);
where
e(k)=Qothermeasured(k)−{circumflex over (Q)}otherdeterministic(k)
where k is a time step, ê(k+1) is a residual for a next time step, a1is a constant, Qothermeasured(k) is measurement of the heat disturbance for time step k, and {circumflex over (Q)}otherdeterministic(k) is the deterministic heat disturbance for time step k. As Qothermeasured(k) cannot be directly obtained, Qothermeasured(k)=d1(k) as estimated in Eq. P and Eq. R. As a result, a residual can be determined for a current time step and estimated for future time steps by the AR model. In some embodiments, the AR model trained to predict residuals is similar to and/or the same as an autoregressive model to predict residuals as described in U.S. patent application Ser. No. 14/717,593, filed May 20, 2015, incorporated by reference herein in its entirety. In some embodiments, step2006is performed by controller212.

Process2000is shown to include collecting current environmental data (step2008), according to some embodiments. The current environmental data can include, for example, temperature measurements, humidity measurements, external weather conditions, air quality, etc. The current environmental data can be collected by various sensors in a building, a weather service, and/or any device/system capable of collecting and communicating environmental data to controller212. For example, temperature sensor214and temperature sensor216may collect and communicate an indoor air temperature and an outdoor air temperature to controller212respectively. The current environmental data may also include information that can result in a heat disturbance such as, for example, occupancy in a building, current electrical consumption, time of day, etc. In some embodiments, step2008is performed by controller212and various devices/services capable of collecting and communicating environmental data (e.g., temperature sensor214).

Process2000is shown to include performing an online state estimation using the system matrices and the current environmental data (step2010), according to some embodiments. In step2010, the five state model described in Eq. P and Eq. R can be converted to discrete time. Further, a state estimation can be conducted using the Kalman gain identified for the model. For each time step k, a state vector of the following form can be calculated:

[Ti⁢a⁡(k)Tm⁡(k)I⁡(k)d1⁡(k)d2⁡(k)]
where all variables are the same as described in Eq. P and Eq. R for each time step k. In some embodiments, step2010is performed by controller212.

Process2000is shown to include estimating a deterministic heat disturbance value for a current time step using the disturbance model and the current environmental data (step2012), according to some embodiments. Specifically, the current environmental data can be applied to the disturbance model to obtain the estimation of the deterministic heat disturbance value. In some embodiments, step2012is performed by controller212.

Process2000is shown to include calculating a current residual based on the online state estimation and the estimated deterministic heat disturbance for the current time step (step2014), according to some embodiments. In general, the current residual can be calculated by the following equation:
e(k)=d1(k)−{circumflex over (Q)}otherdeterministic(k)  (Eq. W)
where e(k) is the current residual, d1(k) is a value of heat disturbance for the current time step determined based on the online state estimation performed in step2010, and {circumflex over (Q)}otherdeterministic(k) is the deterministic heat disturbance estimated in step2012. The current residual can illustrate the inaccuracy of {circumflex over (Q)}otherdeterministic(k) due to the stochastic heat disturbance. In general, a large value (negative or positive) of the current residual indicates {circumflex over (Q)}otherdeterministic(k) is inaccurate while values of the current residual close to 0 indicate {circumflex over (Q)}otherdeterministic(k) is more accurate. If the current residual is equal to 0, the stochastic heat disturbance may have no effect on heat disturbance. In some embodiments, step2014is performed by controller212.

Process2000is shown to include predicting a stochastic heat disturbance for subsequent time steps using the autoregressive (AR) model and the current residual (step2016), according to some embodiments. As the current residual can indicate a stochastic heat disturbance for a current time step, the current residual can be used as input to the AR model for predicting the stochastic heat disturbance for the subsequent time steps (i.e., the model trained at step2006). Based on a magnitude of the stochastic heat disturbance, the AR model can determine an accuracy of the estimated deterministic heat disturbance for the current time step. As the magnitude of the stochastic heat disturbance grows, additional correction may be necessary by the AR model to predict an accurate stochastic heat disturbance for the subsequent time steps. In some embodiments, the stochastic heat disturbance for a next time step is described as {circumflex over (Q)}otherstochastic(k+1) where k+1 is the next time step (i.e., the time step after time step k). Similarly, the stochastic heat disturbance for any subsequent time steps can be described as {circumflex over (Q)}otherstochastic(k+n) for some nthsubsequent time step. In some embodiments, step2016is performed by controller212.

Process2000is shown to include predicting a deterministic heat disturbance for the subsequent time steps using the disturbance model and the current environmental data (step2018), according to some embodiments. The disturbance model can be configured to estimate future deterministic heat disturbances based on current data. As such, the current environmental data can be used as input to the disturbance model such that the disturbance model can output the predicted deterministic heat disturbance for the subsequent time steps. When predicting the deterministic heat disturbance for the subsequent time steps, the disturbance model may not use the stochastic heat disturbance for the subsequent time steps predicted in step2016. As such, step2018may be similar to step2012for the subsequent time steps. In general, the deterministic heat disturbance for a next time step is described as {circumflex over (Q)}otherdeterministic(k+1) where k+1 is the next time step. Similarly, the deterministic heat disturbance for any subsequent time step can be described as {circumflex over (Q)}otherdeterministic(k+n) for some nthsubsequent time step. In some embodiments, step2018is performed by controller212.

Process2000is shown to include predicting a forecasted heat disturbance for the subsequent time steps as a sum of the deterministic heat disturbance and the stochastic heat disturbance (step2020), according to some embodiments. The forecasted heat disturbance {circumflex over (Q)}otherforecast(k+1) is an estimated value of total heat disturbance for the next time step k+1. In general {circumflex over (Q)}otherforecast(k+1) can be predicted by the following equation:
{circumflex over (Q)}otherforecast(k+1)={circumflex over (Q)}otherdeterministic(k+1)+{circumflex over (Q)}otherstochastic(k+1)  (Eq. X);
or as:
{circumflex over (Q)}otherforecast(k+1)={circumflex over (Q)}otherdeterministic(k+1)+ê(k+1)  (Eq. Y);
where all variables are as defined above. If {circumflex over (Q)}otherforecast(k+1) is a positive value, the predicted heat disturbance for the next time step may be adding additional heat into a space (e.g., building10). If {circumflex over (Q)}otherforecast(k+1) is a negative value, the predicted heat disturbance for the next time step may be taking heat out of the space. {circumflex over (Q)}otherforecastcan similarly be determined for other subsequent time steps based on {circumflex over (Q)}otherdeterministicand {circumflex over (Q)}otherstochasticfor each subsequent-time step. In some embodiments, step2020is performed by controller212.

Process2000is shown to include performing a model predictive control process to control building equipment using the forecasted heat disturbance and the identified model (step2022), according to some embodiments. By taking into account the forecasted heat disturbance, the model predictive control process can further optimize (e.g., reduce) costs related to operation of building equipment. For example, if the forecasted heat disturbance is positive and a building zone requires heating to maintain occupant comfort, the model predictive control process may determine that a heater is not required to be operated as the heat disturbance will increase a temperature of the building zone regardless. Without accounting for the forecasted heat disturbance, the model predictive control process may otherwise make control decisions that do not maintain occupant comfort and/or do not optimize costs. In some embodiments, step2022is performed by controller212.

Referring now toFIG. 21, a process2100for determining a deterministic heat disturbance and performing a model predictive control process based on the deterministic heat disturbance is shown, according to some embodiments. Process2100utilizes multistep system identification to determine a Kalman gain and a stochastic disturbance model that is part of a state space system. In process2100, the heat disturbance predictions from the deterministic heat disturbance model can be used as input for determining overall heat disturbance. Further, the Kalman gain and the stochastic disturbance model can be identified to account for inaccuracy in prediction of the deterministic heat disturbance due to inherent randomness, thereby allowing a stochastic heat disturbance to be calculated. In some embodiments, some and/or all steps of process2000are performed by controller212.

Process2100is shown to include estimating historical heat disturbance (step2102), according to some embodiments. In some embodiments, step2102is accomplished by performing process1900described with reference toFIG. 19. In some embodiments, step2102is performed by controller212.

Process2100is shown to include training a deterministic heat disturbance model for online predictions of the deterministic portion of heat disturbances (step2104), according to some embodiments. In some embodiments, the deterministic heat disturbance model uses a linear regression and data fitting and/or pattern recognition to estimate the deterministic heat disturbance value similar to and/or the same as a deterministic load as described in U.S. patent application Ser. No. 14/717,593 filed May 20, 2015, incorporated by reference herein in its entirety. The disturbance model can be trained using the training data of process1900and the estimated historical heat disturbance so that the disturbance model can properly estimate deterministic heat disturbance values based on environmental conditions. In some embodiments, step2104is performed by controller212.

Process2100is shown to include receiving a system model and augmenting the system model with a stochastic state space model and a Kalman gain (step2106), according to some embodiments. Specifically, the stochastic state space model used to augment the system model can be a new disturbance model Adethat captures the stochastic piece of the heat disturbance. In general, the augmented system model can have the following form:

Process2100is shown to include determining an identified model by performing a system identification to determine values of the Kalman gain and stochastic disturbance model parameter(s) based on the historical heat disturbances and the training data (step2108), according to some embodiments. Based on the historical heat disturbances, the Kalman gain can be determined such that the Kalman gain accounts for previous heat disturbances during previous time steps. By accounting for the historical heat disturbances, the Kalman gain (i.e., an adjustment for how inaccurate the deterministic heat disturbance is due to the stochastic heat disturbance) can be more precisely estimated to account for actual dynamics of the system. In step2108, values of Ade(i.e., values of the stochastic disturbance model parameter(s)) can also be determined as part of the performed system identification. In some embodiments, multi-step ahead prediction is utilized in identification of Adeand the Kalman gain for improved estimation of the stochastic heat disturbance. In some embodiments, step2108is performed by controller212.

Process2100is shown to include collecting an environmental condition forecast and/or a weather forecast (step2110), according to some embodiments. The environmental condition forecast can include current and/or predicted environmental conditions for future times. The environmental condition forecast can include, for example, current/predicted temperature measurements, current/predicted humidity measurements, current/predicted external weather conditions, current/predicted air quality, etc. The environmental condition forecast can be collected by various sensors in a building, a weather service, and/or any device/system capable of collecting and communicating environmental data to controller212. For example, temperature sensor214and temperature sensor216may collect and communicate an indoor air temperature and an outdoor air temperature to controller212respectively. The environmental condition forecast may also include information regarding factors that can result in a heat disturbance such as, for example, occupancy in a building, current electrical consumption, time of day, etc. The weather forecast collected in step2110can include predictions of external weather conditions at future times. The weather forecast can be obtained by, for example, requesting the weather forecast from an application programming interfaces (APIs) that provides weather forecasts to requesting services. The weather forecast can be utilized for estimating heat disturbances due to external conditions. In some embodiments, the environmental condition forecast and the weather forecast are part of a single forecast. In some embodiments, step2110is performed by controller212and various devices/services capable of collecting and communicating environmental data (e.g., temperature sensor214).

Process2100is shown to include predicting deterministic heat disturbances using the disturbance model along with the environmental condition forecast and/or the weather forecast (step2112), according to some embodiments. The deterministic heat disturbances may need to be estimated as the deterministic heat disturbance is used as input to the disturbance model for predicting the stochastic heat disturbances. In some embodiments, step2112is similar to and/or the same as step2012described with reference toFIG. 20. In general, the disturbance model can utilize the environmental condition forecast and/or the weather forecast as input and output a deterministic heat disturbance prediction. In some embodiments, the disturbance model utilizes the entirety of the environmental condition forecast and/or the weather forecast as input. In some embodiments, the disturbance model utilizes a portion of the environmental condition forecast and/or the weather forecast as input. In some embodiments, step2112is performed by controller212.

Process2100is shown to include performing a model predictive control process to control building equipment using the deterministic heat disturbance predictions and the identified model (step2114), according to some embodiments. Based on results of step2112, the deterministic heat disturbance predictions can be used as input to the identified model to predict the stochastic heat disturbances. The deterministic heat disturbance predictions and the stochastic heat disturbance predictions can be used to determine a total heat disturbance to be considered during the model predictive control process. By taking into account the total heat disturbance, the model predictive control process can further optimize (e.g., reduce) costs related to operation of building equipment operable to heat or cool a building. For example, if the total heat disturbance is positive and a building zone requires heating to maintain occupant comfort, the model predictive control process may determine that a heater is not required to be operated as the heat disturbance will increase a temperature of the building zone regardless. Without accounting for the total heat disturbance, the model predictive control process may otherwise make control decisions that do not maintain occupant comfort and/or do not optimize costs. In some embodiments, step2114is performed by controller212.

Heat Disturbance Estimation and Prediction

Overview

Referring generally toFIGS. 22-30, methods for estimating a heat disturbance {dot over (Q)}otherfor use in performing a control process (e.g., MPC) is shown, according to some embodiments. As described above, an inaccurate disturbance model, quantized sensor measurements, and/or high frequency input characteristics can result in inaccurate heat disturbance estimations if not accounted for. As such,FIGS. 22-30illustrate how these issues can be addressed to ensure accurate heat disturbance estimations can be obtained for use in control processes such as MPC.

Processes for Estimating Heat Disturbances

Referring now toFIG. 22, a process2200to training a deterministic model of heat disturbance is shown, according to some embodiments. In some embodiments, some and/or all steps of process2200are performed by controller212as described with reference toFIG. 2.

Process2200is shown to include obtaining zone group data describing thermal dynamics of a zone group (step2202). To accurately train a deterministic model of heat disturbance, a minimum amount of data can be required. For example, step2202may include gathering a minimum of two weeks' worth of training data. During data gathering, step2202may include ensuring that specific data is obtained for use in training the deterministic model. For example, step2202can include gathering data such as an outdoor relative humidity (e.g., as a percentage), a dry bulb weather (e.g., in degrees Kelvin), a temperature of a zone (e.g., in degrees Kelvin), measurements of {dot over (Q)}HVAC(e.g., in kilowatts), a time stamp associated with each data point, and a sampling rate of the data (i.e., how often data is collected).

In step2202, the sampling rate of data can be set to an appropriate value to ensure a sufficient amount of data is gathered over a data gathering period (also referred to herein as a learning period). For example, the sampling rate can be set to sample data every five minutes. The zone group data gathered over the learning period can be sampled at a constant rate (e.g., every five minutes) to ensure that the data does not include variations due to differences in the sampling rate. Further, to ensure data representative of system dynamics is gathered, the learning period may be set to start and end at a same time on different days (e.g., from midnight to midnight). In some embodiments, step2202is performed by controller212.

Process2200is shown to include performing a linear regression to estimate a thermal resistance between indoor air and outdoor air (step2204). The thermal resistance between the indoor air and the outdoor air can be defined as

1Ro⁢i.
The linear regression performed in step2204can help in estimating heat transfer effects from outside air temperature. In performing the linear regression, {dot over (Q)}othercan be modeled as a constant offset. Further, the following zone temperature steady state equation can be utilized in performing the linear regression:

0=1Ro⁢i⁢(To⁢a-Ti⁢a)+Q.H⁢V⁢A⁢C+Q.o⁢t⁢h⁢e⁢r
In this example of the zone temperature steady state equation, effects of building mass dynamics are not included.

To perform the linear regression, values of {dot over (Q)}HVACcan be extracted from the zone group data to create a first column vector containing the values. The first column vector can have the following form:

[Q.H⁢V⁢A⁢C1Q.H⁢V⁢A⁢C2⋮Q.H⁢V⁢A⁢Cn]
where {dot over (Q)}HVACxindicates a value of {dot over (Q)}HVACat a time x and where n indicates a last recorded value of {dot over (Q)}HVAC. Likewise, a second column vector containing values of (Toa−Tia) of the zone temperature steady state equation can be generated to have the following form:

[To⁢a1-Ti⁢a1To⁢a2-Ti⁢a2⋮To⁢an-Ti⁢an]
where Toax−Tiaxindicates a value of Toa−Tiaat time x.

To ensure accurate relationships between operation of HVAC equipment and the heat disturbance are maintained, any data points in the first column vector where {dot over (Q)}HVACequals zero can be removed as well as any corresponding data points from the second column vector (i.e., data points with a same timestamp as those removed from the first column vector). Removal of said data points can result in a modified first column vector of the following form:

[Q.H⁢V⁢A⁢C1Q.H⁢V⁢A⁢C2⋮Q.H⁢V⁢A⁢Cm]
where m indicates a new last value of {dot over (Q)}HVACsuch that m≤n. It should be understood that {dot over (Q)}HVAC1and/or {dot over (Q)}HVAC2of the first column vector may or may not be the same as {dot over (Q)}HVAC1and/or {dot over (Q)}HVAC2of the modified first column vector depending on whether {dot over (Q)}HVAC1and/or {dot over (Q)}HVAC2of the first column vector are equal to zero. Likewise, a modified second column vector can be determined to have the form:

[To⁢a1-Ti⁢a1To⁢a2-Ti⁢a2⋮To⁢am-Ti⁢am]
where values of Toax−Tiaxfrom the second column vector are removed if {dot over (Q)}H{dot over (Q)}HVACx=0 for some time x.

Based on the modified first column vector without data points where {dot over (Q)}HVACequals zero, a third column vector can be generated. The third column vector can be generated such that the third column vector is of the same length as the modified first column vector and includes only the value 1. In other terms, the third column vector can be represented as:

[11⋮1]
where a length of the third column vector is the same as the length of the modified first column vector corresponding to the modified first column vector.

Further, a two column matrix A can be generated where A is a combination of the modified second column vector and the third column vector. In particular, a first column of A can include values of (Toa−Tia) excluding those removed that have the same timestamp as data points where {dot over (Q)}HVACequals zero. Likewise, a second column of A can be of the same length as the third column vector and include all 1's. As such, A can have the following form:

A=[To⁢a1-Ti⁢a111To⁢a2-Ti⁢a212⋮⋮To⁢am-Ti⁢am1m]
Further, in performing the linear regression, a vector b can be generated such that b=−{dot over (Q)}HVAC. As such, b can have the following form:

Based on A and b, an equation Ax=b can be solved by using a pseudoinverse computed as the following:
x=(ATA)−1ATb
where the T superscript indicates a transposed matrix. In this way, a first element of the x vector can define

1Ro⁢i.
In some embodiments, step2204is performed by controller212.

Process2200is shown to include computing an initial profile of heat disturbance base on the thermal resistance between the indoor air and the outdoor air (step2206). As {dot over (Q)}otheris modeled as a constant offset in the linear regression performed in step2204,

1Ro⁢i
may be an approximation. Based on the approximate value of

1Ro⁢i,
the initial profile of {dot over (Q)}other(referred to as {dot over (Q)}otherBase)can be calculated using the zone temperature steady state equations not including the building mass temperature effects. In general, {dot over (Q)}otherBasecan be calculated by the following equation:

Q.o⁢t⁢h⁢e⁢rB⁢a⁢s⁢e=-1Ro⁢i⁢(To⁢a-Ti⁢a)-Q.H⁢V⁢A⁢C
In some embodiments, step2206is performed by controller212.

Process2200is shown to include filtering the initial profile of heat disturbance (step2208). If generated {dot over (Q)}otherBasemay include various noise and spikes that can lead to inaccuracies of predictions and/or other determinations made based on {dot over (Q)}otherBase. As such, to ensure {dot over (Q)}otherBaseincludes accurate estimations of a magnitude of {dot over (Q)}otherfor each day {dot over (Q)}otherBasecan be filtered to reduce an amount of noise and spikes. To filter {dot over (Q)}otherBase, {dot over (Q)}otherBasecan be passed through an anti-spike filter to remove spikes in {dot over (Q)}otherestimations. Examples of anti-spike filters can include median filters, clamp filters, etc. {dot over (Q)}otherBasecan also be passed through a smoothing filter to reduce an impact of noise on estimations of {dot over (Q)}otherExamples of smoothing filters can include an exponentially-weighted moving average filter, a Savitzky-Golay filter, etc. In some embodiments, {dot over (Q)}otherBaseis not passed through any filters if a determination is made that {dot over (Q)}otherBasedoes not include excessive noise and/or spikes. In some embodiments, {dot over (Q)}otherBaseis only passed through one of the anti-spike filter or the smoothing filter. In some embodiments, {dot over (Q)}otherBaseis passed through additional filters and/or is passed through the anti-spike filter and/or the smoothing filter multiple times. In other words, step2208can include passing {dot over (Q)}otherBasethrough any necessary filters to get {dot over (Q)}otherBaseto include good (i.e., accurate) estimates of {dot over (Q)}other. In some embodiments, step2208is performed by controller212.

Process2200is shown to include training a deterministic model of heat disturbance based on the filtered initial profile of heat disturbance (step2210). In some embodiments, {dot over (Q)}otherBaseis used to train an approximate heat disturbance profile {dot over (Q)}otherApprox. If {dot over (Q)}otherBaseis used to train {dot over (Q)}otherApprox, {dot over (Q)}otherApproxcan be used to train the deterministic model. Generation of {dot over (Q)}otherApproxis described in greater detail below with reference toFIG. 23. As such, step2210may include performing some and/or all of a process2300described with reference toFIG. 23.

Alternatively, the deterministic model can be trained directly based on filtered {dot over (Q)}otherBase. If {dot over (Q)}otherBaseis used to directly train the deterministic model, processing efficiency of performing process2200may increase as fewer approximations and other calculations related to the heat disturbance may be needed as compared to generating {dot over (Q)}otherApproxin order to train the deterministic model. However, {dot over (Q)}otherApproxmay provide a more accurate deterministic model, thereby resulting in more accurate heat disturbance predictions and better control decisions that more precisely operate building equipment to account for the heat disturbance. As such, step2210may include determining whether to train the deterministic model directly based on {dot over (Q)}otherBaseor if to first generate {dot over (Q)}otherApproxand use {dot over (Q)}otherApproxto train the deterministic model. In some embodiments, step2210is performed by controller212.

Referring now toFIG. 23, process2300for training a deterministic model of heat disturbance based on an initial profile of the heat disturbance is shown, according to some embodiments. In some embodiments, process2300is performed based on a determination in step2210of process2200as described with reference toFIG. 22that the deterministic model is to be trained based on an approximate profile of {dot over (Q)}other, referred to as {dot over (Q)}otherApprox. In some embodiments, some and/or all steps of process2300are performed by controller212.

Process2300is shown to include obtaining an initial profile of heat disturbance (step2302). As described above with reference toFIG. 22, the initial profile of heat disturbance can be referred to {dot over (Q)}otherBase. In some embodiments, step2302is completed by performing some and/or all of process2200as described with reference toFIG. 22. In some embodiments, step2302is performed by controller212.

Process2300is shown to include fitting each day's worth of the initial profile to a predefined function form to obtain an approximate profile of the heat disturbance (step2304). {dot over (Q)}otherBasecan include accurate or relatively accurate estimates of a magnitude of {dot over (Q)}otherin each day. However, {dot over (Q)}otherBasemay not include a correct shape of {dot over (Q)}otherif a high frequency {dot over (Q)}HVACinput signal is present. The shape of {dot over (Q)}othercan be defined by how estimates of {dot over (Q)}otherdetermined based on {dot over (Q)}otherBasechange over each day. The high frequency {dot over (Q)}HVACinput signal can result in difficulty in distinguishing high frequency noise from an actual signal, thereby leading to a loss of signal contents even if {dot over (Q)}otherBaseis filtered. Advantageously, effects of the high frequency {dot over (Q)}HVACinput signal can be mitigated by fitting each day's worth of {dot over (Q)}otherBasedata to a predefined function form. The total of each fitted profile can be used to define {dot over (Q)}otherApprox.

The function forms used to fit each day's worth of {dot over (Q)}otherBasedata can be chosen based on an expected heat load profile of an associated building (e.g., building10). Different types of buildings may have varying expected heat load profiles based on a purpose of a building. For example, a residential building may have a different expected heat load profile than a commercial building as times when occupants are present may be different, usage of electronics may differ, etc. To obtain {dot over (Q)}otherApprox, various predefined function forms can be used. In particular, the function forms can include, for example, a Gaussian function, a sine function, and a user defined function, all described in greater detail below. In some embodiments, a different function form is used rather than the Gaussian function, the sine function, or the user defined function to obtain {dot over (Q)}otherApprox. The Gaussian function, the sine function, and the user defined function are provided purely for sake of example.

The Gaussian function can be defined by a target daily fitting function shown as:

f⁡(x)=a·e-(x-b)22⁢⁢c2+d
where a is an amplitude factor where a∈[−∞, ∞], b is a mean or horizontal shift where b∈[0, samples in one day], c is a standard deviation where c∈[−∞, ∞], d is a vertical shift where d∈[−∞, ∞], and x is a sample index defined as a consecutive integer array starting from one to a total number of samples per day. In some embodiments, c can change shape to be wide or narrow as needed.

For each day, a separate f (x) fit can be found by determining optimal values of a, b, c, and d that minimize a two-norm error between f(x) and {dot over (Q)}otherBasefor a particular day. In general, a two-normal measure can be defined as the square root of the squares of each data point. For example, the two-norm of {dot over (Q)}otherBasecan be shown as:

Q.otherB⁢a⁢s⁢ex
is an estimation of {dot over (Q)}otherbased on {dot over (Q)}otherBasefor a time step x and n is a total number of time steps. Based on each f(x), a total of all f(x) can be determined to make up {dot over (Q)}otherApprox. Possible shapes of the Gaussian function are described in greater detail below with reference toFIGS. 25A-C.

As for the sine function, the target daily fitting function can be shown as:
f(x)=a·sin(b·(x−c))+d
where a is an amplitude factor where a∈[−∞, ∞], b is a period where b∈

[0.25⁢πsamples⁢⁢in⁢⁢one⁢⁢day,2⁢πsamples⁢⁢in⁢⁢one⁢⁢day],
c is a phase shift where c∈[0, samples in one clay], d is a vertical shift where d∈[−∞, ∞], and x is a sample index defined as a consecutive integer array starting from one to a total number of samples per day.

Similar to the Gaussian function, a separate f (x) can be found for each day by determining optimal values of a, b, c, and d that minimize the two-norm error between f(x) and {dot over (Q)}otherBasefor a particular day. Based on each f(x), a total of all f(x) can be determined to make up {dot over (Q)}otherApprox.

As for the user defined function, a user can set their own function to generate {dot over (Q)}otherApprox. An example of a target daily fitting function can be determined by generating a normalized {dot over (Q)}otherprofile that has an expected shape with a peak around noon. An example of the normalized{dot over (Q)}otherprofile is shown below in graph2600as described with reference toFIG. 26. Based on (filtered) {dot over (Q)}otherBaseprofile obtained in step2302, the normalized {dot over (Q)}otherprofile can be set such that a peak of the normalized {dot over (Q)}otherprofile can be set equal to a maximum value of {dot over (Q)}otherBasefor a particular day. In other words, the maximum value of {dot over (Q)}otherBasecan be multiplied by the normalized{dot over (Q)}otherprofile for each day. Said multiplication can result in a scaled profile (i.e., {dot over (Q)}otherApprox) that has a unique scale for each day.

As mentioned above, the Gaussian function, the sine function, and the user defined function are given purely for sake of example. The present disclosure contemplates using various other functions to generate {dot over (Q)}otherApprox. In some embodiments, step2304is performed by controller212.

Process2300is shown to include training a deterministic model of the heat disturbance based on the approximate profile, a relative humidity, and an outdoor air temperature (step2306). The relative humidity and the outdoor air temperature can be obtain through various sensors, services, etc. For example, the relative humidity and/or the outdoor air temperature can be obtained by outdoor air temperature sensor216. As another example, the relative humidity and/or the outdoor air temperature can be obtained from a weather service that provides data regarding various environmental conditions. The relative humidity and the outdoor air temperature can be particularly helpful in training the deterministic model as the relative humidity and the outdoor air temperature can potentially impact the heat disturbance affecting a building. In some embodiments, other various environmental conditions (e.g., cloud cover) are used to train the deterministic model. In some embodiments, the deterministic heat disturbance model is obtain by performing a linear regression and data fitting process as described in U.S. patent application Ser. No. 14/717,593 filed May 20, 2015, incorporated by reference herein in its entirety. In some embodiments, step2306is performed by controller212.

Referring now toFIG. 24, a process2400for identifying a state space system for use in a model predictive control (MPC) process is shown, according to some embodiments. Process2400can be useful for obtaining accurate predictions of heat disturbance based on both a deterministic model and a stochastic model. Based on said predictions, a control process such as MPC can determine how to operate building equipment (e.g., HVAC equipment210) over the duration of an optimization period in order to optimize (e.g., reduce) costs and/or maintain occupant comfort. By accounting for the heat disturbance, control decisions related to operation of the building equipment can result in desired changes to temperature and/or other environmental conditions in a building (e.g., building10). In some embodiments, some and/or all steps of process2400are performed by controller212.

Process2400is shown to include obtaining a deterministic model of heat disturbance (step2402). In some embodiments, step2402includes performing process2300described with reference toFIG. 23and/or process2400described with reference toFIG. 24. As such, the deterministic model obtained in step2402may be trained based on either {dot over (Q)}otherBaseor {dot over (Q)}otherApprox. In some embodiments, step2402is performed by controller212.

Process2400is shown to include predicting deterministic portion profiles of the heat disturbance based on the deterministic model (step2404). The deterministic portion profiles can define predictions of the deterministic portion of {dot over (Q)}otherover a time period. As the deterministic model is trained to represent how {dot over (Q)}otherchanges over time in a building based on varying conditions (e.g., relative humidity, outdoor air temperature, etc.), conditions can be supplied as input to the model to generate the deterministic portion profiles. In some embodiments, step2404is performed by controller212.

Process2400is shown to include performing a system identification process to identify a state space system of a building thermal model, a Kalman gain, a stochastic heat disturbance model, and a scaling parameter of the heat disturbance deterministic portion profiles (step2406). The building thermal model can describe various thermal dynamics of a building (e.g., building10). In some embodiments, the building thermal model identified in the system identification has the following form:

Process2400is shown to include augmenting the state space system with a stochastic model of the heat disturbance and a Kalman gain identified by multi-step system identification (step2408). In some embodiments, the stochastic model and/or the Kalman gain are identified in the system identification process performed in step2406. In some embodiments, the stochastic model and/or the Kalman gain are identified in a separate multi-step system identification process performed in step2408. The augmented state space model can be reflected by the following:

Process2400is shown to include performing a model predictive control (MPC) process based on the deterministic model, the stochastic model, and the building thermal model to determine an optimal control strategy (step2410). Effectively, by applying the state space model identified above, the heat disturbance can be determined. Based on the heat disturbance, MPC can determine how building equipment should be operated to optimize (e.g., reduce) costs and/or maintain occupant comfort. For example, if a space in a building is currently too hot for occupants and a heat disturbance due to solar radiation is present, MPC may determine that additional operation of building equipment is required to cool the space and offset effects of the heat disturbance. In general, by utilizing each model, MPC can generate modified control decisions that obtain desired adjustments to temperature of a space given heat disturbances affecting the space. The building equipment can be thereby controlled in accordance with an output of processes2200, process2300, and/or process2400. In some embodiments, step2410is performed by controller212.

Example Target Daily Fitting Functions

Referring generally toFIGS. 25A-26, example fitting functions that can be utilized by controller212to generate an approximate heat disturbance profile {dot over (Q)}otherApproxare shown, according to some embodiments. In some embodiments, some and/or all of the functions ofFIGS. 25A-26are utilized in step2304of process2300as described with reference toFIG. 23.

Referring now toFIG. 25A, a graph2500of an estimated heat disturbance profile of a large zone group of a commercial building is shown, according to some embodiments. Graph2500is shown to include a series2502. Series2502can illustrate how values a normalized {dot over (Q)}otherprofile are expected to change over time for the large zone group of the commercial building. Series2502is shown to reach a maximum at approximately 12 hours (i.e., noon) of a day. As graph2500is relative to the commercial building, the heat disturbance can be expected to be large towards the middle of the day when most employees are expected to be working, solar radiation is near a maximum, etc. For example, since employees are typically expected to be at work near the middle of the day, a heat disturbance due to body heat and electronic equipment used by the occupants (e.g., computers, phones, etc.) may expected to be highest around noon.

Referring now toFIG. 25B, a graph2504of an estimated heat disturbance profile of a small zone group of a commercial building is shown, according to some embodiments. In some embodiments, graph2504is similar to graph2500as described with reference toFIG. 25A. Graph2504is shown to include a series2506illustrating values of a normalized {dot over (Q)}otherprofile over a day for the small zone group. A maximum value of series2506is shown to be lower than a maximum value of series2502of graph2500. Similarly, values of series2506are shown to change at a more gradual pace as compared to values of series2502. As graph2504illustrates the small zone group, fewer effects of heat disturbance may impact the small zone group as compared to those impacting the large zone group illustrated by graph2500. Similarly, as series2506illustrates normalized values, the heat disturbance may change more gradually as a maximum and a minimum heat disturbance are closer than as compared to the large zone group of graph2500.

Referring now toFIG. 25C, a graph2508of an estimated heat disturbance profile of a small zone group of a residential building is shown, according to some embodiments. In some embodiments, graph2508is similar to graph2500and/or graph2504ofFIGS. 25A and 25Brespectively. Graph2508is shown to include a series2510. Series2510can illustrate values of a normalized {dot over (Q)}otherprofile over a day for the small zone group of the residential building. Unlike series2502or2506, series2510is shown to have a minimum around noon rather than a maximum. As series2510represents a residential building rather than a commercial building, the expected profile of heat disturbance may decrease towards the middle of the day when occupants are expected to be gone (e.g., at work). Further, as occupants may expected to be away, the expected heat disturbance due to operation of electronic equipment may decrease as well. Also, as series2510applies to a residential building, expected variations in the heat disturbance are shown to be less than both series2502and2506due to, for example, fewer people occupying a residential building as compared to a commercial building.

Referring now toFIG. 26, a graph2600of a user defined target daily fitting function used for fitting an initial profile of heat disturbance is shown, according to some embodiments. Graph2600is shown to include a series2602. Similar to series2502, series2506, and series2510ofFIGS. 25A-C, series2602can be used to generate an approximate heat disturbance profile {dot over (Q)}otherApprox. In particular, for each day's worth of data, the maximum value of {dot over (Q)}otherBasecan be multiplied by series2602to generate {dot over (Q)}otherApprox. As series2602can be determined by a user, series2602may include various adjustments that reflect a building managed by the occupant. For example, as shown by series2602, series2602has a maximum around 12 hours (i.e., noon), but also increases temporarily around 14 hours (i.e., 2:00 p.m.). The temporary increase around 14 hours may be defined by the user if, for example, the user expects a meeting in a zone defined by series2602to occur around 2:00 p.m. daily. In this way, {dot over (Q)}otherApproxcan be obtained such that {dot over (Q)}otherApproxaccounts for an additional heat disturbance around 2:00 p.m. due to the meeting. Of course, series2602is shown purely for sake of example. A user can define the normalized {dot over (Q)}otherprofile to follow any various shape.

Experimental Results of Applying Heat Disturbance Estimations and Prediction to MPC

To illustrate the advantages of the processes described above, the following experimental results are included and shown inFIGS. 27-30. Particularly, the experimental results shown below are based on a two week data collection period.

Referring now toFIG. 27, a graph2700of a regression performed on data representing a correspondence between heat transfer from HVAC equipment and a temperature difference is shown, according to an example experiment. Graph2700is shown to include a regression line2702. Regression line2702can be calculated by performing a linear regression based on various data points indicating a correspondence between a measurements of {dot over (Q)}HVACin kilowatts and a difference between a zone temperature and an outside air temperature in degrees Celsius represented by the various circles in graph2700. In the example experiment illustrated inFIG. 27, a slope of regression line2702is equal to 0.824, thereby indicating

Based on

1Ro⁢i,
{dot over (Q)}otherBasecan be calculated by the following:

{dot over (Q)}otherBasecan be passed through an anti-spike filter. In this example experiment, {dot over (Q)}otherBasecan be particularly passed through a median filter with a horizon looking forward and backward two samples as given by the following equation:
y(t)=median(u(t+2)+u(t+1)+u(t)+u(t−1)+u(t−2))

{dot over (Q)}otherBasecan be also passed through a smoothing filter. In the example experiment, {dot over (Q)}otherBasecan be particularly passed through a Savitzky-Golay filter with a first order polynomial and a horizon looking forward and backward 3 samples. In this way, a {dot over (Q)}otherBaseprofile without excessive noise and/or spikes can be obtained.

Referring now toFIG. 28, a graph2800of an approximate heat disturbance profile over time is shown, according to an example experiment. Graph2800is shown to include a series2802illustrating values of {dot over (Q)}otherApproxover time. {dot over (Q)}otherApproxcan be calculated based on the {dot over (Q)}otherBaseprofile determined above with reference toFIG. 27by taking a maximum value of the filtered {dot over (Q)}otherBasefor each day and scaling a normalized {dot over (Q)}otherprofile to have peaks of the normalized {dot over (Q)}otherprofile equal to the maximum value for each day. Series2802illustrates two weeks of {dot over (Q)}otherApproxprofile that can be used in conjunction with relative humidity and outside air temperature measurements to obtain a deterministic model of heat disturbance using linear regression and data fitting as described above with reference toFIG. 22.

Referring now toFIG. 29, a graph2900of deterministic heat disturbance predictions scaled by a scaling parameter is shown, according to some embodiments. Graph2900is shown to include a series2902illustrating deterministic predictions of the heat disturbance over a four day period determined based on the deterministic model described above with reference toFIG. 28. However, as {dot over (Q)}otherBasemay be obtained using steady state assumptions, series2902may include a correct shape of heat disturbances over time, but may not be scaled correctly. As such, graph2900is also shown to include a series2904and a series2906. Series2904illustrates a situation in which series2902is scaled up to increase values of the deterministic predictions indicated by series2902. Alternatively, series2906illustrates a situation in which series2902is scaled down to reduce values of the deterministic predictions.

To determine a proper scaling factor by which to scale series2902, a system identification process can be performed to identify the scaling factor and a building thermal model. In the context of the example experiment illustrated byFIGS. 27-30, a scaling factor of 0.637 can be calculated to be an optimal scaling value by which to scale values of series2902. As such, series2906can illustrate the scaled values of series2902. It should be noted that series2902-2904may not be to scale.

Referring now toFIG. 30, a graph3000, a graph3002, a graph3004, and a graph3006illustrating results of performing MPC to determine an optimal heating profile that maximizes cost savings while maintaining temperatures comfortable to occupants are shown, according to some embodiments. Based on the building thermal model and scaling parameter as described above with reference toFIG. 29, the building thermal model can be augmented with a stochastic disturbance model and a Kalman gain which can be identified through a multi-step system identification. In regards to the example experiment described inFIGS. 27-30, the multi-step horizon can be set to 72 steps (i.e., 6 hours for a sampling rate of 5 minutes). In this case, the stochastic disturbance model can be calculated to be Ade=−θ6=−0.974.

By using the deterministic, stochastic, and building thermal models, MPC can be performed to generate control decisions that control building equipment in accordance with the estimated heat disturbance. In MPC, the control decisions can be generated to optimize (e.g., reduce) costs and/or maintain occupant comfort. In the context of the example experiment described inFIGS. 27-30, a comfortable range can be set such that a temperature of a space (e.g., a zone of a building) is to stay between 79° F. and 84° F.

Graph3000is shown to include a series3008illustrating outside air temperature over time. Tracking the outside air temperature may be crucial for accurately estimating the heat load disturbance. As the outside air temperature increases, a larger heat disturbance due to the outside air and/or solar irradiance may affect the space, thereby impacting how building equipment should be operated.

Graph3002is shown to include a series3010, a series3012, and a series3014. Series3010illustrates a predicted building mass temperature of a building associated with the experiment. Series3012illustrates a predicted zone temperature for a zone of the building over time. Series3014illustrates a temperature setpoint generated based on MPC decisions.

Graph3004is shown to include a series3016illustrating an electricity rate profile over time. Series3016is shown to include two time periods at which a cost of electricity (as shown in dollars per kilowatt) is higher. The two time periods can illustrate peak consumption times at which an electricity supplier charges more for electricity. In order to generate cost-effective decisions, MPC may determine values of decision variables that reduce an amount of electricity consumed during the peak consumption times. For example, during the first peak consumption time, the MPC temperature setpoint as shown by series3014is shown to decrease to reduce an amount of electricity consumed by HVAC equipment.

Graph3006is shown to include a series3018illustrating power consumption (in kilowatts) of HVAC equipment over time. Similar to series3014, series3018is shown to decrease at the start of the first peak consumption time in order to reduce an amount of electricity used during the peak consumption time. By pre-heating the zone before the peak consumption time, the predicted zone temperature shown by series3012is maintained above 78° F. (i.e., within comfortable bounds).

The control decisions shownFIG. 30can be determined based on the deterministic model, the stochastic model, and the building thermal model determined from the data/calculations described above with reference toFIGS. 27-29. In this way, MPC can determine control decisions for the building/zone that optimize (e.g., reduce) costs and maintain occupant comfort even with heat disturbances affecting the building/zone.

Configuration of Exemplary Embodiments