#  How to Secure Valid Quantizations

John R<EMAIL_ADDRESS>
Department of Physics and Department of Mathematics
University of Florida, Gainesville, FL 32611-8440

###### Abstract

Classical mechanics involves position and momentum variables that must be
special coordinates chosen to promote to suitable quantum operators. Since
classical variables may be broadly chosen, only unique variables should be
chosen. We will outline how the favored variables and their suitable quantum
operators is guaranteed to assure a truly valid quantization. Invalid
quantizations may be mistaken for valid ones, which then leads to incorrect
physics. Besides particle examples, there is also a brief run-through for
fields and gravity.

##  1 The Special Rules for a Valid
Canonical Quantization

In our search for the proper classical variables, $p\;\&\;q$, and thus the
proper variables to promote to quantum operators, $P\;\&\;Q$, we introduce
suitable coherent states. They are defined as

$\displaystyle|p,q\rangle\equiv
e^{-iq\,P/\hbar}\,e^{ip\,Q/\hbar}\,|\omega\rangle\;,$ (1)

where we choose the normalized fiducial vector, $|\omega\rangle$, so that
$[Q+iP/\omega]\,|\omega\rangle=0$.

We first examine two important Hamiltonian functions given by

$\displaystyle H(p,q)=\langle p,q|{\cal
H}(P,Q))|p,q\rangle=\langle\omega|{\cal H}(P+p,Q+q)|\omega\rangle$
$\displaystyle\hskip 33.99998pt={\cal H}(p,q)+{\cal{O}}(\hbar;p,q)\>.$ (2)

Although $p\times q$ has the dimensions of $\hbar$ it follows that $|QP-
PQ|=\hbar$ as well, and we may assume that $Q\;\&\;P$ can each be proportional
to $\sqrt{\hbar}$. Adopting a general polynomial style, with $|{\cal
H}(p,q)|<\infty$, then the term ${\cal H}(p,q$ has no $\hbar$, while the last
term in (2) contains all the factors with $\hbar$. We now let
$\hbar\rightarrow 0$, and it follows that $H(p,q)={\cal
H}(p,q)\Rightarrow{\cal H}(P,Q)$, where we add what becomes the proper quantum
Hamiltonian. That completes rule number one for CQ.222In what follows, we
often use CQ for canonical quantization, and we often use AQ for the second
procedure, called affine quantization.

Rule number two, following Dirac [1], will fix what kind of classical
variables to choose for CQ. For that purpose, we again use the coherent
states, as shown by

$\displaystyle d\sigma^{2}_{CQ}=2\hbar[\>|\\!|\;d|p,q\rangle|\\!|^{2}-|\langle
p,q|\;d|p,q\rangle|^{2}\,]=\omega^{-1}\,dp^{2}+\omega\,dq^{2}\;,$ (3)

which requires a flat, Cartesian phase space of these two different variables,
which may also be called a constant zero curvature. Observe, that this
relation relies on $\hbar>0$, and is therefore a semi-classical relation. The
second factor of the coherent states is included to ensure that phase factors,
such as $|p,q:f\rangle=e^{if(p,q)}\;|p,q\rangle$, do not appear in the final
metric. This safety precaution is known as the Fubini-Study metric, which was
introduced in the early 19 hudreds [2].

Clearly, using coherent states, can help provide some tools to create a valid
quantization while using a canonical quantization. In using CQ, it is
important to remember that $-\infty<P=P^{\dagger}\;\&\;Q=Q^{\dagger}<\infty$,
which ensures that the classical variables, $-\infty<p\;\&\;q<\infty$, are
also respected.

###  1.1 Finding the primitive ground state

Solving a simple ground state is elementary for canonical quantization, and we
include this because our next quantization procedures will require a new
formulation.

Let us choose at first an elementary Hilbert space that only involves a finite
coordinate space which is that of a circle of radius $R$. Net we introduce the
operator $P$ and seek a wave function for which $P\;|\alpha\rangle=0$.
Choosing Schrödinger’s representation, the solution is $\langle
x|\alpha\rangle=c$, which is simply a constant. Now, for a Hamiltonian that is
simply $P^{2}\geq 0$, it follows that $P^{2}\;c=0$, and since
$\langle\alpha|P^{2}|\alpha\rangle\geq 0$, the ‘ground state’ for this tiny
problem is the constant $c$. We can also choose $c=1$, in which case
$1=\Pi_{y}\,1$. Briefly stated, the ground state for the kinetic operator,
$P^{2}$, is simply $1$.

Of course, this state will change when we introduce a potential, $V(Q)$, and
perhaps we can even expand the radius $R\rightarrow\infty$ to become the usual
coordinate space.

Next, we introduce a newer quantization procedure, namely. AQ, that is a very
useful tool for many different problems. The next formulation of a ‘simple
ground state’ will then be something, very likely, quite new.

##  2 The Special Rules for a Valid
Affine Quantization

The fundamental difference between AQ and CQ is that, while $-\infty<p<\infty$
remains, we remove a coordinate point, say $q=0$, which then leaves $q<0$ and
$q>0$. We discard $q<0$ and keep only $q>0$. That implies that
$Q=Q^{\dagger}>0$, but that $P^{\dagger}\neq P$, which can lead to
complications. We can fix that by introducing the classical ‘dilation
variable’, $d=pq\Rightarrow D=[P^{\dagger}\,Q+Q\,P]/2=D^{\dagger}$. Note that
where $Q$ is sitting it forces any wave function to vanish at the point
removed; this property also permits $P^{\dagger}\,Q=P\,Q$, which can be very
useful for any calculations.

In addition, we choose to treat $q\;\&\;Q$ as dimensionless to simplify some
equations. That means that both $p\;\&\;P$ and $d\;\&\;D$ have the dimensions
of $\hbar$.

Comparing AQ and CQ, we find that CQ has $[Q,P]=i\hbar 1\\!\\!1$, while AQ has
$[Q,D]=i\hbar\,Q$. This latter equation resembles the Lie algebra for the
affine group, from which our procedures have adopted ‘affine’ for its name. We
are now prepared to introduce an AQ paragraph to exhibit the differences from
those of CQ.

In our search for the proper classical variables, $d\;\&\;q$, which lead to
the proper quantum variables, $D\;\&\;Q$, we introduce suitable coherent
states. They are defined as (Note: $p,q\rightarrow p;g$ in $|p:q\rangle$)

$\displaystyle|p;q\rangle\equiv
e^{ip\,Q/\hbar}\,e^{-i\,\ln(q)\,D/\hbar}\,|\beta\rangle\;,$ (4)

where we choose the normalized fiducial vector, $|\beta\rangle$, so that
$[(Q-1\\!\\!1)+iD/\beta\hbar]\,|\beta\rangle=0$, while $\beta>1/2$. Note that
$-\infty<\ln(q)<\infty$, which involves the natural logarithm so that
$e^{\ln(q)}=q>0$.

We first examine two Hamiltonian-like expressions, denoted by a prime, e.g.,
$H^{\prime}\;\&\;{\cal H}^{\prime}$, given by

$\displaystyle H^{\prime}(pq,q)=\langle p;q|{\cal
H}^{\prime}(D,Q)|p;q\rangle=\langle\beta|{\cal
H}^{\prime}(D+pqQ,qQ)|\beta\rangle$ $\displaystyle\hskip 41.00012pt={\cal
H}^{\prime}(pq,q)+{\cal{O}}^{\prime}(\hbar;pq,q)\;.$ (5)

In addition, we observe that any wave function that is even, i.e.,
$\psi(-x)=\psi(x)$, it follows that
${\textstyle\int}_{-\infty}^{\infty}(1-x)^{p}\,|\psi(x)|^{2}\;dx=1+{\cal{O}}_{p}(\hbar)$,
while the result is only $1$ after $\hbar\rightarrow 0$ for a conventual
classical Hamiltonian is $H=\textstyle{\frac{1}{2}}\,p^{2}+V(|q|)$.

Adopting a general polynomial style, with $|{\cal H}(pq,q)|<\infty$, and the
final term ${\cal{O}}^{\prime}(\hbar;pq,q)\rightarrow 0$ when
$\hbar\rightarrow 0$, it follows that $H^{\prime}(pq,q)={\cal
H}^{\prime}(pq,q)\Rightarrow{\cal H}^{\prime}(D,Q)$, where again, the quantum
Hamiltonian is the correct choice. That is rule number one for AQ.

Rule number two will fix what kind of classical variables to choose to promote
to the valid AQ operators. For that purpose, we again use the coherent states,
as shown here

$\displaystyle
d\sigma^{2}_{AQ}=2\hbar[\>|\\!|\;d|p;q\rangle\,|\\!|^{2}-|\langle
p;q|\;d|p;q\rangle|^{2}\,]=(\beta\hbar)^{-1}\,q^{2}\,dp^{2}+(\beta\hbar)\,q^{-2}\,dq^{2}\;,$
(6)

which has led us to a constant negative curvature surface [3], a very special
phase space of two variables.333The surface of a simple round ball acts as a
constant positive curvature, which is related to spin operators, such as
SO(3). With that, we have fulfilled the three surfaces – positive, zero, and
negative – while the negative curvature cannot be seen like the other two. The
spin coherent state stony is covered in [4]. Observe, this metric relation
also relies on $\hbar>0$, and is not a purely classical result. The second
coherent state factor is included to ensure that no phase factor, such as
$|p;q:g\rangle=e^{ig(p,q)}\;|p;q\rangle$, appear in the metric. As stated
before, this safety precaution is known as the Fubini-Study metric; see [2].

###  2.1 Finding the new primitive ground state

For a system for which $q>0$ we have adopted affine variables for which the
kinetic term becomes $p^{2}=d^{2}/q^{2}\Rightarrow D(Q^{-2})D$. Now we seek
the ground state for this operator. While $P\;1=0$, we find that
$D\;x^{-1/2}=0$, again using Schrödinger’s representation. It follows that an
acceptable wave function, near $x=0$, is like $\gamma(x)=x^{3/2}\,(remainder)$
so that it is continuous and remains so after the first derivative. This
particular expression will be explained further toward the end of the next
section.

Now, still using Schrödinger’s representation, it follows that our ground
state for this kinetic operator is $\Pi_{y}\;y^{-1/2}$, and
$D\;\Pi_{y}\;y^{-1/2}=0$, with an eigenvalue of $0$, as before. As standard
(positive) potential terms are now added, other eigenfunctions will emerge,
and all others have larger eigenvalues. Now, a general wave function may be
written as $\gamma(x)=w(x)\,x^{-1/2}$, where $w(x)=x^{2}\,(remaider)$.
Normalization leads to

$\displaystyle
1={\textstyle\int}\;|\gamma(x)|^{2}\;dx={\textstyle\int}\;|w(x)|^{2}\;\;dx/|x|\;.$
(7)

Clearly, relying on coherent states, has provided fine tools to create valid
quantizations when using CQ or AQ.

##  3 Two Toy Models: One for CQ, One for AQ

The two toy models are first a full-harmonic oscillator with classical
Hamiltonian $H=(p^{2}+q^{2})/2$, while $-\infty<p\;\&\;q<\infty$. The second
toy model is the half-harmonic oscillator with the same Hamiltonian, but while
$-\infty<p<\infty$, now $0<q$. Physically, in the second case, the particles
bounce back off the ‘wall’ at $q=0$ and immediately changes direction, meaning
that $p(t)$ changes sign immediately. That behavior is similar to someone
throwing a tennis ball at a brick wall.

An important item of the harmonic oscillator quantum Hamiltonian, ${\cal
H}=(P^{2}+Q^{2})/2$, which, using CQ, is the eigenvalues ,
$E_{n}=\hbar(n+1/2)$, where $n=0,1,2,3,...$. Note that the eigenvalues are
equally spaced by $\hbar$!

The second model, using AQ and $q>0$, first replaces the classical Hamiltonian
with $H^{\prime}=(d^{2}/q^{2}+q^{2})/2$ which then points to

$\displaystyle{\cal H}^{\prime}=[D(Q^{-2})D+Q^{2}]/2$ $\displaystyle\hskip
15.00002pt=[P^{2}+(3/4)\hbar^{2}/Q^{2}+Q^{2}]/2\>,$ (8)

for which the eigenvalues now have the form $E^{\prime}_{n}=2\hbar(n+1)$, with
$n$ just as before. This implies that it is still equally spaced, now by
$2\hbar$, with a ground state value of $2\hbar$ [5]. Clearly, both of these
two different quantizations are valid.

Although the operators $P^{\dagger}\neq P$, because $Q>0$, the $\hbar$-term in
(8) renders both $P^{\dagger}$ and $(P^{\dagger})^{2}$ to act like $P$ and
$P^{2}$. Now every eigenfunction is of the form
$\psi_{n}(x)=x^{3/2}(polynomial_{n})e^{-x^{2}/2\hbar}$, as they were fully
displayed in [5]. To see what the second derivative does, take a break and
‘solve’ $[-\hbar^{2}\,(d^{2}/dx^{2})+(3/4)\hbar^{2}/x^{2}]\;x^{3/2}=\,?$, for
the interval $0<x<\infty$.

Each toy model fully followed the rules for their quantization procedures to
achieve valid results. Interested readers might try to solve the half-harmonic
oscillator using CQ, and try to prove that it followed the rules; one example
of that effort is presented in [6]. Otherwise, these particular results are
most surely invalid, and will offer only incorrect physics.

This concludes our simple examples with the fact that CQ can truly solve
certain problems, while AQ can truly solve other problems. Difficult problems
in one procedure may find success in the other process. To build many valid
examples just trade the potential $q^{2}/2$ for a general potential
$V(q)\Rightarrow V(Q)$, providing that it obeys $|V(q)|<\infty$.

By just adding a different potential, you are certain to have valid results.
The focus on getting things right is on the kinetic factor, i.e.,
$p^{2}\Rightarrow?$, and much less on the potential, $V(q)$.

##  4 The Benefits of an Affine Quantization for Field Theories

Next, we study some common classical Hamiltonians, such as

$\displaystyle
H(\pi,\varphi)={\textstyle\int}\\{\textstyle{\frac{1}{2}}[\pi(x)^{2}+(\overrightarrow{\nabla}\varphi(x))^{2}+m^{2}\varphi(x)^{2}]+g\,\varphi(x)^{p}\\}\;d^{s}\\!x$
(9)

In a path integration of these models, there are cases where $|\pi(x)|$ and/or
$|\varphi(y)|$ reach $\infty$, while still having $H(\pi,\varphi)<\infty$.

###  4.1 A Brief Overview of
Quantum Field Theory

To begin with, we focus on a feature of mathematics that impacts physics. As
one example, consider the function $f(x)=1/|x|^{1/3}$ in the interval
$-1<x<1$, where $f(0)=\infty$. It follows that
${\textstyle\int}^{1}_{-1}f(x)^{2}\;dx<\infty$, while
${\textstyle\int}_{-1}^{1}f(x)^{4}\;dx=\infty$. Later, we will refer to this
situation as an ‘f-issue’ which involves a field that reaches infinity, and
also is part of an integration that is finite.

Physics is frequently engaged in studying nature, and then having any field
with an f-issue, as in the last paragraph, seems impossible. Stated bluntly,
the strength of any field of nature can not reach infinity. The author accepts
that no field of nature reaches infinity, and we now endeavors to make that
happen in our analysis. For example, a classical Hamiltonian density, $H(x)$,
that describes a part of nature, should not reach infinity for any $x$, and
that fact needs to be part of the mathematics involved.

Presently, the mathematical focus only requires that $\int
H(x)\;d^{s}\\!x<\infty$, and this approach can lead to a nonrenormalizable
behavior if the interaction power of terms, such as for $\varphi^{p}_{n}$ in
(9), when $p\geq 2n/(n-2)$. Let us examine a procedure in which we can
actually favor nature.

###  4.2 A typical model of a covariant scalar field

Our models of interest have the classical Hamiltonians, already introduced
above, and once again is

$\displaystyle
H(\pi,\varphi)={\textstyle\int}\\{\textstyle{\frac{1}{2}}[\pi(x)^{2}+(\overrightarrow{\nabla}\varphi(x))^{2}+m^{2}\varphi(x)^{2}]+g\,\varphi(x)^{p}\\}\;d^{s}\\!x\;.$
(10)

When such a model is quantized, say by a path integration procedure, the
classical Hamiltonian is evaluated by a vast number of functions for both
$\pi(x)$ and $\varphi(x)$. In so doing, fields that can diverge but still
offer finite integrations – like the example of an f-issue in a previous
paragraph – are conventionally introduced along with fields without any
divergencies.

How can we limit the classical fields so that f-issues do not arise? The
answer to that question appears in the next section, and it is much easier
than could have been expected.

Following the simple rules of Sec. 2, we introduce the dilation field
$\kappa(x)=\pi(x)\,\varphi(x)$, along with $\varphi(x)=0$, because if
$\varphi(x)=0$, then $\kappa(x)=0$ and $\pi(x)$ cannot help. But is it
acceptable to remove $\varphi(x)=0$?

Comment: Consider, a rain storm that has a rain quantity=(qu-r) per hour that
is $0<(qu-r)<\infty$. If $(qu-r)=0$, there simply is no rain. If instead, it
was a snow storm, a similar story could be a snow quantity=(qu-s), where
$0<(qu-s)<\infty$. In fact, the physics of both is identical when
$(qu-r)=(qu-s)=0$. So we ignore rain and snow when they are absent.

If $\varphi(x)$ represents a physical particle, say an electron, the field can
remove $\varphi(x)=0$ if the particle is not present. Maybe math may be
concerned, but physics can certainly accept removing any portion of a field
representation when that field represents absolutely nothing.

When using these new classical variables, the classical Hamiltonian, from
(10), becomes

$\displaystyle
H^{\prime}(\kappa,\varphi)={\textstyle\int}\\{\textstyle{\frac{1}{2}}[\kappa(x)^{2}/\varphi(x)^{2}+(\overrightarrow{\nabla}\varphi(x))^{2}+m^{2}\varphi(x)^{2}]+g\,\varphi(x)^{p}\\}\>d^{s}\\!x\;,$
(11)

Already we see that $0<|\varphi(x)|<\infty$ for if $\varphi(x)^{2}=0$ or
$\varphi(x)^{-2}=0$, it means that $\kappa(x)$ fails to offer any value for
$\pi(x)$. It already follows that $|\kappa(x)|<\infty$ and
$0<|\varphi(x)|^{p}<\infty$, and therefore, nonrenormalizability vanishes!
Such field models clearly keeps both sides of $\varphi(x)\neq 0$, i.e., both
$\varphi(x)>0$ and $\varphi(x)<0$. Since there is the gradient term, the field
will appear to be continuous and integrations should not be affected.

To ensure our quantization follows the rules, we introduce suitable coherent
states

$\displaystyle|\pi;\varphi\rangle=e^{i{\textstyle\int}\,\pi(x)\,\hat{\varphi}(x)\,d^{s}\\!x/\hbar}\;e^{-i{\textstyle\int}\,\ln(|\varphi(x)|)\,\hat{\kappa}(x)\,d^{s}\\!x/\hbar}\,|\beta\rangle\;,$
(12)

while the new fiducial vector obeys
$[(|\hat{\varphi}(x)|-1\\!\\!1)+i\,\hat{\kappa}(x)/\beta]\,|\beta\rangle=0$,
and the same label now referees to a new fiducial vector.

It follows that

$\displaystyle H^{\prime}(\kappa,\varphi)=\langle\pi;\varphi|{\cal
H}^{\prime}(\hat{\kappa},\hat{\varphi})|\pi;\varphi\rangle$
$\displaystyle\hskip 40.00006pt=\langle\beta|{\cal
H}^{\prime}(\hat{\kappa}+\pi|\varphi|\,\hat{\varphi},|\varphi|\,\hat{\varphi})|\beta\rangle$
$\displaystyle\hskip 40.00006pt=\langle\beta|{\cal
H}^{\prime}(\hat{\kappa}+\pi\,\varphi\,|\hat{\varphi}|,\varphi\,|\hat{\varphi}|)|\beta\rangle$
$\displaystyle\hskip 40.00006pt={\cal
H}^{\prime}(\kappa,\varphi)+{\cal{O}}(\hbar;\kappa,\varphi)\>,$ (13)

which, as was like the case for particles, we find, if $\hbar\rightarrow 0$,
then $H^{\prime\prime}(\kappa,\varphi)={\cal
H}^{\prime\prime}(\kappa,\varphi)\Rightarrow{\cal
H}^{\prime\prime}(\hat{\kappa},{\hat{g}})$.

To offer an affine quantization for this example, we first introduce the
dilation operator
$\hat{\kappa}(x)=[{\hat{\pi}}(x)^{\dagger}\hat{\varphi}(x)+\hat{\varphi}(x)\hat{\kappa}(x)]/2$
and $\hat{\varphi}(x)\neq 0$. Next, adopting a Schrödinger representation for
the quantum Hamiltonian, we are led to

$\displaystyle{\cal
H}^{\prime}(\hat{\kappa},\varphi)={\textstyle\int}\\{\textstyle{\frac{1}{2}}[\hat{\kappa}(x)\,(\varphi(x)^{-2})\,\hat{\kappa}(x)+(\overrightarrow{\nabla}\varphi(x))^{2}+m^{2}\varphi(x)^{2}]$
$\displaystyle\hskip 140.00021pt+g\,\varphi(x)^{p}\\}\;d^{s}\\!x\>,$ (14)

where it is evident that $\varphi(x)^{2}>0$ and $1/\varphi(x)^{2}>0$, and
therefore $0<|\varphi(x)|^{p}<\infty$, which eliminates nonrenormalizability,
as desired. Finally, we offer Schrödinger’s equation as

$\displaystyle i\hbar\,\partial\,\Psi(\varphi,t)/\partial t={\cal
H}^{\prime}(\hat{\kappa},\varphi)\,\Psi(\varphi,t)\;.$ (15)

As usual, it may be necessary to introduce some version of a regularization
for these equations, but these same equations should point the way to proceed.
To offer some support, we note that although
${\hat{\pi}}(x)^{\dagger}\neq{\hat{\pi}}(x)$ it can be helpful to know that
${\hat{\pi}}(x)^{\dagger}\varphi(x)={\hat{\pi}}(x)\varphi(x)$; see footnote 1.

If the reader has already accepted the expressions for the half-harmonic
oscillator in (8), they may be willing to accept (14) and (15) for this
quantized covariant scalar model as well.444We note that Monte Carlo studies
of the scalar fields $\varphi^{4}_{4}$ and $\varphi^{12}_{3}$ using canonical
quantization have led to “free-results”, as if the interaction term was absent
when it was not. However, using affine quantization has led to “non-free-
results”, when the interaction term makes different results when the coupling
constant changes; see [8 - 12 ].

##  5 Applying Affine Quantization to
Einstein’s Gravity

In order to quantize gravity it is important to render a valid quantization of
the ADM classical Hamiltonian [7]. We first choose our new classical variables
that includes what we will call the dilation field
$\pi^{a}_{b}(x)\equiv\pi^{ac}(x)\,g_{bc}(x)$ (summed by $c$) along with the
metric field $g_{ab}(x)$. We don’t need to impose any restriction on the
metric field because physics already requires that
$ds(x)^{2}=g_{ab}(x)\;dx^{a}\;dx^{b}>0$ provided that $\\{dx^{a}\\}\neq 0$; we
choose that $g_{ab}(x)$ is dimensionless because that will simplify further
equations. The metric field can also be diagonalized by non-physical,
orthogonal matrices, and then it includes only
$g_{11}(x),\;g_{22}(x),\;\&\;g_{33}(x)$, or what we will call
$g_{[ab]}(x)\equiv g_{ab\>(diagonal)}(x)$ each element of which must be
strictly positive as required by physics.555The reader should compare the
three diagonalized positive metric variables with $q>0$ in the half-harmonic
oscillator story, which then required an affine quantization. Now we introduce
two diagonal matrices, $A\equiv g_{[ab]}(x)$ and $B\equiv\eta_{[d]}^{[c]}(x)$.
Next we let $A=e^{B}$, and use $O$ and $O^{T}$, one field element from $SO(3)$
and its transpose, and where $O^{T}(x)\,O(x)=O(x)\,O^{T}(x)=I$. Now, we are
led to

$\displaystyle\bar{A}\equiv O^{T}A\,O=O^{T}e^{B}\,O=e^{O^{T}B\,O}\equiv
e^{\bar{B}}\;,$ (16)

where we now let $\bar{B}(x)=[\eta(x)]$ and
$[e^{[\eta(x)]}]_{ab}=[g(x)]_{ab}\equiv g_{ab}(x)>0$.

Moreover, we now have

$\displaystyle\hskip-8.99994ptH^{\prime\prime}(\pi,g)=\langle\pi;g|{\cal
H}^{\prime\prime}(\hat{\kappa},{\hat{g}})|\pi;g\rangle$ $\displaystyle\hskip
30.00005pt=\langle\alpha|{\cal
H}^{\prime\prime}(\hat{\kappa}^{a}_{b}+\pi^{ac}g_{bc}\;{\hat{g}}_{bc},\,g_{mn}\;{\hat{g}}_{mn})|\alpha\rangle$
$\displaystyle\hskip 30.00005pt={\cal
H}^{\prime\prime}(\pi,g)+{\cal{O}}(\hbar;\pi,g)\;,$ (17)

where here $\pi$ stands for $\\{\pi^{a}_{b}(\cdot)\\}$ and $g$ stands for
$\\{g_{cd}(\cdot)\\}$, and use the facts that
$\langle\alpha|{\hat{g}}_{de}(x)|\alpha\rangle=\delta_{de}$ and
$\langle\alpha|{\hat{\pi}}^{c}_{d}(x)|\alpha\rangle=0$. Since
${\hat{\pi}}^{a}_{b}\propto\hbar$, it follows that if $\hbar\rightarrow 0$,
then $H^{\prime\prime}(\pi,g)={\cal H}^{\prime\prime}(\pi,q)\Rightarrow{\cal
H}^{\prime\prime}({\hat{\pi}},{\hat{g}})$, pointing the way to a valid
quantization of gravity.

Next we present the ADM classical Hamiltonian in our chosen affine variables,
which, now introducing $g(x)\equiv\det[g_{ab}(x)]>0$, leads to

$\displaystyle
H^{\prime\prime}(\pi,g)={\textstyle\int}\\{g(x)^{-1/2}[\pi^{a}_{b}(x)\pi^{b}_{a}(x)-\textstyle{\frac{1}{2}}\pi^{a}_{a}(x)\pi^{b}_{b}(x)]$
$\displaystyle\hskip
100.00015pt+g(x)^{1/2}\,{}^{(3)}\\!\\!R(x)\\}\;d^{3}\\!x\;,$ (18)

where ${}^{(3)}\\!\\!R(x)$ is the Ricci scalar for three spatial coordinates,
and which contains all of the spatial derivatives of the metric field. Already
this version of the classical Hamiltonian contains reasons that restrict
$g(x)$ to $0<g(x)<\infty$, plus $|\pi^{a}_{b}(x)|<\infty$, and
$|^{(3)}\\!\\!R(x)|<\infty$. Just like the field theory example using affine
variables, our gravity fields have no f-issues for the gravity story.

We next introduce the coherent states for gravity. For that we use the
dilation gravity operator
${\hat{\pi}}^{a}_{b}(x)=[{\hat{\pi}}^{ac}(x)^{\dagger}\,{\hat{g}}_{bc}(x)+{\hat{g}}_{bc}(x)\,{\hat{\pi}}^{ac}(x)]/2$
along with ${\hat{g}}_{ab}(x)>0$. Our gravity coherent states are

$\displaystyle|\pi;g\rangle=e^{i{\textstyle\int}\,\pi^{ab}\,{\hat{g}}_{ab}(x)\;d^{s}\\!x/\hbar}\;e^{-i{\textstyle\int}\;\eta^{a}_{b}(x)\,{\hat{\pi}}^{b}_{a}(x)\;d^{s}\\!x/\hbar}\;|\alpha\rangle\;,$
(19)

where we choose $|\alpha\rangle$ by
$[({\hat{g}}_{ab}(x)-\delta_{ab}1\\!\\!1)+i({\hat{\pi}}^{c}_{a}(x)+{\hat{\pi}}^{c}_{b}(x))/2\alpha(x)\hbar]\;|\alpha\rangle=0$.

Using Schrödinger‘s representation, we are led to

$\displaystyle{\cal
H}^{\prime}({\hat{\pi}},g)={\textstyle\int}\\{\;[\,{\hat{\pi}}^{a}_{b}(x)\,g(x)^{-1/2}\,{\hat{\pi}}^{b}_{a}(x)-\textstyle{\frac{1}{2}}{\hat{\pi}}^{a}_{a}(x)\,g(x)^{-1/2}\,{\hat{\pi}}^{b}_{b}(x)\;]$
$\displaystyle\hskip
100.00015pt+\,g(x)^{1/2}\;{}^{(3)}\\!\\!R(x)\;\\}\;d^{3}\\!x\>.$ (20)

Finally, we close with Schrödinger’s equation

$\displaystyle i\hbar\,\partial\>\Psi(g,t)/\partial t={\cal
H}^{\prime}({\hat{\pi}},g)\;\Psi(g,t)\>,$ (21)

which offers the necessary ingredients for the foundation of a valid
quantization of the classical Hamiltonian, and is an important part of the
full story. An unexpected property is that
${\hat{\pi}}^{a}_{b}(x)\;g(x)^{-1/2}=0$, which would simplify matters
considerably!

As before, it may be necessary to introduce some version of regularization for
these equations, but these same equations point the way to proceed. In that
effort, note that although
${\hat{\pi}}^{ac}(x)^{\dagger}\neq{\hat{\pi}}^{ac}(x)$ it can be helpful to
know that
${\hat{\pi}}^{ac}(x)^{\dagger}\,{\hat{g}}_{bc}(x)={\hat{\pi}}^{ac}(x)\,{\hat{g}}_{bc}(x)$;
see footnote 1.

A full quantization of gravity must deal with first and likely second order
constraints, which are designed to reduce the overall Hilbert space to secure
a final quantization. This paper is not the proper place to finalize a
quantization of gravity.

For those who like path integration, a recent paper has found that affine
quantization provides a royal path to a valid path integration of Einstein’s
gravity [13], in which dealing with constraints is included, along with
additional references therein.

##  6 Summary

We have stressed the definition, procedures, and advantages of affine
quantization in offering to secure valid quantizations of several different
examples. We first used simple models with different coordinate spaces to
quantize models that fit their coordinate space in order to insure a valid
result. Following that path, we examined field models by letting affine
variables remove all f-issues as unphysical for any of nature’s fields.
Finally, we were able to put affine procedures to work on an f-issue-free
version of an affine quantization of gravity, a contribution that has been
needed for a long tine. It is hoped that readers can use affine quantization
procedures to help solve some of their problems.

## References

  * [1] P.A.M. Dirac, The Principles of Quantum Mechanics, (Claredon Press, Oxford, 1958), page 114, in a footnote.
  * [2] Wikipedia; “Fubini-Study metric”.
  * [3] Scholarpedia: “negative curvature 2d”.
  * [4] J. Klauder, “The Favored Classical Variables to Promote to Quantum Operators”, Journal of High Energy Physics, Gravitation and Cosmology, 6, 828-832 (2020); doi: 10.4236/jhepgc.2020.64055.
  * [5] L. Gouba, “Affine Quantization on the Half Line”, Journal of High Energy Physics, Gravitation and Cosmology 7, 352-365 (2021); arXiv:2005.08696.
  * [6] J. Klauder, “Using a Toy Model to Improve the Quantization of Gravity and Field Theories”, Journal of High Energy Physics, Gravitation and Cosmology, 8, 303-308 (2022); doi: 10.4236/jhepgc.2022.82024.
  * [7] R. Arnowitt, S. Deser, and C. Misner, in Gravitation: An Introduction to Current Research, Ed. L. Witten, (Wiley & Sons, New York, 1962), p. 227; arXiv:gr-qc/0405109.
  * [8] R. Fantoni and J. Klauder, “Affine Quantization of $(\varphi^{4})_{4}$ Succeeds While Canonical Quantization Fails”, Phys. Rev. D 103 076013 (2021); arXiv:2012.09991.
  * [9] B. Freedman, P. Smolensky, D. Weingarten, “Monte Carlo of the Continuum Limit of $\varphi^{4}_{4}$ and $\varphi^{4}_{3}$”, Phys. Lett. B 113, 481 (1982).
  * [10] M. Aizenman, “Proof of the Triviality. of $\varphi^{4}_{d}$ Field Theory and Some Mean-Field Features of Ising Models for $d>4$”, Phys. Rev. Lett.47, 1-4, E-886 (1981).
  * [11] J. Fröhlich, “On the Triviality of $\lambda\varphi^{4}_{d}$ Theories and the Approach to the Critical Point in $d\geq 4$ Dimensions”, Nuclear Physics B 200, 281-296 (1982)
  * [12] M. Aizenman and H. Duminil-Copin, “Marginal triviality of the scaling limits of critical 4D Ising and $\varphi^{4}$ models”, Annals of Mathematics 194, 163 (2021).
  * [13] J. Klauder, “A Straight Forward Path to a Path Integration of Einstein’s Gravity”; arXiv:2203.15141.


# Distributed Control of Multi-zone HVAC Systems Considering Indoor Air
Quality

Yu Yang, Seshadhri Srinivasan,
Guoqiang Hu, and Costas J. Spanos This work was supported by the Republic of
Singapore’s National Research Foundation through a grant to the Berkeley
Education Alliance for Research in Singapore (BEARS) for the Singapore-
Berkeley Building Efficiency and Sustainability in the Tropics (SinBerBEST)
Program. BEARS has been established by the University of California, Berkeley
as a center for intellectual excellence in research and education in
Singapore.Yu Yang and Seshadhri Srinivasan are with SinBerBEST, Berkeley
Education Alliance for Research in Singapore, Singapore, 138602. E-mail: (
{yu.yang, seshadhri.srinivasan}@bears-berkeley.sg).Guoqiang Hu is with the
School of Electrical and Electronic Engineering, Nanyang Technological
University, Singapore, 639798. E-mail: (gqhu@ntu.edu.sg).Costas J. Spanos is
with the Department of Electrical Engineering and Computer Sciences,
University of California, Berkeley, CA, 94720 USA. Email:
(spanos@berkeley.edu).

###### Abstract

This paper studies a scalable control method for multi-zone heating,
ventilation and air-conditioning (HVAC) systems to optimize the energy cost
for maintaining thermal comfort and indoor air quality (IAQ) (represented by
CO2) simultaneously. This problem is computationally challenging due to the
complex system dynamics, various spatial and temporal couplings as well as
multiple control variables to be coordinated. To address the challenges, we
propose a two-level distributed method (TLDM) with a _upper_ level and _lower_
level control integrated. The _upper_ level computes zone mass flow rates for
maintaining zone thermal comfort with minimal energy cost, and then the
_lower_ level strategically regulates zone mass flow rates and the ventilation
rate to achieve IAQ while preserving the _near_ energy saving performance of
_upper_ level. As both the _upper_ and _lower_ level computation are deployed
in a distributed manner, the proposed method is scalable and computationally
efficient. The near-optimal performance of the method in energy cost saving is
demonstrated through comparison with the centralized method. In addition, the
comparisons with the existing distributed method show that our method can
provide IAQ with only little increase of energy cost while the latter fails.
Moreover, we demonstrate our method outperforms the demand controlled
ventilation strategies (DCVs) for IAQ management with about $8$-$10\%$ energy
cost reduction.

###### Note to Practitioners

The high portion of building energy consumption has motivated the energy
saving for HVAC systems. Concurrently, the living standards for indoor
environment are rising among the occupants. Nevertheless, the status quo on
improving building energy efficiency has mostly focused on maintaining thermal
comfort (such as temperature), and the indoor air quality (IAQ) (usually
represented by CO2 level) has been seldom incorporated.

In our previous work with the similar setting, we observed that the CO2 levels
will surge beyond tolerance during the high occupancy periods if only thermal
comfort is considered for HVAC control. This deduces the IAQ and thermal
comfort should be jointly considered while pursuing the energy cost saving
target and thus studied in this paper. This task is computationally cumbersome
due to the complex system dynamics (thermal and CO2) and tight correlations
among the different control components (variable air volume and fresh air
damper). To cope with these challenges, this work develops a two-level
distributed computation paradigm for HVAC systems based on problem structures.
Specifically, the _upper_ level control first calculates zone mass flow rates
for maintaining comfortable zone temperature with minimal energy cost and then
the _lower_ level strategically regulates the computed zone mass flow rates as
well as ventilation rate to satisfy IAQ while preserving the _near_ energy
saving performance of the _upper_ level control. As both the _upper_ and
_lower_ level calculation can be implemented in a distributed manner, the
proposed method is scalable to large multi-zone deployment. The method’s
performance both in maintaining comfort (i.e., thermal comfort and IAQ) and
energy cost saving is demonstrated via simulations in comparisons with the
centralized method, the distributed token-based scheduling strategy and the
demand controlled ventilation strategies.

###### Index Terms:

multi-zone HVAC system, two-level, distributed approach, IAQ, CO2.

##  I Introduction

The heating, ventilation and air-conditioning (HVAC) systems account for
40%-50% of building energy consumption for maintaining comfortable indoor
environment [1]. Aware of the high portion of energy consumption, the
proposition to improve energy efficiency for HVAC systems has stimulated
widespread attention [2, 3]. A plenty of works have demonstrated substantial
energy can be saved by deploying advanced HVAC control strategies (see [4, 5],
for examples). However, regarding the human sensation, the main focus has been
placed on thermal comfort which is usually indicated by temperature and
humidity [6, 7], the indoor air quality (IAQ), such as the carbon dioxide
(CO2) concentration has been seldom studied.

IAQ is closely related to the mechanical ventilation rate (i.e., fresh air
infusion) of HVAC systems, especially nowadays where most buildings are
constructed in closed envelope for energy saving concerns. IAQ has emerged as
a critical issue along developing building automation. On one hand, the living
standards and awareness of health are rising among the occupants; on the other
hand, IAQ is closely related to occupant well-being and working productivity
[8, 9]. A study conducted by the National Institute of Environmental Health
Science (NIEHS) in 2015 has demonstrated the impact of IAQ on cognitive
abilities [10]. Significantly, the findings show that the strategic skills of
occupants reduce to 20% with an indoor CO2 concentration 1,400 ppm compared to
the normal outdoor level of 400 ppm. Straightforwardly, for energy-efficient
HVAC control, as the thermal comfort only relates to zone mas flow rates, the
fresh air ventilation of HVAC systems is usually suppressed to the lower
boundary to minimize the cooling load. However, for IAQ management, sufficient
ventilation is required to dilute the indoor air. This implies considering
thermal comfort alone is not enough for developing energy-efficient HVAC
control strategies. As in such setting, the energy saving target may be
achieved at the expense of awful IAQ. Specifically, we can imagine that the
CO2 concentration will accumulate and surge beyond tolerance during the time
periods with high occupancy though the temperature is maintained comfortable
if sufficient fresh air ventilation is not provided by the HVAC systems.

Notably, increasing focus has been shifted towards IAQ management for HVAC
systems both from industries and research communities in recent years. Some
industries have informed that the thermostats will integrate CO2 sensors
inside for monitoring IAQ in foreseeable future [11]. Within the research
communities, some pioneering works have jointly considered IAQ and thermal
comfort while studying energy-efficient HVAC control (see [12, 13, 14], for
examples). Nevertheless, this problem hasn’t been well addressed yet as _i)_
the related works are still fairly limited (see the references therein); and
_ii)_ most of them are for single-zone cases [12, 13] or not scalable to
multi-zone commercial buildings due to the centralized computation framework
[14]. The computational challenges of the complex problem remain to be
addressed:

  * i)

Multiple control variables need to be coordinated for achieving the multiple
objectives. Both the ventilation rate and zone mass flow rates require to be
coordinated as to optimize energy efficiency for achieving thermal comfort and
IAQ simultaneously. Particularly, a good IAQ requires high ventilation rate
and zone mass flow rates. However, this may induce lower temperature bound
violations and high energy cost.

  * ii)

The problem is non-linear and non-convex due to the complex system dynamics.
Both the energy cost and system dynamics (i.e., temperature and CO2) are non-
linear w.r.t the control inputs: zone mass flow rates and ventilation rate.

  * iii)

There exist various spatial and temporal couplings. The optimal operation of
HVAC systems should consider zone temperature and CO2 inertia, which
correspond to substantial temporal constraints. Besides, the inter-zone heat
transfer and the recirculated air to AHU induce tight spatial couplings.

###  I-A Contributions

Motivated by the literature, this paper studies the control of multi-zone
commercial HVAC systems to optimize the energy cost for maintaining thermal
comfort and IAQ simultaneously. We make the following contributions to
overcome the computational challenges.

  * (C1)

We propose a two-level control method integrated with the _upper_ level
control (ULC) and the _lower_ level control (LLC) by exploiting the problem
structures.

  * (C2)

While the ULC adopts an existing distributed method, we develop a distributed
method for the LLC to achieve scalability and computation efficiency.

  * (C3)

We demonstrate the method’s performance in energy cost saving, human comfort
(i.e., temperature and CO2) as well as computation efficiency through
simulations.

Our two-level structure is motivated by the independent zone temperature and
zone CO2 dynamics (but both subject to the control inputs), which makes it
possible to tackle the two comfort indexes successively. Therefore, we
decompose the problem into two levels: the ULC first computes the optimal zone
mass flow rates to for maintaining zone thermal comfort with minimal energy
cost. Successively, the LLC strategically regulates the computed zone mass
flow rates from ULC and ventilation rate to achieve the desirable IAQ. Such
two-level paradigm makes it computationally tractable to achieve the two
comfort indexes simultaneously while preserving the _near-optimal_ energy
cost. Moreover, the two-level structure favors computation efficiency as the
LLC is only invoked when the CO2 upper bounds are to be violated.
Particularly, as the ULC on thermal comfort has been comprehensively studied
in our previous work [15] and thus adopted, we place our main focus on
bringing in the two-level structure and the LLC.

The remainder of this paper is structured as below. In Section II, we review
the related works. In Section III, we present the problem formulation. In
Section IV, we discuss the two-level distributed method. In Section V, we
evaluate the performance of the method via simulations. In Section VI, we
briefly conclude this paper.

##  II Related Works

Developing advanced control for HVAC systems to improve energy efficiency has
been an edge issue along building automation over the past decades. Various
methods have been studied and some comprehensive review can refer to [2, 3].
Significantly, these existing works have proved that substantial energy can be
saved by deploying advanced HVAC control while not compromising human comfort.
From the computation structure standpoint, these methods can be categorized
into centralized methods [4] and decentralized methods [5]. Centralized
approaches are usually developed for single-room/zone cases and not amenable
to commercial buildings due to the substantial computation burden (see [6, 7],
for examples). By contrast, a number of decentralized approaches have been
developed to address the computation challenges of commercial HVAC systems [5,
15, 16].

Nevertheless, one critical issue to be noticed is that most of these works
have focused on thermal comfort (i.e., temperature, humidity, etc.) and lack
IAQ management while pursuing energy cost savings (see [17, 18] and the
references therein). In this backdrop, the IAQ may be awful and beyond
tolerance though the temperature is maintained comfortable. For example, the
pollutants and particles, especially CO2 concentration accumulate over the
time if sufficient ventilation (i.e, outdoor fresh air) is not provided. This
is easy to understand as high proportion of recirculated air is preferred to
reduce the cooling load. In the literature, the IAQ management is mostly
addressed by some simple ventilation rules referred to the demand controlled
ventilation strategies (DCVs). Such methods could be CO2-based [19, 20, 21] or
occupancy-based [22, 23, 24]. The main ideas are adjusting the amounts of
fresh air infusion based on the detected instantaneous CO2 concentration or
occupancy. Whereas for multi-zone commercial buildings, the DCVs tend to cause
over-ventilation or under-ventilation due to zone CO2 or occupancy variations
[25]. Another drawback is that the ventilation regulation for IAQ management
are not coordinated with the thermal comfort, which will jeopardize the energy
saving objective of HVAC systems.

It’s imperative to jointly consider both thermal comfort and IAQ
simultaneously both for achieving HVAC energy cost savings and maintaining
human comfort. Such awareness has motivated the joint management of thermal
comfort and IAQ for single-zone cases [12, 26, 13]. Particularly, most of
these methods depend on simplified linear models to capture system dynamics.
Therefore, the methods and modelings are generally not amenable to multi-zone
commercial buildings. As a scarce exception, [14] studied both thermal comfort
and IAQ management for commercial HVAC systems based on Lyapunov optimization.
The performance of the method both in maintaining human comfort (i.e., thermal
comfort and IAQ) and saving energy cost was demonstrated on a $4$-zone case
study. However, the computation tractability for larger scale applications
remains to be addressed and motivates this work.

##  III Problem Formulation

###  III-A HVAC Systems for Commercial Buildings

The configuration for a typical commercial HVAC system is shown in Fig. 1,
which mainly consists of an Air Handling Unit (AHU), Variable Air Volume (VAV)
boxes and a chiller system (not shown here). The AHU is basically integrated
with a damper, a cooling/heating coil and a supply fan. The heating/cooling
coil cools down/heats up the mixed air (the mixture of outside fresh air and
inside recirculated air) to the set-point temperature. The damper regulates
the ventilation rates: the fraction of recirculated air $d_{r}$ or the
fraction of outdoor fresh air $1-d_{r}$. Therefore, a smaller $d_{r}$
specifies higher proportion of fresh air infusion but induces higher cooling
demand for the HVAC system, and vice versa. The supply fan is for driving the
circulation of air within the building duct network. As another main
component, the VAV boxes are connected to the zones, each of which consists of
a damper and a heating coil. The damper regulates the zone mass flow rate and
the heating coil reheats the supply air if necessary (not discussed in this
paper). The chiller system is usually constituted by a chiller pump, water
tank and the chiller, which provides continuous chilled water to the cooling
coils of AHU. Except for the chiller, the chiller pump is also partially
responsible for the HVAC’s energy consumption for circulating the water
between the water tank and the chiller. This paper studies HVAC systems with
constant water flow system [27] where the chiller pump’s energy consumption
can be regarded as fixed and thus not explicitly discussed. Besides, without
loss of generality, this paper studies the cooling mode. More details of the
HVAC systems can refer to [4, 28].

As aforementioned, we seek to minimize the energy cost for maintaining both
zone thermal comfort and IAQ simultaneously. In such setting, both zone air
flow rates and ventilation rate ($d_{r}$) need to be coordinated. The problem
is studied in a discrete-time framework on a daily basis with $\Delta_{k}=30$
$\mathrm{min}$’s sampling and computing epoch. To account for the multi-stage
control under uncertainties (e.g., weather, occupancy, etc.), we deploy the
model predictive control (MPC) framework: the control inputs at each executed
epoch are computed based on the predicted information over a look-ahead
planning horizon $H=10$ (5$\mathrm{h}$). This process is repeated with the
time evolving until the end of the optimization horizon $\mathcal{T}$.

Figure 1: The schematic of the HVAC system for multi-zone buildings.

###  III-B Zone Thermal Dynamics

We consider a commercial building with $I$ zones indexed by
$\mathcal{I}=\\{1,2,\cdots,I\\}$. At each computing epoch, the zone thermal
dynamics over the planning horizon $\mathcal{H}=\\{0,1,\cdots,H-1\\}$ are
captured by a Resistance-Capacitance (RC) network [29, 30]:

$\begin{split}&C^{\rho}_{i}(T_{i}(k+1)\\!-\\!T_{i}(k))\\!=\\!\sum_{\mathclap{j\in\mathcal{N}_{i}}}\\!\frac{T_{j}(k)\\!-\\!T_{i}(k)}{R_{ij}}\Delta_{k}\\!+\\!\frac{T_{o}(k)\\!-\\!T_{i}(k)}{R_{oi}}\\!\Delta_{k}\\\
&\quad+c_{p}m_{i}^{z}(k)(T_{c}-T_{i}(k))\Delta_{k}+Q_{i}(k)\Delta_{k},~{}\forall
i\in\mathcal{I},k\in\mathcal{H}.\end{split}$ (1)

where $k\in\mathcal{H}$ and $i,j\in\mathcal{N}$ denote the time and zone
indices. $C^{\rho}_{i}$ is the zone air heat capacity. $T_{i}(k)$, $T_{o}(k)$
and $T_{c}$ denote the zone temperature, outdoor temperature and the set-point
temperature of AHU, respectively. $R_{ij}$ ($R_{ji}$) denotes the adjacent
zone thermal resistance. Particularly, $R_{oi}$ denotes the thermal resistance
between zone $i$ and the outside, We use $\mathcal{N}_{i}$ to indicate the set
of spatially adjacent zones to zone $i$. $c_{p}$ is the specific heat of the
air. $m^{z}_{i}(k)$ indicates the zone mass flow rate. $Q_{i}(k)$ captures
zone internal heat gains from the occupants and electrical equipment, which
can estimated by the occupancy [31].

We translate (1) into a standard form:

$\begin{split}&T_{i}(k+1)=A_{ii}T_{i}(k)+\sum_{j\in\mathcal{N}_{i}}A_{ij}T_{j}(k)\\\
&\quad\quad+C_{ii}m^{z}_{i}(k)(T_{i}(k)-T_{c})+D_{i}(k),~{}\forall
i\in\mathcal{I},k\in\mathcal{H}.\end{split}$ (2)

where
$A_{ii}=1\\!-\\!(\sum_{j\in\mathcal{N}_{i}}\frac{\Delta_{k}}{R_{ij}C^{\rho}_{i}}+\frac{\Delta_{k}}{C^{\rho}_{i}R_{oi}})$,
$A_{ij}\\!=\\!\frac{\Delta_{k}}{C^{\rho}_{i}R_{ij}}$,
$C_{ii}\\!=\\!-\frac{\Delta_{k}\cdot c_{p}}{C^{\rho}_{i}}$, and
$D_{i}(k)=\frac{\Delta_{k}T_{o}(k)}{C^{\rho}_{i}R_{oi}}+\frac{\Delta_{k}\cdot
Q_{i}(k)}{C^{\rho}_{i}}$.

Similar to [12, 14], this paper uses CO2 concentration as an IAQ indicator.
The zone CO2 dynamics for multi-zone commercial buildings can be described by
[14]

$\begin{split}m_{i}(C_{i}(k+1)&-C_{i}(k))=N_{i}(k)C_{g}\Delta_{k}\\\
&+m_{i}^{z}(k)(C_{z}(k)-C_{i}(k))\Delta_{k}~{}\forall
i\in\mathcal{I},k\in\mathcal{H}\\\ \end{split}$ (3)

where $m_{i}$ denotes zone air mass. $C_{i}(k)$ (in ppm) denotes zone CO2
concentration. We suppose the occupants are the main source of CO2 generation
and thus the CO2 accumulation can be estimated by the average CO2 generation
rate per person $C_{g}$ ($\mathrm{g}\text{\,}{\mathrm{h}}^{-1}$) multiplied by
the occupancy $N_{i}(k)$ as indicated in the first term of right-hand side.
$C_{z}(k)$ denotes the CO2 concentration of supply air, which can be estimated
by

$\begin{split}&C_{z}(k)=\big{(}1-d_{r}(k)\big{)}C_{o}(k)+d_{r}(k)C_{\textrm{m}}(k),\\\
&\textrm{with}~{}C_{\textrm{m}}(k)=\frac{\sum_{i\in\mathcal{I}}m_{i}^{z}(k)C_{i}(k)}{\sum_{i\in\mathcal{I}}m_{i}^{z}(k)},~{}\forall
k\in\mathcal{H}\end{split}$ (4)

where $d_{r}(k)\in[0,1]$ denotes the fraction of return air delivered to AHU
(i.e., the ventilation rate to be controlled). $C_{m}(k)$ characterizes the
average CO2 concentration of the return air from all zones.

From (3) and (4), we note that the zone CO2 dynamics are non-linear and fully
coupled through the recirculated air.

###  III-C Decision Variables

Our objective is to minimize the HVAC energy cost for maintaining thermal
comfort and IAQ. Therefore, our decision variables can be divided into control
variables and state variables. The control variables are the control inputs of
the HVAC system that affect indoor condition, which include the ventilation
rate $d_{r}(k)$ and zone mass flow rate $m_{i}^{z}(k)$. Our state variables
are zone temperature $T_{i}(k)$ and zone CO2 concentration $C_{i}(k)$, which
indicate human comfort and requires to be maintained.

###  III-D Objective Function

The energy consumption of the HVAC system is mainly incurred by the cooling
coil $P_{f}(k)$ and supply fan $P_{f}(k)$ within AHU, i.e.,

$\begin{split}&P_{c}(k)=c_{p}\eta(1-d_{r}(k))\sum_{i\in\mathcal{I}}m_{i}^{z}(k)(T_{o}(k)-T_{c})\\!\\\
&\quad\quad\quad+c_{p}\eta
d_{r}(k)\sum_{i\in\mathcal{I}}m_{i}^{z}(k)(T_{i}(k)-T_{c})\\\
&P_{f}(k)=\kappa_{f}(\sum_{i\in\mathcal{I}}m_{i}^{z}(k))^{2}\end{split}$ (5)

where $\eta$ is the reciprocal of the coefficient of performance (COP) of the
chiller, which captures the ratio of provided cooling to the total consumed
electrical power.

Considering the energy consumption is not easy to inspect in practice, we
selected the total energy cost charged by the electricity price $c_{k}$
(s$$\mathrm{/}{\mathrm{k}\mathrm{W}}$) on a daily basis as our objective:

$J=\sum_{k\in\mathcal{H}}c_{k}\big{(}P_{c}(k)+P_{f}(k)\big{)}\Delta_{k}$ (6)

###  III-E System Constraints

The HVAC operation should respect the thermal comfort and IAQ requirements,
which are represented by zone temperature bounds [5] and CO2 bounds [14] in
this paper:

$\begin{split}{T}^{\min}_{i}&\leq T_{i}(k)\leq{T}^{\max}_{i},\quad\forall
i\in\mathcal{I},~{}k\in\mathcal{H}.\end{split}$ (7) $\begin{split}C_{i}(k)\leq
C^{\max}_{i},\quad\forall i\in\mathcal{I},~{}k\in\mathcal{H}.\end{split}$ (8)

where ${T}^{\min}_{i}$ and ${T}^{\max}_{i}$ characterize the comfortable zone
temperature ranges. $C_{i}^{\max}$ denotes the desirable zone CO2 upper bound.
Particularly, the formulation can account for personalized human comfort by
differentiating ${T}^{\min}_{i}$ , ${T}^{\max}_{i}$ and $C_{i}^{\max}$ w.r.t
the zones.

Additionally, the HVAC operation must abide by the physical limits of VAV
boxes and AHU:

$\begin{split}&{m}^{z,\min}_{i}\leq m^{z}_{i}(k)\leq
m_{i}^{z,\max},\quad\forall i\in\mathcal{I},~{}k\in\mathcal{H}\\\ \end{split}$
(9)
$\begin{split}&\sum_{i\in\mathcal{I}}m_{i}^{z}(k)\leq{m}^{\max},~{}~{}~{}\forall
k\in\mathcal{H}\quad\quad\quad\quad~{}~{}\\\ \end{split}$ (10)

where $m^{z,\min}_{i}$ and $m^{z,\max}_{i}$ denote the operation range of VAV
boxes. $m^{\max}$ denotes the supply capacity of AHU.

Similarly, the damper of AHU for regulating ventilation rate is usually
restricted to the operation range characterized by $d^{\min}_{r}$ and
${d}^{\max}_{r}$ to reduce wear and tear:

$\begin{split}d^{\min}_{r}\leq d_{r}(k)\leq d^{\max}_{r},~{}~{}\forall
k\in\mathcal{H}.\end{split}$ (11)

###  III-F The Problem

Overall, the optimization problem at each computing epoch can be summarized as
($\mathcal{P}$).

$\displaystyle\min_{\bm{m}^{z}_{i},\bm{T}_{i},\bm{C}_{i},i\in\mathcal{I}.\bm{d}_{r}}J$
($\mathcal{P}$ ) $\displaystyle s.t.~{}(\ref{zone temperature dynamics2})\\!$
$\displaystyle-\\!(\ref{CO2 dynamics of supply air}),(\ref{temperature
ranges})\\!-\\!(\ref{CO2 ranges}),~{}\eqref{zone air flow rate
bound}\\!-\\!\eqref{total air flow rate bound},~{}\eqref{ventilation
rate}.\quad\quad$

where we have $\bm{m}^{z}_{i}=[m^{z}_{i}(k)]_{k\in\mathcal{H}}$,
$\bm{T}_{i}=[T_{i}(k)]_{k\in\mathcal{H}}$,
$\bm{C}_{i}=[C_{i}(k)]_{k\in\mathcal{H}}$ ($\forall i\in\mathcal{I}$) and
$\bm{d}_{r}=[d_{r}(k)]_{k\in\mathcal{H}}$ denoting the concatenating decision
variables over the computing epoch.

It’s clear to see that problem ($\mathcal{P}$) is non-linear and non-convex.
The non-linearity and non-convexity both arise from the objective function and
the constraints. Moreover, the substantial spatially and temporally coupled
constraints imposed by zone temperature and zone CO2 makes it computationally
intractable even for moderate commercial HVAC systems with centralized
methods. Indeed, this motivates our two-level distributed method to be
discussed.

##  IV Two-Level Distributed Method

To address the computational challenges of problem ($\mathcal{P}$), we propose
a two-level distributed method (TLDM) based on the problem structures: the
zone temperature and CO2 dynamics are independent (but both subject to the
control inputs). We divide the problem into two levels: the _upper_ level
control (ULC) and the _lower_ level control (LLC). The ULC is mainly
responsible for thermal comfort and the LLC corresponds to IAQ management. The
underlying idea of the two-level structure is to address the two comfort
indexes successively while preserving the _near-optimal_ energy saving
performance. Moreover, such two-level structure favors computation efficiency
as the LLC requires to be invoked only if the CO2 is to be violated. Both the
ULC and LLC use distributed implementation to achieve scalability and
computation efficiency. We present the holistic framework of the TLDM in Fig.
2. Particularly, the ULC first computes the zone flow rates $\bm{m}^{z,U}$ for
maintaining thermal comfort with minimal energy cost and then the LLC
strategically regulates the ventilation rate $\bm{d}^{L}_{r}$ to achieve IAQ
while preserving the _near_ energy cost of ULC. The ULC and LLC alternate and
communicate with each other. The implementation of ULC and LLC are as below.

Figure 2: The framework of two-level distributed method.

###  IV-A The Upper Level Control (ULC)

The ULC focuses on maintaining zone thermal comfort with minimal energy cost.
Considering that _i)_ zone temperature is only affected by zone mass flow rate
$m^{z}_{i}(k)$ and _ii)_ the HVAC’s energy cost is non-decreasing w.r.t. the
ventilation rate $d_{r}(k)$, we define the ULC problem as

$\displaystyle\quad\quad\min_{\bm{m}^{z}_{i},\bm{T}_{i},i\in\mathcal{I}.}J^{{U}}$
$\displaystyle s.t.~{}~{}(\ref{zone temperature dynamics2}),(\ref{temperature
ranges}),(\ref{zone air flow rate bound})\\!\\!-\\!\\!(\ref{total air flow
rate bound}).\quad\quad$ ($\mathcal{P}_{U}$) $\displaystyle\quad\quad
d_{r}(k)=d^{L}_{r}(k),k\in\mathcal{H}.$

where we have $J^{{U}}=J$ and
$\bm{d}^{L}_{r}=[d^{L}_{r}(k)]_{k\in\mathcal{K}}$. The intuitive
interpretation of ULC is to identify the zone mass flow rates for maintaining
thermal comfort with minimal energy cost under the specified ventilation rate
$d^{L}_{r}(k)$ from LLC (we can start with ${d}^{\max}_{r}$
($k\in\mathcal{H}$)). For problem ($\mathcal{P}_{U}$), there already exists a
number of distributed methods [15, 5]. We adopt our previously proposed
distributed method [15] and place our main attention on the implementation of
LLC.

###  IV-B The Lower Level Control (LLC)

As indicated in Fig. 2, the LLC only needs to be invoked when the zone CO2
concentrations violate the upper bounds (i.e., $C_{i}(k)\geq C^{\max}_{i}$).
Clearly, both increasing zone mass flow rates and ventilation rate (decrease
$d_{r}(k)$) can dilute zone CO2 thus improve IAQ. However, the latter is
generally more expensive as it increase the cooling load of all zones. With
the concern to preserve the _near_ energy cost computed in ULC, our LLC adopts
a _two-phase_ method to achieve IAQ. We first seek to regulate zone mass flow
rates to satisfy the user-defined zone CO2 bounds. To preserve the _near_
energy cost of ULC, we select the deviation of zone mass flow rates relative
to the ULC calculations as the objective. We have the following problem for
the first phase of LLC:

$\displaystyle\min_{\bm{m}^{z}_{i},\bm{C}_{i},i\in\mathcal{I}.}J^{L}=\sum_{k\in\mathcal{H}}\sum_{i\in\mathcal{I}}\big{(}m^{z}_{i}(k)-m^{z,U}_{i}(k)\big{)}^{2}$
($\mathcal{P}_{L}$) $\displaystyle~{}~{}~{}s.t.~{}~{}m^{z,U}_{i}(k)\leq
m^{z}_{i}(k)\leq{m}^{z,\max}_{i},\forall i\in\mathcal{I}.$ (12)
$\displaystyle\quad\quad\quad\quad\quad\quad~{}(\ref{CO2
dynamics})\\!-\\!(\ref{CO2 dynamics of supply air}),(\ref{CO2
ranges}),(\ref{total air flow rate bound}).\quad\quad$

where constraints (3)-(4) and (8) capture zone CO2 dynamics and bounds.
Constraints (10) denotes the supply capacity of AHU. Particularly, constraints
(12) model the zone mass flow rate lower bounds for maintaining thermal
comfort, which is obtained from the ULC.

We suppose to achieve thermal comfort and IAQ simultaneously by solving
problem ($\mathcal{P}_{L}$), however things may not always happen in that way.
We note that only the upper comfortable zone temperature bounds can be
maintained but not necessarily the lower ones by (12). Therefore, it may occur
that the desirable CO2 is achieved by solving problem ($\mathcal{P}_{L}$) but
the zone temperature drops beyond the lower bounds due to the substantial
increase of zone mass flow rates. Actually, this means regulating zone mass
flow rates alone is not viable to achieve the two comfort indexes
simultaneously, and the second phase of LLC should be invoked to increase the
ventilation rate (decrease $\bm{d}_{r}$). As illustrated in Fig. 2, such two
phases alternate until we achieve the two comfort indexes. Suppose the _two-
phase_ structure gets well across, the remaining problem is about the
implementation.

_First phase of LLC_ : Clearly, for the first phase, we require to solve the
non-linear and non-convex problem ($\mathcal{P}_{L}$), which is a nontrivial
task. We note that the zone CO2 dynamics are fully coupled though the
recirculated air, posing the primary challenge to develop a distributed
method. To handle this, we introduce a learning framework to estimate CO2
concentration for the supply air $C_{z}(k)$. Therefore, the procedure to solve
problem ($\mathcal{P}_{L}$) contains two successive steps: first estimate
$C_{z}(k)$ and then solve problem of ($\mathcal{P}_{L}$) with the estimated
$C_{z}(k)$. Specifically, with the estimated $\bm{C}_{z}(k)$, we have the
decoupled zone CO2 dynamics:

$\begin{split}C_{i}(k+1)=C_{i}(k)+&E_{i}(k)m^{z}_{i}(k)\\\
&+F_{i}(k)m^{z}_{i}(k)C_{i}(k)+G_{i}(k)\end{split}$ (13)

where $E_{i}(k)={C_{z}(k)\Delta_{k}}/{m_{i}}$,
$F_{i}(k)={-\Delta_{k}}/{m_{i}}$ and
$G_{i}(k)={N_{i}(k)C_{g}\Delta_{k}}/{m_{i}}$ are now constant parameters.

We note that the zone CO2 dynamics (13) are bilinear w.r.t the zone mass flow
rate $m^{z}_{i}(k)$ and zone CO2 $C_{i}(k)$. To address the non-linearity, we
introduce some auxiliary decision variables $Z_{i}(k)=m^{z}_{i}(k)C_{i}(k)$
and use the _McCormick envelopes_ [32] to relax the bilinear terms:

$\displaystyle Z_{i}(0)=m^{z}_{i}(0)C_{i}(0).$ (13a) $\displaystyle
Z_{i}(k)\geq
m^{z,\min}_{i}C_{i}(k)+m^{z}_{i}(k)C^{\min}_{i}-m^{z,\min}_{i}C_{i}^{\min},$
(13b) $\displaystyle
Z_{i}(k)\\!\geq\\!m^{z,\max}_{i}C_{i}(k)\\!+\\!m^{z}_{i}(k){C}^{\max}_{i}\\!-\\!{m}^{z,\max}_{i}{C}^{\max}_{i},$
(13c) $\displaystyle Z_{i}(k)\leq
m^{z}_{i}(k){C}^{\max}_{i}+{m}^{z,\min}_{i}C_{i}(k)-{m}^{z,\min}_{i}{C}^{\max}_{i},$
(13d) $\displaystyle
Z_{i}(k)\leq{m}^{z,\max}_{i}C_{i}(k)+m^{z}_{i}(k){C}^{\min}_{i}-{m}^{z,\max}_{i}{C}^{\min}_{i},$
(13e) $\displaystyle\quad\quad\quad\quad\quad\forall
k\in\mathcal{H}\setminus\\{0\\}.$ (14)

where $C_{i}(0)$ denotes the measured zone CO2 at the beginning of current
computing epoch.

By invoking the auxiliary decision variables $Z_{i}(k)$, we have the following
relaxed convex problem for problem ($\mathcal{P}_{L}$):

$\displaystyle\min_{\bm{m}^{z}_{i},\bm{C}_{i},\bm{Z}_{i},\forall
i\in\mathcal{I}.}J^{L}=\sum_{k\in\mathcal{H}}\sum_{i\in\mathcal{I}}(m^{z}_{i}(k)-m^{z,U}_{i}(k))^{2}$
($\mathcal{P}^{{}^{\prime}}_{L}$) $\displaystyle
s.t.~{}~{}C_{i}(k+1)=C_{i}(k)+E_{i}(k)m^{z}_{i}(k)+F_{i}(k)Z_{i}(k)$
$\displaystyle\quad\quad\quad\quad\quad\quad+G_{i}(k),~{}~{}\forall
i\in\mathcal{I},k\in\mathcal{H}.$ (14a)
$\displaystyle\quad\quad\quad\quad\quad(\ref{CO2 ranges}),(\ref{total air flow
rate bound}),\eqref{eq:11},(\ref{13a})-(\ref{13e}).\quad\quad$

We note problem ($\mathcal{P}^{{}^{\prime}}_{L}$) is characterized by _i)_ a
decomposable objective function w.r.t. the zones and _ii)_ coupled linear
constraints. This problem can be efficiently tackled by the Accelerated
Distributed Augmented Lagrangian (ADAL) method [33], which may be more clearly
seen by recasting it into a standard form:

$\displaystyle\min_{\bm{x}_{i},i\in\mathcal{I}}$ $\displaystyle
J^{L}=\sum_{i\in\mathcal{I}}J^{L}_{i}(\bm{x}_{i})$ $\displaystyle s.t.\quad$
$\displaystyle\sum_{i=0}^{I}\bm{A}^{c}_{i}\bm{x}_{i}=\bm{b}^{c}.$
($\mathcal{P}^{{}^{\prime\prime}}_{L}$)
$\displaystyle\bm{x}_{i}\in\mathcal{\bm{X}}_{i},\forall
i\in\mathcal{I}\cup\\{0\\}.$

where $\bm{x}_{i}\\!=\\![(\bm{x}_{i}(k))^{T}]_{k\in\mathcal{H}}^{T}$ ,
$\bm{x}_{i}(k)\\!\\!=\\!\\!\big{(}C_{i}(k),m^{z}_{i}(k),Z_{i}(k)\big{)}^{T}$
concatenating the decision variables of zone $i$.
$J^{L}_{i}(\bm{x}_{i})=\sum_{k\in\mathcal{H}}\big{(}m^{z}_{i}(k)-m^{z,U}_{i}(k)\big{)}^{2}$
represents the local objective function of zone $i$.
$\sum_{i=0}^{I}\bm{A}^{c}_{i}\bm{x}_{i}=\bm{b}^{c}$ accounts for the coupled
linear constraints (10) with an additional slack variables $\bm{x}_{0}(k)\geq
0$ introduced at each stage (transform the inequality constraints to equality
constraints). $\mathcal{\bm{X}}_{i}$ ($\forall i\in\mathcal{I}$) indicates the
local constraints (8), (12), (13a)-(13e) corresponding to zone $i$ and
particularly $\mathcal{X}_{0}=\\{\bm{x}_{0}|\bm{x}_{0}\geq\bm{0}\\}$. We have
$\bm{A}^{c}_{i}=\small\left(\begin{array}[]{ccccccc}0&1&0&0&0&0&\cdots\\\
0&0&0&0&1&0&\cdots\\\ \cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\\
\end{array}\small\right)\in\mathbb{R}^{H\times 3H}$ ($\forall
i\in\mathcal{I}$), and
$\bm{b}^{c}=\big{(}{m}^{\max},{m}^{\max},\cdots,{m}^{\max}\big{)}^{T}\in\mathbb{R}^{H}$.

Following the standard procedure of ADAL [33], we have the augmented
Lagrangian function:

$\displaystyle\begin{split}\mathbb{L}_{\rho}\big{(}\bm{x}_{0},\\!\bm{x}_{1},\\!\cdots,\bm{x}_{I},\bm{\alpha}\big{)}\\!=&\\!\\!\sum_{i\in\mathcal{I}}J^{L}_{i}\\!+\\!\bm{\alpha}^{T}\big{(}\sum_{i=0}^{I}\bm{A}^{c}_{i}\bm{x}_{i}\\!-\\!\bm{b}^{c}\big{)}\\\
&\\!+\\!\frac{\rho}{2}\|\sum_{i=0}^{I}\bm{A}^{c}_{i}\bm{x}^{i}\\!-\bm{b}^{c}\|^{2}\\\
\end{split}$ (16)

where
$\bm{\alpha}=\big{(}\alpha_{0},{\alpha}_{1},\cdots,{\alpha}_{H-1}\big{)}^{T}$
are Lagrangian multipliers. $\rho(\rho>0)$ is penalty parameter.

Therefore, we have the following primal problem with given Lagrangian
multipliers $\bm{\alpha}$:

$\begin{split}\min_{\bm{x}_{0}.\bm{x}^{i},\forall
i\in\mathcal{I}}&\mathbb{L}_{\rho}\big{(}\bm{x}_{0},\bm{x}_{1},\cdots,\bm{x}_{I},\bm{\alpha}\big{)}\\\
s.t.&\quad\bm{x}_{i}\in\mathcal{X}_{i},~{}~{}\forall i\in\mathcal{I}.\\\
&\quad\bm{x}_{0}\geq\bm{0}.\end{split}$ (17)

The main procedures of using ADAL to solve problem
($\mathcal{P}^{{}^{\prime\prime}}_{L}$) generally contain three steps: _i)_
solving the primal problem (17), _ii)_ updating the Lagrangian multipliers
$\bm{\alpha}$, and _iii)_ updating the penalty factor $\rho$. While the last
two procedures are standard, we illustrate how to solve the primal problem
($\mathcal{P}^{{}^{\prime\prime}}_{L}$) in a distributed fashion.
Specifically, we define $I+1$ agents, where Agent $1\sim I$ correspond to the
$I$ zones and Agent $0$ is a virtual agent for managing the slack decision
variable $\bm{x}_{0}$. At each iteration $q$, we have the local objective
functions the agents:

$\begin{split}\mathbb{L}_{\rho}^{0}(\bm{x}_{0},&\bm{x}^{q}_{-0},\bm{\alpha})=\bm{\alpha}^{T}\bm{A}^{c}_{0}\bm{x}_{0}+\frac{\rho}{2}\|\bm{A}^{c}_{0}\bm{x}_{0}+\sum_{i=1}^{I}\bm{A}^{c}_{i}\bm{x}_{i}^{q}-\bm{b}^{c}\|^{2}\\\
\mathbb{L}_{\rho}^{i}(\bm{x}_{i},&\bm{x}_{-i}^{q},\bm{\alpha})=J^{L}_{i}+\bm{\alpha}^{T}\bm{A}^{c}_{i}\bm{x}_{i}+\frac{\rho}{2}\|\bm{A}^{c}_{i}\bm{x}_{i}+\sum_{\mathclap{j\in\mathcal{I}\cup\\{0\\},j\neq
i}}\bm{A}^{c}_{j}\bm{x}_{j}^{q}-\bm{b}^{c}\|^{2},\\\
&\quad\quad\quad\quad\forall i\in\mathcal{I}.\\\ \end{split}$

It’ clear that we can obtain the estimated zone mass flow rates $m^{z}_{i}(k)$
and zone CO2 $C_{i}(k)$ by solving problem
($\mathcal{P}^{{}^{\prime\prime}}_{L}$). Afterward, we invoke the procedure to
update the estimation of $C_{z}(k)$ according to (4).

Overall, we present the details to solve problem ($\mathcal{P}_{L}$) based on
ADAL in Algorithm 1. We use the superscripts $p$ and $q$ to denote the
iteration corresponding to solving the relaxed problem
($\mathcal{P}^{{}^{\prime\prime}}_{L}$) and updating the estimation of
$C_{z}(k)$, respectively. For notation,
$\bm{x}^{q}=[(\bm{x}_{j}^{q})^{T}]^{T}_{j\in\mathcal{I}\cup\\{0\\}}$ represent
the augmented control and state trajectories for all agents while
$\bm{x}_{-i}^{q}=[(\bm{x}_{j}^{q})^{T}]^{T}_{j\in\mathcal{I}\\{0\\},j\neq i}$
excluding Agent $i$. For the ADAL method, we uses the residual error of
coupled constraints:
$r_{q}(\bm{x})=\|\sum_{i=0}^{I}\bm{A}^{c}_{i}\bm{x}_{i}^{q}\\!-\\!\bm{b}^{c}\|\leq\epsilon^{\text{in}}$
(C1) as the stopping criterion. For the estimation of CO2 concentration for
the supply air, we evaluate the deviations of any two successive estimation:
${\big{\|}\bm{C}^{p+1}_{z}-\bm{C}^{p}_{z}\big{\|}}\leq\epsilon^{\text{out}}$
(C2), where $\epsilon^{\text{in}}$ and $\epsilon^{\text{out}}$ are small
positive thresholds.

Algorithm 1 Solve problem ($\mathcal{P}^{{}^{\prime\prime}}_{L}$) based on
ADAL

1:Initialize $p\leftarrow 0$, $\bm{C}^{0}_{z}$ and $\bm{d}_{r}$.

2: Initialize $q\\!\leftarrow\\!0$, $\\!\bm{\alpha}^{0}$, and $\bm{x}^{0}_{i}$
($\forall i\\!\in\\!\mathcal{I}\\!\cup\\!\\{0\\}$),
$\bm{C}_{z}=\bm{C}^{p}_{z}$.

3: for $i\in\mathcal{I}\cup\\{0\\}$, do
$\begin{split}&\bm{x}_{i}^{q+1}\\!\\!=\\!\arg\min_{\bm{x}_{i}}\mathbb{L}^{i}_{\rho}(\bm{x}_{i},\\!\bm{x}_{-i}^{q},\bm{\alpha}^{q})~{}~{}s.t.~{}~{}\bm{x}_{i}\in\mathcal{X}_{i}.\\\
\end{split}$ (18)

4: end for

5: Update the Lagrangian multipliers:
$\bm{\alpha}^{q+1}\\!=\\!\bm{\alpha}^{q}\\!+\\!\rho\big{(}\sum_{i=0}^{I}\bm{A}^{c}_{i}\bm{x}_{i}^{q+1}\\!-\\!\bm{b}^{c}\Big{)}$

6: If (C1) is satisfied, stop with $\bm{x}^{p}=\bm{x}^{q+1}$, otherwise
$q\rightarrow q+1$ and go to Step 3.

7:Estimate $\bm{C}^{p+1}_{z}$ according to
$\begin{split}\bm{C}^{p+1}_{z}=(1-\bm{d}_{r})\bm{C}_{o}+\bm{d}_{r}\frac{\sum_{i\in\mathcal{I}}\bm{m}^{z,p}_{i}\bm{C}_{i}^{p}}{\sum_{i\in\mathcal{I}}\bm{m}^{z,p}_{i}}.\\\
\end{split}$ (19)

8:If (C2) is satisfied, stop with $\bm{x}^{p}$, otherwise set $p\rightarrow
p+1$ and go to Step 2.

9:$\bm{x}^{*}\\!\\!=\\!\\![(\bm{x}^{p}_{i}(k))^{T}]_{k\in\mathcal{H}}^{T}$ ,
$\bm{x}^{p}_{i}(k)\\!\\!=\\!\\!\big{(}C^{p}_{i}(k),m^{z,p}_{i}(k),Z^{p}_{i}(k)\big{)}^{T}$.

_Recursive feasibility_ : Recall that the bilinear terms
$Z_{i}(k)=m^{z}_{i}(k)C_{i}(k)$ are relaxed in problem
($\mathcal{P}^{{}^{\prime}}_{L}$) (or ($\mathcal{P}^{{}^{\prime\prime}}_{L}$)
), therefore the recursive feasibility remains to be addressed. We propose a
heuristic method (see Algorithm 2) to recover a control input
$\hat{\bm{x}}=[(\hat{\bm{x}}_{i}(k))^{T}]_{k\in\mathcal{H}}^{T}$ from the
optimal solution of problem ($\mathcal{P}^{{}^{\prime\prime}}_{L}$), which can
ensure the recursive feasibility and _near_ energy saving performance. The
main idea is to preserve the computed zone mass flow rates:
$\hat{m}^{z}_{i}(k)=m^{z,*}_{i}(k)$ ($\forall
i\in\mathcal{I},k\in\mathcal{H}$) as they determine the HVAC’s energy cost. In
such setting, the desirable CO2 at each executed epoch $t$ will be sustained.
As the method is deployed in MPC framework, the desirable CO2 over the day
(optimization horizon $\mathcal{T}$) will be achieved along the computing
epoch. We illustrate the _feasibility issue_ by induction. For notation, we
indicate the real zone zone CO2 trajectories as $[C_{i,t}]_{t\in\mathcal{T}}$
($\forall i\in{\cal I}$), which suppose to be obtained at the end of
optimization horizon $\mathcal{T}$. At each execution instant $t$, we can
obtain the current zone CO2 measurements $C_{i}(0)=C_{i,t}$ ($\forall
i\in\mathcal{I}$). Suppose at time $t$ we have $C_{i,t}\leq C^{\max}_{i}$
($\forall i\in\mathcal{I}$), it suffices to address the _feasibility issue_ by
proving $C_{i,t+1}\leq C^{\max}_{i}$ ($\forall i\in\mathcal{I}$) as below.

$\begin{split}C_{i,t+1}&=C_{i}(1)\\\
&=C_{i}(0)+E_{i}(0)\hat{m}^{z}_{i}(0)+F_{i}(0)\hat{m}^{z}_{i}(0)C_{i}(0)+G_{i}(0)\\\
&=C_{i}(0)+E_{i}(0)m^{z,*}_{i}(0)+F_{i}(0)m^{z,*}_{i}(0)C_{i}(0)+G_{i}(0)\\\
&=C_{i}(0)+E_{i}(0)m^{z,*}_{i}(0)+F_{i}(0)Z^{*}_{i}(0)+G_{i}(0)\\\ &\leq
C^{\max}_{i},\forall i\in\mathcal{I}.\end{split}$

where the first two equalities are directly from zone CO2 dynamics (13). The
third and fourth equality are derived from $\hat{m}^{z}_{i}(k)=m^{z,*}_{i}(k)$
and (13a). The inequality is deduced from constraint (8).

Algorithm 2 Recover Recursive Feasibility

1:$\bm{x}^{*}\\!\\!=\\!\\![(\bm{x}^{p}_{i}(k))^{T}]_{k\in\mathcal{H}}^{T}$ ,
$\bm{x}^{p}_{i}(k)\\!\\!=\\!\\!\big{(}C^{p}_{i}(k),m^{z,p}_{i}(k),Z^{p}_{i}(k)\big{)}^{T}$
(from Algorithm 1).

2:Set $\hat{C}_{i}(0)=C_{i}(0)$ and $\hat{T}_{i}(0)=T_{i}(0)$ ($\forall
i\in\mathcal{I}$).

3:for $k\in\mathcal{H}$ do

4: for $i\in\mathcal{I}$ do

5: Set $\hat{m}^{z}_{i}(k)=m^{z,*}_{i}(k)$ and
$\hat{Z}_{i}(k)=\hat{m}^{z}_{i}(k)\hat{C}_{i}(k)$.

6: Estimate $\hat{C}_{i}(k+1)$ and $\hat{T}_{i}(k+1)$ by
$\begin{split}&\quad\quad\hat{C}_{i}(k+1)=\hat{C}_{i}(k)+\hat{E}_{i}(k)\hat{m}^{z}_{i}(k)+F_{i}(k)\hat{Z}_{i}(k)+G_{i}(k),\\\
&\quad\quad\hat{T}_{i}(k+1)=A_{ii}\hat{T}_{i}(k)\\!+\\!\\!\sum_{j\in\mathcal{N}_{i}}A_{ij}\hat{T}_{j}(k)\\\
&\quad\quad\quad\quad\quad\quad\quad\quad+C_{ii}\hat{m}^{z}_{i}(k)(\hat{T}_{i}(k)-T_{c})+\\!D_{i}(k),\\\
&\quad\quad\textrm{with}~{}\hat{E}^{z}_{i}(k)={\hat{C}_{z}(k)\Delta_{t}}/{m_{i}},\textrm{and}~{}\hat{C}^{z}(k)=(1-d^{*}_{r}(k))C_{o}(k)\\\
&\quad\quad+d^{*}_{r}(k)\frac{\sum_{i\in\mathcal{I}}\hat{m}^{z}_{i}(k)\hat{C}_{i}(k)}{\sum_{i\in\mathcal{I}}\hat{m}^{z}_{i}(k)},\forall
i\in\mathcal{I},k\in\mathcal{H}.\\\ \end{split}$

7: end for

8:end for

9:$\bm{\hat{m}}^{z}_{i}\\!=\\![\hat{m}^{z}_{i}(k)]_{k\in\mathcal{H}}$,
$\bm{\hat{C}}_{i}=[\hat{C}_{i}(k)]_{k\in\mathcal{H}}$, and
$\bm{\hat{T}}_{i}=[\hat{T}_{i}(k)]_{k\in\mathcal{H}}$.

_Second phase of LLC_ : As discussed previously, it may occur that the lower
comfortable zone temperature bounds may be violated for solving problem
($\mathcal{P}_{L}$). In such situation, we are required to invoke the second
_phase_ of LLC to increase ventilation rate (decrease $d_{r}$). In principle,
the regulating step-size should be determined by the temperature violation
amplitude, however it’s difficult to coordinate such two different dimensions
and we thus use a small constant step-size $\Delta d_{r}$. The overall
framework of the proposed TLDM which embeds the ULC and LLC is presented in
Algorithm 3 . We use $l$ to denote the iteration. Function $\mathbb{I}(A)$ is
an indicator function, where we have $\mathbb{I}(A)=1$ with true condition
$A$, otherwise $\mathbb{I}(A)=0$.

Algorithm 3 Two-level Distributed Method (TLDM)

1:Initialize $l\leftarrow 0$ and $\bm{d}^{0}_{r}$.

2:$[\bm{m}^{z,L}_{i}]_{i\in\mathcal{I}}=ULC(\bm{d}^{l}_{r})$ [15].

3:$[\bm{\hat{m}}^{z}_{i},\bm{\hat{C}}_{i},\bm{\hat{T}}_{i}]_{i\in\mathcal{I}}=LLC([\bm{m}^{z,L}_{i}]_{i\in\mathcal{I}},\bm{d}^{l}_{r})$(Algorithm
1-2).

4:If $\bm{\hat{T}}^{i}\geq{T}^{\min}_{i}$, then stop, otherwise continue.

5:Update ventilation rate $\bm{d}_{r}$:
$\bm{d}_{r}^{l+1}(k)=\bm{d}^{l}_{r}(k)-\Delta
d_{r}\mathbb{I}(\bm{\hat{T}}_{i}(k+1)<{T}^{\min}_{i}),~{}\forall
k\in\mathcal{H}.$

6:Set $l\rightarrow l+1$ and go to Step 2\.

##  V Application

We evaluate the performance of the TLDM on multi-zone commercial HVAC system
via simulations. We first study the energy saving performance and
computational advantage of the method via a _benchmark_ ($5$ zones). After
that the capability and scalability of the method to _medium_ scale (10,20
zones) and _large_ scale (50,100 zones) are illustrated.

###  V-A Benchmark

This part considers a $5$-zone _benchmark_. We select the general comfortable
zone temperature range $[24,26]^{\circ}$C and zone CO2 range $[0,800]$ ppm.
The set-point temperature of AHU is $T_{c}=15^{\circ}$C. We assume the zones
are spatially next to each other $1\leftrightarrow 2\leftrightarrow
3\leftrightarrow 4\leftrightarrow 5\leftrightarrow 1$. The initial zone
temperature is set as $[29,30,31,30,29]^{\circ}$C (zone 1-5). The predicted
outdoor temperature and zone occupancy are shown in Fig. 3. The HVAC’s energy
cost is calculated according to the time-of-use (TOU) price in Singapore [7].
The other parameters refer to TABLE I.

TABLE I: Simulation Parameters Param. |  Value |  Units
---|---|---
$C_{i}(i\in\mathcal{I})$ |  $1.5\times 10^{3}$ |  $\mathrm{kJ}\text{\,}{\mathrm{K}}^{-1}$
$c_{p}$ |  $1.012$ |  $\mathrm{k}\mathrm{J}\mathrm{/}\mathrm{k}\mathrm{g}\mathrm{-}\mathrm{K}$
$R_{oi}$ |  $50$ |  $\mathrm{kW}\text{\,}{\mathrm{K}}^{-1}$
$R_{ij}(i,j\in\mathcal{I})$ |  $14$ |  $\mathrm{kW}\text{\,}{\mathrm{K}}^{-1}$
$\kappa_{f}$ |  $0.08$ |  -
$\eta$ |  $1$ |  -
$C_{g}$ |  $40$ |  $\mathrm{g}\text{\,}{\mathrm{h}}^{-1}$
$\Delta d_{r}$ |  $0.05$ |  -
${m}^{z,\min}_{i}$ |  $0$ |  $\mathrm{kg}\text{\,}{\mathrm{h}}^{-1}$
${m}^{z,\max}_{i}$ |  $0.5$ |  $\mathrm{kg}\text{\,}{\mathrm{h}}^{-1}$

Figure 3: (a) Outdoor temperature. (b) Zone occupancy.

In the _benchmark_ , we compare the proposed method with _i)_ distributed
Token-based scheduling strategy (DTBSS) [5], ii) centralized method, iii) the
commonly-used DCV strategies [34, 25] and iv) sequential quadratic programming
(SQP) [35]. Akin to most existing works, the DTBSS lacks IAQ management [5]
and we therefore fix the ventilation rate in DTBSS as $d_{r}(k)=d_{r}^{\max}$
($\bm{k}\in\mathcal{H}$) to achieve the energy cost saving target. For the
centralized method, we obtain the optimal solution by solving the non-linear
optimization problem ($\mathcal{P}$) using the IPOPT solver embedded in MATLAB
[36]. The DCVs calculate the amount of fresh air by the zone occupancy and
zone area [34, 25], i.e.,

$\begin{split}m^{z,\textrm{fresh}}_{i}(k)=N_{i}(k)R_{p}+A_{i}R_{a},~{}\forall
k\in\mathcal{H}.\end{split}$ (20)

where $A_{i}$ denotes the area of zone $i$. $R_{p}$ and $R_{a}$ denote the
average occupancy and space ventilation rate. This paper considers two DCVs:
_i)_ DCV I: calculating zone fresh air flow rates based on zone occupancy
($R_{p}\geq 0$, $R_{a}=0$), and _ii)_ DCV II: computing zone fresh air flow
rates by zone occupancy and space ($R_{p}\geq 0$, $R_{a}\geq 0$).

For single-zone case, the ventilation rate ($d_{r}$) of DCVs can be
straightforwardly determined by the fresh air infusion and mass flow rate.
Nevertheless, for multi-zone case, the calculation of ventilation rate needs
to account for zone diversity [25]:

$\displaystyle\begin{split}&m^{z,\textrm{fresh}}(k)=\sum_{i\in\mathcal{I}}m_{i}^{z,\textrm{fresh}}(k),~{}m^{z}(k)=\sum_{i\in\mathcal{I}}m^{z}_{i}(k),\\\
&Z(k)=\max_{i\in\mathcal{I}}\Big{\\{}\frac{m^{z,\textrm{fresh}}_{i}(k)}{m^{z}_{i}(k)}\Big{\\}},\quad
X(k)=\frac{m^{z,\textrm{fresh}}(k)}{m^{z}(k)},\\\
&Y(k)=\frac{X(k)}{1+X(k)-Z(k)},\quad d_{r}(k)=1-Y(k),~{}\forall
k\in\mathcal{H}.\\\ \end{split}$ (21)

One critical issue to be noticed is that the DCVs only account for IAQ by
determining the amounts of fresh air infusion, the thermal comfort which
corresponds to the zone mass flow rates is not considered. Therefore, we use
the ULC of TLDM to calculate zone mass flow rates $m^{z}_{i}(k)$ in DCVs for
fair comparisons. Therefore, similar to our TLDM, the DCVs correspond to two
alternative procedures: obtain the zone mass flow rates by solving problem
($\mathcal{P}_{U}$) (start with $d_{r}(k)=d^{\max}_{r}$) and amending the
ventilation rate $d_{r}(k)$ according to (21).

TABLE II: Performance comparisons

Method |  $R_{p}$ |  $R_{a}$ |  Cost |  Time |  TC |  IAQ
---|---|---|---|---|---|---
|  (L/p) |  ($\mathrm{L}\mathrm{/}\mathrm{m}^{2}$) |  (s)
DTBSS |  - |  - |  $245.03$ |  $2.11$ |  Y |  N
Centralized |  - |  - |  $247.15$ |  $575.07$ |  Y |  Y
SQP |  - |  - |  $275.91$ |  $151.26$ |  Y |  Y
DCV I |  $21$ |  $0$ |  $276.23$ |  - |  Y |  Y
DCV II |  $16$ |  $0.04$ |  $274.99$ |  - |  Y |  Y
TLDM |  - |  - |  $257.02$ |  $5.21$ |  Y |  Y

N=No, Y=Yes. |
---|---

Figure 4: Zone mass flow rate (_Benchmark_).
Figure 5: Zone temperature (_Benchmark_).
Figure 6: Zone CO2 (_Benchmark_).
Figure 7: Ventilation rate $\bm{d}_{r}$ (_benchmark_).

We present the energy cost, average computation time of each computing epoch
111MATLAB R2016a on PC with Intel(R) Core(TM) i7-5500U CPU @2.40GHZ
processor., satisfaction of thermal comfort (TC) and IAQ induced by the
different methods in TABLE II. In particular, we exclude the DCVs while
studying the computation time as they depend on _off-line_ regulation for the
occupancy and space ventilation rates, which is time-consuming. We display the
zone mass flow rates, zone temperature, zone CO2 and ventilation rate in Fig.
4-7. First of all, we find the zone thermal comfort is ensured by each method
(Fig. 5) but not necessary for IAQ (Fig. 6) as we observe the desirable zone
CO2 are out of range $[0,800]$ ppm for DTBSS. This is not surprising as the
DTBSS lacks IAQ management. Therefore, the method favors low ventilation rate
(see Fig. 7) to save energy cost and fails to ensure IAQ. Further, we study
the energy cost and computational efficiency of different methods. From TABLE
II, we find the TLDM provides about $7.0\%$ lower energy cost compared with
the DCVs and the SQP. Since we observe close peaks over the zone temperature
(Fig. 5) and zone CO2 (Fig. 6) trajectories for these methods, we conclude
that the TLDM can maintain the same thermal comfort and IAQ but with less
energy cost. As for the computational efficiency, the average computation time
of TLDM for each executive epoch is about 5.21$\mathrm{s}$ (in parallel) with
a slight increase over the DTBSS (2.11$\mathrm{s}$). This is mainly attributed
to the LLC to achieve IAQ. However, the TLDM obviously outperforms SQP in
computation efficiency. While compared with centralized method, we imply the
sub-optimality of TLDM in energy cost is around $4\%$ in the _benchmark_.
Nevertheless, the computational benefit of TLDM is significant as the average
computation time is reduced from 575.07s to 5.21s.

###  V-B Scalability

In this part, the proposed TLDM is applied to _medium_ (10,20 zones) and
_large_ (50,100 zones) scale cases. Considering that centralized method and
SQP are now computationally intractable, we compare the TLDM with the other
three methods (i.e., DTBSS, and DCV I, II). For each case, we randomly
generate a network to represent the spatial connectivity of the zones (the
maximum number of adjacent zones zone is set as $4$). In particular, the space
and occupancy ventilation rates for the DCVs to achieve the IAQ are documented
in TABLE III, which are obtained from _off-line_ regulation. The other
parameters refer to the _benchmark_ in Section IV-A.

TABLE III: Occupancy and Space Ventilation Rates in DCV I, II  #zones |  DCV I | DCV II
---|---|---
$R_{p}$ ($\mathrm{L}$/p) |  $R_{p}$ ($\mathrm{L}$/p) |  $R_{a}$ ($\mathrm{L}\mathrm{/}\mathrm{m}^{2}$)
$10$ |  $19$ |  $15$ |  $0.03$
$20$ |  $20$ |  $19$ |  $0.03$
$50$ |  $21$ |  $19$ |  $0.03$
$100$ |  $23$ |  $21$ |  $0.03$

Similarly, we investigate the energy cost and computational efficiency of
different methods, with the results shown in TABLE IV. Compared with the
DTBSS, we observe a minor increase in energy cost and computation time with
the TLDM to achieve IAQ. However, the TLDM is still computationally efficient
and scalable by inspecting the average computation time (e.g.,
$21.68$$\mathrm{s}$ for $100$-zone case) _versus_ the decision epoch
$30$$\mathrm{min}$s. Besides, we observe similar gains of the TLDM over the
DCVs in energy cost saving: $8.0$-$9.8\%$ (DCV I) and $8.1$-$10.2\%$ (DCV II).
Moreover, except for the energy cost reduction, the TLDM is expected to be
more applicable over the DCVs as it doesn’t depend on time-consuming _off-
line_ regulation for the occupancy and space ventilation rate.

TABLE IV: Performance comparisons

Method | _Medium_ |  TC |  IAQ
---|---|---|---
$10$ | $20$
Cost (s$) |  +/- (%) |  Time (s) |  Cost (s$) |  +/- (%) |  Time (s)
TLDM |  $407.32$ |  - |  $6.55$ |  $939.83$ |  - |  $8.03$ |  Y |  Y
DTBSS |  $387.07$ |  -$4.97$ |  $2.35$ |  $887.12$ |  -$5.61$ |  $2.71$ |  Y |  N
DCV I |  $447.42$ |  +$9.84$ |  - |  $1015.50$ |  +$8.05$ |  - |  Y |  Y
DCV II |  $440.42$ |  +$8.13$ |  - |  $1020.60$ |  +$8.59$ |  - |  Y |  Y

Method | _Large_ |  TC |  IAQ
---|---|---|---
$50$ | $100$
Cost (s$)$\times 10^{3}$ |  +/- (%) |  Time (s) |  Cost (s$)$\times 10^{3}$ |  +/- (%) |  Time (s)
TLDM |  $2.65$ |  - |  $13.28$ |  $5.80$ |  - |  $21.68$ |  Y |  Y
DTBSS |  $2.54$ |  -$4.15$ |  $4.47$ |  $5.49$ |  -$5.34$ |  $6.90$ |  Y |  N
DCV I |  $2.91$ |  +$9.81$ |  - |  $6.31$ |  +$8.79$ |  - |  Y |  Y
DCV II |  $2.92$ |  +$10.19$ |  - |  $6.36$ |  +$9.66$ |  - |  Y |  Y

N=No, Y=Yes. |
---|---

Figure 8: The dynamic occupancy for the three zones among $I=50$ zones.
Figure 9: The zone air flow rates for the three zones among $I=50$ zones.
Figure 10: The temperature for the four zones among $I=50$ zones.
Figure 11: The CO2 concentration for the four zones among $I=50$ zones
Figure 12: The ventilation rate ($\bm{d}_{r}$) of the HVAC system

As an instance, we further investigate the 50-zone case study. We display the
zone occupancy, zone mass flow rates, zone temperature, and zone CO2 for three
randomly selected zones in Fig. 8-11. From Fig. 10 -11, we see both the TLDM
and the DCVs maintain the comfortable zone temperature range
($[24,26]^{\circ}$C) and zone CO2 range ($[0,800]$ ppm). Besides, we observe
close peak (800 ppm) over the zone CO2 trajectories under the TLDM and the
DCVs, which implies the close IAQ maintained by these methods. Besides, from
Fig. 10, we can have some insights in the characteristics of the TLDM to
achieve thermal comfort and IAQ. Specifically, we see the zone temperature
approaches the lower bound during the working hours with high occupancy. This
phenomenon is caused by the _two-phase_ structure of LLC in TLDM where the
zone mass flow rates is first adjusted to achieve IAQ and then the ventilation
rate is regulated if necessary.

We also compare the ventilation rates ($\bm{d}_{r}$) of different methods in
Fig. 12, which exhibit some interesting phenomenon to be interpreted. First of
all, we see the ventilation rate of TLDM corresponds well to the occupancy.
Specifically, we see relatively low ventilation rate (larger $\bm{d}_{r}$)
during the off-working hours but the opposite over the working hours (smaller
$\bm{d}_{r}$). This is reasonable as the LLC is only invoked to regulate the
ventilation rate if the CO2 is to be violated, which generally results from
high occupancy. Surprisingly, the results for DCV I and DCV II are almost
opposite and quite different from the TLDM. Specifically, we see that the
ventilation rate ($\bm{d}_{r}$) of DCV I corresponds well to the occupancy
whereas the DCV II presents the opposite. These phenomenon are attributed to
the rules that used to determine zone fresh air flow rates in the DCVs as
discussed before. For DCV I, the zone fresh air totally depends on zone
occupancy and thus we observe a synchronous pace of ventilation rate
($\bm{d}_{r}$) with the occupancy. However, for DCV II, the zone mass flow
rates are jointly determined by the occupancy and space. Therefore, during the
off-working hours with low occupancy, we can observe higher ventilation rate
as the proportions of zone fresh air flow rates dominate in the total zone
mass flows rates (higher proportion) for maintaining thermal comfort.

##  VI Conclusion

This paper studied scalable control of multi-zone HVAC systems with the
objective to reduce the energy cost for maintaining thermal comfort and IAQ
simultaneously. This problem is computationally challenging due to the complex
system dynamics. To cope with the difficulties, we proposed a two-level
distributed method (TLDM) which integrates the _upper_ and _lower_ level
control by exploiting the problem structures. Specifically, the _upper_ level
first optimizes zone mass flow rates to satisfy thermal comfort with minimal
energy cost and the _lower_ level regulates zone mass flow rates and the
ventilation rate to achieve IAQ while preserving the _near_ energy saving
performance of ULC. As both the _upper_ and _lower_ level control use
distributed computation, the proposed method is computationally efficient and
scalable. The method’s performance in energy saving and scalability was
demonstrated by numeric studies. The sub-optimality of the proposed method in
energy cost is around 4% but with significant computational benefits over the
centralized method. Compared to the distributed Token-based scheduling
strategy (DTBSS), the proposed method induces a marginal increase of energy
costs but provides IAQ. In addition, the proposed method provided 8-10% energy
savings over the demand controlled ventilation strategies (DCVs).

## References

  * [1] K. Ku, J. Liaw, M. Tsai, and T. Liu, “Automatic control system for thermal comfort based on predicted mean vote and energy saving,” IEEE Transactions on Automation Science and Engineering, vol. 12, no. 1, pp. 378–383, 2015.
  * [2] H. Mirinejad, K. C. Welch, and L. Spicer, “A review of intelligent control techniques in HVAC systems,” in 2012 IEEE Energytech, pp. 1–5, IEEE, 2012.
  * [3] A. Afram and F. Janabi-Sharifi, “Theory and applications of HVAC control systems–a review of model predictive control (MPC),” Building and Environment, vol. 72, pp. 343–355, 2014.
  * [4] A. Kelman and F. Borrelli, “Bilinear model predictive control of a HVAC system using sequential quadratic programming,” in Ifac world congress, vol. 18, pp. 9869–9874, 2011.
  * [5] N. Radhakrishnan, Y. Su, R. Su, and K. Poolla, “Token based scheduling for energy management in building HVAC systems,” Applied energy, vol. 173, pp. 67–79, 2016.
  * [6] Y. Yang, G. Hu, and C. J. Spanos, “Stochastic optimal control of HVAC system for energy-efficient buildings,” arXiv preprint arXiv:1911.00840, 2019\.
  * [7] Z. Xu, G. Hu, C. J. Spanos, and S. Schiavon, “Pmv-based event-triggered mechanism for building energy management under uncertainties,” Energy and Buildings, vol. 152, pp. 73–85, 2017.
  * [8] H. Dasi, F. Xiaowei, and C. Daisheng, “On-line control strategy of fresh air to meet the requirement of IAQ in office buildings,” in 2010 5th IEEE Conference on Industrial Electronics and Applications, pp. 845–848, IEEE, 2010\.
  * [9] W. Li, C. Koo, S. H. Cha, T. Hong, and J. Oh, “A novel real-time method for HVAC system operation to improve indoor environmental quality in meeting rooms,” Building and Environment, vol. 144, pp. 365–385, 2018.
  * [10] “Good indoor air quality leads to good decisions.” https://www.vaisala.com/sites/default/files/documents/VIM-G-HVAC-Good-Indoor-Air-Quality-Application-note-B211681EN.pdf, 2015.  Accessed Dec. 13, 2020.
  * [11] “Why measure CO2 in HVAC applications?.” https://www.co2meter.com/blogs/news/why-measure-co2-hvac-applications, 2020.  Accessed Dec. 13, 2020.
  * [12] A. Parisio, M. Molinari, D. Varagnolo, and K. H. Johansson, “A scenario-based predictive control approach to building HVAC management systems,” in 2013 IEEE International Conference on Automation Science and Engineering (CASE), pp. 428–435, IEEE, 2013.
  * [13] A. Parisio, D. Varagnolo, M. Molinari, G. Pattarello, L. Fabietti, and K. H. Johansson, “Implementation of a scenario-based MPC for HVAC systems: an experimental case study,” IFAC Proceedings Volumes, vol. 47, no. 3, pp. 599–605, 2014.
  * [14] L. Yu, D. Xie, C. Huang, T. Jiang, and Y. Zou, “Energy optimization of hvac systems in commercial buildings considering indoor air quality management,” IEEE Transactions on Smart Grid, vol. 10, no. 5, pp. 5103–5113, 2018.
  * [15] Y. Yang, G. Hu, and C. J. Spanos, “HVAC Energy cost optimization for a multi-zone building via a decentralized approach.” https://arxiv.org/abs/1905.10934, 2019.  Accessed Feb. 17, 2020.
  * [16] J. Mei and X. Xia, “Distributed control for a multi-evaporator air conditioning system,” Control Engineering Practice, vol. 90, pp. 85–100, 2019.
  * [17] C. C. Okaeme, S. Mishra, and J. T.-Y. Wen, “Passivity-based thermohygrometric control in buildings,” IEEE Transactions on Control Systems Technology, vol. 26, no. 5, pp. 1661–1672, 2017.
  * [18] C. C. Okaeme, S. Mishra, and J. T. Wen, “A comfort zone set-based approach for coupled temperature and humidity control in buildings,” in 2016 IEEE International Conference on Automation Science and Engineering (CASE), pp. 456–461, IEEE, 2016.
  * [19] N. Nassif, S. Kajl, and R. Sabourin, “Ventilation control strategy using the supply CO2 concentration setpoint,” Hvac&R Research, vol. 11, no. 2, pp. 239–262, 2005.
  * [20] T. R. Nielsen and C. Drivsholm, “Energy efficient demand controlled ventilation in single family houses,” Energy and buildings, vol. 42, no. 11, pp. 1995–1998, 2010.
  * [21] A. Hesaraki, J. A. Myhren, and S. Holmberg, “Influence of different ventilation levels on indoor air quality and energy savings: A case study of a single-family house,” Sustainable cities and society, vol. 19, pp. 165–172, 2015.
  * [22] Z. Wang and L. Wang, “Intelligent control of ventilation system for energy-efficient buildings with co2 predictive model,” IEEE Transactions on Smart Grid, vol. 4, no. 2, pp. 686–693, 2013.
  * [23] N. Nassif, “A robust CO2-based demand-controlled ventilation control strategy for multi-zone HVAC systems,” Energy and buildings, vol. 45, pp. 72–81, 2012.
  * [24] M. Marinov, T. Djamiykov, B. Ganev, and V. Zerbe, “Sensor-based multi-zone demand-controlled ventilation,”
  * [25] K. Shan, Y. Sun, S. Wang, and C. Yan, “Development and in-situ validation of a multi-zone demand-controlled ventilation strategy using a limited number of sensors,” Building and environment, vol. 57, pp. 28–37, 2012.
  * [26] A. Parisio, D. Varagnolo, D. Risberg, G. Pattarello, M. Molinari, and K. H. Johansson, “Randomized model predictive control for HVAC systems,” in Proceedings of the 5th ACM Workshop on Embedded Systems For Energy-Efficient Buildings, pp. 1–8, ACM, 2013.
  * [27] Hattersley, “Variable flow vs constant flow.” https://www.hattersley.com/page/hvac/variable-flow-systems.  Accessed Dec. 20, 2019.
  * [28] X. Zhang, W. Shi, B. Yan, A. Malkawi, and N. Li, “Decentralized and distributed temperature control via HVAC systems in energy efficient buildings,” arXiv preprint arXiv:1702.03308, 2017.
  * [29] Y. Lin, T. Middelkoop, and P. Barooah, “Issues in identification of control-oriented thermal models of zones in multi-zone buildings,” in Decision and Control (CDC), 2012 IEEE 51st Annual Conference on, pp. 6932–6937, IEEE, 2012.
  * [30] M. Maasoumy, A. Pinto, and A. Sangiovanni-Vincentelli, “Model-based hierarchical optimal control design for HVAC systems,” in ASME 2011 Dynamic Systems and Control Conference and Bath/ASME Symposium on Fluid Power and Motion Control, pp. 271–278, American Society of Mechanical Engineers, 2011\.
  * [31] I. Korolija, L. Marjanovic-Halburd, Y. Zhang, and V. I. Hanby, “Influence of building parameters and hvac systems coupling on building energy performance,” Energy and Buildings, vol. 43, no. 6, pp. 1247–1253, 2011\.
  * [32] G. P. McCormick, “Computability of global solutions to factorable nonconvex programs: Part I—convex underestimating problems,” Mathematical programming, vol. 10, no. 1, pp. 147–175, 1976.
  * [33] N. Chatzipanagiotis, D. Dentcheva, and M. M. Zavlanos, “An augmented lagrangian method for distributed optimization,” Mathematical Programming, vol. 152, no. 1-2, pp. 405–434, 2015.
  * [34] Z. Sun, S. Wang, and Z. Ma, “In-situ implementation and validation of a CO2-based adaptive demand-controlled ventilation strategy in a multi-zone office building,” Building and Environment, vol. 46, no. 1, pp. 124–133, 2011.
  * [35] J. Nocedal and S. J. Wright, “Sequential quadratic programming,” Numerical optimization, pp. 529–562, 2006.
  * [36] Y. Kawajir, C. Laird, and A. Wachter, “Introduction to IPOPT: A tutorial for downloading, installing, and using IPOPT,” 2006.