## Abstract

An autonomous adaptive model predictive control (MPC) architecture is presented for control of heating, ventilation, and air condition (HVAC) systems to maintain indoor temperature while reducing energy use. Although equipment use and occupant changes with time, existing MPC methods are not capable of automatically relearning models and computing control decisions reliably for extended periods without intervention from a human expert. We seek to address this weakness. Two major features are embedded in the proposed architecture to enable autonomy: (i) a system identification algorithm from our prior work that periodically re-learns building dynamics and unmeasured internal heat loads from data without requiring re-tuning by experts. The estimated model is guaranteed to be stable and has desirable physical properties irrespective of the data; (ii) an MPC planner with a convex approximation of the original nonconvex problem. The planner uses a descent and convergent method, with the underlying optimization problem being feasible and convex. A yearlong simulation with a realistic plant shows that both of the features of the proposed architecture—periodic model and disturbance update and convexification of the planning problem—are essential to get performance improvement over a commonly used baseline controller. Without these features, long-term energy savings from MPC can be small while with them, the savings from MPC become substantial.

## 1 Introduction

Heating, ventilation, and air conditioning (HVAC) systems are responsible for approximately 40% of the total energy consumption of buildings in USA [1]. It has been recognized by many researchers that instead of the traditional rule-based control systems, an optimization-based controller—especially model predictive control (MPC)—is a highly promising approach to reduce energy use; see, for instance, the review paper [2].

In spite of extensive studies and even successful demonstration projects, e.g., Refs. [35], MPC has not been widely adopted in practice. The bottlenecks—which have been discussed extensively as well—can be summarized into lack of autonomy of existing control architectures that use MPC. By autonomous MPC, we mean an MPC scheme capable of reliably computing high-quality control decisions at all times without the need for human intervention. A building’s and its equipment’s behavior are quite complex and uncertain, so the models needed by MPC need to be learned from data. Since the building’s behavior also changes with time—albeit slowly—the models need to be updated over time. The overall architecture thus needs to be adaptive.

Although there is an extensive literature on identification of HVAC system models from data, the vast majority of the existing methods cannot be used for autonomous adaptation. These algorithms fit model parameters by solving a non-convex optimization problem, e.g., Refs. [69]. Depending on the type and quality of data used, they require re-tuning of hyper-parameters by a human expert. Clearly such an approach cannot lead to an autonomous control system. Another issue is that although the unmeasurable internal heat gains from occupants are substantial, most identification methods still ignore them which can lead to poor model quality. Works on model identification in the presence of large unknown disturbance in a principled manner are limited [810].

The planning problem that MPC solves at every decision instant to compute control commands should be feasible and convex. With a nonconvex problem, the planner can fail to converge to a local minimum within the allowed computation time. Infeasibility has the same effect. In either case, a rule-based controller must be used as back up when the non-convex planner cannot provide a control command. Switching between controllers can cause poor performance. The MPC planning problem is usually non-convex due to bilinearities in models and cost functions [1115]. Most works on HVAC MPC ignore the issue of reliability of the decisions computed by a non-convex planner, especially over long periods of operation.

In this paper, we propose an adaptive MPC architecture for HVAC systems, shown in Fig. 1, that is capable of operating autonomously for long periods of time without intervention of a human expert.

Fig. 1
Fig. 1
Close modal

The “ID + prediction” block uses an algorithm proposed in our prior work [10] to identify the plant model and the unmeasured internal disturbance from easily measured input-output data. This algorithm involves solving an optimization problem that is always feasible and convex, and the model it identifies ($M^$) is guaranteed to be stable and possess properties that are consistent with properties of a building HVAC system. The algorithm has one hyper parameter that needs to be tuned only once. In short, the identification algorithm does not need any human intervention when new data is fed into it periodically, say, every week, to update the model. The past disturbance identified by the algorithm is used to forecast the future disturbance ($w¯^$) that is in turn used by the MPC planner.

The “MPC planner” block of the system uses the model and disturbance forecasts to decide control commands so as to maintain indoor climate while reducing energy use. We provide a convex approximation algorithm to approximately solve the nominal non-convex planning problem. We show that the algorithm is a feasible, descent, and convergent algorithm. Thus, the MPC planner block can compute decisions autonomously without human expert. We also show that among the many ways of convexifying these types of non-convex problems, the proposed approach is the only one applicable to our specific problem structure.

The proposed convex planner and the associated analysis is the first novel contribution of the paper. The second contribution is the performance assessment of the closed-loop system for a yearlong period. Numerical results show that the proposed MPC scheme is not only more energy efficient and better at indoor climate control than a conventional baseline controller. More importantly, these simulations show that both features of the proposed design—periodic update of the model and disturbance and convexified MPC planner—are necessary to get the performance improvement over the baseline controller. These discoveries are made possible only due to the long duration for which simulations are conducted. While that is perhaps not surprising for the role of periodic model update, the discovery on the role of convexification also required the yearlong simulation. In particular, the non-convex controller was seen to perform as well as the convex one in all but a few rare instances. In these few rare cases, however, the performance of the non-convex planner was catastrophic.

This study makes three additional contributions over the preliminary version [16]. (i) We provide new analysis regarding the appropriate convexification method for the MPC planner, while Ref. [16] did not have any such results. (ii) This study includes closed-loop simulations for a yearlong period while the preliminary version only included 3 weeks of simulation. (iii) Comparison with three additional architectures, each obtained by removing model update or convexification or both, is provided. These comparisons reveal which of these features of the proposal are useful or necessary (or not).

The rest of this paper is organized as follows. Section 1.1 provides a review of relevant work on HVAC MPC. Section 2 describes the HVAC control problem. Sections 3 and 4 present the components of “MPC planner” and “ID + prediction,” respectively. The simulation is set up in Sec. 5, and the results are presented in Sec. 6. Finally Sec. 7 concludes this work.

### 1.1 Literature Review.

Dynamic models of building HVAC systems are typically nonlinear, which makes the planning problem in MPC a non-convex optimization problem. The nonlinearity comes from the existence of the bilinear terms—products between a state, temperature, and a control command, airflow rate; we will see this in Eqs. (4) and (6) later in this paper. Sometimes, the dynamic model is linearized to obtain a convex problem. Among works adopting this approach, some assume the value for a certain control command is known so that the product in which that command appears becomes linear in the remaining decision variables [14,17]. This reduces a degree-of-freedom that MPC can use. Others linearize around a trajectory, which requires an optimal (or at least a near optimal) trajectory first [13,18]. The quality of a linearized model is sensitive to the choice of the trajectory, and determining such a trajectory is challenging. After all, if it were easy, there would be little need for MPC. Identifying a linear black box model directly from data is also not straightforward (we discuss this later in Sec. 3). Recent progress in this direction is made in Refs. [10,19], which identify a linear model in which the input is the heat gain due to the HVAC system. However, the model is still not linear with respect to the control commands such as air flowrate. Convex relaxation of the MPC planning problem is thus far from trivial.

There has been recent attempts at convexification of the non-convex planning problems encountered in HVAC MPC [12,15,20]. Reference [20] does not require constraints to be satisfied at all time, but only with a pre-defined probability. Therefore, the resulting solution may not satisfy actuator constraints. In Ref. [12], values of the Lagrange multipliers are required for reformulation of the problem. The convexification approach using a McCormick envelope considered in Ref. [15] requires feasibility of the original problem (without slack) for all time. The original problem is likely to be infeasible when disturbance is large, since the actuator limits will prevent them from maintaining state constraints.

Another particular challenge is that the internal disturbance is also a large part of the heat load and hence a large part of the energy consumption in buildings. Therefore, internal disturbance prediction is needed to achieve the promised performance of MPC. Some works ignore the effects of internal disturbance in the MPC formulation altogether, e.g., Ref. [21]. Many leave the disturbance forecast question aside, assuming that future disturbance is somehow known to MPC, e.g., Refs. [12,15,22,23]. Some use stochastic optimization to address uncertainty in disturbance forecasts, e.g., Ref. [13]. A few works forecast internal disturbance from its estimate, which is obtained from measuring its various surrogates and modeling the relationship, such as plug loads [5] and occupancy and CO2 [4].

## 2 Architecture

### 2.1 Problem Description.

The focus of this study is the indoor climate control of a single-zone HVAC system shown in Fig. 2.

Fig. 2
Fig. 2
Close modal

In such a system, part of the air exhausted from the zone is recirculated and then mixed with outdoor air (OA) at a specified ratio. This mixed air (MA) is usually warm and humid, especially for hot-humid climates, and is therefore cooled and dehumidified by passing through a cooling coil. Dehumidification requires that the air is cooled enough for the water vapor to condense out of the air stream, so the conditioned air (CA) temperature (after the cooling coil), is usually too cold for a comfortable indoor climate. It is reheated by the reheat coil up to supply air (SA) temperature and then delivered into the zone.

The goal of the control system designed in this study is to decide the control commands to maintain the zone temperature (Tz) within time-varying pre-determined bounds, while keeping the energy use as small as possible. The control commands are the setpoints for total airflow rate ($m˙$) and supply air temperature (Tsa). Lower level PI controllers will maintain these setpoints by varying fan speed and reheat valve position.

Although conditioned air temperature $Tca$ and the outdoor air ratio α (ratio of outdoor air flowrate to supply air flowrate) can also be varied, in this study we assume they are fixed. The conditioned air temperature is typically set to 12.8 °C in order to maintain zone humidity, which is an important aspect of thermal comfort [24]. Similarly, the outdoor air ratio α is pre-specified at a constant value, and the minimum allowed value for the supply air flowrate $m˙$ is computed so that the OA flowrate $αm˙$ meets ventilation requirements [25].

### 2.2 Control System Architecture.

The control architecture proposed in this study is shown in Fig. 1. It involves two main components: (i) the ID and prediction block and (ii) MPC planner block that uses the models and forecasts to compute control commands.

Model predictive control of a system xk+1 = fk(xk, uk, vk), yk = hk(xk, uk, vk), with x being the state, u being the control command and v being the uncontrollable inputs, involves minimization of a cost function $Ji=∑k=ii+N−1ck(x^k+1,uk,v^k)$ over the planning horizon N with ck(·) being the energy used during the interval between k and k + 1. At time index i, an optimization problem of minimizing Ji subject to the system model and other constraints is posed based on the current estimate $x^i$ of the state xi and forecasts $v^$ of uncontrollable inputs v. The solution to this problem yields optimal commands ui, ui+1, …, ui+N−1. The first entry, ui, is implemented. At the next time index i + 1, the procedure is repeated.

In the adaptive architecture proposed here, the model is updated periodically by the system identifier, although at a much slower time scale than that of the control command update. In the numerical studies later reported, the model is updated every week while control commands are updated every 15 min.

## 3 (Block II) Model Predictive Control Planner

We start with the second block (MPC planner) of Fig. 1, since it is in charge of computing control commands, and the other blocks are there merely to support it. The MPC planner needs models to describe the energy consumption and the temperature dynamics, which appear as part of the cost and constraints in the planning problem. Energy is the integral of power, and the power consumption of the HVAC system consists of fan power, cooling power, and reheat power. The fan power is modeled as [26]
$Pkfan=afm˙k2$
(1)
with $m˙$ being the total airflow rate (kg/s). The cooling power $Pkcc$ is the electrical power consumed by the chiller to cool down the warm mixed air as it passes through the cooling coil:
$Pkcc={Cpam˙k(Tkma−Tca)COPTkma>Tca0otherwise$
(2)
where Cpa is the specific heat of air at constant pressure, Tca (°C) is the conditioned air temperature, COP is the chiller performance coefficient, and the mixed air temperature $Tkma$ (°C) is given by
$Tkma=αTkoa+(1−α)Tkz$
(3)
where $α(=m˙koam˙k)$ is the outdoor air ratio, $Tkoa$ is the outdoor air temperature (°C), $m˙koa$ is the outdoor airflow rate (kg/s), and $Tkz$ (°C) is the indoor zone temperature. The power consumed by reheat coil is modeled as the heat it adds to the conditioned air:
$Pkrh=Cpam˙k(Tksa−Tca)$
(4)
where $Tksa$ is the supply air temperature (°C).
Dynamics of zone temperature are modeled by the following second-order linear system with one output (indoor zone temperature $Tz$), and four inputs (heat injected by the HVAC system $qhvac$, outdoor temperature Toa (°C), solar irradiance ηsol (kW/m2), transformed disturbance $w¯$):
$xk+1=Axk+Bqkhvac+FvkTkz=Cxk+Dqkhvac+Gvk$
(5)
where $x∈R2$ is the state and $A∈R2×2$, $B∈R2×1$, $C∈R1×2$, $D∈R$, $F∈R2×3$, and $G∈R1×3$ are appropriate system matrices. The four inputs are separated into the single controllable input $qkhvac$ and three uncontrollable inputs $vk:=[Tkoa,ηksol,w¯k]T∈R3$, where $w¯$ is the transformed version of the internal heat load $qkint$ (kW); see Ref. [10] for details. It captures the effect of $qkint$ on the zone temperature $Tkz$. The quantity $qkhvac$ is related to the two actuation commands (supply airflow rate $m˙k$ and the deviation of supply air temperature $Tksa$) via the bilinear relation:
$qkhvac=Cpam˙k(Tksa−Tkz)$
(6)
Remark 1

Although qhvac is considered the controllable input in Eq. (5), it cannot be commanded directly. Only $m˙$ and Tsa can be commanded. Treating qhvac as the controllable input helps in two ways. First, it makes the model (5) linear, which aids model identification (discussed in Sec. 4). Second, the linear model is a convex constraint in the optimization problem the planner has to solve. We emphasize that a linear model structure with $m˙$ and Tsa as inputs, even though conceptually possible, is not useful for eventual use in MPC. The reason is that the sign of the DC gain (from $m˙$ to Tz) depends on whether the control commands are having a cooling or heating effect on the zone. If the supply air temperature Tsa is higher than the zone temperature Tz, increasing $m˙$ will increase the zone temperature. So, the DC gain is positive in such a scenario. The opposite happens when Tsa is lower than Tz. Now the DC gain has to be negative. However, a-priori knowledge of whether the control inputs will lead to heating or cooling is not available since that depends on both the state an control command.

### 3.1 Nominal Non-Convex Planner.

The goal of the MPC planner is to compute the control commands over the planning horizon, supply airflow rate $m˙$, and supply air temperature Tsa, to maintain thermal comfort while reducing energy use over that horizon. A direct translation of this goal into an optimization problem will be a non-convex problem, partly due to the bilinearity in Eq. (6). We first present this problem below, and then use it as a stepping stone to formulate a convex approximation that is actually used in the proposed MPC planner.

For notational simplicity, the current time index i is assumed to be 0 in this section. Define the decision variables as $zk:=[m˙k,Tksa,Tkma,Tkz,qkhvac,xk+1T,ϵkmin,ϵkmax]T∈R9$, in which $x∈R2$ is the state of the thermal model (5), $m˙$ and Tsa are the control commands, and N is the planning horizon. Let $x^0$ be the estimate of the current state obtained from a state estimator, and let $v^k$ ($:=[T^koa,η^ksol,w¯^k]T$) be the prediction of the uncontrollable inputs, for $k=0,…,N−1$. Specifically, $T^oa$ and $η^sol$ are from weather forecast, and $w¯^$ is provided by a disturbance predictor which will be discussed later in Sec. 4.3.

The nominal non-convex planning problem is
$minzk|k=0N−1J,J:=∑k=0N−1(Δt(Pkrh+Pkcc+Pkfan)+ρ(ϵkmin+ϵkmax))$
(7)
$s.t.−qkhvac+12zkTPczk=0$
(7a)
$−xk+1+Axk+Bqkhvac+Fv^k=02×1,x0=x^0$
(7b)
$−Tkz+Cxk+Dqkhvac+Gv^k=0$
(7c)
$−Tkma+(1−α)Tkz+αTkoa=0$
(7d)
$−m˙k≤−m˙min,m˙k≤m˙max$
(7e)
$−m˙k+1+m˙k≤m˙rateΔt$
(7f)
$m˙k+1−m˙k≤m˙rateΔt$
(7g)
$−Tksa≤−Tsa,min,Tksa≤Tsa,max$
(7h)
$−Tk+1sa+Tksa≤Tsa,rateΔt$
(7i)
$Tk+1sa−Tksa≤Tsa,rateΔt$
(7j)
$−Tkz−ϵkmin≤−Tz,min$
(7k)
$Tkz−ϵkmax≤Tz,max$
(7l)
$−ϵkmin≤0,−ϵkmax≤0k=0,…,N−1$
(7m)
where
$Pc=[0Cpa0−Cpa01×5Cpa008×8−Cpa05×1]$
(8)
and is obtained by rewriting Eq. (6).

Actuator constraints $[m˙min,m˙max]$ and $[Tsa,min,Tsa,max]$ represent the lower and upper bounds of the airflow rate and the supply air temperature, respectively. The minimum supply airflow rate, $m˙min$, is computed based on ventilation requirements [25]. To ensure reheat coil can only add heat, we require $Tsa,min=Tca$. Thermal comfort bounds are $[Tz,min,Tz,max]$. Slack variables εmin, εmax are used to relax the thermal comfort bounds from a fixed range $[Tz,min,Tz,max]$ to a variable range $[Tz,min−ϵmin,Tz,max+$$ϵmax]$. These slack variables help ensure that the problem is feasible. A high penalty parameter ρ encourages the slacks variables to be small so that temperature violation—when it occurs—is small.

For later convenience, we note that J can be compactly expressed as
$J=∑k=0N−1(12zkTPzk+qTzk)$
(9)
where
$P=[2αfCpaCpaCOP01×6CpaCpaCOP08×806×1]Δt$
(10)
$q=[−CpaTca1+COPCOPΔt06×1ρρ]$
(11)
by rewriting Eqs. (1)(2) and (4).
Proposition 1

Problem (7) is feasible.

The proof of Proposition 1 is provided in the  Appendix.

### 3.2 Proposed Convex Planner.

The optimization problem (7) is non-convex since the equality constraint (7b) is bilinear, and the quadratic term in the cost (9) involves the indefinite matrix P. The goal now is to approximate the problem (7) with a convex problem, so that the approximation is easy to solve and the obtained solution provides good approximation to that of problem (7).

The algorithm we propose to this end is described in Algorithm 1. It uses the convex-concave procedure (CCP) [27]. In Algorithm 1, the following terminology is used. Let P = Q+ + Λ)QT be the eigen-decomposition of the real symmetric matrix P from Eq. (10), where $Λ+≽0$ is the positive semi definite part and $Λ−≺0$ is the negative definite part. Define $P+:=QΛ+QT$ and $P−:=QΛ−QT$.

#### Convex planner

Algorithm 1

Input: Initial guess ζ(0).

n ← 0.

repeat

Convexify: Form:
$J^=∑k=0N−1(12zkTP+zk+(P−ζk(n)+q)Tzk−12ζk(n)TP−ζk(n)),h^1,k:−qhvac+ζk(n)TPczk−12ζk(n)TPcζk(n)=0.$
(12)
Solve forz*:
$z*=argminzk|k=0N−1J^$
(13)

s. t. equality constraints (12), (7b)(7d)

inequality constraints (7e)(7m)

k = 0, … N − 1

Update iteration: Set $n←n+1,ζ(n)←z*$.

until$‖ζ(n)−ζ(n−1)‖≤δ$;

Output:$z*←ζ(n)$

Proposition 2

[16] Problem (13) is feasible and convex, and Algorithm 1 is a descent and convergent algorithm.

Remark 2

Proposition 2 guarantees reliable performance of Algorithm 1. Since problem (13) is feasible and convex, if the algorithm converges within the allowable time, it converges to a local minimum of the original non-convex problem. If the algorithm must be stopped before convergence due to inadequate time, the solution obtained has a lower cost than solutions from previous iterates since it is a descent algorithm.

#### 3.2.1 Choice of Convex Approximation Method.

Apart from the convex-concave procedure we used, there are many approximation methods for non-convex optimization problems that involve bilinearities. The commonly used methods are (i) Branch-and-Bound (BnB) [28] and (ii) Alternate Convex Search (ACS) [29]. Next, we show that these methods are not applicable to our problem (7), leaving CCP as the only candidate. The following two propositions will be needed for that discussion.

Proposition 3

[16] The dual of problem (7) is unbounded from below.

Proposition 4

Every solution of Problem (7) is a boundary solution.

The proof of Proposition 4 is provided in the  Appendix.

##### 3.2.1.1 Inapplicability of branch-and-bound (BnB).

BnB requires construction of a tight convex under-estimator of the NLP within any given region of the space of the variables [28]. The most widely used under-estimators are Lagrangian relaxation [30] or convex relaxation. However, Proposition 3 shows the dual of our problem (7) is unbounded from below. Therefore, Lagrangian relaxation cannot be applied. For convex relaxation, common options are McCormick envelope [28] and reformulation linearization technique (RLT) [31]. Both of them reformulate a problem via the addition of certain nonlinear constraints that are generated by using the products of the bounding constraints. However, constructing such products require knowledge of bounds on variables that are involved. In our problem, thermal comfort limits do not have known bounds because of the introduction of slack variables. Hence, convex relaxation is also not applicable for our case.

##### 3.2.1.2 Inapplicability of alternate convex search (ACS).

ACS [29] divides variable set into disjoint blocks and in every step, only the variables of an active block are optimized while those of the other blocks are fixed. Analyses and examples from Refs. [32,33] show that this method will most likely fail to find a local optimum for problem with boundary solutions (our case). Only initial guesses that belong to a particular set will converge to a local optimum. Because there is no guarantee on convergence to local minima, we do not use ACS.

## 4 (Block I) Identification and Prediction

### 4.1 Identification.

The job of the identification block of Fig. 1 is to use data to identify parameters of the zone temperature dynamics model (5), along with the unmeasurable occupant-induced disturbance. We rewrite the model in a different form to describe the identification problem precisely:
$xk+1=Axk+Bidukid+Fidw¯kyk(=Tkz)=Cxk+Didukid+Gidw¯k$
(14)
Here, the state $xk∈R2$, the output $yk∈R$, and the matrices A, C are the same as in Eq. (5). But while the four inputs in Eq. (14) were divided into controllable and not controllable; here, they are divided into measurable and non-measurable. In particular, $ukid:=[qhvac,Toa,ηsol]k∈R3$ consists of the measurable inputs to the thermal dynamics and the transformed disturbance $w¯k∈R$ is the non-measurable input. Other than the regrouping, the two models are identical. Among the three components of $ukid$, qhvac is computed from measurements of $m˙$ and Tsa using Eq. (6), and the remaining two inputs can be obtained from a weather station. The output $Tz$ is measured with a sensor.

The system identification algorithm used here is the SPDIR method proposed in our earlier work [10]. Fix i as the current time when system identification is to be carried out. Define τi: = {iN, iN + 1, …, i − 1} and (uid, y)j, jτi be the measured input–output data for the model (14) over that time interval. The algorithm SPDIR takes this data and produces an estimate of the model parameters $M:=(A,Bid,Fid,C,Did,Gid)$ and an estimate of the transformed disturbance $w¯j,j∈τi$. We denote these estimates $M^i$ and $w¯^j,j∈τi$ since they depend on i. The SPDIR algorithm is executed at time instants i, i + Nad, i + 2Nad, …, with Nad large so that enough time has after the previous identification to warrant updating the estimates of the model and disturbance.

The SPDIR algorithm comes with the following guarantees [10]:

1. The computation involved in obtaining the estimates (model and disturbance signal) is a feasible and convex optimization problem with a strictly convex cost.

2. The model $M^i$ is BIBO stable and has a positive DC gain from each of the three measurable inputs (outdoor temperature, solar irradiance, and HVAC heat injection) to indoor temperature.

3. There is exactly one parameter that requires tuning by a human expert. This tuning can be done once (one data set). The two properties mentioned above hold irrespective of the value of this parameter.

Remark 3

The first property ensures that the system identification algorithm can be executed periodically without any human intervention, i.e., autonomously. Autonomy is also helped by the third feature. The second feature helps in two ways. One, it ensures that the model identified is consistent with the physics of HVAC systems. Two, it helps in state estimation. At every decision instant i, a Kalman filter is used to estimate the state of the thermal model (5), which is then used as the initial state by the MPC planner : $x^0$ in Eq. (7b). The stability guarantee of the model mentioned above ensures that the Kalman filter is stable [34].

### 4.2 Forecasts of Uncontrollable Inputs.

Two types of uncontrollable inputs appear in the thermal model (5), and thus their forecasts over the planning horizon are needed by the MPC planner: weather variables and transformed disturbance $w¯$. These forecasts are obtained as follows:

1. Weather variables: Obtain forecast of $[Toa,ηsol]kT$ over the next planning horizon from an online weather service.

2. Transformed disturbance $w¯$: If the prediction horizon does not contain a holiday, assign the disturbance for the same time interval from the previous week estimated by the system identifier, as the forecast. If the prediction horizon contains a holiday, use the disturbance estimate from the same time interval of the previous Saturday as the forecast. This method is similar to the one used in Ref. [5], except for the holiday corrections.

### 4.3 Putting Them All Together.

The components described so far are now combined to form the proposed controller. The architecture is described in Algorithm 2. Recall that Fig. 1 shows the complete closed loop system.

#### Proposed MPC architecture

Algorithm 2

Input: Planning horizon $N∈Z+$, control horizon $Nc∈Z+$, and model updating interval $Nad∈Z+$.

Setup:$Sc:={Nc,2Nc,…}$, $Sad:={Nad,2Nad,…}$

for$i=1,2,…$do

if$i∈Sad$then

Measure: Input uid and output y of the model (14), over the time interval [i − Nad : i − 1].

System ID: Estimate model $M^i$ and disturbance $w¯^[i−Nad:i−1]$ using the SPDIR algorithm from [10].

end

Estimate state: Estimate current state $x^[i]$ of thermal model (5) using a Kalman filter.

Predict disturbance: As described in Section 4.2

Optimize: Compute control decisions $m˙[i:i+N−1]$ and $Tsa[i:i+N−1]$ using Algorithm 1.

Implement: Apply $m˙[i:i+Nc−1]$ and $Tsa[i:i+Nc−1]$ to the plant.

end

### 4.4 Baseline Controller for Comparison.

The baseline controller is chosen to be the single-maximum controller which is widely used in practice [35]. The single-maximum controller operates the HVAC system in three modes depending on where the zone temperature $Tz$ is compared with the deadband $[Tz,min,Tz,max]$. When $Tz$ exceeds the upper bound $Tz,max$, reheat is turned off and the supply airflow rate $m˙$ is increased with the help of a PI controller. When the zone temperature is below the lower bound $Tz,min$, the airflow rate $m˙$ is kept at the minimum allowed value but the supply air temperature is increased with the help of a PI controller. When the zone temperature is in the deadband $[Tz,min,Tz,max]$, the supply air temperature is kept at $Tca$ and the flowrate are both kept at the minimum allowed value.

## 5 Simulation Setup

To assess performance of the proposed control system, we perform closed loop simulations for nearly a yearlong period with a realistic time varying plant. Simulations with the baseline controller are also performed on the same plant for comparison. The plant model on which the controller acts is calibrated to mimic a large auditorium in a building in the University of Florida campus (Pugh Hall). The auditorium in Pugh Hall is served by an air handling unit, and it has the same HVAC system configuration as shown in Fig. 2.

### 5.1 Plant Description.

The plant is a time varying non-linear ordinary differential equation, and with a large internal heat load qint (kW).
$Cz(t)T˙z(t)=Tw(t)−Tz(t)Rz(t)+qhvac(t)+Ae(t)ηsol(t)+qint(t)Cw(t)T˙w(t)=Toa(t)−Tw(t)Rw(t)+Tz(t)−Tw(t)Rz(t)$
(15)
where Tw (°C) is the wall temperature, Cz(t), Cw(t), Rz(t), and Rw(t) are the time-varying thermal capacitances and resistances of the zone and wall, respectively, and $Ae(t)$ is the effective area of the building for incident solar radiation. One can view this model as a time-varying version of the commonly used RC-network models of building thermal dynamics.

The time-varying plant parameters are shown in Fig. 3, which are chosen as follows. The average values of the time-varying parameters are chosen to be the same as the values given in Ref. [9], which contains the plant parameters estimated using data from an auditorium in Pugh Hall located in the University of Florida.

Fig. 3
Fig. 3
Close modal

### 5.2 Closed Loop Parameters.

The planning horizon for MPC is 1 day and the control horizon is 15 min, with a sampling time Δt = 5 min, so $N=288$ and $Nc=3$. These choices are inspired by the study presented in Ref. [36]. The total time span for MPC is 50 weeks. The number of decision variables for problems (7) and (13) is 2592($=9N$). The plant was simulated in matlab by discretizing the differential equation (13). Future work will explore using publicly available matlab-based simulators such as Ref. [37].

Thermal comfort and flowrate constraints depend on whether the building is in occupied or unoccupied mode [24]. The maximum occupancy for Pugh Hall auditorium is 229 persons, and its occupied mode (occ) is scheduled from 6:30 AM to 10:30 PM while the remaining time is deemed unoccupied (unocc). We used these parameters for the simulation. The thermal comfort bounds are [21.9, 23.6]°C for occupied mode and [21.1, 24.4]°C for unoccupied mode. The minimum allowed value for the supply airflow rate $m˙min$ is computed based on the ventilation requirements specified in ASHRAE 62.1 [25]. More specifically, $m˙min,unocc$ is computed assuming 0 occupancy for the unoccupied period with 31% occupancy for the occupied period. Note for the baseline controller, $m˙min,occ$ is kept as high as 1.90 kg/s; otherwise, the baseline controller fails to maintain the zone temperature comfort satisfactorily. The remaining parameters are listed in Table 1.

Table 1

Parameters for baseline and MPC controllers

 Unoccupied Occupied $Tsa,min$ 12.8 °C $Tca$ 12.8 °C $Tz,min$ 21.1 21.9 °C $Tsa,max$ 37.8 °C $COP$ 3.5 N/A $Tz,max$ 24.4 23.6 °C $Tsa,rate$ 0.56 °C/min α 0.3 N/A $m˙min$ 0.95 1.47, $1.90a$ kg/s $m˙max$ 4.72 kg/s af 417.5 W/(kg/s)2 $m˙rate$ 0.2 kg/s/min
 Unoccupied Occupied $Tsa,min$ 12.8 °C $Tca$ 12.8 °C $Tz,min$ 21.1 21.9 °C $Tsa,max$ 37.8 °C $COP$ 3.5 N/A $Tz,max$ 24.4 23.6 °C $Tsa,rate$ 0.56 °C/min α 0.3 N/A $m˙min$ 0.95 1.47, $1.90a$ kg/s $m˙max$ 4.72 kg/s af 417.5 W/(kg/s)2 $m˙rate$ 0.2 kg/s/min

$am˙min,occ=1.47$ is used for the MPC controllers and $m˙min,occ=1.90$ is used for the baseline controller.

The uncontrollable input signals are chosen as follows: solar irradiance data ηsol is taken from NSRDB [38], and ambient temperature Toa is taken from online,1 both for Gainesville, FL, from the year of 2013. The internal heat load (qint) is chosen by scaling CO2 data collected from the auditorium in Pugh Hall during the same year, which is shown in Fig. 4. The high-resolution and long-term data collection was made possible by using a custom made data logger [39]. The rationale is that occupancy is correlated to the CO2 level. Note that the heat load is by design a large, time-varying, and aperiodic signal.

Fig. 4
Fig. 4
Close modal

All numerical results presented in this work are obtained through matlab. Specifically, the plant is simulated in simulink©. The system identification problem from Ref. [10] for estimating model and disturbance is solved using cvx© [40] package. For control computation, the nominal non-convex problem (7) is solved using ipopt© [41] package, and the proposed convex problem (13) is solved using cvx© [40] package. We used a desktop computer with a 3.60GHz × 8 CPU and 16 GB RAM, running Linux, for the closed loop simulations.

## 6 Simulation Results

A total of five distinct controllers are tested through simulations on the same plant:

• Baseline: the single-max controller described in Sec. 4.4.

• Proposed (Adapt-CVX): the proposed controller (Algorithm 2), with both model update and convex planner for control computation.

• NAdapt-CVX: the proposed controller (Algorithm 2), but without updating the dynamic model and the disturbance estimates.

• Adapt-NCVX: the proposed controller (Algorithm 2), but using the non-convex problem (7) instead of the convex problem from Algorithm 1 to compute commands.

• NAdapt-NCVX: MPC with the nominal non-convex problem (7) for computing control commands, and without updating the dynamic model and the disturbance estimates. Note this the MPC architecture generally described in the literature.

In all the controllers that uses a non-convex optimization, if the NLP solver is unable to converge before the control update interval is over, decisions computed by the baseline controller are sent to the actuators.

### 6.1 Comparison With the Baseline Controller.

The proposed MPC scheme outperforms the baseline controller in both maintaining zone temperature and reducing energy use, see Table 2. Data on uncontrollable inputs, control command, and the output (zone temperature) are shown in Fig. 5 for the full 50 weeks. Figure 6 zooms into one week: Aug. 26, 2013, to Sept. 1, 2013.

Fig. 5
Fig. 5
Close modal
Fig. 6
Fig. 6
Close modal
Table 2

Performance comparison among various controllers

ControllerSite EUI (kBtu/ (ft2 ·yr))Planner failure (%)RMSE of Tz violation (°C)Max Tz violation (°C)
Baseline72.9N/A0.452.3
ControllerSite EUI (kBtu/ (ft2 ·yr))Planner failure (%)RMSE of Tz violation (°C)Max Tz violation (°C)
Baseline72.9N/A0.452.3

In particular, the proposed controller (Algorithm 2) reduces energy use by $26.8%$ over the baseline controller, to EUI = 53.5 kBtu/(ft2 ·year), see Table 2. The baseline controller is already more efficient than the average controller in the field: its site EUI for the tested period is 72.9 kBtu/(ft2 ·year), which is lower than the median site EUI = 84.3 kBtu/(ft2 ·year) for college buildings in the United States [42].

The improvement in performance over the baseline controller is consistent with results in the literature that have compared MPC with baseline controllers. MPC’s ability to use disturbance forecasts and prediction from the model allows it to make better decisions than a purely output feedback controller.

### 6.2 Benefit and Necessity of the Design Features

#### 6.2.1 Need for Model and Disturbance Ppdate.

We tested the role and/or value of adaptation by turning off the adaptation block. A model and disturbance (for a week) are estimated from data from the first week of 2013. They are used by the controller for every week of the year. The resulting MPC controller is referred to by the “NAdapt-” prefix, e.g., in Table 2. We see from the table that adaptation reduces energy use by about 16% and reduces zone temperature violation over the non-adaptation case.

Thus, adaptation—periodically updating models and disturbances from data—is both necessary (for indoor comfort) and beneficial (improves energy use) for an MPC-based controller for HVAC systems.

#### 6.2.2 Need for Convexification of the MPC Planner.

NLP solvers such as ipout [41] are quite powerful. Thus, solving the non-convex MPC planning problem (7) is usually not an issue. On average it takes 2.7 s for ipopt to find a local minimum of the non-convex problem, failing to do so with the available 15 min only $0.1%$ of the time. When this happens, decision from the baseline controller is used as control commands. The resulting switching control action can lead to large violation in the indoor temperature. See Fig. 7 for an example of this phenomenon. The zone temperature exceeds the upper bound by 3.2 °C for an extended period of time. Thus, though a non-convex planner rarely fails, when it does it leads a catastrophic loss of performance that will render the control system unacceptable to the user.

Fig. 7
Fig. 7
Close modal

In contrast to MPC with a non-convex planner, the proposed MPC scheme with a convexified planner always finds a minimum within the available 15 min, taking 1.7 s on average to compute the control decisions. Partially as a result of that, it is able to provide the best performance in maintaining zone temperature among all five controllers tested.

Therefore, even though solving the nominal non-convex problem is rarely an issue, in those rare occasions the controller can cause serious disruption to occupant’s thermal comfort. It is unlikely such a control system will be acceptable to building owners and occupants. In short, the convex approximation of the MPC planner is necessary.

It should be noted that the NAdapt-NCVX controller is the MPC scheme generally used in the literature, e.g., Refs. [20,22]. Without the benefits from both of the designed features, this controller has a maximum zone temperature violation of 4.0 °C, even though it occurs rarely and does not perform as well as the proposed controller in terms of energy use.

Remark 4

We remark here the performance delivered by the proposed MPC scheme is obtained under strong plant-model mismatch in the following aspects: (i) The plant is time-varying and nonlinear, while the MPC planner uses a linear model. (ii) The proposed MPC scheme assumes both the plant and the disturbance are the same as that from the previous week, but the plant and the disturbance do not satisfy those properties.

## 7 Conclusion

This paper takes a first stab at designing an MPC-based control system for HVAC systems that can operate autonomously for long periods, without requiring intervention of human experts. Autonomy is made possible by two features: (i) automated periodic update of thermal model and internal disturbance signals and (ii) a convex approximation of the MPC planner’s optimization problem. The yearlong simulations shows that both of the features are essential to get the performance improvement over the simple baseline controller over a long period. The need for periodic re-learning the model and disturbance is easy to see in the context of buildings. The need for convexity in the planning problem is less obvious at the design stage, but was discovered from the simulation results. Even though the nominal non-convex planning problem can be used effectively nearly 100% of times, the rare instances it fails to converge causes dramatic fluctuations in the indoor temperature rendering the control system an unlikely contender for real-life application. Without these features, though MPC can outperform the baseline controller in certain scenarios, the benefits may not be substantial enough to defray the additional cost of implementing MPC.

At the current stage, the proposed MPC architecture uses arguably one of the simplest schemes for forecasting of the internal disturbance. It is envisioned that a more accurate prediction scheme, possibly with the aid of technologies such as occupancy recognition or CO2 level sensing, should further improve performance of the MPC controller.

Many extensions of this work are possible. The most immediate next step is extending the proposed control scheme to include humidity dynamics and ventilation requirements, which will require including as part of the control commands the outdoor airflow and conditioned air temperature (downstream of the cooling/dehumidification coil; see Fig. 2). These two have been assumed fixed in this paper but in fact can be commanded through the building automation system. It should reduce energy use even more and provide better thermal comfort by including outdoor airflow and conditioned air temperature into the list of control commands. The challenge is to incorporate the nonlinear humidity dynamics in zone thermal models and the nonlinear process models of the cooling/dehumidification coil [22]. The autonomy achieved by the control system proposed here is due to the use of linear dynamic models. Other useful directions include extension to multi-zone buildings, improvements in the forecasting methodology for the internal disturbance, etc.

## Acknowledgment

This research is partially supported by NSF through Grant Nos. 1463316 and 1934322.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. Data provided by a third party are listed in Acknowledgment.

## References

1.
2018
. Annual Energy Outlook, April, https://www.eia.gov/outlooks/archive/aeo18/pdf/AEO2018.pdf
2.
Serale
,
G.
,
Fiorentini
,
M.
,
Capozzoli
,
A.
,
Bernardini
,
D.
, and
,
A.
,
2018
, “
Model Predictive Control (MPC) for Enhancing Building and HVAC System Energy Efficiency: Problem Formulation, Applications and Opportunities
,”
Energies
,
11
(
3
), p.
631
.
3.
Castilla
,
M.
,
Álvarez
,
J.
,
Normey-Rico
,
J.
, and
Rodríguez
,
F.
,
2014
, “
Thermal Comfort Control Using a Non-linear MPC Strategy: A Real Case of Study in a Bioclimatic Building
,”
J. Process. Control.
,
24
(
6
), pp.
703
713
.
4.
Bengea
,
S. C.
,
Kelman
,
A. D.
,
Borrelli
,
F.
,
Taylor
,
R.
, and
Narayanan
,
S.
,
2014
, “
Implementation of Model Predictive Control for An HVAC System in a Mid-Size Commercial Building
,”
HVAC&R Res.
,
20
(
1
), pp.
121
135
.
5.
De Coninck
,
R.
, and
Helsen
,
L.
,
2016
, “
Practical Implementation and Evaluation of Model Predictive Control for An Office Building in Brussels
,”
Energy Buil.
,
111
, pp.
290
298
.
6.
,
H.
, and
Holst
,
J.
,
1995
, “
Estimation of Continuous-Time Models for the Heat Dynamics of a Building
,”
Energy Buil.
,
22
(
1
), pp.
67
79
.
7.
Berthou
,
T.
,
Stabat
,
P.
,
Salvazet
,
R.
, and
Marchio
,
D.
,
2014
, “
Development and Validation of a Gray Box Model to Predict Thermal Behavior of Occupied Office Buildings
,”
Energy Buil.
,
74
, pp.
91
100
.
8.
Kim
,
D.
,
Cai
,
J.
,
Ariyur
,
K. B.
, and
Braun
,
J. E.
,
2016
, “
System Identification for Building Thermal Systems Under the Presence of Unmeasured Disturbances in Closed Loop Operation: Lumped Disturbance Modeling Approach
,”
Build. Environ.
,
107
, pp.
169
180
.
9.
Coffman
,
A. R.
, and
Barooah
,
P.
,
2018
, “
Simultaneous Identification of Dynamic Model and Occupant-Induced Disturbance for Commercial Buildings
,”
Buil. Environ.
,
128
, pp.
153
160
.
10.
Zeng
,
T.
,
Brooks
,
J.
, and
Barooah
,
P.
,
2021
, “
Simultaneous Identification of Linear Building Dynamic Model and Disturbance Using Sparsity-Promoting Optimization
,”
Automatica
,
129
, p.
109631
.
11.
Smarra
,
F.
,
Jain
,
A.
,
de Rubeis
,
T.
,
Ambrosini
,
D.
,
D’Innocenzo
,
A.
, and
Mangharam
,
R.
,
2018
, “
Data-Driven Model Predictive Control Using Random Forests for Building Energy Optimization and Climate Control
,”
Appl. Energy
,
226
, pp.
1252
1272
.
12.
Kelman
,
A.
, and
Borrelli
,
F.
,
2011
, “
Bilinear Model Predictive Control of a HVAC System Using Sequential Quadratic Programming
,”
IFAC Proc. Vol.
,
44
(
1
), pp.
9869
9874
.
13.
Ma
,
Y.
,
Matuško
,
J.
, and
Borrelli
,
F.
,
2015
, “
Stochastic Model Predictive Control for Building HVAC Systems: Complexity and Conservatism
,”
IEEE Trans. Control Syst. Tech.
,
23
(
1
), pp.
101
116
.
14.
Parisio
,
A.
,
Fabietti
,
L.
,
Molinari
,
M.
,
Varagnolo
,
D.
, and
Johansson
,
K. H.
,
2014
, “
Control of HVAC Systems Via Scenario-Based Explicit MPC
,”
53rd IEEE Conference on Decision and Control
,
Los Angeles, CA
,
IEEE, pp. 5201–5207
.
15.
Atam
,
E.
, and
Helsen
,
L.
,
2015
, “
A Convex Approach to a Class of Non-Convex Building HVAC Control Problems: Illustration by Two Case Studies
,”
Energy Build.
,
93
, pp.
269
281
.
16.
Zeng
,
T.
, and
Barooah
,
P.
,
July 2020
, “
An Autonomous MPC Scheme for Energy-Efficient Control of Building HVAC Systems
,”
2020 American Control Conference (ACC)
,
Denver, CO
,
IEEE, pp. 4213–4218
.
17.
Razmara
,
M.
,
Maasoumy
,
M.
,
Shahbakhti
,
M.
, and
Robinett III
,
R.
,
2015
, “
Optimal Exergy Control of Building HVAC System
,”
Appl. Energy
,
156
, pp.
555
565
.
18.
Koehler
,
S.
, and
Borrelli
,
F.
,
2013
, “
Building Temperature Distributed Control Via Explicit MPC and Trim and Respond Methods
,”
2013 European Control Conference (ECC)
,
Zurich, Switzerland
,
IEEE, pp. 4334–4339
.
19.
Zeng
,
T.
, and
Barooah
,
P.
,
2020
, “
Identification of Network Dynamics and Disturbance for a Multi-Zone Building
,”
IEEE Trans. Control Syst. Tech.
,
28
(
5
), pp.
2061
2068
.
20.
Oldewurtel
,
F.
,
Parisio
,
A.
,
Jones
,
C. N.
,
Gyalistras
,
D.
,
Gwerder
,
M.
,
Stauch
,
V.
,
Lehmann
,
B.
, and
Morari
,
M.
,
2012
, “
Use of Model Predictive Control and Weather Forecasts for Energy Efficient Building Climate Control
,”
Energy Buil.
,
45
, pp.
15
27
.
21.
Prívara
,
S.
,
Širokỳ
,
J.
,
Ferkl
,
L.
, and
Cigler
,
J.
,
2011
, “
Model Predictive Control of a Building Heating System: The First Experience
,”
Energy Buil.
,
43
(
2–3
), pp.
564
572
.
22.
Raman
,
N. S.
,
,
K.
,
Chen
,
B.
,
Ingley
,
H. A.
, and
Barooah
,
P.
,
2020
, “
Model Predictive Control for Energy-Efficient HVAC Operation with Humidity and Latent Heat Considerations
,”
Appl. Energy
,
279
, p.
115765
.
23.
Hou
,
X.
,
Xiao
,
Y.
,
Cai
,
J.
,
Hu
,
J.
, and
Braun
,
J. E.
,
2017
, “
Distributed Model Predictive Control Via Proximal Jacobian ADMM for Building Control Applications
,”
2017 American Control Conference (ACC)
,
Seattle, WA
,
IEEE, pp. 37–43
.
24.
American Society of Heating, Refrigerating and Air-Conditioning Engineers
,
2017
. ANSI/ASHRAE Standard 55-2017, Thermal Environmental Conditions for Human Occupancy.
25.
ASHRAE
,
2016
. ANSI/ASHRAE Standard 62.1-2016, Ventilation for Acceptable Air Quality.
26.
Raman
,
N. S.
, and
Barooah
,
P.
,
2020
, “
On the Round-Trip Efficiency of An HVAC-Based Virtual Battery
,”
IEEE Trans. Smart Grid
,
11
(
1
), pp.
403
410
.
27.
Lipp
,
T.
, and
Boyd
,
S.
,
2016
, “
Variations and Extension of the Convex-Concave Procedure
,”
Optim. Eng.
,
17
(
2
), pp.
263
287
.
28.
McCormick
,
G. P.
,
1976
, “
Computability of Global Solutions to Factorable Nonconvex Programs: Part I–Convex Underestimating Problems
,”
Math. Program.
,
10
(
1
), pp.
147
175
.
29.
Wendell
,
R. E.
, and
Hurter Jr
,
A. P.
,
1976
, “
Minimization of a Non-Separable Objective Function Subject to Disjoint Constraints
,”
Operat. Res.
,
24
(
4
), pp.
643
657
.
30.
d’Aspremont
,
A.
, and
Boyd
,
S.
,
2003
, “
Relaxations and Randomized Methods for Nonconvex QCQPs
,”
EE392o Class Notes, Stanford Univ.
,
1
, pp.
1
16
.
31.
Sherali
,
H. D.
, and
Alameddine
,
A.
,
1992
, “
A New Reformulation-Linearization Technique for Bilinear Programming Problems
,”
J. Global Optim.
,
2
(
4
), pp.
379
410
.
32.
Gorski
,
J.
,
Pfeuffer
,
F.
, and
Klamroth
,
K.
,
2007
, “
Biconvex Sets and Optimization With Biconvex Functions: A Survey and Extensions
,”
Math. Methods Operat. Res.
,
66
(
3
), pp.
373
407
.
33.
Floudas
,
C. A.
, and
Visweswaran
,
V.
,
1990
, “
A Global Optimization Algorithm (gop) for Certain Classes of Nonconvex Nlps–i. Theory
,”
Comput. Chem. Eng.
,
14
(
12
), pp.
1397
1417
.
34.
Rhodes
,
I.
,
1971
, “
A Tutorial Introduction to Estimation and Filtering
,”
IEEE. Trans. Automat. Contr.
,
16
(
6
), pp.
688
706
.
35.
ASHRAE
,
2009
. The ASHRAE Handbook Fundamentals (SI Edition).
36.
Huang
,
S.
,
Lin
,
Y.
,
Chinde
,
V.
,
Ma
,
X.
, and
Lian
,
J.
,
2021
, “
Simulation-Based Performance Evaluation of Model Predictive Control for Building Energy Systems
,”
Appl. Energy
,
281
, pp.
116027
.
37.
Kircher
,
K. J.
,
Schaefer
,
W.
, and
Max Zhang
,
K.
,
2020
, “
A Computationally Efficient, High-Fidelity Testbed for Building Climate Control
,”
ASME J. Eng. Sustainable Buil. Cities
,
2
(
1
), p.
011002
.
38.
Sengupta
,
M.
,
Xie
,
Y.
,
Lopez
,
A.
,
Habte
,
A.
,
Maclaurin
,
G.
, and
Shelby
,
J.
,
2018
, “
The National Solar Radiation Data Base (NSRDB)
,”
Renewable. Sustainable. Energy. Rev.
,
89
, pp.
51
60
.
39.
Middelkoop
,
T.
,
2020
, An Asynchronous BACnet Logger Written in Python.
40.
Grant
,
M.
, and
Boyd
,
S.
,
2011
, CVX: Matlab Software for Disciplined Convex Programming, version 1.21. http://cvxr.com/cvx
41.
Wächter
,
A.
, and
Biegler
,
L. T.
,
2006
, “
On the Implementation of an Interior-Point Filter Line-Search Algorithm for Large-Scale Nonlinear Programming
,”
Math. Program.
,
106
(
1
), pp.
25
57
.
42.
Star
,
E.
,
2018
, “
Portfolio Manager
.”
US Energy Use Intensity by Property Type
.