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

1 MB·

by  Takeshi Hatanaka

An Integrated Design of Optimization and Physical Dynamics for Energy
Efficient Buildings: A Passivity Approach
Takeshi Hatanaka, Xuan Zhang, Wenbo Shi, Minghui Zhu and Na Li
Abstract—In thispaper, weaddressenergy management for thendemonstratedthroughsimulation.Theauthorsof[8],[9]
heating, ventilation, and air-conditioning (HVAC) systems in incorporate the grid dynamics into the optimization process
buildings,andpresentanovel combinedoptimizationandcon-
by identifying the physical dynamics with a subprocess
7 trol approach. We first formulate a thermal dynamics and an
of seeking the optimal solution. A scheme to eliminate
1 associated optimization problem. An optimization dynamics is
0 thendesignedbased onastandard primal-dualalgorithm, and structural constraints required in [8], [9] is presented by
2 its strict passivity is proved. We then design a local controller Zhang et al. [10] while instead assuming measurements of
and prove that the physical dynamics with the controller is disturbances.Asimilarapproachisalsotakenforpowergrid
n
ensured to be passivity-short. Based on these passivity results,
a control in Stegink et al. [11].
we interconnect the optimization and physical dynamics, and
J Inthispaper,weaddressintegrateddesignofoptimization
proveconvergenceoftheroomtemperaturestotheoptimalones
9 definedforunmeasurabledisturbances.Finally,wedemonstrate andphysicaldynamicsforHVACcontrolbasedonpassivity,
1 the present algorithms through simulation. whereweregardtheoptimizationprocessasadynamicalsys-

tem similarly to [7]–[11]. Interconnectionsof such dynamic
] I. INTRODUCTION
Y HVAC optimization with a building dynamics are partially
Stimulatedby strongneedsfor reducingenergyconsump-
studiedin[12],wheretemperaturedataforallzonesandtheir
S
tion of buildings, smart building energy management algo-
. derivatives in the physics side are fed back to the dynamic
s rithms have been developed both in industry and academia.
optimizationprocessto recoverthe disturbanceterms. How-
c
Inparticular,abouthalfofthecurrentconsumptionisknown
[ ever,suchdata arenotalwaysavailableinpracticalsystems.
to be occupied by heating, ventilation, and air-conditioning
We thus present an architecture relying only on temperature
1
(HVAC) systems, and a great deal of works have been
data of a subgroup of zones with HVAC systems.
v
devotedtoHVACoptimizationandcontrol[1].Inthispaper,
7 The contents of this paper are as follows. A thermal
7 we address the issue based on a novel approach combining dynamic model with unmeasurable disturbances and an as-
5 optimization and physical dynamics. sociated optimization problem are first presented. We then
5 Interplays between optimization and physical dynamics
formulate an optimization dynamics based on the primal-
0
have been most actively studied in the field of Model
dualgradientalgorithm[13]. Thedesigneddynamicsis then
.
1 Predictive Control (MPC), which has also been applied to
provedto be strictly passive from a transformeddisturbance
0 building HVAC control [1]–[6]. While the MPC approach
estimate to an estimated optimalroomtemperature.We next
7
regards the optimization process as a static map from phys-
1 designacontrollersothattheactualroomtemperaturetracks
ical states to optimal inputs, another approach to integrating
: agivenreference,andproducesadisturbanceestimate.Then,
v optimization and physical dynamicsis presented in [7]–[11]
the physical dynamics is proved to be passivity-short from
i
X mainly motivated by power grid control. There, the solution the reference to the disturbance estimate. From these two
processoftheoptimizationisviewedasadynamicalsystem,
r passivity-related results, we then interconnect the optimiza-
a and the combination of optimization and physical dynamics
tion and physical dynamics, and prove convergence of the
is regardedas an interconnectionofdynamicalsystems. The
actual room temperature to the optimal solution defined for
benefits of the approach relative to MPC are as follows.
the unmeasurable actual disturbance. Finally, the presented
First,theapproachallowsonetoavoidcomplicatedmodeling
algorithm is demonstrated through simulation.
and prediction of factors hard to know in advance, while
MPCneedstheirmodelstopredictfuturesystemevolutions. II. PROBLEM SETTINGS
Second, since the entire system is a dynamical system, its A. Preliminary
stability and performance are analyzed based on unifying
In this section, we introduce the concept of passivity.
dynamical system theory.
Consider a system with a state-space representation
In this line of works, Shiltz et al. [7] addressessmartgrid
control, and interconnects a dynamic optimization process x˙ =φ(x,u), y =ϕ(x,u), (1)
anda locallycontrolledgriddynamics.The entireprocessis where x(t) ∈ Rn is the state, u(t) ∈ Rp is the input and
y(t)∈Rp is the output. Then, passivity is defined as below.
T.HatanakaiswithSchoolofEngineering,TokyoInstituteofTechnology,
Tokyo 152-8552, JAPAN ([email protected]). X. Zhang, W. Shi
andN.LiarewithElectrical Engineering andAppliedMathematics ofthe
School of Engineering and Applied Sciences, Harvard Univ., 33 Oxford
Definition 1 The system (1) is said to be passive if there
St,Cambridge,MA02138,USA.M.ZhuiswithDepartment ofElectrical
Engineering, Pennsylvania StateUniv.,University Park,PA16802,USA. exists a positive semi-definite function S : Rn → R+ :=
[0,∞), called storage function, such that T ) (i = 1,2,...,n ), 1 is the n-dimensional real vector
i 1
t whose elements are all 1, and B =[In1 0]⊤ ∈Rn×n1.
S(x(t))−S(x(0))≤ yT(τ)u(τ)dτ (2) We next linearize the model at around an equilibrium as
Z
0
holds for all inputs u : [0,t] → Rp, all initial states x(0) ∈ Cδ˙T =RδTa1−RδT −LδT +BG(T¯)δm−U¯δT +Bδq,
Rnandallt∈R+.Inthecaseofthestaticsystemy =ϕ(u), where δT, δm, δTa and δq describe the errors from the
it is passive if yTu = ϕT(u)u ≥ 0 for all u ∈ Rp. The equilibrium states and inputs and U¯ ∈ Rn×n is a diagonal
system (1) is also said to be output feedback passive with
matrix whose diagonal elements are m¯ ,...,m¯ ,0,...,0,
index ε>0 if (2) is replaced by 1 n1
where m¯ is the i-th element of the equilibrium input m¯.
i
t Using the variable transformations
S(x(t))−S(x(0))≤ yT(τ)u(τ)−εky(τ)k2dτ (3)
Z0 x := C1/2δT, u:=B⊤C−1/2BG(T¯)δm,
If the right-hand side of (3) is changed as w := C−1/2RδTa1, w :=B⊤C−1/2Bδq, (8)
a q
t
S(x(t))−S(x(0))≤ yT(τ)u(τ)+εku(τ)k2dτ (4) (6) is rewritten as
Z
0
x
withε>0,thesystemissaidtobepassivity-short,andthen x˙ =−Ax+Bu+Bw +w , x:= 1 (9)
ε is called impact coefficient [18]. q a (cid:20)x2(cid:21)
B. System Description where A:=C−1/2(R+L+U¯)C−1/2, x1 ∈Rn1 and x2 ∈
Rn2. Remark that the matrix A is positive definite [19].
In this paper, we consider a building with multiple zones
From control engineering point of view, x is the system
i = 1,2,...,n. The zones i = 1,2,...,n are divided into
state,uisthecontrolinput,andw andw aredisturbances.
q a
two groups:The first groupconsists of zones equippedwith
Wesupposethatw ismeasurableaswellasx .Meanwhile,
a 1
VAV (Variable Air Volume) HVAC systems whose thermal
it is in general hard to measure the heat gain w .
q
dynamicsisassumedtobemodeledbytheRCcircuitmodel
[2]–[6], [12] as C. Optimization Problem
C T˙ = Ta−Ti + Tj −Ti +a (Ts−T )m +q , (5) Regarding the above system, we formulate the optimiza-
i i R R i i i i i tion problem to be solved as follows:
i jX∈Ni ij
wrahteereatTzioinsethie, tTema piserathtuereamofbzieonntetie,mmpeiriastutrhee, mTasssisfltohwe z=[zx⊤ zmu⊤]i⊤n∈Rn+n1kzx1−hk2+f(zu) (10a)
i subject to: g(z )≤0 (10b)
air temperature supplied to zone i, which is treated as a u
constant throughout this paper, qi is the heat gain at zone −Azx+Bzu+Bdq+da =0 (10c)
i from external sources like occupants, C is the thermal
i Thevariablesz :=[z⊤ z⊤]⊤ (z ∈Rn1,z ∈Rn2),and
capacitance,R isthethermalresistanceofthewall/window, x x1 x2 x1 x2
R is the thermi alresistance betweenzone i and j and a is zu ∈ Rn1 correspond to the zone temperature x and mass
ij i flowrateuafterthetransformation(8).Theparametersd ∈
the specific heat of the air. q
The second group is composed of other spaces such as Rn1 andda ∈Rn areDCcomponentsofthedisturbanceswq
and w , respectively. These variables are coupled by (10c)
walls and windows whose dynamics is modeled as a
which describes the stationary equation of (9).
C T˙ = Ta−Ti + Tj −Ti. (6) Throughoutthis paper, we assume the following assump-
i i R R tion.
i jX∈Ni ij
Roomsnotinusecanbecategorizedintothisgroup.Without
Assumption 1 The problem (10) satisfies the following
loss of generality,we assume that i=1,2,...,n1 belongto properties: (i) f : Rn1 → R is convex and its gradient is
the first groupand i=n +1,...,n (n :=n−n ) belong
1 2 1 locallyLipschitz,(ii)everyelementoftheconstraintfunction
to the second. Remark that the system parameters Ci, Ri, g :Rn1 →Rc isconvexanditsgradientislocallyLipschitz,
andThReijcoclalenctbiveeiddeynntaifimeidcsusoifng(5t)haentdoo(6lb)oixsidnes[c1r4ib].ed as and (iii) there exists zu ∈Rn1 such that g(zu)<0.
CT˙ =RTa1−RT −LT +BG(T)m+Bq, (7) The first term of (10a) evaluates the human comfort,
where h ∈ Rn1 is the collection of the most comfortable
where T, q and m are collections of T (i = 1,2,...,n), temperatures for occupants in each zone, which might be
i
q (i=1,2,...,n ), and m (i=1,2,...,n ) respectively. determined directly by occupants in the same way as the
i 1 i 1
The matrices C and R are diagonal matrices with diagonal current systems, or computed using human comfort metrics
elements C (i = 1,2,...,n) and 1 (i = 1,2,...,n), likePMV(PredictedMeanVote).Thequadraticfunctionfor
i Ri
respectively. The matrix L describes the weighted graph the error is commonly employed in the MPC papers [3],
Laplacian with elements 1 , G(T) ∈ Rn1×n1 is a block [22]–[24] and [12]. Note that it is common to put weights
Rij
diagonal matrix with diagonal elements equal to a (Ts − on each element of z −h to give priority to each zone,
i i x1
but this can be done by appropriately scaling each element
of z in the function f. This is why we take (10a).
u
The function f is introduced to reduce power consump-
tion. The papers [3], [22]–[24] simply take a linear or
quadratic function of control efforts as such a function
and then Assumption 1(i) is trivially satisfied. Also, as
mentioned in [5], [12], the power consumption of supply
fansisapproximatedbythecubeofthesumofthemassflow
rates, which also satisfies Assumption 1(i). A simple model
ofconsumptionatthecoolingcoilisgivenbytheproductof Fig. 1. Block diagram of optimization dynamics, which is passive from
themassflowratemi and|Tis−Ta|[25],whichalsobelongs v˜o=vo−vo∗ y˜o=yˆo−yo∗ withvo∗:=Mdqandyo∗:=zx∗1(Lemma2).
totheintendedclass1.Theconstraintfunctiong :Rn1 →Rc
reflects hardware constraints and/or an upper bound of the
primal-dualgradientalgorithm[13]asoneofsuchsolutions2.
powerconsumption.For example,the constraintsin [12] are
However,it is hard to obtain d since w is notmeasurable.
reduced to the above form. q q
We thus need to estimate d from the measurements of
The objective of this paper is to design a controller so q
physical quantities. This motivates us to interconnect the
as to ensure convergenceof the actual room temperature x
1 physical dynamics with the optimization dynamics.
to the optimal room temperature z∗ , the solution to (10),
x1 Taking account of the above issues, we present
without direct measurements of w .
q
zˆ˙ = −α{M2(zˆ +dˆ )+N(w −¯h)+∇f(zˆ )+p},(14a)
u u q a u
III. OPTIMIZATION DYNAMICS
λˆ˙ = [g(zˆ )]+, p=∇g(zˆ )λˆ, (14b)
A. Optimization Dynamics u λˆ u
In this subsection, we present a dynamics to solve the where zˆ and λˆ are estimates of z∗ andλ∗ respectively,and
u u
above optimization problem. Before that, we eliminate z α>0. The notation [b]+ for real vectors a,b with the same
x a
from (10) using (10c). Then, the problem is rewritten as dimensionprovidesa vectorwhose l-thelement,denotedby
([b]+) , is given by
min kB⊤A−1(Bz +Bd +d −h¯)k2+f(z ) (11a) a l
u q a u
zu∈Rn1 ([b]+) = 0, if al =0 and bl <0 , (15)
subject to: g(zu)≤0 (11b) a l (cid:26) bl, otherwise
with h¯ =ABh. It is easy to confirm from positive definite- where al,bl are the l-th element of a,b, respectively. Note
ness of A that the cost function of (11) is strongly convex. that(14)isdifferentfromtheprimal-dualgradientalgorithm
In this case, (11) has the unique optimal solution, denoted for (12) in that the term da is replaced by the measurement
by zu∗, and it satisfies the following KKT conditions [26]. wa, and dq is replaced by its estimate dˆq whose production
will be mentioned later. The system is illustrated in Fig. 1.
M2(z∗+d )+N(d −h¯)+∇f(z∗)+∇g(z∗)λ∗ =0,
u q a u u
(12a) B. Passivity Analysis for Optimization Dynamics
g(z∗)≤0, λ∗ ≥0, λ∗◦g(z∗)=0, (12b) Hereafter, we analyze passivity of the above optimization
u u
process assuming that w is constant. In this case, w ≡d
a a a
where M = B⊤A−1B, N := B⊤A−1BB⊤A−1, and ◦ holds. In practice, the disturbance w , namely the ambient
a
represents the Hadamard product. Since (11) is essentially temperature Ta, is time-varying but the following results
equivalent to (10), the solution to (11) also provides a solu- are applied to the practical case if w is approximated by
a
tionforthe problem(10).Precisely,if we definezx∗ ∈Rn as a piecewise constant signal, which is fully expected since
zx∗ :=A−1(Bzu∗+Bdq+da),thepairz∗ :=[(zx∗)⊤ (zu∗)⊤]⊤ Ta usually varies slowly.
is a solution to (10). In the sequel, we also use the notation Under the above assumption, we define the output ν :=
zx∗ := [(zx∗1)⊤ (zx∗2)⊤]⊤ (zx∗1 ∈ Rn1 and zx∗2 ∈ Rn2). It is −M2(zˆu+dˆq) for(14).We thenhavethe followinglemma.
then easy to confirm
z∗ =M(z∗+d )+B⊤A−1d . (13)
x1 u q a Lemma 1 Consider the system (14) with w ≡ d and
a a
Givend andd ,itwouldbeeasytosolve(12).However, λˆ(0) ≥ 0. Then, under Assumption 1, it is passive from
q a
in the practical applications, it is desired that dq and da are d˜q := dˆq − dq to −ν˜, where ν˜ := ν − ν∗ and ν∗ :=
updatedinrealtimeaccordingtothechangesofdisturbances. −M2(zu∗+dq).
In this regard, it is convenient to take a dynamic solution
Proof: See Appendix I.
processofoptimizationsinceittriviallyallowsonetoupdate
the parameters in real time. In particular, we employ the
2Another benefit of using the dynamic solution is that it provides a
distributedsolutionwhenthepresentresultsareextendedtoamoreglobal
1Forsimplicity,weskipthedependenceofdaonf,butsubsequentresults problem, although it exceeds the scope of this paper. Please refer to [20]
areeasily extended tothecasethatf depends onda. formoredetails ontheissue.
Let us next transform the output ν to
y :=−M−1ν+B⊤A−1d =M(zˆ +dˆ )+B⊤A−1d ,
o a u q a
y∗ :=−M−1ν∗+B⊤A−1d .
o a
Comparing(13)andtheabovedefinitionofy , thesignaly
o o
isregardedasanestimateofz∗ .Wealsodefinev :=Mdˆ
x1 o q
and v∗ :=Md . Then, we can prove the following lemma.
o q
Lemma 2 Consider the system (14) with w ≡ d and Fig.2. Blockdiagramofthephysicaldynamicswiththelocalcontroller,
pλˆa(s0s)iv≥e f0r.oTmhev˜n,:=undver−Avs∗sutmopy˜tio:n=1y, it−isy∗ouwtapituhtifnedeaedxba1c.k wanhdicyhp∗is:=pas−siMvitdy-qsh(oLretmfrmoma5v˜)p. =vp−vp∗toy˜p=yp−yp∗withvp∗ :=zx∗1
o o o o o o
Proof: It is easy to see from ν∗ =−M2(z∗+d ) and
u q Substituting (17) into (9) yields
(13) that y∗ = z∗ . From (32), y˜ = −M−1ν˜ and v˜ =
o x1 o o
Md˜q, we have the following inequality. x˙ = −Ax+kPB(r−x1)+κBr+Bξ
+Bw +(BF +I )w , (19a)
D+S ≤y˜⊤v˜ −ky˜ k2 (16) q n a
o o o o ξ˙ = k (r−x ). (19b)
I 1
Integrating this in time completes the proof.
For the system, we have the following lemma.
IV. PHYSICAL DYNAMICS
Lemma 3 The steady states x∗ and ξ∗ of (19) for r ≡ r∗,
In this section, we design a physical dynamics and prove
w ≡d and w ≡d are given as follows.
a a q q
its passivity. A passivity-based design for the model (9) is
porredseernttoedinitner[c2o1n]n,ebcuttitwweimthotdhiefyotphteimcioznattrioonl adrycnhaitmecitcusr.e in x∗ = −F⊤r∗+(cid:20)A−310Bc⊤(cid:21)da, Bc :=(cid:20)I0n2(cid:21)∈Rn×n2
ξ∗ = (M−1−κI )r∗−d . (20)
A. Controller Design n2 q
Equation (19) is now rewritten as
In this subsection, we design a controller to determine
the input u so that x1 tracks a reference signal r. Here we x˙ = −A¯x+BM−1ζ+Bwq+F¯wa (21a)
assume the following assumption, where ξ˙ = k (r−x ), ζ =k¯ M(r−x )+Mξ, (21b)
I 1 P 1
A A⊤
A= 1 2 , A ∈Rn1×n1, A ∈Rn2×n2. where
(cid:20)A2 A3(cid:21) 1 3 A¯:=A−κBB⊤, k¯ :=k +κ, F¯ := −A⊤2A−31 B⊤.
P P (cid:20) In2 (cid:21) c
Assumption 2 ThematrixMA +A M ispositivedefinite. Remark that, at the steady state, the variable ζ is equal to
1 1
ζ∗ :=Mξ∗ =Kr∗−Md , K :=I −κM. (22)
This property does not always hold for any positive definite q n1
matrices A and M, but it is expected to be true in many B. Passivity Analysis for Physical Dynamics
1
practical cases since the diagonal elements tend to be dom-
Inthissubsection,weanalyzepassivityofthesystem(21)
inant both for A and M = (A −A⊤A−1A )−1 in this
1 1 2 3 2 assuming that wq and wa are constant. The case of the time
application [12]. Note that this assumption holds in the full
varying w will be treated in the end of the next section.
actuation case (A =A,M =A−1). q
1 Choose ζ as the output and prove passivity as follows.
Inspired by the fact that many existing systems employ
Proportional-Integral (PI) controllers (with logics) as the Lemma 4 Consider the system (21) with r ≡ r∗, w ≡ d
a a
local controller, we design the following controller adding and w ≡ d . Then, under Assumption 2, the system is
q q
reference and disturbance feedforward terms. passive from r˜:=r−r∗ to ζ˜:=ζ−ζ∗.
ξ˙ = kI(r−x1) (17a) Proof: See Appendix II.
u = k (r−x )+ξ+κr+Fw (17b) Remark that this lemma holds regardless of the value of r∗.
P 1 a
To extract the term Md from (22), we define the output
where k >0, k >0 and F :=[−I A⊤A−1] and I is q
the n1-bPy-n1 idenItitymatrix.Thefeednf1orw2ard3gainκ>n01 is yp :=ζ−Kr, yp∗ :=ζ∗−Kr∗ =−Mdq, (23)
selected so that and v :=r and v∗ :=z∗ . Then, we have the following.
p p x1
P :=MA +A M −2κM >0. (18)
1 1 Lemma 5 Consider the system (21) with r ≡ r∗, w ≡ d
a a
Such a κ exists under Assumption 2. and w ≡ d . Then, under Assumption 2, the system from
q q
v˜ := v − v∗ to y˜ := y − y∗ is passivity-short with
p p p p p p
the impact coefficient 1−κσ, where σ > 0 is the minimal
eigenvalue of M.
Proof: If we take r∗ = z∗ , we have y˜ = ζ˜−Kv˜ .
x1 p p
Substituting this into (39) yields
S˙ ≤y˜⊤v˜ +v˜⊤Kv˜ −k (v˜ −x˜ )⊤M(v˜ −x˜ ) (24)
p p p p p P p 1 p 1
≤y˜⊤v˜ +(1−κσ)kv˜ k2−k σkv˜ −x˜ k2 (25)
p p p P p 1
This completes the proof.
V. INTERCONNECTIONOFOPTIMIZATIONANDPHYSICAL
DYNAMICS
Let us interconnect the optimization dynamics (14) and
physical dynamics (21). Remark that the stationary value of
y , y∗ = z∗ , is equivalent to that of v , v∗ = z∗ . Also, Fig.3. Hierarchical controlarchitecture.
o o x1 p p x1
v∗ = −y∗ holds. Inspired by these facts, we interconnect
o p
these systems via the negativefeedback as v =−y , v =
o p p
We give some remarks on the present architecture.
y . We then have the following main result of this paper.
o
Figs. 1 and 2 are oriented by theoretical analysis, but
Theorem 1 Suppose that λˆ(0)≥0, w ≡d and w ≡d . the implementation does not need to follow the information
a a q q
processing in the figures. Indeed, the interconnected system
Then, if Assumptions 1 and 2 hold, the interconnection of
(14)and(21)viav =−y , v =y ensuresthatx →z∗ . is equivalently transformed into Fig. 3. If we let the oper-
o p p o 1 x1
ations shaded by dark gray be executed in the high-level
Proof: Define S := S +S . Then, combining (16), controller, the low-level controller can be implemented in a
o p
(25) and v =−y , v =y yields decentralized fashion similarly to the existing systems.
o p p o
In Fig. 3, both of the high-level and low-level controller
D+S ≤ −κσky˜ k2−k σkv˜ −x˜ k2. (26)
o P p 1 withthephysicaldynamicsarebiproperandhenceaproblem
Thismeansthatbothofy˜ =y −z∗ andy˜ −x˜ =y −x ofalgebraicloopscanoccur.Thishoweverdoesnotmatterin
o o x1 o 1 o 1
belong to class L . Since S is positive definite, all of the practice since the information transmissions between high-
2
statevariableszˆu,λˆ,xandξ belongtoL∞.From(14),y˜˙o = and low-level processes usually suffer from possibly small
−M−1ν˙ = Mzˆ˙ is bounded. Also, (19) means that x˜˙ = x˙ delays. Although the high-level controller itself contains an
u
is bounded and hence y˙ −x˙ is bounded. Thus, invoking algebraic loop, it is easily confirmed that the loop can be
o 1
Barbalat’slemma,wecanprovey −z∗ →0, y −x →0, solved by direct calculations of the algebraic constraint.
o x1 o 1
which means x →z∗ . This completes the proof. The transfer functionfrom w to the disturbanceestimate
1 x1 q
Lyapunov stability of the desirable equilibrium, tuple of dˆ , roughly speaking, is almost the same as the comple-
q
z∗,λ∗,x∗ and ξ∗, is also proved in the above proof. mentary transfer function and hence only the low frequency
u
It is to be emphasized that the optimal solution is depen- components are provided by the physical dynamics. This is
dentonthe unmeasurabledisturbance.Nevertheless,conver- why dˆ is regarded as an estimate of the DC component
q
gence to the solution is guaranteed owing to the feedback of w . The cutoff frequency of the disturbance can be
q
path from physics to optimization. in principle tuned by k and k , but, once a closed-loop
P I
The above results are obtained assuming that both of w system is designed, the cutoff is also automatically decided.
q
and w are constant. This is likely valid for w since the It is however not always a drawback at least qualitatively.
a a
ambient temperature is in general slowly varying. However, Actually,evenifoptimalsolutionsreflectingmuchfasterdis-
the heat gain w may contain high frequency components. turbance variations are provided, the physical states cannot
q
To address the issue, we decompose the signal w into the respond to variations faster than the bandwidth.
q
DCcomponentsdq andothersw˜q aswq =dq+w˜q.Wealso Itis a consequenceof the internalmodelcontrolandcon-
assume that w˜q belongs to an extended L2 space [16]. The stant disturbancesthat the actual disturbance wq is correctly
followingcorollarythenholds,whichisprovedfollowingthe estimated. A controlarchitecture based on a similar concept
proof procedure of the well-known passivity theorem [16], is presented in Section VII of Stegink et al. [11]. However,
[17] and using the fact that the right-hand side of (26) is it is clear that the problem (10) does not meet the structural
upper bounded by −κκ+kPkσPkx˜1k2. constraints assumed in [11] and hence the architecture in
[11] cannot be directly applied to our problem.
Corollary 1 Suppose that λˆ(0) ≥ 0, w ≡ d and w = Zhang et al. [12] present another kind of interconnec-
a a q
d +w˜ . Then, if Assumptions 1 and 2 hold, the intercon- tion between physical and optimization dynamics based on
q q
nected system (14), (21) and v = −y , v = y from w˜ a quasi-disturbance feedforward, where the disturbance is
o p p o q
to x˜ =B⊤x˜=x −z∗ has a finite L gain. computed by state measurements and their derivatives, and
1 1 x1 2
then fed back to the optimization dynamics.The differences ee]35 0.7
of the present scheme from [12] are listed as follows: egr 0.6
d
The approach of [12] requires the measurements of state p.[c-30 kW]0.5
vanardi/aobrlewsa.llIsf,itthsetecshtantoesloginiccaluldfeeastiebmiliptyermatuayrebseopfrowblienmdoawtics ent tem25 heat[0.4
bi 0.3
or at least increases the system cost. On the other hand, our m
a
20 0.2
approach needs only x1 which is usually measurable. 0 100 200 300 400 500 0 100 200 300 400 500
time[min] time[min]
Since there is no sensor to measure x˙, it has to be com-
puted using the difference approximation, which provides Fig.4. Ambienttemperature (left)andexternal heats(right).
approximationerrors.Meanwhile,thepresentapproachdoes
0.7
not need such an approximation.
0.6
0.6
stuernIbsnaon[r1c2ens]o,isatherese dadinifrdfeecrhtelinygchesefanrpetpqrtuooexnitcmhyeatcoioopmntipmeorirnzoearntsitotsongoedftyhtnhearemwdiicitssh-, heat[kW]00..45 heat[kW]00..5568
which may cause fluctuations for the output and internal 0.3 0.54
variables in the optimization process unless it is carefully 0.2 0.52
designed in the sense of the noise reduction. Adding a 0 100 200 300 400 500 140 160 180 200 220
time[min] time[min]
low-pass filter to the computed disturbance might eliminate
e]
these undesirablefactors. However,the filter is notdesigned gre28 ee]28
taiohnnfedduenqtphcueeenarsfidtiael-tinfenetrtleiyidessfooifinnrscwtxla˙uabrdadineliddtsyyitnhsottefoesmtthyhesebteeelnmoctooirmpme.oesMdsyesealtaesnmifwneceihnedi,lbteiha,ncetkhtpherisesnysocseatinesscemees, nce temp. [c-de2246 m temp. [c-degr2246
are automatically rejected by the physical dynamics in our efere22 roo22
r
algorithm. 0 100 200 300 400 500 0 100 200 300 400 500
time[min] time[min]
VI. SIMULATION Fig. 5. Time responses of the estimated heat gain (top 2 figures), the
estimated optimal room temperatures (bottom-left) and actual room tem-
In this section, we demonstrate the presented control peratures(bottom-right). Inallfigures,thesolidcurvesshowtheresponses
deliveredbytheproposedmethod,andthedottedoneswithlightcolorsare
architecture through simulation. For this purpose, we build
thosebythedisturbancefeedforwardscheme.Inthetopfigures,thedotted
a building on 3D modeling software SketchUp (Trimble lines coincide withtheactual heatgain.
Inc.), which contains three rooms (n = 3) and other 46
1
zones (n = 46) including walls, ceilings, and windows.
2
The buildingmodelis then installed into EnergyPlus[27] in connection from optimization to physics. It is to be noted
order to simulate the evolution of zone temperatures. Then, that it is hard to implement this in practice.
usingtheacquireddata,we identifythemodelparametersin The trajectories of the estimated heat gains, the estimated
(5) and (6) via BRCM toolbox [14]. optimal room temperatures and temperatures T ,T ,T are
1 2 3
We next specify the optimization problem (10). All the illustrated by solid curves in Fig. 5. The dotted lines with
elements of h are set to 22C1/2◦C and we take f(z ) = lightcolorsshowthetrajectoriesusingtheabovedisturbance
u
150kz k2. The constraints are also selected as |z | ≤ feedforwardscheme.Weseefromthetopleftfiguresthatthe
u ui
0.61 i = 1,2,3 and 3 |z | ≤ 1.25, where z is the presentedalgorithmalmostcorrectlyestimatethedisturbance
i=1 ui ui
i-th element of zu. CPollecting these constraints, we define wq. It is also observed from the right fine-scale figure that
thefunctiong.However,sinceitturnsoutthatdirectlyusing high frequency signal components are filtered out in the
g(z ) ≤ 0 has a response speed problem in penalizing the case of our algorithm. The trajectories in the bottom-left
u
constraintviolation in the primal-dualalgorithm,we instead figure sometimes get far from 22◦C since the constraint
take the constraint θg(z )≤0 with θ =15, which does not 3 |z | ≤ 1.25 gets active during the periods due to
u i=1 ui
essentially change the optimization problem. tPhe high ambient temperature and heat gains. We see from
In the simulation, we take the feedback gains k =6.0× the bottom figures that the response of the present method
P
10−2 and k = 1.0×10−3, and κ = 1.0× 10−3, which to the constraint violations is slower than the disturbance
I
are tuned so that the peak gain of σ-plot from r to x is feedforward scheme because high frequency components of
1
smaller than the well-known criterion. It is then confirmed the disturbance are filtered out by the physical dynamics.
that Assumption 2 is satisfied. Let us next show that the high sensitivity of disturbance
In the following simulation, we use the disturbance data feedforward can cause another problem. Here, small noises
shown in Fig. 4. Here, we compare the results with the are added to the heat gain w at every 20s, whose absolute
q
ideal case that the disturbance w is directly measurable. value is upper bounded by 1.0 × 10−3. We then run the
q
In this case, the feedback path from the physical dynamics abovetwoalgorithms.Theresultingresponsesareillustrated
to optimizationisnotneededandhencewe takethe cascade in Fig. 6. It is observed from the top-left figure that the
0.7 signedaprimal-dualalgorithm-basedoptimizationdynamics
0.6
0.6 and a local physical control system, and proved the system
heat[kW]00..45 heat[kW]00..5568 porpotipmeriztiaetsiornelaanteddphtoyspicaaslsdivyintya.mWices,tahnednpirnotveercdocnonnevcetregdenthcee
of the room temperatures to the optimal ones. We finally
0.3 0.54
demonstrated the present algorithms through simulation.
0.2 0.52
0 100 200 300 400 500 140 160 180 200 220
time[min] time[min]
e]30 30 REFERENCES
gre28 ee]
nce temp. [c-de22220246 m temp. [c-degr222468 [[21]] EsSFAy.tn.asvOAutiecrlfmdohrena,swmmBeue,.ArnFLtte.,erlJ,ehvavmonAilaea..bwnP7in-2aoSr,afhinspamidproio.,fMd,3Ce“.4l.T3MNp–h.ro3eerJo5daor5irync,i,eta2is“vn0,Ued1Ds4cae.o.pnGoptflyriocmalaliot(sidMotrenaPlssCp,or)Mfe,”dH.iBcVGutAiwivlCedeirccndoogennra,ttrrnVoodll.
refere1168 roo22 EannderwgyeaatnhderBfuoirledcinasgtss,fvoorl.e4n5e,rgpyp.e1ffi5c–i2e7n,t2b0u1il2d.ing climate control,”
0 100 200 300 400 500 0 100 200 300 400 500
time[min] time[min] [3] A.Aswani,N.Master,J.Taneja,D.CullerandC.Tomlin,“Reducing
transient and steady state electricity consumption in HVAC using
Fig. 6. Time responses in the presence of noises on the heat gain wq, learning-based model-predictive control,” Proc. IEEE, vol. 100, no.
whereeverylinehasthesamemeaningasFig.5. 1,pp.240–253, 2012.
[4] Y. Ma, A. Kelman, A. Daly and F. Borrelli, “Predictive control for
e] energyefficientbuildingswiththermalstorage:modeling,stimulation,
gre28 ee]28 andexperiments,”IEEEControlSyst.vol.32,no.1,pp.44–64,2012.
nce temp. [c-de2246 m temp. [c-degr2246 [[56]] fTbYSoar..rasMGenbdsoau.,yoiCJland.lo,iMonnHctgarc.otuuHlIpsnSVakgynoAlscetCayyenmdiasnsynFfsoTd.terBemcPmoh.arnstBr:ioeoallCnolriog,ofoym“oa,Srhpvte,oloencl“x.ehZir2atgoys3ynt,iecane-nfolmfied.vco1eicd,eloenpnlctpsop.benrur1etvri0dlao1idtcl–iitsn1iamvgl1egs,6”,o”c,ro2iIPtnE0hrt1Emoro5cEsl..
efere22 roo22 AmericanControlConf.,pp.3063–3068, 2012.
r [7] D.J. Shiltz, M. Cvetkovic and A.M. Annaswamy, “An integrated
0 100 200 300 400 500 0 100 200 300 400 500
time[min] time[min] dynamic market mechanism for real-time markets and frequency
regulation,”IEEETrans.Sust.Ener.,vol.7,no.2,pp.875–885,2016.
Fig.7. Timeresponsesofthereference(left)androomtemperatures(right) [8] C.Zhao,U.Topcu,N.LiandS.Low“Designandstabilityofload-side
withthescalingfactor θ=50. primaryfrequencycontrolinpowersystems,”IEEETrans.Automatic
Control, vol.59,no.5,pp.1177–1189, 2014.
[9] E. Mallada, C. Zhao and S. Low, “Optimal load-side control for
noise is filtered out in the present algorithm, and its effects frequencyregulationinsmartgrids,”Proc.52ndAnnualAllertonConf.
do not appear on the estimates. Accordingly,the trajectories Communication, Control, andComputing, pp.731–738,2014.
[10] X. Zhang, A. Papachristodoulou and N. Li, “Distributed optimal
of the estimated optimal and actual room temperatures are
steady-state control using reverse- and forward-engineering,” Proc.
almost the same as Fig. 5. Meanwhile, the disturbance 54thIEEEConf.DecisionandControl,pp.5257–5264,2015.
feedforward approach suffers significant effects from the [11] T. Stegink, C.D. Persis and A. van de Shaft, “A unifying
energy-based approach to optimal frequency and market regula-
noise. The trajectories of the bottom-left get smaller than
tion in power grids,” arXiv:1510.05420v1, 2015 (downloadable at
22◦C.Namely,thefluctuationsarecausedbytheundesirable https://arxiv.org/abs/1510.05420v1).
over-andundershoots.Althoughthetrajectoriesoftheactual [12] X.Zhang,W.Shi,X.Li,B.Yan,A.MalkawiandN.Li,“Decentralized
and distributed temperature control via HVAC systems in energy
temperatures get smooth, this behavior of the reference is
efficient buildings,” Automatica, submitted, 2016.
not desirable from an engineering point of view. Adding [13] A.Cherukuri,E.Mallada,andJ.Corte´s,“Asymptoticconvergence of
a low-pass filter or reducing the gain of the optimization constrained primal-dual dynamics,” System and Control Letters, vol.
87,pp.10–15,2016.
dynamics would eliminate the oscillations but it spoils the
[14] D. Sturzenegger, D. Gyalistras, V. Semeraro, M. Morari and R.S.
advantage, namely response speed. It is to be noted that
Smith, “BRCM Matlab toolbox: Model generation formodel predic-
if the disturbance feedforward is implemented using the tivebuildingcontrol,”Proc.2014AmericanControlConf.,pp.1063–
recoverytechniquein[12],thelow-passfilterisnotdesigned 1069,2014.
[15] T. Hatanaka, N. Chopra, T. Ishizaki and N. Li, “Passivity-based
independently of system stability as stated in Section V.
distributed optimization with communication delays using PI con-
If the response speed of the present algorithm in Fig. sensus estimator,” IEEE Trans. Automatic Control, submitted, 2016
5 is still problematic, it can be accelerated by tuning the (downloadable athttps://arxiv.org/abs/1609.04666).
[16] A.J. van der Schaft, L2-Gain and Passivity Techniques in Nonlinear
scaling factor θ. The results for θ = 50 are shown in Fig.
Control, 2nd edn. Communications and Control Engineering Series.
7, where it is observed that almost the same speed as the Springer, London,2000.
disturbance feedforwardin Fig. 5 is achieved by the present [17] T.Hatanaka, N.Chopra,M.Fujita andM.W.Spong,Passivity-Based
ControlandEstimationinNetworked Robotics,Communications and
algorithm. Simulation for a larger-scale system with more
ControlEngineering Series,Springer-Verlag, 2015.
practical settings is left as a future work of this paper.
[18] Z.QuandM.A.Simaan,“Modularizeddesignforcooperativecontrol
and plug-and-play operation of networked heterogeneous systems,”
VII. CONCLUSION Automatica, vol.50,no.9,pp.2405–2414, 2014.
[19] Y. Hong, J. Hu, and L. Gao, “Tracking control for multi-agent
Inthispaper,wepresentedanovelcombinedoptimization
consensus with an active leader and variable topology,” Automatica,
andcontrolalgorithmforHVACcontrolofbuildings.Wede- vol.42,no.7,pp.1177–1182, 2006.
[20] T.Hatanaka,X.Zhang,W.Shi,M.Zhu,andN.Li,“PhysicsIntegrated Since ν˜=−M2(z˜u+d˜q), it follows
Hierarchical/Distributed HVAC Optimization for Multiple Buildings
withRobustnessagainstTimeDelays,”arXiv,2017(downloadable at D+So ≤−ν˜⊤d˜q−(z˜u+d˜q)⊤M2(z˜u+d˜q)
https://arxiv.org/abs/1510.05420v1). −(zˆ −z∗)⊤(∇f(zˆ )−∇f(z∗))≤−ν˜⊤d˜ . (32)
[21] J.T.Wen,S.MishraS.Mukherjee,N.TantisujjathamandN.Minakais, u u u u q
“Building Temperature Control with Adaptive Feedforward,” Proc.
This completes the proof.
52ndIEEEConf.Decision andControl,pp.4827–4832, 2013.
[22] S. Privara, J. Siroky, L. Ferkl, and J. Cigler, “Model predictive APPENDIX II
control of a building heating system: the first experience,” Energy
andBuildings, vol.43,pp.564–572, 2011. PROOF OF LEMMA4
[23] P.-D. Morosan, R. Bourdais, D. Dumur, and J. Buisson, “Building Lemma 8 Under Assumption 2 and (18), the system
temperature regulation using a distributed model predictive control,” (A¯ ,A P−1/2) is stabilizable and (P−1/2MA⊤,A¯ ) is de-
EnergyandBuildings,vol.42,pp.1445–1452, 2010. 3 2 2 3
[24] S. Yuan and R. Perez, “Multiple-zone ventilation and temperature tectable, where A¯3 :=−A3+A2MP−1A⊤2.
controlofasingle-ductVAVsystemusingmodelpredictivestrategy,”
EnergyandBuildings,vol.38,pp.1248-1261, 2006. Proof: Define Φ := P1/2MP−1A⊤ and Φ :=
s 2 d
[25] M. Maasoumya, M. Razmara, M. Shahbakhti, and A.S. Vincentelli, A MP−1M−1P1/2. Then,
“Handling model uncertainty in model predictive control for energy 2
efficientbuildings,”EnergyandBuildings,vol.77,pp.377–392,2014. A¯ −A P−1/2Φ =−A , A¯ −Φ P−1/2MA⊤ =−A
[26] S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge 3 2 s 3 3 d 2 3
University Press,2004. hold and −A is stable. This completes the proof.
3
[27] USDept.ofEnergy’s:EnergyPlus,https://energyplus.net/ Using Lemma 8, we next prove the following result.
APPENDIX I Lemma 9 Consider the system (21a) with r ≡r∗, w ≡d
a a
PROOF OF LEMMA1 and w ≡ d . Then, under Assumption 2, the system is
q q
Lemma 6 Consider the system (14b) with λˆ(0)≥0. Then, passive from ζ˜to x˜1 :=B⊤x˜ with x˜:=x−x∗.
under Assumption 1, it is passive from z˜ = zˆ − z∗ to
u u u Proof: We first formulate the error system
p˜=p−p∗ with p∗ :=∇g(z∗)λ∗.
u x˜˙ = −A¯x˜+BM−1ζ˜ (33a)
Proof: Define the energy function U := 21kλˆ−λ∗k2. ξ˜˙ = kI(r˜−Bx˜), ζ˜=k¯PM(r˜−B⊤x˜)+Mξ˜ (33b)
Then, following the same procedure as [15], we have
where ξ˜ := ξ − ξ∗. Take a positive definite matrix Ψ ∈
D+U ≤(p−p∗)⊤(zˆu−zu∗)=p˜⊤z˜u, (27) Rn2×n2 and define Ψ¯ := M 0 ∈ Rn×n. Then, by
(cid:20)0 Ψ(cid:21)
where the notation D+ represents the upper Dini derivative. calculation, we have
Integrating this in time completes the proof. P MA⊤+A⊤Ψ
We next consider (14a). Now, replace −M2(zˆu+dˆq)−p Ψ¯A¯+A¯Ψ¯ =(cid:20)ΨA2+A2M ΨA23+A32Ψ (cid:21). (34)
by an external input µ and consider the system
FromSchur complement,underAssumption2, Ψ¯A¯+A¯Ψ¯ >
zˆ˙ =−α{N(w −h¯)+∇f(zˆ )−µ}. (28) 0 is equivalent to the following Riccati inequality.
u a u
−A¯ Ψ−ΨA¯⊤+ΨA P−1A⊤Ψ+A MP−1MA⊤ <0 (35)
Then, we have the following lemma. 3 3 2 2 2 2
A positive semi-definite solutionΨ to (35) is shown to exist
Lemma 7 Supposewa ≡da.Then,underAssumption1,the from Lemma 8. Now, define an energy function Sx :=
system (28) is passive from µ˜ := µ−µ∗ to z˜ = zˆ −z∗, 1x˜⊤Ψ¯x˜ forthe solutionΨ to (35).Then,the timederivative
u u u 2
where µ∗ :=−M2(zu∗+dq)−p∗. of Sx along the trajectories of (33a) is given by
1
Proof: Subtracting (12a) from (28) under yields S˙x = − x˜(Ψ¯A¯+A¯Ψ¯)x˜+x˜⊤Ψ¯BM−1ζ˜
2
z˜˙ =−α(∇f(zˆ )−∇f(z∗))+αµ˜. (29) ≤ x˜⊤Bζ˜=(B⊤x˜)ζ˜=x˜⊤1ζ˜. (36)
u u u
This completes the proof.
Now,defineV := 21αkz˜uk2 = 21αkzˆu−zu∗k2. Then,thetime We are now ready to prove Lemma 4. Replace r˜−B⊤x˜
derivative of V along the trajectories of (29) is given by
in (33b) by e˜as
V˙ =−(zˆu−zu∗)⊤(∇f(zˆu)−∇f(zu∗))+z˜u⊤µ˜. (30) ξ˜˙=kIe˜, ζ˜=k¯PMe˜+Mξ˜. (37)
From convexity of f, (zˆu −zu∗)⊤(∇f(zˆu)−∇f(zu∗)) ≥ 0 Define Sξ := 21kIξ˜⊤Mξ˜. Then, the time derivative of Sξ
holds [26]. This completes the proof. along the trajectories of (37) is given as
Thesystem(14)isgivenbyinterconnecting(14b)and(28) S˙ =ξ˜⊤Me˜=(ζ˜−k¯ Me˜)⊤e˜=ζ˜⊤e˜−k e˜⊤Me˜. (38)
ξ P P
via µ=ν−p. It is then easy to confirm that µ˜=ν˜−p˜.
Define S :=S +S . Then, from (36) and (38), we have
WearenowreadytoproveLemma1.DefineS :=V+U. p x ξ
o
Then, combining (27), (30) and µ˜=ν˜−p˜, we have S˙ ≤ ζ˜⊤r˜−k (r˜−x˜ )⊤M(r˜−x˜ )≤ζ˜⊤r˜. (39)
p P 1 1
D+S ≤−(zˆ −z∗)⊤(∇f(zˆ )−∇f(z∗))+z˜⊤ν˜. (31) This completes the proof.
o u u u u u