An Integrated Design of Optimization and Physical Dynamics for Energy Efficient Buildings: A Passivity Approach

TL;DR

This paper presents a passivity-based integrated approach combining optimization and physical dynamics to improve energy efficiency in building HVAC systems, ensuring convergence to optimal temperature setpoints despite disturbances.

Contribution

It introduces a novel passivity-based framework that unifies optimization and physical dynamics for energy-efficient building temperature control.

Findings

Proved passivity of the optimization dynamics.

Ensured convergence of room temperatures to optimal values.

Validated the approach through simulations.

Abstract

In this paper, we address energy management for heating, ventilation, and air-conditioning (HVAC) systems in buildings, and present a novel combined optimization and control approach. We first formulate a thermal dynamics and an associated optimization problem. An optimization dynamics is then designed based on a standard primal-dual algorithm, and its strict passivity is proved. We then design a local controller and prove that the physical dynamics with the controller is ensured to be passivity-short. Based on these passivity results, we interconnect the optimization and physical dynamics, and prove convergence of the room temperatures to the optimal ones defined for unmeasurable disturbances. Finally, we demonstrate the present algorithms through simulation.

Equations

$$\dot{x}=\varphi (x,u),y=\phi (x,u),$$

$$S(x(t))-S(x(0))\le {\int }_0^t{y}^T(\tau )u(\tau )d\tau$$

$$S(x(t))-S(x(0))\le {\int }_0^t{y}^T(\tau )u(\tau )-\epsilon \Vert y(\tau ){\Vert }^{2}d\tau$$

$$S(x(t))-S(x(0))\le {\int }_0^t{y}^T(\tau )u(\tau )+\epsilon \Vert u(\tau ){\Vert }^{2}d\tau$$

$$C_i{\dot{T}}_i=\frac{T^{\mathrm{a}}-T_i}{R_i}+\sum _{j\in {\mathcal{N}}_i}\frac{T_j-T_i}{R_{ij}}+a_i(T_i^{\mathrm{s}}-T_i)m_i+q_i,$$

$$C_i{\dot{T}}_i=\frac{T^{\mathrm{a}}-T_i}{R_i}+\sum _{j\in {\mathcal{N}}_i}\frac{T_j-T_i}{R_{ij}}.$$

$$C\dot{T}=R{T}^{\mathrm{a}}1-RT-LT+BG(T)m+Bq,$$

$$C\dot{\delta T}=R\delta T^{\mathrm{a}}1-R\delta T-L\delta T+BG(\bar{T})\delta m-\bar{U}\delta T+B\delta q,$$

$$x$$

$$w_{\mathrm{a}}$$

$$\dot{x}=-Ax+Bu+B{w}_{\mathrm{q}}+w_{\mathrm{a}},x:=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]$$

$$\underset{z=[z_x^{\top }{z}_u^{\top }{]}^{\top }\in {\mathbb{R}}^{n+n_1}}{\min }\Vert z_{x1}-h{\Vert }^2+f(z_u)$$

$$\text{mbox}subjectto:g(z_u)\le 0$$

$$-A{z}_x+B{z}_u+B{d}_{\mathrm{q}}+d_{\mathrm{a}}=0$$

$$\underset{z_u\in {\mathbb{R}}^{n_1}}{\min }\Vert B^{\top }{A}^{-1}(B{z}_u+B{d}_{\mathrm{q}}+d_{\mathrm{a}}-\bar{h}){\Vert }^2+f(z_u)$$

$$\text{mbox}subjectto:g(z_u)\le 0$$

$$M^2(z_u^{\ast }+d_{\mathrm{q}})+N(d_{\mathrm{a}}-\bar{h})+\nabla f(z_u^{\ast })+\nabla g(z_u^{\ast }){\lambda }^{\ast }=0,$$

$$g(z_u^{\ast })\le 0,{\lambda }^{\ast }\ge 0,{\lambda }^{\ast }\circ g(z_u^{\ast })=0,$$

$$z_{x1}^{\ast }=M(z_u^{\ast }+d_{\mathrm{q}})+B^{\top }{A}^{-1}{d}_{\mathrm{a}}.$$

$${\dot{\hat{z}}}_u$$

$$\dot{\hat{\lambda }}$$

$$\displaystyle([b]^{+}_{a})_{l}=\left\{\begin{array}[]{l}0,\mbox{ if }a_{l}=0\mbox{ and }b_{l}<0\\ b_{l},\mbox{ otherwise}\end{array}\right.,$$

$$y_{\mathrm{o}}$$

$$y_{\mathrm{o}}^{\ast }$$

$$D^{+}{S}_{\mathrm{o}}\le {\tilde{y}}_{\mathrm{o}}^{\top }{\tilde{v}}_{\mathrm{o}}-\Vert {\tilde{y}}_{\mathrm{o}}{\Vert }^2$$

$$A=\left[\begin{array}{cc}{A}_1& A_2^{\top }\\ A_2& A_3\end{array}\right],A_1\in {\mathbb{R}}^{n_1\times n_1},A_3\in {\mathbb{R}}^{n_2\times n_2}.$$

$$\dot{\xi }$$

$$u$$

$$P:=M{A}_1+A_{1}M-2\kappa M>0.$$

$$\dot{x}$$

$$\dot{\xi }$$

$$x^{\ast }$$

$${\xi }^{\ast }$$

$$\dot{x}$$

$$\dot{\xi }$$

$$\bar{A}:=A-\kappa B{B}^{\top },{\bar{k}}_{\mathrm{P}}:=k_{\mathrm{P}}+\kappa ,\bar{F}:=\left[\begin{array}{c}-A_2^{\top }{A}_3^{-1}\\ I_{n_2}\end{array}\right]B_c^{\top }.$$

Full text

I Introduction

Stimulated by strong needs for reducing energy consumption of buildings, smart building energy management algorithms have been developed both in industry and academia. In particular, about half of the current consumption is known to be occupied by heating, ventilation, and air-conditioning (HVAC) systems, and a great deal of works have been devoted to HVAC optimization and control [1]. In this paper, we address the issue based on a novel approach combining optimization and physical dynamics.

Interplays between optimization and physical dynamics have been most actively studied in the field of Model Predictive Control (MPC), which has also been applied to building HVAC control [1]–[6]. While the MPC approach regards the optimization process as a static map from physical states to optimal inputs, another approach to integrating optimization and physical dynamics is presented in [7]–[11] mainly motivated by power grid control. There, the solution process of the optimization is viewed as a dynamical system, and the combination of optimization and physical dynamics is regarded as an interconnection of dynamical systems. The benefits of the approach relative to MPC are as follows. First, the approach allows one to avoid complicated modeling and prediction of factors hard to know in advance, while MPC needs their models to predict future system evolutions. Second, since the entire system is a dynamical system, its stability and performance are analyzed based on unifying dynamical system theory.

In this line of works, Shiltz et al. [7] addresses smart grid control, and interconnects a dynamic optimization process and a locally controlled grid dynamics. The entire process is then demonstrated through simulation. The authors of [8, 9] incorporate the grid dynamics into the optimization process by identifying the physical dynamics with a subprocess of seeking the optimal solution. A scheme to eliminate structural constraints required in [8, 9] is presented by Zhang et al. [10] while instead assuming measurements of disturbances. A similar approach is also taken for power grid control in Stegink et al. [11].

In this paper, we address integrated design of optimization and physical dynamics for HVAC control based on passivity, where we regard the optimization process as a dynamical system similarly to [7]–[11]. Interconnections of such dynamic HVAC optimization with a building dynamics are partially studied in [12], where temperature data for all zones and their derivatives in the physics side are fed back to the dynamic optimization process to recover the disturbance terms. However, such data are not always available in practical systems. We thus present an architecture relying only on temperature data of a subgroup of zones with HVAC systems.

The contents of this paper are as follows. A thermal dynamic model with unmeasurable disturbances and an associated optimization problem are first presented. We then formulate an optimization dynamics based on the primal-dual gradient algorithm [13]. The designed dynamics is then proved to be strictly passive from a transformed disturbance estimate to an estimated optimal room temperature. We next design a controller so that the actual room temperature tracks a given reference, and produces a disturbance estimate. Then, the physical dynamics is proved to be passivity-short from the reference to the disturbance estimate. From these two passivity-related results, we then interconnect the optimization and physical dynamics, and prove convergence of the actual room temperature to the optimal solution defined for the unmeasurable actual disturbance. Finally, the presented algorithm is demonstrated through simulation.

II Problem Settings

II-A Preliminary

In this section, we introduce the concept of passivity. Consider a system with a state-space representation

[TABLE]

where $x(t)\in{\mathbb{R}}^{n}$ is the state, $u(t)\in{\mathbb{R}}^{p}$ is the input and $y(t)\in{\mathbb{R}}^{p}$ is the output. Then, passivity is defined as below.

Definition 1

The system (1) is said to be passive if there exists a positive semi-definite function $S:{\mathbb{R}}^{n}\to{\mathbb{R}}_{+}:=[0,\infty)$, called storage function, such that

[TABLE]

holds for all inputs $u:[0,t]\to{\mathbb{R}}^{p}$, all initial states $x(0)\in{\mathbb{R}}^{n}$ and all $t\in{\mathbb{R}}^{+}$. In the case of the static system $y=\varphi(u)$, it is passive if $y^{T}u=\varphi^{T}(u)u\geq 0$ for all $u\in{\mathbb{R}}^{p}$. The system (1) is also said to be output feedback passive with index $\varepsilon>0$ if (2) is replaced by

[TABLE]

If the right-hand side of (3) is changed as

[TABLE]

with $\varepsilon>0$, the system is said to be passivity-short, and then $\varepsilon$ is called impact coefficient [18].

II-B System Description

In this paper, we consider a building with multiple zones $i=1,2,\dots,n$. The zones $i=1,2,\dots,n$ are divided into two groups: The first group consists of zones equipped with VAV (Variable Air Volume) HVAC systems whose thermal dynamics is assumed to be modeled by the RC circuit model [2]–[6], [12] as

[TABLE]

where $T_{i}$ is the temperature of zone $i$, $m_{i}$ is the mass flow rate at zone $i$, $T^{\rm a}$ is the ambient temperature, $T^{\rm s}_{i}$ is the air temperature supplied to zone $i$, which is treated as a constant throughout this paper, $q_{i}$ is the heat gain at zone $i$ from external sources like occupants, $C_{i}$ is the thermal capacitance, $R_{i}$ is the thermal resistance of the wall/window, $R_{ij}$ is the thermal resistance between zone $i$ and $j$ and $a_{i}$ is the specific heat of the air.

The second group is composed of other spaces such as walls and windows whose dynamics is modeled as

[TABLE]

Rooms not in use can be categorized into this group. Without loss of generality, we assume that $i=1,2,\dots,n_{1}$ belong to the first group and $i=n_{1}+1,\dots,n\ (n_{2}:=n-n_{1})$ belong to the second. Remark that the system parameters $C_{i}$, $R_{i}$, and $R_{ij}$ can be identified using the toolbox in [14].

The collective dynamics of (5) and (6) is described as

[TABLE]

where $T$, $q$ and $m$ are collections of $T_{i}\ (i=1,2,\dots,n)$, $q_{i}\ (i=1,2,\dots,n_{1})$, and $m_{i}\ (i=1,2,\dots,n_{1})$ respectively. The matrices $C$ and $R$ are diagonal matrices with diagonal elements $C_{i}\ (i=1,2,\dots,n)$ and $\frac{1}{R_{i}}\ (i=1,2,\dots,n)$, respectively. The matrix $L$ describes the weighted graph Laplacian with elements $\frac{1}{R_{ij}}$, $G(T)\in{\mathbb{R}}^{n_{1}\times n_{1}}$ is a block diagonal matrix with diagonal elements equal to $a_{i}(T^{\rm s}_{i}-T_{i})\ (i=1,2,\dots,n_{1})$, ${\bf 1}$ is the $n$-dimensional real vector whose elements are all 1, and $B=[I_{n_{1}}\ 0]^{\top}\in{\mathbb{R}}^{n\times n_{1}}$.

We next linearize the model at around an equilibrium as

[TABLE]

where $\delta T$, $\delta m$, $\delta T^{\rm a}$ and $\delta q$ describe the errors from the equilibrium states and inputs and $\bar{U}\in{\mathbb{R}}^{n\times n}$ is a diagonal matrix whose diagonal elements are $\bar{m}_{1},\dots,\bar{m}_{n_{1}},0,\dots,0$, where $\bar{m}_{i}$ is the $i$-th element of the equilibrium input $\bar{m}$.

Using the variable transformations

[TABLE]

(6) is rewritten as

[TABLE]

where $A:=C^{-1/2}(R+L+\bar{U})C^{-1/2}$, $x_{1}\in{\mathbb{R}}^{n_{1}}$ and $x_{2}\in{\mathbb{R}}^{n_{2}}$. Remark that the matrix $A$ is positive definite [19].

From control engineering point of view, $x$ is the system state, $u$ is the control input, and $w_{\rm q}$ and $w_{\rm a}$ are disturbances. We suppose that $w_{\rm a}$ is measurable as well as $x_{1}$. Meanwhile, it is in general hard to measure the heat gain $w_{\rm q}$.

II-C Optimization Problem

Regarding the above system, we formulate the optimization problem to be solved as follows:

[TABLE]

The variables $z_{x}:=[z_{x1}^{\top}\ z_{x2}^{\top}]^{\top}\ (z_{x1}\in{\mathbb{R}}^{n_{1}},z_{x2}\in{\mathbb{R}}^{n_{2}})$, and $z_{u}\in{\mathbb{R}}^{n_{1}}$ correspond to the zone temperature $x$ and mass flow rate $u$ after the transformation (8). The parameters $d_{\rm q}\in{\mathbb{R}}^{n_{1}}$ and $d_{\rm a}\in{\mathbb{R}}^{n}$ are DC components of the disturbances $w_{\rm q}$ and $w_{\rm a}$, respectively. These variables are coupled by (10c) which describes the stationary equation of (9).

Throughout this paper, we assume the following assumption.

Assumption 1

The problem (10) satisfies the following properties: (i) $f:{\mathbb{R}}^{n_{1}}\to{\mathbb{R}}$ is convex and its gradient is locally Lipschitz, (ii) every element of the constraint function $g:{\mathbb{R}}^{n_{1}}\to{\mathbb{R}}^{c}$ is convex and its gradient is locally Lipschitz, and (iii) there exists $z_{u}\in{\mathbb{R}}^{n_{1}}$ such that $g(z_{u})<0$.

The first term of (10a) evaluates the human comfort, where $h\in{\mathbb{R}}^{n_{1}}$ is the collection of the most comfortable temperatures for occupants in each zone, which might be determined directly by occupants in the same way as the current systems, or computed using human comfort metrics like PMV (Predicted Mean Vote). The quadratic function for the error is commonly employed in the MPC papers [3, 22]–[24] and [12]. Note that it is common to put weights on each element of $z_{x1}-h$ to give priority to each zone, but this can be done by appropriately scaling each element of $z_{u}$ in the function $f$. This is why we take (10a).

The function $f$ is introduced to reduce power consumption. The papers [3, 22]–[24] simply take a linear or quadratic function of control efforts as such a function and then Assumption 1(i) is trivially satisfied. Also, as mentioned in [5, 12], the power consumption of supply fans is approximated by the cube of the sum of the mass flow rates, which also satisfies Assumption 1(i). A simple model of consumption at the cooling coil is given by the product of the mass flow rate $m_{i}$ and $|T_{i}^{\rm s}-T^{\rm a}|$ [25], which also belongs to the intended class 111For simplicity, we skip the dependence of $d_{\rm a}$ on $f$, but subsequent results are easily extended to the case that $f$ depends on $d_{\rm a}$.. The constraint function $g:{\mathbb{R}}^{n_{1}}\to{\mathbb{R}}^{c}$ reflects hardware constraints and/or an upper bound of the power consumption. For example, the constraints in [12] are reduced to the above form.

The objective of this paper is to design a controller so as to ensure convergence of the actual room temperature $x_{1}$ to the optimal room temperature $z_{x1}^{*}$, the solution to (10), without direct measurements of $w_{\rm q}$.

III Optimization Dynamics

III-A Optimization Dynamics

In this subsection, we present a dynamics to solve the above optimization problem. Before that, we eliminate $z_{x}$ from (10) using (10c). Then, the problem is rewritten as

[TABLE]

with $\bar{h}=ABh$. It is easy to confirm from positive definiteness of $A$ that the cost function of (11) is strongly convex. In this case, (11) has the unique optimal solution, denoted by $z^{*}_{u}$, and it satisfies the following KKT conditions [26].

[TABLE]

where $M=B^{\top}A^{-1}B$, $N:=B^{\top}A^{-1}BB^{\top}A^{-1}$, and $\circ$ represents the Hadamard product. Since (11) is essentially equivalent to (10), the solution to (11) also provides a solution for the problem (10). Precisely, if we define $z_{x}^{*}\in{\mathbb{R}}^{n}$ as $z_{x}^{*}:=A^{-1}(Bz_{u}^{*}+Bd_{\rm q}+d_{\rm a})$, the pair $z^{*}:=[(z_{x}^{*})^{\top}\ (z_{u}^{*})^{\top}]^{\top}$ is a solution to (10). In the sequel, we also use the notation $z_{x}^{*}:=[(z_{x1}^{*})^{\top}\ (z_{x2}^{*})^{\top}]^{\top}$ ($z_{x1}^{*}\in{\mathbb{R}}^{n_{1}}$ and $z_{x2}^{*}\in{\mathbb{R}}^{n_{2}}$). It is then easy to confirm

[TABLE]

Given $d_{\rm q}$ and $d_{\rm a}$, it would be easy to solve (12). However, in the practical applications, it is desired that $d_{\rm q}$ and $d_{\rm a}$ are updated in real time according to the changes of disturbances. In this regard, it is convenient to take a dynamic solution process of optimization since it trivially allows one to update the parameters in real time. In particular, we employ the primal-dual gradient algorithm [13] as one of such solutions222Another benefit of using the dynamic solution is that it provides a distributed solution when the present results are extended to a more global problem, although it exceeds the scope of this paper. Please refer to [20] for more details on the issue.. However, it is hard to obtain $d_{\rm q}$ since $w_{\rm q}$ is not measurable. We thus need to estimate $d_{\rm q}$ from the measurements of physical quantities. This motivates us to interconnect the physical dynamics with the optimization dynamics.

Taking account of the above issues, we present

[TABLE]

where $\hat{z}_{u}$ and $\hat{\lambda}$ are estimates of $z_{u}^{*}$ and $\lambda^{*}$ respectively, and $\alpha>0$. The notation $[b]^{+}_{a}$ for real vectors $a,b$ with the same dimension provides a vector whose $l$-th element, denoted by $([b]^{+}_{a})_{l}$, is given by

[TABLE]

where $a_{l},b_{l}$ are the $l$-th element of $a,b$, respectively. Note that (14) is different from the primal-dual gradient algorithm for (12) in that the term $d_{\rm a}$ is replaced by the measurement $w_{\rm a}$, and $d_{\rm q}$ is replaced by its estimate $\hat{d}_{\rm q}$ whose production will be mentioned later. The system is illustrated in Fig. 1.

III-B Passivity Analysis for Optimization Dynamics

Hereafter, we analyze passivity of the above optimization process assuming that $w_{\rm a}$ is constant. In this case, $w_{\rm a}\equiv d_{\rm a}$ holds. In practice, the disturbance $w_{\rm a}$, namely the ambient temperature $T^{\rm a}$, is time-varying but the following results are applied to the practical case if $w_{\rm a}$ is approximated by a piecewise constant signal, which is fully expected since $T^{\rm a}$ usually varies slowly.

Under the above assumption, we define the output $\nu:=-M^{2}(\hat{z}_{u}+\hat{d}_{\rm q})$ for (14). We then have the following lemma.

Lemma 1

Consider the system (14) with $w_{\rm a}\equiv d_{\rm a}$ and $\hat{\lambda}(0)\geq 0$. Then, under Assumption 1, it is passive from $\tilde{d}_{\rm q}:=\hat{d}_{\rm q}-d_{\rm q}$ to $-\tilde{\nu}$, where $\tilde{\nu}:=\nu-\nu^{*}$ and $\nu^{*}:=-M^{2}(z^{*}_{u}+d_{\rm q})$.

Proof.

See Appendix A. ∎

Let us next transform the output $\nu$ to

[TABLE]

Comparing (13) and the above definition of $y_{\rm o}$, the signal $y_{\rm o}$ is regarded as an estimate of $z^{*}_{x1}$. We also define $v_{\rm o}:=M\hat{d}_{\rm q}$ and $v_{\rm o}^{*}:=Md_{\rm q}$. Then, we can prove the following lemma.

Lemma 2

Consider the system (14) with $w_{\rm a}\equiv d_{\rm a}$ and $\hat{\lambda}(0)\geq 0$. Then, under Assumption 1, it is output feedback passive from $\tilde{v}_{\rm o}:=v_{\rm o}-v_{\rm o}^{*}$ to $\tilde{y}_{\rm o}:=y_{\rm o}-y_{\rm o}^{*}$ with index $1$.

Proof.

It is easy to see from $\nu^{*}=-M^{2}(z_{u}^{*}+d_{\rm q})$ and (13) that $y_{\rm o}^{*}=z^{*}_{x1}$. From (34), $\tilde{y}_{\rm o}=-M^{-1}\tilde{\nu}$ and $\tilde{v}_{\rm o}=M\tilde{d}_{\rm q}$, we have the following inequality.

[TABLE]

Integrating this in time completes the proof. ∎

IV Physical Dynamics

In this section, we design a physical dynamics and prove its passivity. A passivity-based design for the model (9) is presented in [21], but we modify the control architecture in order to interconnect it with the optimization dynamics.

IV-A Controller Design

In this subsection, we design a controller to determine the input $u$ so that $x_{1}$ tracks a reference signal $r$. Here we assume the following assumption, where

[TABLE]

Assumption 2

The matrix $MA_{1}+A_{1}M$ is positive definite.

This property does not always hold for any positive definite matrices $A_{1}$ and $M$, but it is expected to be true in many practical cases since the diagonal elements tend to be dominant both for $A_{1}$ and $M=(A_{1}-A_{2}^{\top}A_{3}^{-1}A_{2})^{-1}$ in this application [12]. Note that this assumption holds in the full actuation case ($A_{1}=A,M=A^{-1}$).

Inspired by the fact that many existing systems employ Proportional-Integral (PI) controllers (with logics) as the local controller, we design the following controller adding reference and disturbance feedforward terms.

[TABLE]

where $k_{\rm P}>0,\ k_{\rm I}>0$ and $F:=[-I_{n_{1}}\ A_{2}^{\top}A_{3}^{-1}]$ and $I_{n_{1}}$ is the $n_{1}$-by-$n_{1}$ identity matrix. The feedforward gain $\kappa>0$ is selected so that

[TABLE]

Such a $\kappa$ exists under Assumption 2.

Substituting (19) into (9) yields

[TABLE]

For the system, we have the following lemma.

Lemma 3

The steady states $x^{*}$ and $\xi^{*}$ of (21) for $r\equiv r^{*}$, $w_{\rm a}\equiv d_{\rm a}$ and $w_{\rm q}\equiv d_{\rm q}$ are given as follows.

[TABLE]

Equation (21) is now rewritten as

[TABLE]

where

[TABLE]

Remark that, at the steady state, the variable $\zeta$ is equal to

[TABLE]

IV-B Passivity Analysis for Physical Dynamics

In this subsection, we analyze passivity of the system (23) assuming that $w_{\rm q}$ and $w_{\rm a}$ are constant. The case of the time varying $w_{\rm q}$ will be treated in the end of the next section.

Choose $\zeta$ as the output and prove passivity as follows.

Lemma 4

Consider the system (23) with $r\equiv r^{*}$, $w_{\rm a}\equiv d_{\rm a}$ and $w_{\rm q}\equiv d_{\rm q}$. Then, under Assumption 2, the system is passive from $\tilde{r}:=r-r^{*}$ to $\tilde{\zeta}:=\zeta-\zeta^{*}$.

Proof.

See Appendix B. ∎

Remark that this lemma holds regardless of the value of $r^{*}$.

To extract the term $Md_{\rm q}$ from (24), we define the output

[TABLE]

and $v_{\rm p}:=r$ and $v_{\rm p}^{*}:=z_{x1}^{*}$. Then, we have the following.

Lemma 5

Consider the system (23) with $r\equiv r^{*}$, $w_{\rm a}\equiv d_{\rm a}$ and $w_{\rm q}\equiv d_{\rm q}$. Then, under Assumption 2, the system from $\tilde{v}_{\rm p}:=v_{\rm p}-v_{\rm p}^{*}$ to $\tilde{y}_{\rm p}:=y_{\rm p}-y_{\rm p}^{*}$ is passivity-short with the impact coefficient $1-\kappa\sigma$, where $\sigma>0$ is the minimal eigenvalue of $M$.

Proof.

If we take $r^{*}=z_{x1}^{*}$, we have $\tilde{y}_{\rm p}=\tilde{\zeta}-K\tilde{v}_{\rm p}$. Substituting this into (41) yields

[TABLE]

This completes the proof. ∎

V Interconnection of Optimization and Physical Dynamics

Let us interconnect the optimization dynamics (14) and physical dynamics (23). Remark that the stationary value of $y_{\rm o}$, $y_{\rm o}^{*}=z_{x1}^{*}$, is equivalent to that of $v_{\rm p}$, $v_{\rm p}^{*}=z_{x1}^{*}$. Also, $v_{\rm o}^{*}=-y_{\rm p}^{*}$ holds. Inspired by these facts, we interconnect these systems via the negative feedback as $v_{\rm o}=-y_{\rm p},\ v_{\rm p}=y_{\rm o}$. We then have the following main result of this paper.

Theorem 1

Suppose that $\hat{\lambda}(0)\geq 0$, $w_{\rm a}\equiv d_{\rm a}$ and $w_{\rm q}\equiv d_{\rm q}$. Then, if Assumptions 1 and 2 hold, the interconnection of (14) and (23) via $v_{\rm o}=-y_{\rm p},\ v_{\rm p}=y_{\rm o}$ ensures that $x_{1}\to z_{x1}^{*}$.

Proof.

Define $S:=S_{\rm o}+S_{\rm p}$. Then, combining (18), (27) and $v_{\rm o}=-y_{\rm p},\ v_{\rm p}=y_{\rm o}$ yields

[TABLE]

This means that both of $\tilde{y}_{\rm o}=y_{\rm o}-z_{x1}^{*}$ and $\tilde{y}_{\rm o}-\tilde{x}_{1}=y_{\rm o}-x_{1}$ belong to class ${\mathcal{L}_{2}}$. Since $S$ is positive definite, all of the state variables $\hat{z}_{u},\hat{\lambda},x$ and $\xi$ belong to ${\mathcal{L}}_{\infty}$. From (14), $\dot{\tilde{y}}_{\rm o}=-M^{-1}\dot{\nu}=M\dot{\hat{z}}_{u}$ is bounded. Also, (21) means that $\dot{\tilde{x}}=\dot{x}$ is bounded and hence $\dot{y}_{\rm o}-\dot{x}_{1}$ is bounded. Thus, invoking Barbalat’s lemma, we can prove $y_{\rm o}-z_{x1}^{*}\to 0,\ \ y_{\rm o}-x_{1}\to 0$, which means $x_{1}\to z_{x1}^{*}$. This completes the proof. ∎

Lyapunov stability of the desirable equilibrium, tuple of $z_{u}^{*},\lambda^{*},x^{*}$ and $\xi^{*}$, is also proved in the above proof.

It is to be emphasized that the optimal solution is dependent on the unmeasurable disturbance. Nevertheless, convergence to the solution is guaranteed owing to the feedback path from physics to optimization.

The above results are obtained assuming that both of $w_{\rm q}$ and $w_{\rm a}$ are constant. This is likely valid for $w_{\rm a}$ since the ambient temperature is in general slowly varying. However, the heat gain $w_{\rm q}$ may contain high frequency components. To address the issue, we decompose the signal $w_{\rm q}$ into the DC components $d_{\rm q}$ and others $\tilde{w}_{\rm q}$ as $w_{\rm q}=d_{\rm q}+\tilde{w}_{\rm q}$. We also assume that $\tilde{w}_{\rm q}$ belongs to an extended ${\mathcal{L}}_{2}$ space [16]. The following corollary then holds, which is proved following the proof procedure of the well-known passivity theorem [16, 17] and using the fact that the right-hand side of (28) is upper bounded by $-\frac{\kappa k_{\rm P}\sigma}{\kappa+k_{\rm P}}\|\tilde{x}_{1}\|^{2}$.

Corollary 1

Suppose that $\hat{\lambda}(0)\geq 0$, $w_{\rm a}\equiv d_{\rm a}$ and $w_{\rm q}=d_{\rm q}+\tilde{w}_{\rm q}$. Then, if Assumptions 1 and 2 hold, the interconnected system (14), (23) and $v_{\rm o}=-y_{\rm p},\ v_{\rm p}=y_{\rm o}$ from $\tilde{w}_{\rm q}$ to $\tilde{x}_{1}=B^{\top}\tilde{x}=x_{1}-z^{*}_{x1}$ has a finite ${\mathcal{L}}_{2}$ gain.

We give some remarks on the present architecture.

Figs. 1 and 2 are oriented by theoretical analysis, but the implementation does not need to follow the information processing in the figures. Indeed, the interconnected system is equivalently transformed into Fig. 3. If we let the operations shaded by dark gray be executed in the high-level controller, the low-level controller can be implemented in a decentralized fashion similarly to the existing systems.

In Fig. 3, both of the high-level and low-level controller with the physical dynamics are biproper and hence a problem of algebraic loops can occur. This however does not matter in practice since the information transmissions between high- and low-level processes usually suffer from possibly small delays. Although the high-level controller itself contains an algebraic loop, it is easily confirmed that the loop can be solved by direct calculations of the algebraic constraint.

The transfer function from $w_{\rm q}$ to the disturbance estimate $\hat{d}_{\rm q}$, roughly speaking, is almost the same as the complementary transfer function and hence only the low frequency components are provided by the physical dynamics. This is why $\hat{d}_{\rm q}$ is regarded as an estimate of the DC component of $w_{\rm q}$. The cutoff frequency of the disturbance can be in principle tuned by $k_{\rm P}$ and $k_{\rm I}$, but, once a closed-loop system is designed, the cutoff is also automatically decided. It is however not always a drawback at least qualitatively. Actually, even if optimal solutions reflecting much faster disturbance variations are provided, the physical states cannot respond to variations faster than the bandwidth.

It is a consequence of the internal model control and constant disturbances that the actual disturbance $w_{q}$ is correctly estimated. A control architecture based on a similar concept is presented in Section VII of Stegink et al. [11]. However, it is clear that the problem (10) does not meet the structural constraints assumed in [11] and hence the architecture in [11] cannot be directly applied to our problem.

Zhang et al. [12] present another kind of interconnection between physical and optimization dynamics based on a quasi-disturbance feedforward, where the disturbance is computed by state measurements and their derivatives, and then fed back to the optimization dynamics. The differences of the present scheme from [12] are listed as follows:

The approach of [12] requires the measurements of state variables. If the states include temperatures of windows and/or walls, its technological feasibility may be problematic or at least increases the system cost. On the other hand, our approach needs only $x_{1}$ which is usually measurable.

Since there is no sensor to measure $\dot{x}$, it has to be computed using the difference approximation, which provides approximation errors. Meanwhile, the present approach does not need such an approximation.

In [12], the difference approximation errors together with sensor noises and high frequency components of the disturbances are directly sent to the optimization dynamics, which may cause fluctuations for the output and internal variables in the optimization process unless it is carefully designed in the sense of the noise reduction. Adding a low-pass filter to the computed disturbance might eliminate these undesirable factors. However, the filter is not designed independently of stability of the entire system in the presence of uncertainties in $\dot{x}$ and the system model since, in this case, the quasi-feedforward system becomes a feedback system and the filter is included into the loop. Meanwhile, the noises are automatically rejected by the physical dynamics in our algorithm.

VI Simulation

In this section, we demonstrate the presented control architecture through simulation. For this purpose, we build a building on 3D modeling software SketchUp (Trimble Inc.), which contains three rooms ($n_{1}=3$) and other 46 zones ($n_{2}=46$) including walls, ceilings, and windows. The building model is then installed into EnergyPlus [27] in order to simulate the evolution of zone temperatures. Then, using the acquired data, we identify the model parameters in (5) and (6) via BRCM toolbox [14].

We next specify the optimization problem (10). All the elements of $h$ are set to $22C^{1/2}{}^{\circ}$C and we take $f(z_{u})=150\|z_{u}\|^{2}$. The constraints are also selected as $|z_{ui}|\leq 0.61\ i=1,2,3$ and $\sum_{i=1}^{3}|z_{ui}|\leq 1.25$, where $z_{ui}$ is the $i$-th element of $z_{u}$. Collecting these constraints, we define the function $g$. However, since it turns out that directly using $g(z_{u})\leq 0$ has a response speed problem in penalizing the constraint violation in the primal-dual algorithm, we instead take the constraint $\theta g(z_{u})\leq 0$ with $\theta=15$, which does not essentially change the optimization problem.

In the simulation, we take the feedback gains $k_{\rm P}=6.0\times 10^{-2}$ and $k_{\rm I}=1.0\times 10^{-3}$, and $\kappa=1.0\times 10^{-3}$, which are tuned so that the peak gain of $\sigma$-plot from $r$ to $x_{1}$ is smaller than the well-known criterion. It is then confirmed that Assumption 2 is satisfied.

In the following simulation, we use the disturbance data shown in Fig. 5. Here, we compare the results with the ideal case that the disturbance $w_{\rm q}$ is directly measurable. In this case, the feedback path from the physical dynamics to optimization is not needed and hence we take the cascade connection from optimization to physics. It is to be noted that it is hard to implement this in practice.

The trajectories of the estimated heat gains, the estimated optimal room temperatures and temperatures $T_{1},T_{2},T_{3}$ are illustrated by solid curves in Fig. 5. The dotted lines with light colors show the trajectories using the above disturbance feedforward scheme. We see from the top left figures that the presented algorithm almost correctly estimate the disturbance $w_{q}$. It is also observed from the right fine-scale figure that high frequency signal components are filtered out in the case of our algorithm. The trajectories in the bottom-left figure sometimes get far from 22*∘*C since the constraint $\sum_{i=1}^{3}|z_{ui}|\leq 1.25$ gets active during the periods due to the high ambient temperature and heat gains. We see from the bottom figures that the response of the present method to the constraint violations is slower than the disturbance feedforward scheme because high frequency components of the disturbance are filtered out by the physical dynamics.

Let us next show that the high sensitivity of disturbance feedforward can cause another problem. Here, small noises are added to the heat gain $w_{\rm q}$ at every 20s, whose absolute value is upper bounded by $1.0\times 10^{-3}$. We then run the above two algorithms. The resulting responses are illustrated in Fig. 7. It is observed from the top-left figure that the noise is filtered out in the present algorithm, and its effects do not appear on the estimates. Accordingly, the trajectories of the estimated optimal and actual room temperatures are almost the same as Fig. 5. Meanwhile, the disturbance feedforward approach suffers significant effects from the noise. The trajectories of the bottom-left get smaller than 22*∘*C. Namely, the fluctuations are caused by the undesirable over- and undershoots. Although the trajectories of the actual temperatures get smooth, this behavior of the reference is not desirable from an engineering point of view. Adding a low-pass filter or reducing the gain of the optimization dynamics would eliminate the oscillations but it spoils the advantage, namely response speed. It is to be noted that if the disturbance feedforward is implemented using the recovery technique in [12], the low-pass filter is not designed independently of system stability as stated in Section V.

If the response speed of the present algorithm in Fig. 5 is still problematic, it can be accelerated by tuning the scaling factor $\theta$. The results for $\theta=50$ are shown in Fig. 7, where it is observed that almost the same speed as the disturbance feedforward in Fig. 5 is achieved by the present algorithm. Simulation for a larger-scale system with more practical settings is left as a future work of this paper.

VII Conclusion

In this paper, we presented a novel combined optimization and control algorithm for HVAC control of buildings. We designed a primal-dual algorithm-based optimization dynamics and a local physical control system, and proved the system properties related to passivity. We then interconnected the optimization and physical dynamics, and proved convergence of the room temperatures to the optimal ones. We finally demonstrated the present algorithms through simulation.

Appendix A Proof of Lemma 1

Lemma 6

Consider the system (14b) with $\hat{\lambda}(0)\geq 0$. Then, under Assumption 1, it is passive from $\tilde{z}_{u}=\hat{z}_{u}-z_{u}^{*}$ to $\tilde{p}=p-p^{*}$ with $p^{*}:=\nabla g(z^{*}_{u})\lambda^{*}$.

Proof.

Define the energy function $U:=\frac{1}{2}\|\hat{\lambda}-\lambda^{*}\|^{2}$. Then, following the same procedure as [15], we have

[TABLE]

where the notation $D^{+}$ represents the upper Dini derivative. Integrating this in time completes the proof. ∎

We next consider (14a). Now, replace $-M^{2}(\hat{z}_{u}+\hat{d}_{\rm q})-p$ by an external input $\mu$ and consider the system

[TABLE]

Then, we have the following lemma.

Lemma 7

Suppose $w_{\rm a}\equiv d_{\rm a}$. Then, under Assumption 1, the system (30) is passive from $\tilde{\mu}:=\mu-\mu^{*}$ to $\tilde{z}_{u}=\hat{z}_{u}-z_{u}^{*}$, where $\mu^{*}:=-M^{2}(z^{*}_{u}+d_{\rm q})-p^{*}$.

Proof.

Subtracting (12a) from (30) under yields

[TABLE]

Now, define $V:=\frac{1}{2\alpha}\|\tilde{z}_{u}\|^{2}=\frac{1}{2\alpha}\|\hat{z}_{u}-z_{u}^{*}\|^{2}$. Then, the time derivative of $V$ along the trajectories of (31) is given by

[TABLE]

From convexity of $f$, $(\hat{z}_{u}-z_{u}^{*})^{\top}(\nabla f(\hat{z}_{u})-\nabla f(z_{u}^{*}))\geq 0$ holds [26]. This completes the proof. ∎

The system (14) is given by interconnecting (14b) and (30) via $\mu=\nu-p$. It is then easy to confirm that $\tilde{\mu}=\tilde{\nu}-\tilde{p}$.

We are now ready to prove Lemma 1. Define $S_{\rm o}:=V+U$. Then, combining (29), (32) and $\tilde{\mu}=\tilde{\nu}-\tilde{p}$, we have

[TABLE]

Since $\tilde{\nu}=-M^{2}(\tilde{z}_{u}+\tilde{d}_{\rm q})$, it follows

[TABLE]

This completes the proof.

Appendix B Proof of Lemma 4

Lemma 8

Under Assumption 2 and (20), the system $(\bar{A}_{3},A_{2}P^{-1/2})$ is stabilizable and $(P^{-1/2}MA_{2}^{\top},\bar{A}_{3})$ is detectable, where $\bar{A}_{3}:=-A_{3}+A_{2}MP^{-1}A_{2}^{\top}$.

Proof.

Define $\Phi_{\rm s}:=P^{1/2}MP^{-1}A_{2}^{\top}$ and $\Phi_{\rm d}:=A_{2}MP^{-1}M^{-1}P^{1/2}$. Then,

[TABLE]

hold and $-A_{3}$ is stable. This completes the proof. ∎

Using Lemma 8, we next prove the following result.

Lemma 9

Consider the system (23a) with $r\equiv r^{*}$, $w_{\rm a}\equiv d_{\rm a}$ and $w_{\rm q}\equiv d_{\rm q}$. Then, under Assumption 2, the system is passive from $\tilde{\zeta}$ to $\tilde{x}_{1}:=B^{\top}\tilde{x}$ with $\tilde{x}:=x-x^{*}$.

Proof.

We first formulate the error system

[TABLE]

where $\tilde{\xi}:=\xi-\xi^{*}$. Take a positive definite matrix $\Psi\in{\mathbb{R}}^{n_{2}\times n_{2}}$ and define $\bar{\Psi}:=\begin{bmatrix}M&0\\ 0&\Psi\end{bmatrix}\in{\mathbb{R}}^{n\times n}$. Then, by calculation, we have

[TABLE]

From Schur complement, under Assumption 2, $\bar{\Psi}\bar{A}+\bar{A}\bar{\Psi}>0$ is equivalent to the following Riccati inequality.

[TABLE]

A positive semi-definite solution $\Psi$ to (37) is shown to exist from Lemma 8. Now, define an energy function $S_{x}:=\frac{1}{2}\tilde{x}^{\top}\bar{\Psi}\tilde{x}$ for the solution $\Psi$ to (37). Then, the time derivative of $S_{x}$ along the trajectories of (35a) is given by

[TABLE]

This completes the proof. ∎

We are now ready to prove Lemma 4. Replace $\tilde{r}-B^{\top}\tilde{x}$ in (35b) by $\tilde{e}$ as

[TABLE]

Define $S_{\xi}:=\frac{1}{2k_{\rm I}}\tilde{\xi}^{\top}M\tilde{\xi}$. Then, the time derivative of $S_{\xi}$ along the trajectories of (39) is given as

[TABLE]

Define $S_{\rm p}:=S_{x}+S_{\xi}$. Then, from (38) and (40), we have

[TABLE]

This completes the proof.