- Email: [email protected]
Application of model reduction techniques to building energy simulation
Solar Energy, Vol. 53, No. 3, pp. 289-299, 1994 Copyright 0 1994 Elwier Science Ltd Printed in the USA. All rights reserved 0038-092X/94 $6.00 + .OO
Pergamon 003%092X(
94)00068-9
APPLICATION OF MODEL REDUCTION TECHNIQUES BUILDING ENERGY SIMULATION Lehrstuhl
TO
R. RABENSTEIN * fiir Nachrichtentechnik, University of Erlangen-Niirnberg. Cauerstr. 7, D-91058 Erlangen, Germany
Abstract-The physical description of the energy balance of a building leads to a system of coupled differential or difference equations. It can be written in a unifying form using state space techniques, known from control theory. Depending on the level of detail of the building description and on the chosen physical and mathematical models, the resulting system of equations may be of a high degree. Model reduction techniques attempt to simulate the system performance with a reduced set of equations without deteriorating the quality of the simulation results significantly. This enables a computer simulation with reduced numerical expense. After a short introduction into the theory of model reduction by balanced state space realisation, it will be shown how this method is successfully applied to the energy simulation of buildings.
1. INTRODUCTION
The simulation of the energy balance of a building requires the determination of all relevant energy flows throughout the building. The starting point for a calculation of the quantities of interest (e.g., temperatures and energy demands) is the physical modelling of the various energy conversion and transport mechanisms (e.g., conduction, convection, and absorption). The next step is a numerical modelling of the considered physical effects (e.g., finite difference or finite element methods). The result is a set of coupled differential or difference equations, which have to be solved with a suitable numerical method (e.g., by elimination or iteration). There exists an abundance of methods for physical and numerical modelling, as well as for the solution of the resulting equations. An overview is given by Clarke ( 1985 ). Recently, it has been attempted to represent the equation sets from the different physical and mathematical theories in a unifying neutral model format [ Sahlin and Bring ( 199 1)) Tang ( 199 1)] . Models with spatial discretisation only (e.g., energy balancing at finite control volumes) give a continuous time state space representation of the form i(t)
= Ax(t)
+ Bv(t),
y(t) = Cx(t)
+ Dv(t),
(1)
where t is the continuous time variable. The vector v(t) contains the energy input (radiation, heating), the state vector x(t) the temperatures at the discrete space locations, and y(t) the output quantities. The matrices A, B, C, D describe energy storage and transport. Their elements may be constants or time and temperature dependent variables, according to the physical model. Equation ( 1) may be solved by a numerical integration method. * ISES member.
Models with spatial and temporal discretisation (e.g., response factor methods) lead to a discrete time state space representation x(k+
l)=Ax(k)+Bv(k), y(k) = Cx(k)
+ Dv(k)
(2)
for the diskrete time variable k. Vectors and matrices have the same principal meaning as in eqn ( I ). If the set of equations is linear, then eqn (2) gives an explicit algorithm for a computer simulation, otherwise iterative methods are required. Equations ( 1) and (2) give a unifying description of simple as well as highly detailed building models. Simple models with only few temperature nodes lead to equations of low degree with constant coefficients. Complicated models give equation sets of high degree and with variable coefficients. The choice of the appropriate model and the corresponding set of equations is left to the intuition and experience of the human analyst. However, when considering an integration of building performance prediction into the computerbased architectural design process, methods for the automatic determination of the model order are desirable. There have been proposals of graph theoretical methods [ Buhl et ul. ( 1990)] and-for linear problemsmodal reduction [ Lefebvre et al. ( 1990)], Peuportier and Blanc Sommereux ( 1990)]. Here, a further family of methods is presented, which is based on a certain state space structure: the so-called balanced realisation [ Rabenstein ( 1992)]. After an introduction into the theoretical foundation, the application to building simulation will be demonstrated by various examples. 2.
MODEL
REDUCTION
BY BALANCED STATE
SPACE REALISATION
A mathematical model is a set of equations relating input quantities v(t) to output quantities Y(I). These 289
R.
290
RABHVSTE.II\.
:quations may be algebraic for static systems or differential, or difference, equations for dynamic systems. The latter may be arranged in various forms, one being the state space representation eqns ( 1. 2). For linear jystems, transfer functions are also used to describe Input-output relationships. Model reduction means. generally, to simpiify a given mathematical description In such a way that the output calculated by the reduced model does not differ very much from the original model. Many approaches to model reduction have been proposed. An overview is given by Fortuna et al. ’1992). A couple of requirements to a model reduction technique of practical applicability are obvious: At first. 11should be possible to estimate the necessary math:matical complexity of the reduced equations prior to the actual model reduction. Once the possible degree 3f reduction is known, the reduced equation set should be obtained from the original one in a straightforward lnd numerically efficient manner. Furthermore, the reduced equations should exhibit no inherent or numerical instabilities. A successful method which meets these requirenents is the model reduction based on balanced state ;pace realisations. It is applicable to linear models given n continuous or discrete time formulation eqns ( I. ?) A short description will be given here for continuous ime systems with one input and one output. A more .horough discussion of the theory and the required al;orithms can be found in Fortuna rt cd. ( 1992). Laub 1980), Liu and Anderson (1989). and Moore (1981). We consider a minimal state space description of 1 linear, time-invariant system of the form ( 1) with ,ne input v(t) and one output y(t) and a state vector c(t) with n elements x(t)
+ Bu( 1).
= Ax(t)
y(t) = Cx(t)
+ Dv(t).
x(0
Ti(r)
iX)
with the non-singular (t1 x II) transformation matnx T into eqn (3). It can be easily verified, that all state space representations
= (T-‘AT.
T ‘B, CT, D. T-‘KT-‘,
T’WT)
(9;
exhibit the same input-output behavior. However, the internal representation of the states is different. Thus we have gained additional degrees of freedom in our attempt to group the state variables into irn.. portant ones and redundant ones. We are no longer confined to the elements of the state vector w( 0 from eqn ( 3 ) obtained from the chosen physical description of the problem. When a transformed state representation i(t) gives a clearer distinction of important and redundant states than the original one. then we can use the transformed representation as starting point for a reduction. There remain two questions: How to choose the transformation matrix T and how to decide on the elimination of states. They are answered by the theon of balanced state space realisations [Fortuna (>l trl. (1992), Laub (1980). Liu and Anderson (1989). Moore ( 1981 )]. The key results can be summarized as follows: ?? Calculate the transformation matrix ‘il’from the grammian matrices K and W as
(3)
The matrices A, B, C, D are of the sizes (n X n), (12 < I), ( I X n), ( I X 1). respectively. The main notion s to identify states which are neither easily controllable iom the input nor very well observable at the output md to eliminate these states. A mathematical descrip.ion of controllability and observability of states are .he so-called controllability and observability gramnians K and W [ Callier and Desoer ( 199 I ). Reid 1983)]
s s
We will also need the eigenvalues of the product matrix KW. They are real positive numbers 02 ; z- i n, and their roots g, are called second ordu ~~mh. State space representations are not unique. To illustrate this, we insert the transformed statt,
where L is a lower triangular matrix obtained from K by Cholesky decomposition K = LLT. M’ and I_’ are matrices containing the eigenvalues and eigenvectors of the positive definite matrix L7 WL. i.e., LTWL = UM2U7. ?? The resulting state space representation ( A. B. C. fi) is called a h&n& rc~precPntution. It is characterised by the fact that the transformed grammians k and \icr both are equal to a diagonal matrix consisting of the second order modes
x
K=
0
exp(At)BB7
exp(A7t)dt.
(4)
exp(A’l)C’C
exp(At)dt.
(5)
II
w=
0
They satisfy the Lyapunov
equations
AK + KA’ + BB’ = 0.
(6)
ATW + WA + C7C = 0.
(7)
It can be shown that the second order modes descnbe the contribution of single states to the energy throughput of the whole system. Thus elimination of the ith state with a small value of o, causes only a small change in the output signal. There is an upper bound on the deviation between original and reduced system, which can be calculated in terms of the second order modes of the eliminated states
291
Model reduction techniques Since the squared second order modes u: are eigenvalues of the product KW = T-’ (KW )T, they are invariant under the similarity transformation T and can, therefore, be calculated from the original model in advance. The procedure to obtain a reduced-order model from the original state space description eqn (3) is given by the following steps: 1. Calculate the controllability and observability grammians K and W from eqns (6, 7). 2. Calculate the second order modes 0, from the eigenvalues of the product KW. 3. Inspect the numerical values of the second order modes. Large differences indicate that a balanced representation of (3) contains redundant states with minor influence on the output signal. In this case proceed to the next step. 4. Compute the balanced representation (A, fi, c:, 6) from eqns (9) and ( IO). 5. Permute the states of the balanced system such that the second-order modes are sorted in decreasing magnitude. 6. Partition the balanced system into important (ii) and redundant (&) states.
Y(f) =
I. Obtain
1
%(t) + h(t). [ i,(t)
[C, C:2]
a system of reduced
(12)
order r < n
Du(t),
3. A TUTORIAL
EXAMPLE
3.1 Description cfthe problem As an introductory example we consider a very simple two-node model (Fig. 1) of a homogeneous wall. Its only purpose is to gain some insight into the model reduction technique presented in the previous section. Examples of practical relevance will be treated in the following sections. The temperature distribution in a wall of thickness d = 3A is approximated by two nodes according to Fig. 1. The temperatures at the surfaces are given as v(t) and w(t) = 0. The heat flux y(t) at the surface exposed to the temperature v(t) is to be determined. The relations between the node temperatures .x,(t ) and the interior heat fluxes q,(t) are described by
q,(t) = :(x,_,(t)
k(f) = d%(t) + h(t), y(t) = C%(t) + F
tern can be shown to be stable, provided that the original system is stable. The presentation has been restricted to a continuous time system with a single input and a single output. However, the method is also applicable to systems with multiple inputs and outputs and, with minor modifications, to discrete time systems. The theoretical results compiled so far will be applied in a tutorial manner to a very simple problem in the next section. The following sections show the application of this model reduction method to building simulation by various examples. In every case. the starting point is a standard analysis method leading to a set of equations according to ( 1) or (2). These will be taken as a reference for comparison with reduced order systems.
- x,(t))
i = 1. 2, 3
(15)
i = 1. 2
(16)
(13)
q,(t) - q,+,(t) = ; pcA.Cl(t)
with with x0(t) = u(t) and x,(t) = w(t) = 0 (k thermal conductivity, p density, c specific heat capacity). By eliminating the heat fluxes q,(t) from eqns ( 15, 16), we obtain two coupled first-order differential equations, which can be expressed in matrix notation as
There exist also alternate ways to obtain the reduced order system A, B, C, fi from the partitioned balanced representation ( 12). The form of eqn ( 14) was chosen, because its steady state response is identical to the one of the original system. See Liu and Anderson (1989) for details. This model reduction method meets the requirements previously stated: The degree of the reduced system can be estimated prior to the balancing transformation by inspection of the second order modes. The calculation of the reduced system is straightforward, using standard routines from linear algebra software packages, like matrix diagonalisation, inversion, and Cholesky decomposition. Finally, the reduced sys-
where r = A2.pcfk, y = k/A. The matrices in the first equation can be made dimensionless by intro-
40
-t y(t)
_
xp
q1(t)
i& y
q;(t)
+-----t----t--+-~ 0 A Fig.
I. Two-node
2A
-
w(t)=0
43N t 3A=d
model of a homogeneous
wall.
192
R. RABENSTEIN
lucing the reduced time t + f t/r. For simplicity. no special notation for the reduced time is used. The state space matrices for this problem are now given by
C = $1
01,
D = -y
The matrix A can be expressed ind eigenvectors as
is called the yield, since it quantities the observed sys tern response for t r 0. The next step is to express the cost E,. and the yield I:‘,.by the state matrices A, B. C. D. We start with the cost. The solution of the first equation in ( 3 ) is given by Reid ( 1983) and SchiilJler ( 1991)
(18)
by its eigenvalues with
4 = MAM-‘,
f(t)
A=[-;
_;I.
M=[:
_:].
= exp(At)B
-= M exp(At)M-
‘B.
(19)
The unifying description (3, 18) has reduced the physcal problem given by Fig. 1 to an input-output relaionship as shown in Fig. 2. v(t) is the input to a dynamic system which responds with the output y(t). The inner dynamics are represented by the state vector u(t) = [x,(t) x2(t)] ‘. The specific physical properties Ire described by the state matrices A, B, C, D.
as can be veritied by substitution into ( 3 I. There are infinitely many functions v(t) which can drive the system into a state x0. We want to find the one with minimal cost. Mathematically spoken: We want to minimise E, subject to the constraint
3.2 Stute variabh
We will now carefully evaluate the importance of :ach single state x, and x2 for the input-output be‘laviour of the total system. An inspection of Fig. 1 suggests,that the state x1 should be more important for the quantities at the left surface than the more dis:ant state x2. In order to obtain quantitative results, the following strategy is adopted: For t -c 0 the system is excited by a temperature v(t) with the aim to reach a certain predefined state ~0 = x( t = 0). The exact form of the function v( t ) is yet to be determined. In other words: We want to control the system such that a certain state is reached. At t = 0 the input u(t) is set to zero. We now observe the response y( t ) of the system for t > 0 to an initial state at t = 0. Rather than discussing the time dependence of these ‘unctions in detail, we will characterise their magnitude my scalar quadratic measures, which we call the cost md the yield respectively. The cost for driving the sys.em into the state x0 is specified by
&
=
s 0
H(v. A) =
V’(f Id/
7
where X is a vector of Lagrange-multipliers. The sought function v(t) is obtained by minimizing H( L’.X). that is by setting d/dvH( v. X) = 0. Performing the derivation in (25) and solving for v(t) gives
so
l?(t)&.
3n the other hand
s
- y2(t)dt
0
The Lagrange-multipliers the constraint (24) I x,, = -
are found by insertion
description.
into
^!/
J
f(-l)f’(--f)df*
x.
(27)
I
(20) After a simple change of variables and with (23). the integral is recognized as the controllability grammian K (4). The input signal which reaches the state xfJ at minimal cost follows from (26, 27) as c(t) --f’( -t)K-‘R and the corresponding cost is given by the quadratic form E, = xAKm’xo.
Fig. 2. Input-output
( 16 i
V(I) = $ f’( -f)h.
2
-a
E,=
This is known as an isoperimetrical problem in vat-tational calculus [ Sokolnikoff and Redheffer ( 1966 ) 1. It can be solved by constructing the function
(28)
We now turn to the yield of the observed output .I’(t ) for t > 0. The response to an initial state ~0 with u(t) = 0 is given by y(t) = C exp(At)R, t > 0. Insertion
293
Model reduction techniques into (2 1) gives an expression for the yield which contains an integral that is equal to the observability grammian W (5). Thus the yield is given by Ep = x;WxO.
(29)
We see from (28, 29) that the cost and the yield can be expressed very elegantly by the grammians K and W. They can be calculated either from (4,5) with (23) and elementwise integration or, more straightforward, from the Lyapunov equations (6, 7). For the current example ( 18 ) we obtain
Km’ = 8
W = Y’K
according to (8) can give a better disinction between important and redundant states. A well known example for such a change of the basis is the transformation to the diagonal from with the transformation matrix T = M from ( 19) which leads to a diagonal matrix A = A [see (9 )] The components of the transformed state 8(t) = M-‘x( t) are =
I(%(0
+x2(1))
-f’z(l)
=
f(x,(t)
-x2(t)).
K=_
1
I 1 6
3
48 3
2
W = 4y2K.
f,(t):
E,, = 16,
a,(t): Ev2 = 48,
E,,, = iT2.
Ey2 = ; y2,
I
I
i= 1,2
High values for q, render the state i important, values render it redundant. We obtain in our example x,(t): E,, = 8,
E,, = ;
low
y2,
Thus, x”,is only 9 times more important than Zz. This shows that diagonalisation is not a good approach to model reduction. Instead we consider the balanced representation with T according to ( 10) leading to a state matrix A. The new state vector i(t) is related to the original one by g(t) = T-‘x( t). The numerical values for these matrices are A=
I?;,
=
56,
(35)
(31)
7 41 = 192 Y2 = 3.6. lo-‘y2
x2(t):
y2 = 3.5 - lo-‘y2
T
+fL!& 0,
(34)
For the relation between cost and yield for each state follows
q2 = & qr
(33)
This means that the dynamics of the system are no longer represented by the node temperatures x, (t) directly, but by their sum and difference, respectively. The grammians are obtained from (9) as
(30)
How can we use the cost E, and the yield Ey to assess the importance of the individual states? By setting ~0 to the unit vectors e, = [I 01’ and e2 = [0 llT in (28, 29) we can give quantitative answers to the following questions: What is the minimal cost to set either of the states to 1 while maintaining the other one at O? What is the yield to be expected from these initial states? It is obvious that states which can be excited from the input at little cost and which give a high yield to the output are important for the throughput of the system. On the other hand, states which are hardly controlled from the input, and which cannot be observed well from the output, contribute little to the overall system response and are thus redundant. Both aspects-controllability and observability-can be expressed by a measure for the importance of the state i
f,(f)
EL2 = $ ?‘, 1
92= 192.7
y2=
7.4. lO-4y2
(32)
Obviously, the first state is roughly 50 times more important than the second. 3.3 Transformation of the state space representation We will now turn to the question whether a transformation to a different internal state representation
T-’ =
1’ 1
0.832 1.445
-2.555 0.832
0.957
0.290
-0.290
0.957
(36)
A comparison with (33) shows that again the new states are a weighted average ofthe node temperatures, but the weight factors are different from the diagonal form. The evaluation of the cost and the yield is particularly easy, since the grammians are both diagonal (see 11)
R.
294
RABENSTFJN
and thus
.<, : E,, = L , g,
E,.,=o,,
y,=crf,
i-l.2
(38)
This means that the ratio between yield and cost for the balanced realisation is given by the second order modes which can be calculated as eigenvalues of KW, that is in terms of the original system. With the numerical values from (37), it follows that the state 2, is 370 times more important than T2. Thus the distinction between an important and a redundant state is most pronounced for the balanced realisation. 3.4 Reduced order model The next point to consider is how to obtain a rejuced order system. (71early, we use the balanced realsation i,(r)
= a,,_?-, I) + A,,&(z)
+ &u(t)
k,(t)
= A*,,?, f) + ii&(t)
+ &n(t)
md want to eliminate the second state. However we Jo not dispose of the second equation altogether, but neglect only its dynamics. That is by setting i2 ( t ) = 0. we can solve the right hand side for .qz( t), insert the psult into the first equation and obtain a first order system k(f) = z(t)
+&I(r)
with the scalar coefficients
n=a,, -
‘iiJ&&,
= -1.174
B = & - k,,k;;&
‘1
,!
normolised Fig.
!,(!I z .,l:(w)sin(wZ
t cp(u)i
where the amplitude .-I,,(w) and the phase +CJ(Y Jdepend on the frequency of the input signal a( I). They can be obtained from a transfer function representation [Calher and Desoer ( 199 I ), Kwakemaak and Sivan ( 199 1). Reid ( 1983 ), SchiiMer C199 I )I. Figure 3 shows the difference between the amplitudes A,,(o) for the original second order and the reduced first order system ( y = 1). (Note that the angular frequency w .= ZP/~,, is normal&d in the sense that the period 7” is given in terms of the reduced time.) As expected, the deviation is small for low frequencies. The maximum of the error occurs for high frequencies and is bounded by twice the second order mode of the eliminated state ,
= -0.863
---r---
10-l
Note, that neglecting .Gs( t j does not produce any error for the steady state, where the derivatives are zero anyway. However. what will the error be for dynamic situations? We answer this question for periodic temperature variations of the form c(t) = &sin o! with the amplitude 11, and the angular frequency o. The output will also be periodic of the form
10”
10’
1IO2
frequency
3. Amplitude deviation between second order and first order model for sinusoidal
input.
Here. we consider a homogeneous wall with thickness d, thermal diffusivity a, thermal conductivity k, and interior and exterior surface conductance JI, and h,. respectively. Given are the temperatures at both sides of the wall as time dependent functions _x[(f 1and A>( 0. The quantity to be calculated is the heat flux at the interior surface q,( t ). The physical model for this heat flow problem is obtained by dividing the wall into m t- 1 layers of finite thickness LI = d/m or L / 2 according to Fig. 4. The temperature distribution within each layer is assumed to be homogeneous and characterised by a spatially constant node temperature
x,(t). i = 0.
. 111.
Establishing an energy balance for each layer similar to the previous example gives a system of n = m + t equations according to ( 1) with
295
Model reduction techniques
x(t) = [-~ouo(r)Xl(f) * * - Xm(t)lr, v(t) = [x1(t)
-T4(~)lT,
and the state space matrices
y(t) = 4,(t)
(shown
(39)
for n = 4) 0 0
0.5
- (1 + BiA)
1
1
-2
24
18
time 2
C = [-h,
0
0
0]
D = [h,
0]
0
-J 18
of day
model
1
n=4
--~~----~,__-_~_~-_,
,* \
%
\
;I
(41)
--
l-1 layer
.E
1 12
,
(40)
where
! 6
*’ ‘.#’
i; -1 r"
We will now investigate the following questions: How many equations (n) are required to achieve a reasonable accuracy for the resulting heat flux q,(t)? Can a state space structure different from (40) give a higher accuracy? To this end, the layer model with n = 25 equations is used as a reference. At first, we inspect the second order modes for the numerical values: k = 0.22 Wm-’ K~‘,a=0.55~10-6m2s-‘,d=0.15m,h,=8Wm~2 K-’ , hA = 20 WmL2 K-’ . Calculation of the controllability and observability grammians and determination of the eigenvalues of their product gives (sorted by descending magnitude):
18
time
24
6
time
12
=
I
I
I
6
12
time
of day
18
I 18
Fig. 5. Temperature profile of exterior temperature and calculated heat fluxes from reference and low order model.
. . -)
The dominance of the first few second order modes suggests that only a handful of equations should be sufficient. We will investigate two ways of obtaining a low order model: A layer model of the form (40) with few temperature nodes and a balanced realisation obtained from the reference mode1 according to ( 13, 14). For n = 3 the state matrices AI and A, for the layer model and the reduced model are given by: -0.7289 0.0978 0
0.1956 -0.1956 0.1956
0 0.0978 -1.5289
of day
24
18
2.3, 0.92, 0.42, 0.16, 0.048, 0.013, 0.003,
A, =
18
12
of day
Fig. 6. Deviation between the heat fluxes of the low order systems and the reference system for n = 4 and n = 5.
A, = 18
6
24
-0.2432 0.2887 -0.1655
0.2887 -0.4186 0.2585
1
- 10--3,
-0.1952 0.3489 -0.286
1
- 10-2.
Both models were used to simulate the time response of the heat flux q,(t) in response to a constant interior temperature x,(t) = 20°C and to a variable exterior temperature xA (t) with night-day profile as given on the top of Fig. 5. The resulting heat fluxes are shown below together with the response of the high order (n = 25) reference system. Fig. 6 (top) displays the deviation of the responses of both low order systems from the reference. The layer model has a peak deviation of 1.5 Wm-2 while the deviation of the reduced mode1 of the same order stays below 0.5 Wrn-‘. The difference in accuracy is even more pronounced for n
296
R.
RABENSTHN
= 4,5. Figure 6 (bottom) shows the deviations of both models for n = 4. Compared to n = 3, the deviation of the layer model is about 50% while the deviation of the reduced model is hardly visible in this scale. The peak deviations from the reference model are compiled in Fig. 7 for n = 3, 4, 5. They are normalized to the steady state heat flux for a temperature difference of 10 K. The equation set of the reduced model obviously leads to a more accurate simulation than a layer model with the same number of equations.
5.
RESPONSE-FACTOR MULTI-LAYER
METHOD
FOR
WALLS
The response-factor method [Stephenson and Mitalas ( 197 I)] is a widely-used, discrete time model for multi-layer walls. It is based on a transfer function description of the heat conduction process and represents a wall by a difference equation of the form (k is the discrete time variable) y(k) = -i
,I=,
c,y(k - v) + i h,v(k - u). ,a= 0
(42)
where the coefficients h,, c, depend in a complicated manner on the material constants of the different layers and surfaces and on the time step between two samples (usually 1 hour). We will treat two examples: a twolayer wall from Stephenson and Mitalas ( 197 1) Appendix B, 1st method, 1/B; and a four-layer wall according to Table 1. The coefficients for the two-layer wall are given in Stephenson and Mitalas ( 197 1 ), where a difference equation of degree n = 5 is used. The coefficients for the four-layer wall are listed in Table 2 with n = 6. In order to apply the model reduction technique described above, a state space description is required, which corresponds to the the difference eqn (42). Its matrices A, B, C, D are easily constructed from the coefficients b,, c, [see Kwakernaak and Sivan ( 199 1 ), Schtirjler ( 199 1)]. The reduced model is obtained according to the steps listed in Section 2. Finally the reduced state space description can be converted to a difference equation of the form (42), however with
15% i
Table I. Material
Thermal conductivity
Thickness ~__
constants
(Wm. ’EC-‘ ) _________
Cm)
for the four-layer
Heat capacity per volume -3 K ,) (kJm ,
wail --Surface
conductance (Wm 2 K ‘1 ~~~~_~__~,~~_~_~ i
0.0
15
0.175
11.94 0.035 I.1
0.100 0.1 IS _._~_._
l’O0 1800
0.35
25 7000
_-.-. ~_~._._.-.-.-- -____.-...- -_--
1, -..... ~______.
fewer coefficients. Table 3 lists the coefficients of a hrstorder difference equation for the simulation of the twolayer wall and of a second-order difference equation for the four-layer wall. Figure 8 shows simulation results for a problem similar to the one treated in the previous section. The temperature at the interior side of the wall is held constant at 20°C while the exterior side is exposed to a variable temperature as shown in Fig. 8 (top). The resulting heat fluxes are displayed in Fig. 8 (middle and bottom) for the two-layer wall and the four-layer wall, respectively. The difference equations from the response factor models of order tr = 5 and 11 = 6 are used as reference (solid lines), against which the difference equations from the reduced models are corn-pared (dotted lines). The relative deviations between reference and reduced model are below 3%: for the firstorder model of the two-layer wall and below 0.5% for the second-order model of the four-layer wall (deviations normalised to the maximum heat flux of each data set )
6.
I‘HERMAL
NETWORK
MODEL
The preceding sections covered examples of single walls, where the relation between one input quantity (exterior temperature) and one output quantity (heat flux) was considered. This example shows that model reduction techniques can also be applied to more involved models. We will treat a whole building model (though a relatively simple one), with three input quantities (exterior and ground temperature and solar radiation ) and two output quantities (room air tern.. peratures). The energy transport phenomena are not confined to walls and surfaces only. but include absorption and convective heat transport.
??
I \
layer
model Table 2. coefficients for reference of the four-layer wall
10%
modei
-_____
5% I
i;i
l%L 3
Fig. 7. Relative deviation
L
5
from the reference n = 25.
n
model of order
L
!I,
0 1
-0.0000989 0.0004437 -0.0006997 0.0012087 0.0000965 0.0002 138 --0.0000027
2 3 4
5 6
2.47225’) 2.I43.508 -0.766410 0.10393 1 ~0.0049 I h 0.000034
Model reduction techniques
297
Table 3. Coefficients of the reduced models Four-layer wall
Two-layer wall ”
b,
0
-0.347398 0.679234
1 2
C” -0.860142
Figure 9 shows a fictitious building with three zones. Floor, ceiling, and walls are modelled as storage walls, the south-oriented glass facade AF and the roof as massless walls. The capacities C,, i = 2, 3, 4, 6, 7, 8 characterize the heat capacity of the storage walls, the R, , RA , RBZ, RB6 describe heat conductance from the middle of the walls to the room air nodes or to the of the roof and exterior. RD is the heat conductance RF characterises the combined transmission and ventilation loss of the glass facade. RL describes the energy flow between the room air nodes 1 and 5 due to ventilation. The capacities C, , Cs, C9 are the heat capacity of the air volume in each room. The building is exposed to the exterior and the ground temperatures T, and Ts and to the solar radiation G. A simulation of the building energy balance shall predict the room air temperatures at the nodes 1 and 5 from the given meteorological quantities. Setting up the energy balance for the fluxes at each node gives a state space description of the form of ( I), which is equivalent to an ordinary differential equation of degree n = 9. The nonzero entries into the first line of A and B are
, 1 A,, = i = 2, 3, 4, RiCl ’ B,, =-
1 &Cl
’
B,3 = $.
b,
C”
0.00495 14 -0.0 14695 1 0.0121846
-1.791690 0.799863
7. IMPLEMENTATION
IN A SIMULATION
ENVIRONMENT
The model reduction by balanced state space realisation is primarily restricted to linear models. However, this does not prevent its application in the simulation of complex building models which also incorporate nonlinear effects. Simulation programs with a modular structure permit the replacement of linear submodels by reduced models. For instance, consider the thermal-network example in the previous section. A more accurate model could be obtained by replacing the rather crude RC-approximation of the walls by a reduced model derived from a higher-order responsefactor description of multi-layer walls. At the same time, radiative and convective energy transport could be modelled by nonlinear equations. The result would be a more detailed building description not neglecting nonlinear effects with reduced models for linear subsystems. In order to show the applicability of model reduction methods in complex building models, an early version of the Energy Kernel System [ Charlesworth et al. ( 199 1), Rabenstein ( 1993)] was chosen as a simulation environment. The Energy Kernel System (EKS) is not a single simulation program, rather it allows the user to build a custom-made run-time code
1 AIS = RLC, ’
I
The state variables xi(t) i = 1, . . . , 9 are the temperatures at the nodes. The input and output quantities are
G(Ol’,
v(t) = ITA(
T,(r),
Y(l) = [x,(t)
xs(t)lT.
_J
)-lo1 1
10
15 time in
20
25
30
15 20 in days
25
30
days
.‘~~~ I 1
This state space model of order n = 9 was used as a starting point for a model reduction to a reduced system of order n = 3. Simulation runs for both systems with typical values for Ri and Ci , hourly mean values for the exterior temperature TA and the global radiation G, and constant ground temperature Ts gave two predictions for the time history of the node temperatures 1 and 5. The output of the reduced model did not deviate more than 0.3 K from the reference model at either node.
5
5
10 time
time in doys
Fig. 8. Ambient temperature (top), and heat fluxes for the two-layer wall (middle) and the four-layer wall (bottom). Solid lines: reference model: dotted lines: reduced model.
R. RABENS I EIN
298
F
A
TA
TA
TB
TB
B
Fig. 9. Thermal
B
network model of a building.
From a set of existing program modules. The choice of modules contains different physical models of varying level of detail and various numerical models and solvers For the resulting equation systems. Furthermore, its modular structure, based on an object-oriented approach, facilitates the addition of new program modules without alteration of existing program code. A finite difference model for the calculation of heat zonduction through multi-layer walls, already existing n the EKS, was chosen as a reference. It produces a Cscrete space, continuous time model, similar to the me presented in Section 3. However, for arbitrary n&i-layer walls. A new program module was derived using the principle of inheritance in the object-ori:nted sense), which uses the state space formulation If the finite difference model and replaces it by a balmced reduced model. By choosing either the existing inite difference model or, alternatively, the derived .educed model, it is possible to construct simulation )rograms, which are identical except for the modelling If heat conduction through walls. Both program versions were tested for the layernode1 and the response factor method problem dis:ussed earlier. A variety of multi-layer walls has been nvestigated. In each case a reference model was con;tructed using the finite difference method with 20 lodes distributed over the different layers of the wall. rhen lower-order models were found by either reducng the number of nodes in the finite difference model )r by model reduction. Simulation runs with identical nput temperature data showed the trade-off between Iecreased model order and deviation from the refer:nce model. For a quantitative comparison, a number of cases vere identified, where a low order finite difference node1 and a reduced model gave approximately the ,ame deviation from the reference (typically 0.1 Yo- I YG) t turned out that the required order of the finite diftrence models was in the range of 9-13, while the
reduced order models achieved the same or higher accuracy with a degree of only 4-5. Furthermore, simulation runs for finite difterence and reduced models with comparable accuracy were timed on a workstation (SUN SPARCstation WC). lnput data sets with l-year duration of hourly mean values were used. The run time for the reduced models was in the range of 19%36Yo of the run time of the low-order finite difference models.
X. (‘ONCLI
.SlOh
Theoretical considerations and practical implementations have shown that the model reduction technique by balanced state space representation is a useful tool for building energy simulation. It allows to perform simulations with the same accuracy as standard analysis methods. however with a significantly decreased model order. This method is applicable to linear building models and to linear submodels of general nonlineal building descriptions. The calculation of the reduced order state space description from a higher order reEerence is a straightforward numerical procedure. requiring matrix calculus and standard linear algebra routines. Some higher level languages (e.g. MAT1,.4B 1 provide powerful commands for model reduction. .‘1~,~~oll,/~,~~~~~~~~~~\Ihe author
wishes to thank: M. Gerwien, Motz. and .A. Wegehaupt for the programming of the esamples; D. Heinemann for supplying the meteorological data: and J. Clarke. D. Tang, and D. Mac Randal for an introduction to the Energy Kernel System.
r.
RkXERENCES F. Buhl. E. Erdem,
J. -M. Nataf, F. C. Winkelmann. M. A Moshier. and E. F. Sowell. 7%~~US. EKS, Proc. of/he 3nd
Inl. Confrrenw on S~.wm Sirnulution in BuildinKs, LGge. pp. 107-149 (1990). F. M. Callier. C. A. Desoer. Lvzmr s~:trern rhcw,~. SprmgerVerlag, New York ( 199 I ).
Model reduction techniques J. A. Clarke, Energy simulation in building design, Adam Hilger Ltd, Bristol ( 1985). P. Charlesworth, J. A. Clarke, G. Hammond, A. Irving, K. James, B. Lee, S. Lockley, D. Mac Randal, D. Tang, T. J. Wiltshire, A. J. Wright, The energy kernel system, Proc. of Building Simulation ‘91, Nice, pp. 313-322 (1991). L. Fortuna, G. Nunnari, A. Gallo, Model order reduction techniques with applications in electrical engineering, Springer-Verlag, London ( 1992). H. Kwakernaak, R. Sivan, Modern signals and systems, Prentice Hall, Englewood Cliffs, NJ ( 199 I ). A. Laub, Computation of “‘balancing”transformations, Proc. Joint American Control Conference, I, paper FA8-E (1980). G. Lefebvre, A. Neveu, K. el Khoury, J. J. Saigon, Applying the modal method to thermal modelling, Proc. of the 9th Int. Heat Transfer Conference, pp. 157- 16I, 3, Jerusalem (1990). Y. Liu, B. Anderson, Singular perturbation approximation of balanced systems, Int. I. Control 50(4), 1379-1405 (1989). B. C. Moore, Principal Component Analysis in Linear Systems, IEEE Trans. on Automatic Control, Vol. AC-26( I ), 1732(1981).
299
B. Peuportier, 1. Blanc Sommereux, Simulation tool with its expert interface for the thermal design of multizone buildings, Int. I. Solar Energy 8, 109-120 (1990). R. Rabenstein, Modellreduktion als Werkzeug zur ejektiven Simulation des Gebaudewarmehaushalts. Proc. of 8. Internationales Sonnenforum, Berlin, pp. 152-157 ( 1992). R. Rabenstein: Erste Erfahrungen mit dem Energy Kernel System zur energetischen Gebaudesimulation. Proc. of 3. Nat. Symposium Thermische Solarenergie, Regensburg, Germany, pp. 259-263 ( 1993). J. G. Reid, Linear system,findamentals, McGraw-Hill, New York (1983). P. Sahlin and A. Bring, IDA Solver, Conference Proceedings Building Simulation ‘9 1,Nice, pp. 339-348 ( 199 1). H. W. SchilBler, Netzwerke, Signale und S.ysteme 2, 3. ed., Springer-Verlag, Berlin ( 1991). I. S. Sokolnikoffand R. M. Redheffer, Mathematics ofphysics and modern engineering, 2nd ed., M&raw-Hill, New York (1966). D. G. Stephenson and G. P. Mitalas, Caluculation of heat conduction transfer,finctions for multi-layer slabs, ASHRAE Transactions, 77, Part II, pp. I 17-126 ( I97 I ). D. Tang, The generalised system solution classes m the EKS environment, Conference Proceedings Building Simulation ‘9 I, Nice, pp. 323-327 ( I99 I ).
Pergamon 003%092X(
94)00068-9
APPLICATION OF MODEL REDUCTION TECHNIQUES BUILDING ENERGY SIMULATION Lehrstuhl
TO
R. RABENSTEIN * fiir Nachrichtentechnik, University of Erlangen-Niirnberg. Cauerstr. 7, D-91058 Erlangen, Germany
Abstract-The physical description of the energy balance of a building leads to a system of coupled differential or difference equations. It can be written in a unifying form using state space techniques, known from control theory. Depending on the level of detail of the building description and on the chosen physical and mathematical models, the resulting system of equations may be of a high degree. Model reduction techniques attempt to simulate the system performance with a reduced set of equations without deteriorating the quality of the simulation results significantly. This enables a computer simulation with reduced numerical expense. After a short introduction into the theory of model reduction by balanced state space realisation, it will be shown how this method is successfully applied to the energy simulation of buildings.
1. INTRODUCTION
The simulation of the energy balance of a building requires the determination of all relevant energy flows throughout the building. The starting point for a calculation of the quantities of interest (e.g., temperatures and energy demands) is the physical modelling of the various energy conversion and transport mechanisms (e.g., conduction, convection, and absorption). The next step is a numerical modelling of the considered physical effects (e.g., finite difference or finite element methods). The result is a set of coupled differential or difference equations, which have to be solved with a suitable numerical method (e.g., by elimination or iteration). There exists an abundance of methods for physical and numerical modelling, as well as for the solution of the resulting equations. An overview is given by Clarke ( 1985 ). Recently, it has been attempted to represent the equation sets from the different physical and mathematical theories in a unifying neutral model format [ Sahlin and Bring ( 199 1)) Tang ( 199 1)] . Models with spatial discretisation only (e.g., energy balancing at finite control volumes) give a continuous time state space representation of the form i(t)
= Ax(t)
+ Bv(t),
y(t) = Cx(t)
+ Dv(t),
(1)
where t is the continuous time variable. The vector v(t) contains the energy input (radiation, heating), the state vector x(t) the temperatures at the discrete space locations, and y(t) the output quantities. The matrices A, B, C, D describe energy storage and transport. Their elements may be constants or time and temperature dependent variables, according to the physical model. Equation ( 1) may be solved by a numerical integration method. * ISES member.
Models with spatial and temporal discretisation (e.g., response factor methods) lead to a discrete time state space representation x(k+
l)=Ax(k)+Bv(k), y(k) = Cx(k)
+ Dv(k)
(2)
for the diskrete time variable k. Vectors and matrices have the same principal meaning as in eqn ( I ). If the set of equations is linear, then eqn (2) gives an explicit algorithm for a computer simulation, otherwise iterative methods are required. Equations ( 1) and (2) give a unifying description of simple as well as highly detailed building models. Simple models with only few temperature nodes lead to equations of low degree with constant coefficients. Complicated models give equation sets of high degree and with variable coefficients. The choice of the appropriate model and the corresponding set of equations is left to the intuition and experience of the human analyst. However, when considering an integration of building performance prediction into the computerbased architectural design process, methods for the automatic determination of the model order are desirable. There have been proposals of graph theoretical methods [ Buhl et ul. ( 1990)] and-for linear problemsmodal reduction [ Lefebvre et al. ( 1990)], Peuportier and Blanc Sommereux ( 1990)]. Here, a further family of methods is presented, which is based on a certain state space structure: the so-called balanced realisation [ Rabenstein ( 1992)]. After an introduction into the theoretical foundation, the application to building simulation will be demonstrated by various examples. 2.
MODEL
REDUCTION
BY BALANCED STATE
SPACE REALISATION
A mathematical model is a set of equations relating input quantities v(t) to output quantities Y(I). These 289
R.
290
RABHVSTE.II\.
:quations may be algebraic for static systems or differential, or difference, equations for dynamic systems. The latter may be arranged in various forms, one being the state space representation eqns ( 1. 2). For linear jystems, transfer functions are also used to describe Input-output relationships. Model reduction means. generally, to simpiify a given mathematical description In such a way that the output calculated by the reduced model does not differ very much from the original model. Many approaches to model reduction have been proposed. An overview is given by Fortuna et al. ’1992). A couple of requirements to a model reduction technique of practical applicability are obvious: At first. 11should be possible to estimate the necessary math:matical complexity of the reduced equations prior to the actual model reduction. Once the possible degree 3f reduction is known, the reduced equation set should be obtained from the original one in a straightforward lnd numerically efficient manner. Furthermore, the reduced equations should exhibit no inherent or numerical instabilities. A successful method which meets these requirenents is the model reduction based on balanced state ;pace realisations. It is applicable to linear models given n continuous or discrete time formulation eqns ( I. ?) A short description will be given here for continuous ime systems with one input and one output. A more .horough discussion of the theory and the required al;orithms can be found in Fortuna rt cd. ( 1992). Laub 1980), Liu and Anderson (1989). and Moore (1981). We consider a minimal state space description of 1 linear, time-invariant system of the form ( 1) with ,ne input v(t) and one output y(t) and a state vector c(t) with n elements x(t)
+ Bu( 1).
= Ax(t)
y(t) = Cx(t)
+ Dv(t).
x(0
Ti(r)
iX)
with the non-singular (t1 x II) transformation matnx T into eqn (3). It can be easily verified, that all state space representations
= (T-‘AT.
T ‘B, CT, D. T-‘KT-‘,
T’WT)
(9;
exhibit the same input-output behavior. However, the internal representation of the states is different. Thus we have gained additional degrees of freedom in our attempt to group the state variables into irn.. portant ones and redundant ones. We are no longer confined to the elements of the state vector w( 0 from eqn ( 3 ) obtained from the chosen physical description of the problem. When a transformed state representation i(t) gives a clearer distinction of important and redundant states than the original one. then we can use the transformed representation as starting point for a reduction. There remain two questions: How to choose the transformation matrix T and how to decide on the elimination of states. They are answered by the theon of balanced state space realisations [Fortuna (>l trl. (1992), Laub (1980). Liu and Anderson (1989). Moore ( 1981 )]. The key results can be summarized as follows: ?? Calculate the transformation matrix ‘il’from the grammian matrices K and W as
(3)
The matrices A, B, C, D are of the sizes (n X n), (12 < I), ( I X n), ( I X 1). respectively. The main notion s to identify states which are neither easily controllable iom the input nor very well observable at the output md to eliminate these states. A mathematical descrip.ion of controllability and observability of states are .he so-called controllability and observability gramnians K and W [ Callier and Desoer ( 199 I ). Reid 1983)]
s s
We will also need the eigenvalues of the product matrix KW. They are real positive numbers 02 ; z- i n, and their roots g, are called second ordu ~~mh. State space representations are not unique. To illustrate this, we insert the transformed statt,
where L is a lower triangular matrix obtained from K by Cholesky decomposition K = LLT. M’ and I_’ are matrices containing the eigenvalues and eigenvectors of the positive definite matrix L7 WL. i.e., LTWL = UM2U7. ?? The resulting state space representation ( A. B. C. fi) is called a h&n& rc~precPntution. It is characterised by the fact that the transformed grammians k and \icr both are equal to a diagonal matrix consisting of the second order modes
x
K=
0
exp(At)BB7
exp(A7t)dt.
(4)
exp(A’l)C’C
exp(At)dt.
(5)
II
w=
0
They satisfy the Lyapunov
equations
AK + KA’ + BB’ = 0.
(6)
ATW + WA + C7C = 0.
(7)
It can be shown that the second order modes descnbe the contribution of single states to the energy throughput of the whole system. Thus elimination of the ith state with a small value of o, causes only a small change in the output signal. There is an upper bound on the deviation between original and reduced system, which can be calculated in terms of the second order modes of the eliminated states
291
Model reduction techniques Since the squared second order modes u: are eigenvalues of the product KW = T-’ (KW )T, they are invariant under the similarity transformation T and can, therefore, be calculated from the original model in advance. The procedure to obtain a reduced-order model from the original state space description eqn (3) is given by the following steps: 1. Calculate the controllability and observability grammians K and W from eqns (6, 7). 2. Calculate the second order modes 0, from the eigenvalues of the product KW. 3. Inspect the numerical values of the second order modes. Large differences indicate that a balanced representation of (3) contains redundant states with minor influence on the output signal. In this case proceed to the next step. 4. Compute the balanced representation (A, fi, c:, 6) from eqns (9) and ( IO). 5. Permute the states of the balanced system such that the second-order modes are sorted in decreasing magnitude. 6. Partition the balanced system into important (ii) and redundant (&) states.
Y(f) =
I. Obtain
1
%(t) + h(t). [ i,(t)
[C, C:2]
a system of reduced
(12)
order r < n
Du(t),
3. A TUTORIAL
EXAMPLE
3.1 Description cfthe problem As an introductory example we consider a very simple two-node model (Fig. 1) of a homogeneous wall. Its only purpose is to gain some insight into the model reduction technique presented in the previous section. Examples of practical relevance will be treated in the following sections. The temperature distribution in a wall of thickness d = 3A is approximated by two nodes according to Fig. 1. The temperatures at the surfaces are given as v(t) and w(t) = 0. The heat flux y(t) at the surface exposed to the temperature v(t) is to be determined. The relations between the node temperatures .x,(t ) and the interior heat fluxes q,(t) are described by
q,(t) = :(x,_,(t)
k(f) = d%(t) + h(t), y(t) = C%(t) + F
tern can be shown to be stable, provided that the original system is stable. The presentation has been restricted to a continuous time system with a single input and a single output. However, the method is also applicable to systems with multiple inputs and outputs and, with minor modifications, to discrete time systems. The theoretical results compiled so far will be applied in a tutorial manner to a very simple problem in the next section. The following sections show the application of this model reduction method to building simulation by various examples. In every case. the starting point is a standard analysis method leading to a set of equations according to ( 1) or (2). These will be taken as a reference for comparison with reduced order systems.
- x,(t))
i = 1. 2, 3
(15)
i = 1. 2
(16)
(13)
q,(t) - q,+,(t) = ; pcA.Cl(t)
with with x0(t) = u(t) and x,(t) = w(t) = 0 (k thermal conductivity, p density, c specific heat capacity). By eliminating the heat fluxes q,(t) from eqns ( 15, 16), we obtain two coupled first-order differential equations, which can be expressed in matrix notation as
There exist also alternate ways to obtain the reduced order system A, B, C, fi from the partitioned balanced representation ( 12). The form of eqn ( 14) was chosen, because its steady state response is identical to the one of the original system. See Liu and Anderson (1989) for details. This model reduction method meets the requirements previously stated: The degree of the reduced system can be estimated prior to the balancing transformation by inspection of the second order modes. The calculation of the reduced system is straightforward, using standard routines from linear algebra software packages, like matrix diagonalisation, inversion, and Cholesky decomposition. Finally, the reduced sys-
where r = A2.pcfk, y = k/A. The matrices in the first equation can be made dimensionless by intro-
40
-t y(t)
_
xp
q1(t)
i& y
q;(t)
+-----t----t--+-~ 0 A Fig.
I. Two-node
2A
-
w(t)=0
43N t 3A=d
model of a homogeneous
wall.
192
R. RABENSTEIN
lucing the reduced time t + f t/r. For simplicity. no special notation for the reduced time is used. The state space matrices for this problem are now given by
C = $1
01,
D = -y
The matrix A can be expressed ind eigenvectors as
is called the yield, since it quantities the observed sys tern response for t r 0. The next step is to express the cost E,. and the yield I:‘,.by the state matrices A, B. C. D. We start with the cost. The solution of the first equation in ( 3 ) is given by Reid ( 1983) and SchiilJler ( 1991)
(18)
by its eigenvalues with
4 = MAM-‘,
f(t)
A=[-;
_;I.
M=[:
_:].
= exp(At)B
-= M exp(At)M-
‘B.
(19)
The unifying description (3, 18) has reduced the physcal problem given by Fig. 1 to an input-output relaionship as shown in Fig. 2. v(t) is the input to a dynamic system which responds with the output y(t). The inner dynamics are represented by the state vector u(t) = [x,(t) x2(t)] ‘. The specific physical properties Ire described by the state matrices A, B, C, D.
as can be veritied by substitution into ( 3 I. There are infinitely many functions v(t) which can drive the system into a state x0. We want to find the one with minimal cost. Mathematically spoken: We want to minimise E, subject to the constraint
3.2 Stute variabh
We will now carefully evaluate the importance of :ach single state x, and x2 for the input-output be‘laviour of the total system. An inspection of Fig. 1 suggests,that the state x1 should be more important for the quantities at the left surface than the more dis:ant state x2. In order to obtain quantitative results, the following strategy is adopted: For t -c 0 the system is excited by a temperature v(t) with the aim to reach a certain predefined state ~0 = x( t = 0). The exact form of the function v( t ) is yet to be determined. In other words: We want to control the system such that a certain state is reached. At t = 0 the input u(t) is set to zero. We now observe the response y( t ) of the system for t > 0 to an initial state at t = 0. Rather than discussing the time dependence of these ‘unctions in detail, we will characterise their magnitude my scalar quadratic measures, which we call the cost md the yield respectively. The cost for driving the sys.em into the state x0 is specified by
&
=
s 0
H(v. A) =
V’(f Id/
7
where X is a vector of Lagrange-multipliers. The sought function v(t) is obtained by minimizing H( L’.X). that is by setting d/dvH( v. X) = 0. Performing the derivation in (25) and solving for v(t) gives
so
l?(t)&.
3n the other hand
s
- y2(t)dt
0
The Lagrange-multipliers the constraint (24) I x,, = -
are found by insertion
description.
into
^!/
J
f(-l)f’(--f)df*
x.
(27)
I
(20) After a simple change of variables and with (23). the integral is recognized as the controllability grammian K (4). The input signal which reaches the state xfJ at minimal cost follows from (26, 27) as c(t) --f’( -t)K-‘R and the corresponding cost is given by the quadratic form E, = xAKm’xo.
Fig. 2. Input-output
( 16 i
V(I) = $ f’( -f)h.
2
-a
E,=
This is known as an isoperimetrical problem in vat-tational calculus [ Sokolnikoff and Redheffer ( 1966 ) 1. It can be solved by constructing the function
(28)
We now turn to the yield of the observed output .I’(t ) for t > 0. The response to an initial state ~0 with u(t) = 0 is given by y(t) = C exp(At)R, t > 0. Insertion
293
Model reduction techniques into (2 1) gives an expression for the yield which contains an integral that is equal to the observability grammian W (5). Thus the yield is given by Ep = x;WxO.
(29)
We see from (28, 29) that the cost and the yield can be expressed very elegantly by the grammians K and W. They can be calculated either from (4,5) with (23) and elementwise integration or, more straightforward, from the Lyapunov equations (6, 7). For the current example ( 18 ) we obtain
Km’ = 8
W = Y’K
according to (8) can give a better disinction between important and redundant states. A well known example for such a change of the basis is the transformation to the diagonal from with the transformation matrix T = M from ( 19) which leads to a diagonal matrix A = A [see (9 )] The components of the transformed state 8(t) = M-‘x( t) are =
I(%(0
+x2(1))
-f’z(l)
=
f(x,(t)
-x2(t)).
K=_
1
I 1 6
3
48 3
2
W = 4y2K.
f,(t):
E,, = 16,
a,(t): Ev2 = 48,
E,,, = iT2.
Ey2 = ; y2,
I
I
i= 1,2
High values for q, render the state i important, values render it redundant. We obtain in our example x,(t): E,, = 8,
E,, = ;
low
y2,
Thus, x”,is only 9 times more important than Zz. This shows that diagonalisation is not a good approach to model reduction. Instead we consider the balanced representation with T according to ( 10) leading to a state matrix A. The new state vector i(t) is related to the original one by g(t) = T-‘x( t). The numerical values for these matrices are A=
I?;,
=
56,
(35)
(31)
7 41 = 192 Y2 = 3.6. lo-‘y2
x2(t):
y2 = 3.5 - lo-‘y2
T
+fL!& 0,
(34)
For the relation between cost and yield for each state follows
q2 = & qr
(33)
This means that the dynamics of the system are no longer represented by the node temperatures x, (t) directly, but by their sum and difference, respectively. The grammians are obtained from (9) as
(30)
How can we use the cost E, and the yield Ey to assess the importance of the individual states? By setting ~0 to the unit vectors e, = [I 01’ and e2 = [0 llT in (28, 29) we can give quantitative answers to the following questions: What is the minimal cost to set either of the states to 1 while maintaining the other one at O? What is the yield to be expected from these initial states? It is obvious that states which can be excited from the input at little cost and which give a high yield to the output are important for the throughput of the system. On the other hand, states which are hardly controlled from the input, and which cannot be observed well from the output, contribute little to the overall system response and are thus redundant. Both aspects-controllability and observability-can be expressed by a measure for the importance of the state i
f,(f)
EL2 = $ ?‘, 1
92= 192.7
y2=
7.4. lO-4y2
(32)
Obviously, the first state is roughly 50 times more important than the second. 3.3 Transformation of the state space representation We will now turn to the question whether a transformation to a different internal state representation
T-’ =
1’ 1
0.832 1.445
-2.555 0.832
0.957
0.290
-0.290
0.957
(36)
A comparison with (33) shows that again the new states are a weighted average ofthe node temperatures, but the weight factors are different from the diagonal form. The evaluation of the cost and the yield is particularly easy, since the grammians are both diagonal (see 11)
R.
294
RABENSTFJN
and thus
.<, : E,, = L , g,
E,.,=o,,
y,=crf,
i-l.2
(38)
This means that the ratio between yield and cost for the balanced realisation is given by the second order modes which can be calculated as eigenvalues of KW, that is in terms of the original system. With the numerical values from (37), it follows that the state 2, is 370 times more important than T2. Thus the distinction between an important and a redundant state is most pronounced for the balanced realisation. 3.4 Reduced order model The next point to consider is how to obtain a rejuced order system. (71early, we use the balanced realsation i,(r)
= a,,_?-, I) + A,,&(z)
+ &u(t)
k,(t)
= A*,,?, f) + ii&(t)
+ &n(t)
md want to eliminate the second state. However we Jo not dispose of the second equation altogether, but neglect only its dynamics. That is by setting i2 ( t ) = 0. we can solve the right hand side for .qz( t), insert the psult into the first equation and obtain a first order system k(f) = z(t)
+&I(r)
with the scalar coefficients
n=a,, -
‘iiJ&&,
= -1.174
B = & - k,,k;;&
‘1
,!
normolised Fig.
!,(!I z .,l:(w)sin(wZ
t cp(u)i
where the amplitude .-I,,(w) and the phase +CJ(Y Jdepend on the frequency of the input signal a( I). They can be obtained from a transfer function representation [Calher and Desoer ( 199 I ), Kwakemaak and Sivan ( 199 1). Reid ( 1983 ), SchiiMer C199 I )I. Figure 3 shows the difference between the amplitudes A,,(o) for the original second order and the reduced first order system ( y = 1). (Note that the angular frequency w .= ZP/~,, is normal&d in the sense that the period 7” is given in terms of the reduced time.) As expected, the deviation is small for low frequencies. The maximum of the error occurs for high frequencies and is bounded by twice the second order mode of the eliminated state ,
= -0.863
---r---
10-l
Note, that neglecting .Gs( t j does not produce any error for the steady state, where the derivatives are zero anyway. However. what will the error be for dynamic situations? We answer this question for periodic temperature variations of the form c(t) = &sin o! with the amplitude 11, and the angular frequency o. The output will also be periodic of the form
10”
10’
1IO2
frequency
3. Amplitude deviation between second order and first order model for sinusoidal
input.
Here. we consider a homogeneous wall with thickness d, thermal diffusivity a, thermal conductivity k, and interior and exterior surface conductance JI, and h,. respectively. Given are the temperatures at both sides of the wall as time dependent functions _x[(f 1and A>( 0. The quantity to be calculated is the heat flux at the interior surface q,( t ). The physical model for this heat flow problem is obtained by dividing the wall into m t- 1 layers of finite thickness LI = d/m or L / 2 according to Fig. 4. The temperature distribution within each layer is assumed to be homogeneous and characterised by a spatially constant node temperature
x,(t). i = 0.
. 111.
Establishing an energy balance for each layer similar to the previous example gives a system of n = m + t equations according to ( 1) with
295
Model reduction techniques
x(t) = [-~ouo(r)Xl(f) * * - Xm(t)lr, v(t) = [x1(t)
-T4(~)lT,
and the state space matrices
y(t) = 4,(t)
(shown
(39)
for n = 4) 0 0
0.5
- (1 + BiA)
1
1
-2
24
18
time 2
C = [-h,
0
0
0]
D = [h,
0]
0
-J 18
of day
model
1
n=4
--~~----~,__-_~_~-_,
,* \
%
\
;I
(41)
--
l-1 layer
.E
1 12
,
(40)
where
! 6
*’ ‘.#’
i; -1 r"
We will now investigate the following questions: How many equations (n) are required to achieve a reasonable accuracy for the resulting heat flux q,(t)? Can a state space structure different from (40) give a higher accuracy? To this end, the layer model with n = 25 equations is used as a reference. At first, we inspect the second order modes for the numerical values: k = 0.22 Wm-’ K~‘,a=0.55~10-6m2s-‘,d=0.15m,h,=8Wm~2 K-’ , hA = 20 WmL2 K-’ . Calculation of the controllability and observability grammians and determination of the eigenvalues of their product gives (sorted by descending magnitude):
18
time
24
6
time
12
=
I
I
I
6
12
time
of day
18
I 18
Fig. 5. Temperature profile of exterior temperature and calculated heat fluxes from reference and low order model.
. . -)
The dominance of the first few second order modes suggests that only a handful of equations should be sufficient. We will investigate two ways of obtaining a low order model: A layer model of the form (40) with few temperature nodes and a balanced realisation obtained from the reference mode1 according to ( 13, 14). For n = 3 the state matrices AI and A, for the layer model and the reduced model are given by: -0.7289 0.0978 0
0.1956 -0.1956 0.1956
0 0.0978 -1.5289
of day
24
18
2.3, 0.92, 0.42, 0.16, 0.048, 0.013, 0.003,
A, =
18
12
of day
Fig. 6. Deviation between the heat fluxes of the low order systems and the reference system for n = 4 and n = 5.
A, = 18
6
24
-0.2432 0.2887 -0.1655
0.2887 -0.4186 0.2585
1
- 10--3,
-0.1952 0.3489 -0.286
1
- 10-2.
Both models were used to simulate the time response of the heat flux q,(t) in response to a constant interior temperature x,(t) = 20°C and to a variable exterior temperature xA (t) with night-day profile as given on the top of Fig. 5. The resulting heat fluxes are shown below together with the response of the high order (n = 25) reference system. Fig. 6 (top) displays the deviation of the responses of both low order systems from the reference. The layer model has a peak deviation of 1.5 Wm-2 while the deviation of the reduced mode1 of the same order stays below 0.5 Wrn-‘. The difference in accuracy is even more pronounced for n
296
R.
RABENSTHN
= 4,5. Figure 6 (bottom) shows the deviations of both models for n = 4. Compared to n = 3, the deviation of the layer model is about 50% while the deviation of the reduced model is hardly visible in this scale. The peak deviations from the reference model are compiled in Fig. 7 for n = 3, 4, 5. They are normalized to the steady state heat flux for a temperature difference of 10 K. The equation set of the reduced model obviously leads to a more accurate simulation than a layer model with the same number of equations.
5.
RESPONSE-FACTOR MULTI-LAYER
METHOD
FOR
WALLS
The response-factor method [Stephenson and Mitalas ( 197 I)] is a widely-used, discrete time model for multi-layer walls. It is based on a transfer function description of the heat conduction process and represents a wall by a difference equation of the form (k is the discrete time variable) y(k) = -i
,I=,
c,y(k - v) + i h,v(k - u). ,a= 0
(42)
where the coefficients h,, c, depend in a complicated manner on the material constants of the different layers and surfaces and on the time step between two samples (usually 1 hour). We will treat two examples: a twolayer wall from Stephenson and Mitalas ( 197 1) Appendix B, 1st method, 1/B; and a four-layer wall according to Table 1. The coefficients for the two-layer wall are given in Stephenson and Mitalas ( 197 1 ), where a difference equation of degree n = 5 is used. The coefficients for the four-layer wall are listed in Table 2 with n = 6. In order to apply the model reduction technique described above, a state space description is required, which corresponds to the the difference eqn (42). Its matrices A, B, C, D are easily constructed from the coefficients b,, c, [see Kwakernaak and Sivan ( 199 1 ), Schtirjler ( 199 1)]. The reduced model is obtained according to the steps listed in Section 2. Finally the reduced state space description can be converted to a difference equation of the form (42), however with
15% i
Table I. Material
Thermal conductivity
Thickness ~__
constants
(Wm. ’EC-‘ ) _________
Cm)
for the four-layer
Heat capacity per volume -3 K ,) (kJm ,
wail --Surface
conductance (Wm 2 K ‘1 ~~~~_~__~,~~_~_~ i
0.0
15
0.175
11.94 0.035 I.1
0.100 0.1 IS _._~_._
l’O0 1800
0.35
25 7000
_-.-. ~_~._._.-.-.-- -____.-...- -_--
1, -..... ~______.
fewer coefficients. Table 3 lists the coefficients of a hrstorder difference equation for the simulation of the twolayer wall and of a second-order difference equation for the four-layer wall. Figure 8 shows simulation results for a problem similar to the one treated in the previous section. The temperature at the interior side of the wall is held constant at 20°C while the exterior side is exposed to a variable temperature as shown in Fig. 8 (top). The resulting heat fluxes are displayed in Fig. 8 (middle and bottom) for the two-layer wall and the four-layer wall, respectively. The difference equations from the response factor models of order tr = 5 and 11 = 6 are used as reference (solid lines), against which the difference equations from the reduced models are corn-pared (dotted lines). The relative deviations between reference and reduced model are below 3%: for the firstorder model of the two-layer wall and below 0.5% for the second-order model of the four-layer wall (deviations normalised to the maximum heat flux of each data set )
6.
I‘HERMAL
NETWORK
MODEL
The preceding sections covered examples of single walls, where the relation between one input quantity (exterior temperature) and one output quantity (heat flux) was considered. This example shows that model reduction techniques can also be applied to more involved models. We will treat a whole building model (though a relatively simple one), with three input quantities (exterior and ground temperature and solar radiation ) and two output quantities (room air tern.. peratures). The energy transport phenomena are not confined to walls and surfaces only. but include absorption and convective heat transport.
??
I \
layer
model Table 2. coefficients for reference of the four-layer wall
10%
modei
-_____
5% I
i;i
l%L 3
Fig. 7. Relative deviation
L
5
from the reference n = 25.
n
model of order
L
!I,
0 1
-0.0000989 0.0004437 -0.0006997 0.0012087 0.0000965 0.0002 138 --0.0000027
2 3 4
5 6
2.47225’) 2.I43.508 -0.766410 0.10393 1 ~0.0049 I h 0.000034
Model reduction techniques
297
Table 3. Coefficients of the reduced models Four-layer wall
Two-layer wall ”
b,
0
-0.347398 0.679234
1 2
C” -0.860142
Figure 9 shows a fictitious building with three zones. Floor, ceiling, and walls are modelled as storage walls, the south-oriented glass facade AF and the roof as massless walls. The capacities C,, i = 2, 3, 4, 6, 7, 8 characterize the heat capacity of the storage walls, the R, , RA , RBZ, RB6 describe heat conductance from the middle of the walls to the room air nodes or to the of the roof and exterior. RD is the heat conductance RF characterises the combined transmission and ventilation loss of the glass facade. RL describes the energy flow between the room air nodes 1 and 5 due to ventilation. The capacities C, , Cs, C9 are the heat capacity of the air volume in each room. The building is exposed to the exterior and the ground temperatures T, and Ts and to the solar radiation G. A simulation of the building energy balance shall predict the room air temperatures at the nodes 1 and 5 from the given meteorological quantities. Setting up the energy balance for the fluxes at each node gives a state space description of the form of ( I), which is equivalent to an ordinary differential equation of degree n = 9. The nonzero entries into the first line of A and B are
, 1 A,, = i = 2, 3, 4, RiCl ’ B,, =-
1 &Cl
’
B,3 = $.
b,
C”
0.00495 14 -0.0 14695 1 0.0121846
-1.791690 0.799863
7. IMPLEMENTATION
IN A SIMULATION
ENVIRONMENT
The model reduction by balanced state space realisation is primarily restricted to linear models. However, this does not prevent its application in the simulation of complex building models which also incorporate nonlinear effects. Simulation programs with a modular structure permit the replacement of linear submodels by reduced models. For instance, consider the thermal-network example in the previous section. A more accurate model could be obtained by replacing the rather crude RC-approximation of the walls by a reduced model derived from a higher-order responsefactor description of multi-layer walls. At the same time, radiative and convective energy transport could be modelled by nonlinear equations. The result would be a more detailed building description not neglecting nonlinear effects with reduced models for linear subsystems. In order to show the applicability of model reduction methods in complex building models, an early version of the Energy Kernel System [ Charlesworth et al. ( 199 1), Rabenstein ( 1993)] was chosen as a simulation environment. The Energy Kernel System (EKS) is not a single simulation program, rather it allows the user to build a custom-made run-time code
1 AIS = RLC, ’
I
The state variables xi(t) i = 1, . . . , 9 are the temperatures at the nodes. The input and output quantities are
G(Ol’,
v(t) = ITA(
T,(r),
Y(l) = [x,(t)
xs(t)lT.
_J
)-lo1 1
10
15 time in
20
25
30
15 20 in days
25
30
days
.‘~~~ I 1
This state space model of order n = 9 was used as a starting point for a model reduction to a reduced system of order n = 3. Simulation runs for both systems with typical values for Ri and Ci , hourly mean values for the exterior temperature TA and the global radiation G, and constant ground temperature Ts gave two predictions for the time history of the node temperatures 1 and 5. The output of the reduced model did not deviate more than 0.3 K from the reference model at either node.
5
5
10 time
time in doys
Fig. 8. Ambient temperature (top), and heat fluxes for the two-layer wall (middle) and the four-layer wall (bottom). Solid lines: reference model: dotted lines: reduced model.
R. RABENS I EIN
298
F
A
TA
TA
TB
TB
B
Fig. 9. Thermal
B
network model of a building.
From a set of existing program modules. The choice of modules contains different physical models of varying level of detail and various numerical models and solvers For the resulting equation systems. Furthermore, its modular structure, based on an object-oriented approach, facilitates the addition of new program modules without alteration of existing program code. A finite difference model for the calculation of heat zonduction through multi-layer walls, already existing n the EKS, was chosen as a reference. It produces a Cscrete space, continuous time model, similar to the me presented in Section 3. However, for arbitrary n&i-layer walls. A new program module was derived using the principle of inheritance in the object-ori:nted sense), which uses the state space formulation If the finite difference model and replaces it by a balmced reduced model. By choosing either the existing inite difference model or, alternatively, the derived .educed model, it is possible to construct simulation )rograms, which are identical except for the modelling If heat conduction through walls. Both program versions were tested for the layernode1 and the response factor method problem dis:ussed earlier. A variety of multi-layer walls has been nvestigated. In each case a reference model was con;tructed using the finite difference method with 20 lodes distributed over the different layers of the wall. rhen lower-order models were found by either reducng the number of nodes in the finite difference model )r by model reduction. Simulation runs with identical nput temperature data showed the trade-off between Iecreased model order and deviation from the refer:nce model. For a quantitative comparison, a number of cases vere identified, where a low order finite difference node1 and a reduced model gave approximately the ,ame deviation from the reference (typically 0.1 Yo- I YG) t turned out that the required order of the finite diftrence models was in the range of 9-13, while the
reduced order models achieved the same or higher accuracy with a degree of only 4-5. Furthermore, simulation runs for finite difterence and reduced models with comparable accuracy were timed on a workstation (SUN SPARCstation WC). lnput data sets with l-year duration of hourly mean values were used. The run time for the reduced models was in the range of 19%36Yo of the run time of the low-order finite difference models.
X. (‘ONCLI
.SlOh
Theoretical considerations and practical implementations have shown that the model reduction technique by balanced state space representation is a useful tool for building energy simulation. It allows to perform simulations with the same accuracy as standard analysis methods. however with a significantly decreased model order. This method is applicable to linear building models and to linear submodels of general nonlineal building descriptions. The calculation of the reduced order state space description from a higher order reEerence is a straightforward numerical procedure. requiring matrix calculus and standard linear algebra routines. Some higher level languages (e.g. MAT1,.4B 1 provide powerful commands for model reduction. .‘1~,~~oll,/~,~~~~~~~~~~\Ihe author
wishes to thank: M. Gerwien, Motz. and .A. Wegehaupt for the programming of the esamples; D. Heinemann for supplying the meteorological data: and J. Clarke. D. Tang, and D. Mac Randal for an introduction to the Energy Kernel System.
r.
RkXERENCES F. Buhl. E. Erdem,
J. -M. Nataf, F. C. Winkelmann. M. A Moshier. and E. F. Sowell. 7%~~US. EKS, Proc. of/he 3nd
Inl. Confrrenw on S~.wm Sirnulution in BuildinKs, LGge. pp. 107-149 (1990). F. M. Callier. C. A. Desoer. Lvzmr s~:trern rhcw,~. SprmgerVerlag, New York ( 199 I ).
Model reduction techniques J. A. Clarke, Energy simulation in building design, Adam Hilger Ltd, Bristol ( 1985). P. Charlesworth, J. A. Clarke, G. Hammond, A. Irving, K. James, B. Lee, S. Lockley, D. Mac Randal, D. Tang, T. J. Wiltshire, A. J. Wright, The energy kernel system, Proc. of Building Simulation ‘91, Nice, pp. 313-322 (1991). L. Fortuna, G. Nunnari, A. Gallo, Model order reduction techniques with applications in electrical engineering, Springer-Verlag, London ( 1992). H. Kwakernaak, R. Sivan, Modern signals and systems, Prentice Hall, Englewood Cliffs, NJ ( 199 I ). A. Laub, Computation of “‘balancing”transformations, Proc. Joint American Control Conference, I, paper FA8-E (1980). G. Lefebvre, A. Neveu, K. el Khoury, J. J. Saigon, Applying the modal method to thermal modelling, Proc. of the 9th Int. Heat Transfer Conference, pp. 157- 16I, 3, Jerusalem (1990). Y. Liu, B. Anderson, Singular perturbation approximation of balanced systems, Int. I. Control 50(4), 1379-1405 (1989). B. C. Moore, Principal Component Analysis in Linear Systems, IEEE Trans. on Automatic Control, Vol. AC-26( I ), 1732(1981).
299
B. Peuportier, 1. Blanc Sommereux, Simulation tool with its expert interface for the thermal design of multizone buildings, Int. I. Solar Energy 8, 109-120 (1990). R. Rabenstein, Modellreduktion als Werkzeug zur ejektiven Simulation des Gebaudewarmehaushalts. Proc. of 8. Internationales Sonnenforum, Berlin, pp. 152-157 ( 1992). R. Rabenstein: Erste Erfahrungen mit dem Energy Kernel System zur energetischen Gebaudesimulation. Proc. of 3. Nat. Symposium Thermische Solarenergie, Regensburg, Germany, pp. 259-263 ( 1993). J. G. Reid, Linear system,findamentals, McGraw-Hill, New York (1983). P. Sahlin and A. Bring, IDA Solver, Conference Proceedings Building Simulation ‘9 1,Nice, pp. 339-348 ( 199 1). H. W. SchilBler, Netzwerke, Signale und S.ysteme 2, 3. ed., Springer-Verlag, Berlin ( 1991). I. S. Sokolnikoffand R. M. Redheffer, Mathematics ofphysics and modern engineering, 2nd ed., M&raw-Hill, New York (1966). D. G. Stephenson and G. P. Mitalas, Caluculation of heat conduction transfer,finctions for multi-layer slabs, ASHRAE Transactions, 77, Part II, pp. I 17-126 ( I97 I ). D. Tang, The generalised system solution classes m the EKS environment, Conference Proceedings Building Simulation ‘9 I, Nice, pp. 323-327 ( I99 I ).