On a Markovian approach for modeling passive solar devices

E n e r g y a n d B u i l d i n g s , 17 (1991) 103-116

103

On a Markovian approach for modeling passive solar devices F. B o t t a z z i and Th. M. L i e b l i n g C h a i r e de R e c h e r c h e Opdrationnelle, Ecole P o l y t e c h n i q u e FddZ,r a l e d e L a u s a n n e , CH-1015 L a u s a n n e (Switzerland)

J.-L. Scartezzini

a n d M. N y g ~ t r d - F e r g u s o n

L a b o r a t o i r e d ' E n e r g i e S o l a i r e et de P h y s i q u e d u B 5 t i m e n t , Ecole P o l y t e c h n i q u e Fdddrale de L a u s a n n e , CH-1015 L a u s a n n e ( S w i t z e r l a n d )

(Received November 21, 1989; accepted March 11, 1990; revised paper received August 30, 1990)

Abstract Stochastic models for the analysis of the energy and thermal comfort performances of passive solar devices have been increasingly studied for over a decade. A new approach to thermal building modeling, based on Markov chains, is proposed here to combine both the accuracy of traditional dynamic simulation with the practical advantages of simplified methods. A main difficulty of the Markovian approach is the discretization of the system variables. Efficient procedures have been developed to carry out this discretization and several numerical experiments have been performed to analyze the possibilities and limitations of the Markovian model. Despite its restrictive assumptions, it will be shown that accurate results are indeed obtained by this method. However, due to discretization, computer memory requirements are more than inversely proportional to accuracy.

1. I n t r o d u c t i o n

Energy performances and thermal comfort i n d i c a t o r s of p a s s i v e s o l a r devices are highly d e p e n d e n t on m e t e o r o l o g i c a l c o n d i t i o n s a n d on the influence o f u s e r s ' b e h a v i o r on the s o l a r s y s t e m [1]. As p o i n t e d out in ref. 2, w e a t h e r a n d u s e r - d e p e n d e n t d a t a are c h a r a c t e r i z e d b y significant r a n d o m fluctuations: c o n s e q u e n t l y , t h e r m a l p e r f o r m a n c e i n d i c a t o r s of the simulated s o l a r s y s t e m v a r y significantly f r o m y e a r to year. In o r d e r to avoid long c o m p u t e r r u n s u s i n g n u m e r o u s real s e q u e n c e s of d a t a as input files (20- or 3 0 - y e a r - l o n g w e a t h e r files, s i m u l a t i o n of different u s e r s ' b e h a v i o r ) , d y n a m i c simulation p r o g r a m s [ 3 - 7 ] are g e n e r a l l y u s e d with s h o r t e r s e q u e n c e s of d a t a s u c h as 'typical m e t e o r o l o g i c a l y e a r s ' . T h e latter p o r t r a y the typical b e h a v i o r o f the w e a t h e r c o n d i t i o n s for a p e r i o d of s e v e r a l years, b u t u n f o r t u n a t e l y c o n t a i n no i n f o r m a t i o n a b o u t the variability of the m e t e o r o l o g i c a l quantities.

0378-7788/91/$3.50

C o r r e l a t i o n t e c h n i q u e s [ 8 - 1 0 ] a n d simplified m e t h o d s [ 11, 12 ] h a v e the m a i n a d v a n t a g e of simplifying significantly the c o m p u t a t i o n s . H o w e v e r , as n o t e d in ref. 2, their u s e of a v e r a g e i n p u t d a t a a c c o u n t s for s e r i o u s limitations w h e n • c o n s i d e r i n g the r a n d o m n a t u r e of w e a t h e r a n d u s e r s ' b e h a v i o r [ 13]; • e v a l u a t i n g the c o m f o r t c h a r a c t e r i s t i c s of the s o l a r device a c c o r d i n g to F a n g e r ' s t h e o r y b a s e d on probabilistic n o t i o n s [ 14 ]. The M a r k o v i a n a p p r o a c h p r e s e n t e d h e r e w a s p r o p o s e d b y Scartezzini in ref. 15, who, following earlier investigations o f H a s l e t t [16, 17], Devine [18] a n d L a m e i r o [19, 20], s h o w e d h o w s h o r t - t e r m (hourly or half-hourly) distrib u t i o n s of t h e r m a l p e r f o r m a n c e i n d i c a t o r s c o u l d b e c o m p u t e d f r o m the s t a t i o n a r y dist r i b u t i o n of a M a r k o v c h a i n w h i c h c a n t a k e a c c o u n t of all available d a t a of w e a t h e r a n d user-dependent behavior. D y n a m i c m o d e l s for p a s s i v e s o l a r d e v i c e s are usually b a s e d on a n o d a l d e c o m p o s i t i o n t e c h n i q u e . Basically, o u r M a r k o v i a n s t o c h a s t i c a p p r o a c h u s e s the s a m e kind of t h e r m a l build-

© Elsevier Sequoia/Printed in The Netherlands

104

ing m o d e l a n d is b a s e d on the s a m e m a t h e m a t i c a l e q u a t i o n s ( e n e r g y c o n s e r v a t i o n equations) as a traditional d y n a m i c simulation. H o w e v e r , the f u n d a m e n t a l difference b e t w e e n the two t e c h n i q u e s is t h a t the driving v a r i a b l e s (solar radiation, o u t d o o r t e m p e r a t u r e , free gains and o t h e r v a r i a b l e s r e l a t e d to user-dep e n d e n t b e h a v i o r ) are m o d e l l e d via a discrete Markov chain [21 ]. E a c h driving v a r i a b l e c a n be r e p r e s e n t e d b y a probability m a t r i x (Mark o v matrix), which m a k e s up the input of the M a r k o v i a n s t o c h a s t i c p r o g r a m . F r o m this ass u m p t i o n , it c a n be c o n c l u d e d that the v e c t o r of driving a n d s t a t e v a r i a b l e s ( r o o m a n d building e l e m e n t t e m p e r a t u r e s ) is also a M a r k o v chain. In d y n a m i c simulations, a finite differe n c e t e c h n i q u e is u s e d to solve the evolution e q u a t i o n s of the building, w h e r e a s in the Markovian a p p r o a c h a M a r k o v m a t r i x (monitor m a t r i x ) is u s e d to m o d e l this t h e r m a l behavior. F r o m that matrix, e n e r g y c o n s u m p tion and t h e r m a l c o m f o r t distribution, as defined by F a n g e r in ref. 14, can b e c o m p u t e d directly. Although this a p p r o a c h d o e s not i m p l y a n y kind of simulation, we s o m e t i m e s r e f e r to it, in a s o m e w h a t a b u s i v e sense, as a M a r k o v i a n simulation since it c a n p r o d u c e o u t p u t r a t h e r

similar to t h a t o f traditional s i m u l a t i o n methods. Figure 1 s h o w s a c o m p a r i s o n b e t w e e n deterministic d y n a m i c a n d M a r k o v i a n s t o c h a s t i c simulations. In the following p a r a g r a p h s , w e briefly p r e s ent the b a s i c m a t h e m a t i c a l a s p e c t s of this m e t h o d as well as s o m e issues a b o u t its perf o r m a n c e . P a r t i c u l a r a t t e n t i o n will b e given to the discretization of the m o d e l v a r i a b l e s a n d to its influence on the a c c u r a c y o f the results. A m o r e detailed d e s c r i p t i o n of this m e t h o d is given in refs. 15 a n d 22. A n o t h e r o v e r v i e w of the m e t h o d is given in ref. 2.

2. The M a r k o v i a n s t o c h a s t i c m e t h o d D y n a m i c n o d a l m o d e l s for p a s s i v e s o l a r dev i c e s are m a d e up o f a set o f N n o d e s chara c t e r i z e d b y their m a s s a n d their t h e r m a l capacity. The i n t e r a c t i o n s b e t w e e n n o d e s are r e p r e s e n t e d b y t h e r m a l c o n d u c t o r s t h a t acc o u n t for h e a t e x c h a n g e s b y c o n d u c t i o n , conv e c t i o n a n d radiation. T h e s e p h e n o m e n a are linearized t h r o u g h o u t , a l t h o u g h the a p p r o a c h w o u l d also m a k e it p o s s i b l e to t a k e into a c c o u n t the nonlinearities. STOCHASTICSIMULATION Markovmatrix of rnet¢odata

DETERMINISTICSIMULATION

,T

/~r

I

,~ronctua[ resu',ts statistic of auxiliary needs

t

sy:

] globalresults

"V t

statistic of auxiliary needs

/ =

~

CT(aux)

• <::a u x ::>

"

Fig. 1. C o m p a r i s o n of deterministic dynamic simulation and the Markovian stochastic approach. Inputs of the former are hourly data files, inputs to the latter are probability matrices (Markov matrices). The s a m e output results can be obtained by both methods.

105

The e n e r g y c o n s e r v a t i o n principle leads to the evolution e q u a t i o n s which, after s o m e easy and well-known m a t h e m a t i c a l m a n i p u l a t i o n s (finite difference t e c h n i q u e ) , take the f o r m of a set of linear r e c u r r e n c e e q u a t i o n s Xk+l = F X k + G V k

k = O , 1, 2, ...

(1)

w h e r e X k ~ I R n, V ~ I R p, F E M ~ × , ( I R ) , G ~ Mn×p(IR) and w h e r e • Xk is the s y s t e m state v e c t o r at time tk defined by n < N internal n o d e t e m p e r a t u r e s ( t e m p e r a t u r e s of air and building e l e m e n t s ) (state variables), • Vk is the s y s t e m driving v e c t o r at time tk defined by N - n external node temperatures ( a m b i e n t air t e m p e r a t u r e , adjacent r o o m s ) and by internal and e x t e r n a l heat s o u r c e s ( d r i v i n g variables), • m a t r i c e s F and G a c c o u n t for physical par a m e t e r s of the s y s t e m ( t h e r m a l capacities, c o n d u c t a n c e s , etc.). The p r o c e s s {Vk} is a s s u m e d to be a Markov chain and so {Wk} defined by Wk = (Xk, Vk) is also a d i s c r e t e t i m e Markov p r o c e s s , t h o u g h s o m e f u r t h e r a s s u m p t i o n s m u s t be m a d e in o r d e r to obtain a Markov chain ( d i s c r e t e states). Let E X c I R n and E V c I R p be the d i s c r e t e s p a c e s w h e r e v e c t o r s Xk and Vk take their values r e s p e c t i v e l y so that Wk ~ E X × E v = E w. W e define a discretization f u n c t i o n Ax: IR n

>E x

(2)

that associates s o m e p a r t i c u l a r discrete values of E x to e a c h point of IR n. The transition m a t r i x M w of the Markov p r o c e s s {Wk} ( m o n i t o r matrix) is o b t a i n e d f r o m that of {Vk} by defining the conditional probability of a transition f r o m state wk to state wk+~ by

of Wk is o b t a i n e d by r e p e a t e d l y multiplying an initial arbitrary distribution by this m a t r i x until either a c o n v e r g e n c e criterion is m e t or a given m a x i m u m n u m b e r of iterations is r e a c h e d (see refs. 2 and 15). This distribution is u s e d to c o m p u t e the m e a n , v a r i a n c e and d i s t r i b u t i o n of the main variables of interest s u c h as h o u r l y back-up e n e r g y c o n s u m p t i o n , t h e r m a l comfort, t h e r m a l losses and direct solar gains. In fact, s u c h entities are simply e x p e c tations E ( f ) of s o m e functional f : IR n+~

, IR

(5.1)

computed from ~-by E(f)=

~ ~ E

(5.2)

f(w)1r(w) TM

This m e t h o d has b e e n i m p l e m e n t e d in prog r a m SOLSTIS. Figure 2 shows a flowchart of this p r o g r a m . E a c h b l o c k c o r r e s p o n d s to a particular step of the p r o c e d u r e previously described. One of the major a d v a n t a g e s of this m e t h o d is that input data n e e d e d to c o m p u t e building e n e r g y and c o m f o r t p e r f o r m a n c e s during long periods of time (several d e c a d e s ) are cond e n s e d into a relatively small set of p a r a m e t e r s

IoFST+l CH SOLAR DEVICE

1

A * T MARKOV MATRIXES FOR I [e AND FGvS

TIME EVOLUTION I

*

1

T MARKOV MATRIXES FOR Pe AND Po

I

~vl EQUAEONSOF l~I~--~ . . . . .

I

1I

Pr{Wk + 1 = wk + 11Wk = wk} = Pr{Vk +~ = vk +~ IVk = vk} if xk+ 1

=

Ax(FXk + Gl)k)

STATE VECTOR ~

+

(3) i

= 0 otherwise for e a c h pair (wk+~, w k ) ~ E W × E W, k ~ I N , w h e r e Pr{Vk + ~= vk+ ~IVk = vk} is the conditional probability of a transition f r o m state Vk to state vk+~ given by the transition m a t r i x M v of the Markov p r o c e s s {Vk} r e p r e s e n t i n g the driving variables. The s t a t i o n a r y distribution rr: E w

, [0, 1]

(4)

OPTIONALiTERATION(S) --------------------------FOR PARAMETRIC STUDIES

F

LONG RUN VECTOR OF PROBABILITIES

ENERGy 1

AND CO~AFORT I INDICATORS t

................

"~ MARKOVMATRIX CALCULATEDONLY

ONCE FOR ANY NUMBER OF COMPUTER RUNSOF THE PARAMETRIC STUOIES

Fig. 2. Flowchart diagram of SOLSTIS. A Markov matrix characterizing the evolution of the building state vector is calculated using the input Markov matrices and the energy conservation equations within the building nodal model.

106

( M a r k o v m a t r i c e s ) . T h e s e m a t r i c e s n e e d to be e s t i m a t e d only o n c e for a p a r t i c u l a r g e o g r a p h i c p l a c e a n d simulation p e r i o d a n d for e a c h driving variable, so t h a t p r o c e s s i n g of large input d a t a files of m e t e o r o l o g i c a l r e c o r d s do not n e e d to take p l a c e at e a c h simulation run.

In o r d e r to c a r r y out the c o m p u t a t i o n for the M a r k o v i a n a p p r o a c h , it is n e c e s s a r y to p r o c e e d b y the following o p e r a t i o n s : • discretize the driving variables, • build up the input M a r k o v m a t r i c e s , • build up m a t r i x M v, • discretize the state variables, • build up the m o n i t o r m a t r i x M W. The w a y M a r k o v m a t r i c e s of driving v a r i a b l e s (input of M a r k o v i a n s t o c h a s t i c p r o g r a m ) are c a l c u l a t e d f r o m h o u r l y (or half-hourly) d a t a files is illustrated in Fig. 3. The c o n t i n u o u s p r o c e s s r e p r e s e n t i n g e a c h driving v a r i a b l e is first s a m p l e d during the m e a s u r e m e n t proc e d u r e ( s a m p l e d data). After a suitable discretization is defined, a s e q u e n c e of s t a t e s is o b t a i n e d for e a c h driving v a r i a b l e ( d i s c r e t e data). Taking into a c c o u n t the w o r k of Billingsley [23], e a c h transition p r o b a b i l i t y is s i m p l y d e t e r m i n e d by c o u n t i n g the n u m b e r of j u m p s o c c u r r i n g at e a c h s a m p l i n g t i m e interval f r o m one state to a n o t h e r (input M a r k o v m a trices). A s s u m i n g s t o c h a s t i c i n d e p e n d e n c e bet w e e n the driving variables, m a t r i x M v is s i m p l y

i

~

o

SOLAR RADIAT;ON

~ ~

~Tk

IllI

,,,1[ h,

o;

;

oht,,,ll

IIIl!llllllllIImi

'

i

~5 ,g~,2'~ ab d6 4} ~5° [W~

• < ~

~h~

'i4

,ooo4

~

OUT[X)ORTEMPERATURE

3'o ~'G 4 4 °

ci .;rimI

~

['C ] :0

0

0 200500 093 007 0D0

[

300 O~ 064 0.18 6OO 0O0022078

r~q [%3 4

1

6

8

6 004 092 004 ~ OOO 01C O9O

Fig. 3. Discretization a n d s e t - u p o f i n p u t m a t r i c e s for M a r k o v i a n s t o c h a s t i c s i m u l a t i o n . M o n i t o r e d d a t a of driving v a r i a b l e s (solar r a d i a t i o n a n d o u t d o o r t e m p e r a t u r e ) a r e t r a n s f o r m e d in probability m a t r i c e s by t h e following operations: t h e y a r e first d i s c r e t i z e d a n d t h e n t r a n s i t i o n s f r o m s t a t e to s t a t e are c o u n t e d to yield t r a n s i t i o n probabilities ( M a r k o v m a t r i c e s ) .

Outdoor temperature

[oc) 0,490

[W] 0 300 600 0 F0.93 0 0 7 O.O0] 3000.180.640.18 600 L 0.00 0.22 O.TS

F OLOOOO 1

6 8 4 6 0.04 •092 0.04 8 [-0.08 0.10 0.90

\

/

®

(IW],[°C])

3. Discretization procedures

[ W I mi ] 1000

Solar radlatlon

(0,4) I 0,6) (0,8) (300,4) (300,6) (300,8) (600,4) (600,6)

(

0,4)(0,6)

[--08370 --00372 0.0000 0.1620 00072 0.0000 00000 0.0000

16oo 8~

0.0930 08556 00930 00180 0.1656 0.0180 0.0000 00000

ooooo o.0ooo

(0,8)

(300,4) (300,6) (300,8) (600,4) (600,6) (600,8)

O000O 00372 08370 00000 0.0072 0.1620 OOO00 O.OOO0 OOO00

0.0630 0.0028 00000 0.5760 0.0256 00000 01980 00088 00000

00070 00644 00070 00640 05888 00640 00220 02024 0.0220

OOOOO 00028 0.0630 O0000 0.0256 0.5760 00000 00088 0.1980

0.0000 O.OOO0 O.O000 O.OOO0 00000 0.0000 01620 0.0180 00072 01656 O.O000 0.0180 0.7020 0.0780 00312 0.7176 O.OOOO 0.0780

0.00007 0.0000 O.O00O 0.0000 0.0072 0.1620 0.0000 0.0312 0.7020

Fig. 4. Build-up of m a t r i x M v f r o m t h e i n p u t Markov matrices of solar radiation and outdoor temperature.

defined f r o m the input M a r k o v m a t r i c e s b y Pr{Vk +1 = vk + 11Vk = vk} P

= 1-I Pr{Y~+l = V~+l l Y~ = v~}

(6)

i=l

w h e r e Vk=(Y~, Y~ . . . . . V~)~E V, vk=(v 1, v~ . . . . . v ~ ) ~ E v, k E I N and P r { V ~ + I = V ~•+ 1 I V k i -_ v ~ } is the conditional p r o b a b i l i t y of a transition f r o m state v~ to s t a t e v~+l of driving variable i given by its c o r r e s p o n d i n g input M a r k o v matrix. This is illustrated in Fig. 4 for the s i m p l e c a s e of two driving v a r i a b l e s (solar r a d i a t i o n a n d o u t d o o r t e m p e r a t u r e ) . The p r o p e r definition of the discretization of driving a n d s t a t e v a r i a b l e s i n t r o d u c e s f u r t h e r difficulties w h e n c o m p a r e d to d e t e r m i n i s t i c d y n a m i c simulations. P r o c e d u r e s h a v e b e e n d e v e l o p e d in o r d e r to m a k e this as e a s y as possible.

Driving variables A t w o - s t e p p r o c e d u r e to discretize the driving v a r i a b l e s h a s b e e n set up and v a l i d a t e d b y numerical experiments: • s u b d i v i d e the r a n g e of v a r i a t i o n of e a c h v a r i a b l e into m equal intervals Ii, i ~ { 1 , 2, . . . , m}, • c h o o s e the discrete value v* of e a c h interval as the a v e r a g e of the m o n i t o r e d v a l u e s within it. More specifically, let M={vl, v2 . . . . . VNT} b e the s a m p l e set of the N r available m e a s u r e m e n t s for a given driving variable, say the o u t d o o r t e m p e r a t u r e Te (for e a s e of e x p o s i t i o n , we a s s u m e t h a t v~ ~: Vy w h e n e v e r i ~:j). The m intervals Ii = [4-1, 4] are defined b y

107

• lo the lower bound of the first interval, • l m the upper bound of the last interval, • 8 the length of the intervals, • m the num ber of intervals, that are related by

vL . . . . . .

T'7',

i '".,",", i/,"/"

. . . . . . .

i

!

i"

', . . . . . . . . .

.

S [ 1

12

I; ]

lm -

l0 = 8 . m

(9)

13

Fig. 5. The t w o - s t e p driving variables discretization proc e d u r e ( m = 3). v~ = o b s e r v a t i o n at time i o f driving variable v, 8 = (I a - / 0 ) / 3 , v* = m e a n o f o b s e r v a t i o n s b e l o n g i n g to interval Ii.

is needed to specify a discretization for a particular state variable. The bounds l~ and the values v* for each interval Ii = [li-1, li] are given by

l i = i l'~-l----2 + l o = i S + l o m / o = m i n v,

/re=max v

v~M

i ~ { 0 , 1, . . . , m}

(7.1)

(10.1)

v~M

and

lm l~=-i--

-

lo

m

+lo

i~{1,2,...,m-1}

(7.2)

lm v*=(2i-1)--~m

lo

The values v* associated to the intervals I~ are defined by v*= I

~

v

i e { 1 , 2, . . . , m}

(7.3)

N i vEMnli

where N~= IMnI~l, i ~ { 1 , 2, . . . , m}, is the number of monitored values belonging to 1~. The discretization function A~o: IR

, ET~--v--{v*, V'~, ..., V*}

(8.1)

for the o u td o o r tem pe r at ur e T~ is now simply defined by A~,(v) = v*

(8.2)

where i is such that v ~I~. Figure 5 summarizes the two-step discretization p r o ced u r e for the driving variables. It is apparent that this pr oc e dur e preserves the sample mean of M. Both the o u td o or t em pe r at ur e T~ and the solar radiation F~vs have been discretized within 3 intervals with satisfying results (see ref. 22).

8

+/0=(2i-1)~

+/o

i ~ { 1 , 2, . . . , m}

(10.2)

This is illustrated in Fig. 6. Three different mathematical procedures have been developed to carry out this discretization. They are straightforward and automatic in that, once values for some set of parameters are appropriately selected, they are functions of the discretization of the driving variables and only require specification of the building model. These procedures have been validated via numerical experiments (see ref. 22). They are also two-step procedures. In the first one, being the same in the three procedures, we compute the lower and upper bounds of all state variables. One param et er a controls the range of variation of all variables. In the next step, we give three different approaches to compute the interval lengths of the state variables. According to the selected procedure, S

\ :

'a

State v a r i a b l e s The appropriate discrete values of the state variables are more difficult to define: this implies an a priori knowledge of the possible variations of these variables, which cannot be easily obtained. For each state variable, we define m contiguous intervals of equal length each of which is associated with a value given by its midpoint. A set of four parameters

11 II

i-1

Ii [i

lm. I

lm rn

Fig. 6. The discretization function A for state variables. S = s e t o f all p o s s i b l e values o f a state variable, I ~ = i t h interval with b o u n d s l~_~ and l~, v* = m i d p o i n t of interval I,, 8 = l e n g t h o f t h e intervals, A ( x ) = d i s c r c t e value of x given by A ( x ) = v * w h e r e i is s u c h that x~I~.

108 one further p a r a m e t e r m a y be n e e d e d in o r d e r to globally control these interval lengths. Let Vinf d e n o t e the v e c t o r o f lowest values of the driving variables. T h e n the state variables lowest values are defined b y the following vector Xin f =

(I

--

F) -1 GVinf

(11)

Note that this is the equilibrium state r e a c h e d by the s y s t e m if it is p e r m a n e n t l y excited by a c o n s t a n t input Vine. Similarly, d e n o t e by Vsup the v e c t o r of the driving variables highest values, slightly modified by , v if i ~ F G v s (V~up)i = f v,,~,i a if i = F G v s

That is, ~v is the smallest interval length for driving variable Vj and 5 v is the j u m p b e t w e e n its highest and lowest value. V e c t o r 8x is defined, • for p r o c e d u r e A, by 6 X= 2 . m i n gij 8v

i ~ { 1 , 2, . . . , n}

j~{jlgij~O}

(16.1) • for p r o c e d u r e B, by 5x = 2. rain gij Ay

i ~ { 1 , 2, . . . , n}

j~'lgij¢'O}

(16.2) • for p r o c e d u r e C, by

i ~ { 1 , 2 . . . . ,p}

At t~x=2" ~-~f~i 8P

i ~ { 1 , 2 . . . . , n}

(16.3)

(12) w h e r e a is a p a r a m e t e r and v,~.~*vis the highest value of driving variable V~. Again, we define the highest values of the state variables b y the vector X~up = (I - F ) -1GV~op

(13)

The choice ol ---v ~*V GVS,

(14.1)

FGVS

m a y be theoretically justified in that it guara n t e e s that all r e a c h a b l e stochastic states are t a k e n into account, given the s t r u c t u r e of matrices F and G d e s c r i b e d in [2]. However, a large n u m b e r of these states are c h a r a c t e r i z e d by a null probability, i.e. they are t r a n s i e n t . Numerical e x p e r i m e n t s have s h o w n o t h e r choices are also s o m e t i m e s satisfying. W e have tried the following: a = (FGvs>

(14.2)

a = 2

(14.3)

1

,v

a = ~- V,,~FGvs.Fow

(14.4)

w h e r e <.) d e n o t e s mean. As a n e x t step, we c o m p u t e the v e c t o r Bx of interval lengths for the state variables according to the s e l e c t e d p r o c e d u r e . Let us first define two v e c t o r s 8 v and Av by 8v =

min i--1,

(v'J+1-v *j)

j ~ { 1 , 2, . . . , p}

..., vnj-- 1

and 5V=v~-v7

(15.1) j ~ { 1 , 2, . . . , p}

(15.2)

w h e r e f o and g~j are e l e m e n t s (i, 3) of m a t r i c e s F and G respectively (note that they are all non-negative), C, is the t h e r m a l capacity ass o c i a t e d to n o d e i, At is the sampling time interval (usually half an h o u r or an h o u r ) and 8P is a p a r a m e t e r . The idea b e h i n d the above formulas is best u n d e r s t o o d in t e r m s o f p e r t u r b a t i o n of s o m e equilibrium state X, i.e., one satisfying X = F X + GV

(17)

for a given excitation V. P r o c e d u r e A states that the s m a l l e s t c o n t r i b u t i o n o f a n y j u m p in a n y driving variable m u s t p e r t u r b all state variables after discretization. In p r o c e d u r e B we are less restrictive and require only the s m a l l e s t c o n t r i b u t i o n o f the largest j u m p in a n y driving variable to cause a j u m p in all state variables. For the last p r o c e d u r e , we a s s u m e that a fictitious p o w e r s o u r c e is applied to n o d e i and require that a t h r e s h o l d p e r t u r b a t i o n p o w e r ~P m u s t cause at least a one step j u m p i n s t a t e v a r i a b l e Xi. This r e a s o n i n g is applied i n d e p e n d e n t l y to all state variables. Note that the auxiliary p o w e r heating is not explicitly discretized, so we take

Bv=~

and

AV=~

in eqns. (16.1) and (16.2) whenever j designates auxiliary heating. Thus procedures A and B are parameterless so that the state discretization is determined once one of the driving variables is defined. Note that procedure A yields a finer discretization than procedure B does.

109

Transient states e l i m i n a t i o n One of the m a j o r d r a w b a c k s of the discretization p r o c e d u r e s j u s t d e s c r i b e d is t h a t t o o fine a discretization u n a v o i d a b l y leads to a v e r y large n u m b e r of s t o c h a s t i c s t a t e s ( m o r e t h a n 8 0 0 000). H o w e v e r , only a f r a c t i o n of t h e m (often 10%--30%) h a v e a non-null s t a t i o n a r y probability. Intuitively, this m a y be e x p l a i n e d b y noting t h a t t w o tightly c o u p l e d n o d e s s h o w v e r y similar t e m p e r a t u r e b e h a v i o r . T h e stat i o n a r y p r o b a b i l i t y o f a s t a t e c o r r e s p o n d i n g to v e r y different t e m p e r a t u r e s for t h e s e two n o d e s is likely to b e null or v e r y n e a r zero. Such a state is a s s u m e d to be t r a n s i e n t a n d d o e s n o t n e e d to b e c o n s i d e r e d in the m o d e l . In the following lines we i n t r o d u c e a m e t h o d to r e d u c e the n u m b e r of s u c h states. First of all, d e n o t e • ~.£ij=(Xi--X3.), i , j ~ { 1 , 2, . . . , n}, the m e a n temperature deviations between each couple o f internal n o d e s , a n d • o'~j, i , j ~ { 1 , 2 . . . . , n}, the c o v a r i a n c e o f the t e m p e r a t u r e s for e a c h c o u p l e of internal nodes. Now, i n s t e a d o f discretizing X~ in interval [(X~f)i, (Xsu,)il

(18.1)

w e discretize the l a p s e Oi --~-Xi - X j -

the l a p s e s v e c t o r 0 = [ 0 1 0 2 . . , On] T, T h e y are m o d e l e d b y a d i r e c t e d g r a p h ( d i g r a p h ) / ' a = (V, U) s u c h t h a t the set of v e r t i c e s V c o r r e s p o n d s to the s e t of s t a t e n o d e s plus the fictitious n o d e 0 a n d the set of a r c s U is defined b y (i, j ) e U if Xj depends on X~ by X 3.= Or +X~ + ]*£ij

(2o) If /'a is a tree h a v i n g a root 0 (an arbor e s c e n c e ) t h e n the s t a t e X is uniquely determ i n e d b y 0, t h a t is, the d e p e n d e n c e s are s u c h that: (i) e a c h n o d e b u t the r o o t d e p e n d s on one a n d only o n e o t h e r n o d e , a n d (ii) following the d e p e n d e n c e s in their opp o s i t e sense, we e n d up at fictitious n o d e 0 (the r o o t of the a r b o r e s c e n c e ) . Figure 7 illustrates this. T h e p r a c t i c a l m e a n i n g of all the a b o v e definitions is t h a t n o w we c a n easily s t a t e the p r o b l e m of finding the set of d e p e n d e n c e s t h a t m i n i m i z e s the total n u m b e r of d i s c r e t e s t a t e v a l u e s iEXj. F o r this p u r p o s e , c o n s i d e r the d i g r a p h F = (V, ~ ) of all p o s s i b l e d e p e n d e n c e s b e t w e e n n o d e s , i.e., V = { 1 , 2, . . . , n}

~tij in interval [ - ~sij , ~sij ]

~={(i,j)~V×Vli¢j, (18.2)

w h e r e fl is a p a r a m e t e r a n d sij is the s t a n d a r d deviation of X ~ - X j g i v e n b y

sij "-------~/Var(Xi - X j ) = ~o'ii + (r2. - 2cr~y

(19)

In o t h e r w o r d s , we i n t r o d u c e a d e p e n d e n c e of the c o m p u t a t i o n o f X i on t h a t of Xj. X 3. m a y in t u r n d e p e n d on a n o t h e r n o d e a n d so on until we r e a c h a fictitious n o d e 0 with X 0 = ~ o ~ = 0 , i ~ { 1 , 2 . . . . , n}, to e n d the rec u r r e n c e . T h e s e d e p e n d e n c e s i n d u c e an evalu a t i o n o r d e r of the s t a t e v a r i a b l e s X~ out of 0 Fa = (v, u)

V = {0, 1,2, 3, 4, 5, 6,7, 8} U = {(0.2). (0,4), (0, 8), (4,1), (4,3), (8,51, (8,6), (5.7)}

Fig. 7. An arborescence is a tree, i.e., a connected graph with no cycle, such that there exists a vertex called a root from which any other vertex can be reached by a path going from the root to it. In an arborescence all verticcs but the root have exactly one predecessor.

(21.1) i¢0,j¢0}

U {(i, j ) ~ V × Vii = 0, j ~ 0}

(21.2)

E a c h arc ( i , j ) ~ ~ is a s s o c i a t e d with a w e i g h t bij t h r o u g h

bij =

(21.3)

bij is s i m p l y the l o g a r i t h m of the n u m b e r of v a l u e s b y w h i c h w e discretize the t e m p e r a t u r e l a p s e X j - X i of n o d e s j a n d i. A slight m o d i f i c a t i o n o f / " m a y be r e q u i r e d in o r d e r to t a k e into a c c o u n t a p o s s i b l y r e g u l a t e d n o d e [22 ]. T h e n the w e i g h t ~ ( / ~ ) o f an a r b o r e s c e n c e /~a is the l o g a r i t h m of the total n u m b e r of discrete s t a t e v a l u e s IEXi i n d u c e d b y / ' a . More precisely, we h a v e loglEWi -- qt(F~) + l o g i E r 1

(22)

Minimizing iEWI t h u s r e d u c e s s i m p l y to finding a m i n i m u m w e i g h t a r b o r e s c e n c e / ' a with

110 r o o t 0 s p a n n i n g all v e r t i c e s of F. G o o d alg o r i t h m s exist in o p e r a t i o n s r e s e a r c h literature to solve this p r o b l e m ; see, for instance, refs. 24 a n d 25. The m a j o r d r a w b a c k of this m e t h o d is t h a t it r e q u i r e s n 24- 1 additional p a r a m e t e r s to select. Unless one h a s at h a n d g o o d e s t i m a t e s for/~u and ~ru, this m e t h o d is h a r d to i m p l e m e n t . F o r the p u r p o s e of the p r e s e n t s t u d y we h a v e d e t e r m i n e d /~u and Orij USing a c o n v e n t i o n a l deterministic simulation. This w a y only par a m e t e r s a, fl and t~P had to be v a r i e d in o r d e r to find a g o o d discretization.

4. D i s c r e t i z a t i o n performance

influence

on

Our M a r k o v i a n a p p r o a c h h a s b e e n s u c c e s s fully i m p l e m e n t e d in a c o d e (SOLSTIS) w h i c h has allowed us to c a r r y out s e v e r a l n u m e r i c a l e x p e r i m e n t s in o r d e r to e v a l u a t e the influence of discretization on its p e r f o r m a n c e . T h r e e kinds of c o m p a r i s o n s w e r e m a d e . First, conv e n t i o n a l simulation r u n s h a v e b e e n p e r f o r m e d on the following t h r e e sets of input data: (i) half-hourly collected r e c o r d s of e x t e r n a l air t e m p e r a t u r e a n d s o l a r radiation; (ii) discretized v a l u e s of t h e s e r e c o r d s ; (iii) eleven s a m p l e realizations of s y n t h e t i c d a t a g e n e r a t e d b y a M a r k o v chain identified from these records. A deterministic s i m u l a t i o n (PASSIM) w a s e m p l o y e d for this p u r p o s e [3]. It u s e s the s a m e m o d e l d e s c r i p t i o n as t h a t of SOLSTIS so t h a t we could f o c u s o u r a t t e n t i o n on the influence of the driving v a r i a b l e s discretization and of the a s s u m p t i o n s a b o u t the driving p r o c e s s on the m o d e l behavior. Next, an i n t e r m o d e l c o m p a r i s o n w a s perf o r m e d using SOLSTIS a n d PASSIM, with particular a t t e n t i o n given to the s t a t e v a r i a b l e s ' discretization influence on the differences of their output. Finally, the results of s e v e r a l r u n s of SOLSTIS w e r e c h e c k e d for c o n s i s t e n c y in o r d e r to s t u d y the intrinsic reliability of the m e t h o d as a function of the s t a t e v a r i a b l e s discretization. T h o u g h the results of SOLSTIS are p r o b a bility distributions o f the different v a r i a b l e s of interest during a unit p e r i o d ( h e r e 30 m i n u t e s ) , only m e a n v a l u e s o v e r a m o n t h or a w h o l e winter s e a s o n , w h i c h are easily o b t a i n e d f r o m t h e s e distributions, h a v e b e e n u s e d to c o m p a r e

the p e r f o r m a n c e s o f SOLSTIS a n d PASSIM. An alternative w o u l d h a v e b e e n to c o m p a r e p r o b ability distributions, w h i c h c o u l d h a v e also b e e n p r o v i d e d b y PASSIM, b u t a m e a n v a l u e s c o m p a r i s o n s is m u c h easier a n d is e n o u g h for the p u r p o s e s o f the p r e s e n t work. F i g u r e s 8 a n d 9 and T a b l e s 1 and 2 give a detailed d e s c r i p t i o n of the two m o d e l s u s e d in this work. The values c h a r a c t e r i z i n g the discretization of a m b i e n t air t e m p e r a t u r e Te and of s o l a r radiation FGVS are given in T a b l e 3. This is the b e s t a m o n g several discretization p r o c e d u r e s that are r e p o r t e d in ref. 22. T a b l e 4 and Fig. 10 s h o w losses deviations c o m p a r i n g the continuous, the discretized a n d the g e n e r a t e d inp u t s o b t a i n e d using this discretization f r o m data g a t h e r e d on the winter s e a s o n 1 9 8 0 / 1 9 8 1 . T h e s e deviations are a b o u t 1% w h e n considering only the influence of the driving v a r i a b l e s ' discretization a n d a b o u t 10% in m e a n w h e n u s i n g g e n e r a t e d data. The m o r e sensitive p o i n t of the w h o l e m e t h o d is the state v a r i a b l e s discretization. N u m e r i c a l e x p e r i m e n t s h a v e s h o w n that it h a s considerable influence on the a c c u r a c y of e n e r g y flows a n d b a l a n c e s c o m p u t a t i o n , especially the auxiliary h e a t i n g p o w e r a n d e n e r g y flows bet w e e n h e a v y a n d l i g h t nodes. Attention h a s b e e n f o c u s e d on the e n e r g y b a l a n c e e q u a t i o n s at the s t a t i o n a r y point. Simply, t h e y require t h a t w h a t g o e s i n m u s t g o out, i.e., 467 -

-

FGVS

(1)

Pa

13

(3) Te

I/G 1 2 i

--BRICKS

1/G 23

INSULATION

Fig. 8. Direct gain test cell (STESO test cell). A threenodes model with two state variables (T1= T~ inside air temperature (°C), and T2= Tm.~ wall temperature (°C)), two driving variables (T3 = Te ambient air temperature (°C) and F6vs solar radiation CW/m2)) and one optional regulating driving variable (Pa back-up heating power (W)). Thermal capacities of state nodes are not represented to avoid overloading the picture. Their values are given in Table 1. Glazing index = 0.27, U= 0.4 W/m2K,MCp = 5.3 MJ/K.

111

..................... ~'q

~

H .... J (1)

~/-~

1/G12

(2) | J

~,t.-----/~ I

Fig. 9. Direct gain solar room (LESO solar unit). A fivenodes model with four state variables (7"1 inside air

temperature (°C), 7'2 indoor connected wall temperature (°C), T3 outdoor connected wall temperature (°C), and 7"4 south window (°C)); two driving variables (7"5ambient air temperature (°C) and FGvs solar radiation (W/m2)) and one optional regulating driving variable (P~ back-up heating power (W)). Thermal capacities of state nodes are not rcpresented to avoid overloading the picture. Their values are given in Table 2. Volume = 75 m3, floor area = 30 m2, opening area = 6 me or 9 mz, U-value (ext. walls) = 0.38 W/Km2, U-value (frames)= 1.9 W/Km2, U-value (glazing) = 3.1 W/KIn2, thermal capacity= 14.7 MJ/K, thermal loss coefficient=30 W/K or 36 W/K. TABLE 1. Thermal parameters characterizing the three nodes model of the STESO test cell

Symbol Designation C1

c~ G~2 G~z G23

Value

Thermal capacity (cell air) 26.5 (kJ/K) Thermal capacity (walls) 5.3 (MJ/K) Thermal conductance (node 1 to 2) 299.1 (W/K) Thermal conductance (node 1 to 3) 25.8 (W/K) Thermal conductance (node 2 to 3) 18.9 (W/K) Equivalent solar energy collecting 0.85 (m ~) surface (node 2)

i ~ ÷ S ~ -- ~ _ S ~

(23.1)

and IS~] = Io'~1, V i ~ V

(23.2)

+ UV °

where V is the set of all n o d e s ( V = V + U V - U V °)

of the

model

V ÷ is the set of all state n o d e s to w h i c h a h e a t i n g s o u r c e is applied (auxiliary heating, solar radiation, internal gains) V - is the set of external n o d e s ( c o r r e s p o n d i n g to driving variables) V ° is the set of t h o s e state n o d e s to w h i c h no h e a t i n g s o u r c e is applied directly 4)u is the e n e r g y flow f r o m n o d e i to n o d e j (of c o u r s e ~ j i ~- - ~ i j ) S~ is the s u m o f e n e r g y flows to or f r o m n o d e

o'i is the s u m o f e n e r g y flows c o m i n g f r o m a h e a t i n g s o u r c e to n o d e i (auxiliary heating, solar radiation, internal gains). In o u r model, eqn. (23.1) is always satisfied. However, for e v e r y state n o d e i ~ V + U V °, S~ and o-~ can be c o m p u t e d i n d e p e n d e n t l y , since, a p a r t f r o m the c o n t r i b u t i o n of auxiliary heating, ~ , i ~ V + U V °, d e p e n d s only on the distribution of the driving variables. So eqn. (23.2) m a y n o t hold. T h e n o u t p u t results have b e e n c h e c k e d to verify w h e t h e r this e q u a t i o n holds and, if not, to w h a t e x t e n t the differences b e t w e e n the r i g h t - h a n d side a n d the left-hand side d e p e n d o n the diseretization p a r a m e t e r s . Figure 11 s h o w s typical results o b t a i n e d for the STESO test cell. Here the r a n g e of variation of T2 has b e e n s u b d i v i d e d into a n u m b e r of intervals that varied f r o m 3 to 129 a n d the c o r r e s p o n d i n g results are c o m p a r e d to t h o s e o b t a i n e d b y c o n v e n t i o n a l deterministic simulation (PASSIM). It is interesting to p o i n t out that, in this case, e n e r g y flows to or f r o m external n o d e s (i.e., b e l o n g i n g to V - ) are a l m o s t insensitive to state variables diseretization a n d c o r r e s p o n d to t h o s e c o m p u t e d b y PASSIM. On the o t h e r hand, flows b e t w e e n state variables are highly d e p e n d e n t on the state variables diseretization: the m o r e intervals are used, the m o r e the results c o n v e r g e to t h o s e o f PASSIM. ~ l = Q , u x (auxiliary heating) and cr2=Qdt~ (solar radiation) have b e e n c o m p u t e d independently from Sl=4h3+4h~. and $ 2 = 4~2a-4h2. Since b o t h o'1 and $1 are f u n c t i o n s o f state variables a n d driving variables, t h e y are highly d e p e n d e n t o n the discretization of state variables a n d one c a n see that, in this case, ~ ~S~. On the o t h e r hand, o'2 d e p e n d s only on solar radiation while $2 still d e p e n d s on b o t h state a n d driving variables. That is

112 TABLE 2. Thermal parameters characterizing the five nodes model of the LESO solar unit Symbol

Designation

Value

C1 C2 C3 C4

Thermal capacities

Room air Indoor connected wall Outdoor connected wall South window

G~2 G~3 Gl4 G1~ Gea Ge4 G34 Gap G4s

Thermal conductances

Node Nodc Node Node Node Node Node Node Node

S~, Se2 S~3

Equivalent solar energy collecting surface

Node 1 Node 2 Node 3

1 1 1 1 2 2 3 3 4

to to to to to to to to to

90 (kJ/K) 13.13 (MJ/K) 1.465 (MJ/K) 39 (kJ/K)

2 3 4 5 3 4 4 5 5

109.0 35.9 15.5 14.9 72.9 21.1 3.6 6.2 24.9

(W/K) (W/K) (W/K) (W/K) CvV/K) (W/K) (W/K) (W/K) (W/K)

1.05 (m e) 2.11 (m e) 0.35 (m 2)

TABLE 3. Discrctization for driving variables Te (°C) and FGvs CvV/m2) Period

m

FGvs

H8081" October 1980 December 1980 February 1981

3 3 3 3

T,

H8081 October 1980 December 1980 February 1981

3 3 3 3

l0

v*

Ii

v~

12

v*3

13

0 0 0 0

25.19 20.83 14.51 27.12

304.7 272.0 275.0 289.0

456.6 398.8 421.0 429.8

609.3 544 550 578

702.1 684.7 668.3 688.1

914 816 825 867

-7.6 1.6 - 7.5 -7.3

-1.85 6.26 - 3.91 -4.02

2.13 8.27 - 1.77 -2.8

6.46 11.17 1.01 -0.46

11.87 14.93 3.97 1.7

14.77 17.25 5.52 3.15

21.6 21.6 9.7 6.2

*H8081 means winter season 1980/1981. TABLE 4. Thermal losses deviations for a STESO test cell without auxiliary heating. Eleven realizations of driving variables were generated H8081"

October 1980

December 1980

February 1981

Q~ (MJ)

1537.0

295

287.5

287.8

Qg (MJ) AQ~d (%)

1540.0 0.2

294.8 - 0.07

283.6 - 1.36

290.3 0.87

(Q~} (MJ) hV~g (%) aQ~ (MJ)

1541.0 0.27 78.47

265.9 -9.86 38.66

266.8 -7.21 45.56

285.5 -0.79 33.11

*H8081 means winter season 1980/1981. the reason why the state variables discretization h a s n o i n f l u e n c e o n or2 b u t c o n s i d e r a b l y a f f e c t s $2 a n d s o e q n . ( 2 3 . 2 ) d o e s n o t h o l d b u t f o r sufficiently fine discretizations. O n e c o u l d a r g u e t h a t a l l w e n e e d is t o i n c r e a s e t h e n u m b e r o f d i s c r e t e v a l u e s in o r d e r t o g e t

accurate results. To do so, one would simply use procedure C by .conveniently decreasing tiP ( n o t e t h a t a, /~ a n d tip a r e p o s i t i v e p a r a m e t e r s ) . T h e p o i n t is t h a t f r o m e q n s . ( 2 1 . 3 ) a n d (16.3) one gets [EXl ~

O(tiP- 1)

(24)

113 TABLE 5. Discretization of state variables. Procedure C with a = 0.5v~'~vs = 344.05 W/m e, no transient state elimination, February 1981, LESO solar unit

2O

15 10! [%]

5~

~

0

|

5-

a

-

XI,~

X~.p

~

IE X l

(°C)

(°C)

(°C)

(...)

1 21

32.42

9.12X10-3.& °

2 19.28 38.15

2.69X10-4.~P

3 18.4

2.16X10-a-& p

10Winter 80~81 ! • IL

Discrete Input: Relative deviation [%]

October 80 []

December 80

Generated Input: Mean Relative deviation [%]

[]

February 81

Generated Input: Standard Deviation [%1

Fig. 10. Relative deviations, mean relative deviations and relative standard deviations of thermal losses (STESO test cell without auxiliary heating).

36.13

1252 + 0.5 tiPI + 1 aP [70064_~0.5 &oJ+ 1

[

8202 -+0.5 ~PJ +1 ~

[569 ~ 5

4 10.73 24.04 23.4X10-3-& °

&°]+ 1

iMJI

25OO] [-] 1012

2000

10 II 1010

,,K

with transient states reduction

•

withoutIransJemstates reduction

108 10? 106 Oau×

Od,r

Qaux+Odir

Oloss

Flow 1.>3

Flow 2->3

Flow 1 >2

105

February1981

10

Fig. 11. Influence of the number of discretization intervals on energy flows and balances computation (STESO test cell without auxiliary heating).

100

2oo

3oo

4oo

5oo

6oo

7oo

8O0

p

9oo

~000 [w]

Fig. 12. The number of stochastic states ]E Wl as a function of parameter &P (same conditions as in Table 5; with transient state elimination we used /3= 2).

for every state node j. Hence, I E W I e O( ~ P - " )

(25)

w h i c h m e a n s t h a t t h e n u m b e r IEWI o f s t o chastic states increases as a power of n as p a r a m e t e r ~P t e n d s t o z e r o . F u r t h e r m o r e , I E wI grows as an exponential function of the number of state variables n. The following two examples clearly illustrate this situation. Procedure C with a= 0.5v~=344.05 W/m e (February 1981, LESO solar unit) has been used to discretize t h e s t a t e v a r i a b l e s . T h e first e x a m p l e is w i t h o u t transient state elimination. Table 5 gives the values characterizing the discretization as a f u n c t i o n o f p a r a m e t e r &P. F i g u r e 12 c l e a r l y s h o w s t h e b e h a v i o r o f IEWI a g a i n s t tiP. T h e case tiP=300 W (IEWI =888 300) has been r u n g i v i n g t h e r e s u l t s r e p o r t e d in T a b l e 6. C l e a r l y , t h i s d i s c r e t i z a t i o n is n o t fine e n o u g h b u t o n e c a n n o t i m p r o v e it s i m p l y b y r e d u c i n g tiP s i n c e IEWI is t o o h i g h f o r v a l u e s o f t h i s parameter lower than 300 W.

TABLE 6. Output results of SOLSTIS. The conditions are those of Table 5, 8P=300 W and Se=Se, +S~2+S~3 Qaux Qd~ Q~o,,= 1851 Qaux+ Qdw AQ = Qaux+ Q d l r - - Qlo~ AQ~= 100" IAQI/(Q~,= + Q~)

413.9 (MJ) 796.6 (MJ) 1966.1 (MJ) 1210.5 (ME) - 755.6 (MJ) 62.4%

i

o-~ (MJ)

S~ (MJ)

ASi~ (%)

1 2 3 4

Q,,u,,+Sel/Se'Q,~ =892.7 Se2/Se'Q,~r =238.4 Sea/S,.Q~r = 79.4 0

606.3 819.6 497.3 42.8

-32.1 243.8 526.8 -

It t u r n s o u t t h a t , f o r &o__ 3 0 0 W , a t l e a s t 6 4 4 5 7 7 s t o c h a s t i c s t a t e s (i.e., a b o u t 7 0 % o f a l l s t a t e s ) h a v e n u l l p r o b a b i l i t i e s s o in t h e second example we try the transient states elimination technique. Table 7 gives estimates o f cru a n d /zu p r o v i d e d b y c o n v e n t i o n a l d e t e r m i n i s t i c s i m u l a t i o n ( P A S S I M ) . IEWI has

114 TABLE 7. Estimates of oij , ~£ij and s u for the LESO solar unit i

j

(r(j

~j~

so

1 1 1 1 2 2 2 3 3 4

1 2 3 4 2 3 4 3 4 4

0.223 0.190 0.283 0.288 1.136 1.096 0.618 1.086 0.780 1.316

0 -0.46 - 1.39 -8.41 0 -0.93 - 7.95 0 - 7.02 0

0.989 0.861 0.981 0.173 1.103 0.918 -

TABLE 10. Output results of PASSIM and SOLSTIS for the same period, Feb. 1981, and the same model, LESO solar unit. Discretization conditions for SOLSTIS are those of Table 8

Qaux (MJ) Qd~ (MJ) S, (MJ) $2 (MJ) $3 (MJ) $4 (MJ) S~ (MJ) CPU (s)

TABLE 8. Discretization of state variables. Procedure C with transient state elimination: a = 0.5v~,Y~vs = 344.05 W/ m e, ~ = 1.5 5P=50 W, Feb. 1981, LESO solar unit Arc (i, j)

(0, (3, (1, (3,

~j~ (°C)

1) 0 2) 0.93 3) - 1 . 3 9 4) - 7 . 0 2

su

e~

e~u,

~

(°C)

(°C)

(°C)

(°C)

0.173 0.861 0.918

21 -0.26 -1.292 -0.612

32.42 0.456 26 0.26 0.01345 40 1.292 0.108 25 0.612 1.169 3

IEXl

TABLE 9. Output results of SOLSTIS. The conditions are those of Table 8 Q~x Qd~ Q,o~ = ]$51 Q~,~+ Qd~ AQ = Q~u~+ Qd~ -- Qlo~ AQ ~= 100- [AQ [/( Q ~ + Qd~)

1053.5 (MJ) 796.4 (MJ) 1855.6 (MJ) 1849.9 (MJ) - 5.7 (M J) 0.3%

i

a~ (MJ)

S~ (MJ)

~

(%)

1 2 3 4

Qau~+S,i/Se. Qdir=1291.8 Se2/Se'Qdlr = 478.7 S~3/S,'Qdir = 79.4 0

1289.3 444.7 117.7 3.9

-0.2 -7.1 48.2 -

b e e n r e d u c e d b y a f a c t o r o f a b o u t l 0 3 , as s h o w n in F i g . 12, s o w e c a n r u n t h e p r o g r a m w i t h a l o w e r v a l u e f o r 6P a n d w i t h a b o u t t h e same number of stochastic states (a= 344.05 W / m 2, f l = l . 5 , 6 P = 5 0 W, [E W]=702000). Table 8 shows the characterizing values of the discretization corresponding to the minimum weight arborescence found with these parame t e r s . T h e r e s u l t s a r e r e p o r t e d o n T a b l e 9. In T a b l e 10 t h e y a r e c o m p a r e d t o t h o s e o f P A S S I M . They are m u c h m o r e satisfying than t h o s e o b t a i n e d in t h e p r e v i o u s c a s e . N o w , a p a r t f r o m d e v i a t i o n s o n n o d e s 3 a n d 4, w h o s e s u m s o f

PASSIM

SOLSTIS

AQ

AQ~ (%)

1089.0 796.6 1328.0 447.3 77.3 0.06 - 1852.0

1053.5 796.4 1289.3 444.7 117.7 3.9 - 1855.6

35.5 0.2 38.7 2.6 -40.4 - 3.84 3.6

- 3.26 - 0.03 -2.92 -0.59 52.31

30.3

0.19

8298.7

e n e r g y flows are negligible with r e s p e c t to o v e r a l l t h e r m a l l o s s e s , all e n e r g y o u t p u t s o f SOLSTIS are v e r y close to the c o r r e s p o n d i n g quantities yielded by PASSIM.

5. S o m e practical aspects The Markovian approach has two interesting advantages over traditional deterministic simulation: (i) a s p r e v i o u s l y m e n t i o n e d , i n p u t d a t a ( w e a t h e r a n d u s e r b e h a v i o r d a t a ) a r e significantly reduced; (ii) e x e c u t i o n t i m e ( C P U t i m e ) is i n d e p e n d e n t of the length of the simulation period. H o w e v e r , c o m p u t e r r e s o u r c e s s u c h as m e m ory requirements and execution time are the major drawbacks of the method. The present implementation of SOLSTIS runs on a VAX8600 with the VAX/VMS virtual memory operating system. The p r o g r a m n e e d s about 12. IEWI b y t e s o f v i r t u a l m e m o r y p l u s a d d i t i o n a l 4 . ] E W] bytes of scratch disk space so that overall memory requirements for the examples pres e n t e d in t h i s p a p e r w e r e a b o u t 15 M b y t e s . O n t h e o t h e r h a n d , C P U t i m e is a l i n e a r f u n c t i o n o f m e m o r y s p a c e , p r o v i d e d it c a n b e k e p t in central m e m o r y . By e x t r a p o l a t i n g the CPU time r e s u l t s o f T a b l e 10, it c a n b e s e e n t h a t e x e c u t i o n t i m e o f o n e r u n o f S O L S T I S is e q u i v a l e n t t o that of 40 simulations of PASSIM during 'winter' s e a s o n s (7 m o n t h s f r o m b e g i n n i n g o f O c t o b e r till e n d o f A p r i l ) . As a c o n s e q u e n c e , t h e s e c o m p u t e r r e q u i r e ments significantly outweigh the advantages m e n t i o n e d above. H o w e v e r , the basic idea, w h i c h c o n s i s t s in c o n s i d e r i n g d r i v i n g v a r i a b l e s

115 in a s t o c h a s t i c f r a m e w o r k , is still a valid a n d i n t e r e s t i n g one. Figure 10, for e x a m p l e , clearly s u g g e s t s t h a t a h y b r i d a p p r o a c h c o n s i s t i n g in g e n e r a t i n g the t i m e e v o l u t i o n of the driving v a r i a b l e s as realizations o f an u n d e r l y i n g stochastic p r o c e s s a n d t h e n s u b m i t t i n g it to a traditional d e t e r m i n i s t i c s i m u l a t i o n p r o g r a m , c a n yield reliable results. This w o u l d drastically s a v e c o m p u t e r r e s o u r c e s since no discretization of state v a r i a b l e s is n o w n e e d e d . The s t o c h a s t i c m o d e l of driving v a r i a b l e s m a y be f u r t h e r refined a n d the m o d e l s b a s e d on M a r k o v chains and a u t o r e g r e s s i v e p r o c e s s e s (AR p r o c e s s e s ) t h a t are p r o p o s e d in ref. 22 for s o l a r radiation, a m b i e n t t e m p e r a t u r e , artificial lighting a n d r o o m o c c u p a n c y c o u l d b e used.

~ij

the set of all state n o d e s to w h i c h a h e a t i n g s o u r c e is a p p l i e d (auxiliary heating, s o l a r radiation, internal gains) the set of e x t e r n a l n o d e s ( c o r r e s p o n d ing to driving v a r i a b l e s ) the set of t h o s e s t a t e n o d e s to w h i c h no h e a t i n g s o u r c e is a p p l i e d directly e n e r g y flow f r o m n o d e i to n o d e j (of

Si

s u m of e n e r g y flows to or f r o m n o d e

Qp

s u m of e n e r g y flows c o m i n g f r o m a h e a t i n g s o u r c e to n o d e i (auxiliary heating, s o l a r radiation, internal gains) t h e r m a l l o s s e s c o r r e s p o n d i n g to continuous v a l u e d input t h e r m a l losses o b t a i n e d f r o m discrete v a l u e d input

V÷

V_ Vo

course

Qd

6. C o n c l u s i o n s

The n u m e r i c a l p r o p e r t i e s of the M a r k o v i a n a p p r o a c h h a v e b e e n i n v e s t i g a t e d in depth, t h u s i m p r o v i n g the u n d e r s t a n d i n g of this m e t h o d . The effect of the d i s c r e t i z a t i o n s i n v o l v e d h a s b e e n e x a m i n e d in detail t h r o u g h different numerical experiments. It h a s b e e n s h o w n t h a t reliable r e s u l t s c a n be o b t a i n e d p r o v i d e d the discretization of the m o d e l v a r i a b l e s is sufficiently fine. Efficient p r o c e d u r e s h a v e b e e n d e v i s e d to p e r f o r m this o p e r a t i o n . H o w e v e r , the a m o u n t of m e m o r y r e q u i r e d b y this m e t h o d to yield reliable r e s u l t s is v e r y large e v e n for small m o d e l s (15 Mbytes for a f o u r - i n t e r n a l - n o d e m o d e l ) a n d quickly g e t s larger for m o r e c o m p l e x m o d e l s a n d for increasing accuracy. The p r e s e n t w o r k s u g g e s t s t h a t o t h e r m e t h o d s a v o i d i n g s t a t e v a r i a b l e s ' discretization, s u c h as h y b r i d d e t e r m i n i s t i c - s t o c h a s t i c simulation or simplified s t o c h a s t i c t e c h n i q u e s , s h o u l d be p r e f e r a b l e to a M a r k o v i a n a p p r o a c h since t h e i r m e m o r y r e q u i r e m e n t s will n o t g r o w in an e x p l o s i v e m a n n e r a c c o r d i n g to a c c u r a c y or to the c o m p l e x i t y of the model. The perf o r m a n c e s of t h e s e a l t e r n a t i v e s will b e s t u d i e d in the future.

Nomenclature

V

the set of all n o d e s ( V = V ÷ U V - U V °)

of the

model

AQ'pd

~ji = - ~ i j )

_ IQ~-Qp x l 0 0 Q~

relative deviation of t h e r m a l l o s s e s Q~ with r e s p e c t to Q~ m e a n t h e r m a l l o s s e s o b t a i n e d f r o m realizations of g e n e r a t e d i n p u t AQ~g

= (Q~}-Q~

x l00

V~

oV.*

relative Q~ with standard obtained input

deviation of t h e r m a l l o s s e s r e s p e c t to Qp deviation of thermal losses f r o m realizations of g e n e r a t e d

References

1 E. Baumann, T. Baumgartner, J. Gass, P. Hartmann and H. Miihlebach, Users' influence on the energy consumption of 60 identical singlc family houses, Proc. Second Swiss Status-Seminar Energy Conservation in Buildings, Zurich, Switzerland, 1982,

EMPA, Diibendorf, p. 31 (in German). 2 J.-L. Scartezzini, A. Faist and Th. M. Liebling, Using Markovian stochastic modelling to predict energy performances and thermal comfort of passive solar systems, Energy Build., 10 (1987) 135-150. 3 N. Morel, PASSIM Version 4, Manuel d'utilisation, Laboratoire d'Energie Solaixe et de Physique du Bfitimcnt, EPFL, Lausanne, 1987. 4 Building Energy Analysis Group, DOE-2, BDL Summary, User's Guide, Lawrence Berkeley Laboratory, Berkeley, CA, 1979.

116

5 DEROB III, The Derob System, User's Manual, University of Texas, Austin, TX, 1979. 6 PASOLE:A GeneralSimulationProgramforPassive Solar Energy, Los Alamos, NM, 1978. 7 SERI-RES Version 1.0, Solar Energy Research Institute, Golden, CO, 1982. 8 J. D. Balcomb and R. D. Mac Farlan, A simple empirical method for estimating the performance of a passive solar heated building of the thermal storage wall type, Proc. Second National Passive Solar Conf., Philadelphia, PA, 1978, ASES, p. 377. 9 M. Lumsdaine and E. Lumsdainc, Simple design calculation procedure for passive solar houses, Proc. Third National ISES-AS /DOE Passive Solar Conf., San Jos~, CA, 1979, ASES. 10 D. Straub, Design prediction of performance of passive solar home types in Western Washington: analysis using a computer model, Proc. Solar 79-Northwest Conf., Seattle, WA, 1979, ASES, p. 65. 11 A. Faist and J. P. Eggimann, LESOSAI computer program: optimization of passive solar buildings, Proc. Fifth Swiss Symp. Research and Development in Solar Energy, Lausanne, Switzerland, 1985, EPFL, p. 497 (in French). 12 W. L. Glennie, Hand-held calculator aids for passive design, Proc. Third National ISES-AS /DOE Passive Solar CoW~., San Josd, CA, 1979, ASES, p. 433. 13 M.W. Funning and J. Duffy, Comparison of measured and predicted thermal performance of passive solar rcsidences, Proc. Sixth National ISES-AS / E T AL Passive Solar Conf., Oregon, 1981, ASES, p. 362. 14 P.O. Fanger, Thermal Comfort, R. E. Kricger, Malabar, U.S.A., 1982.

15 J.-L. Scartezzini, Application des mSthodes stochastiques ~ l'~tude des syst~mes de captage de l'6nergie solalre, ThAse de Doctorat No. 627, D~partement de Physique, EPFL, Lausanne, 1986. 16 J. Haslett, The analysis by stochastic modelling ot solar energy systems for space and water heating, Proc. ISES World Congress, New Delhi, India, 1976, ISES. 17 J. Haslett and F. P. Hand, Models in solar water heating -- a critical comparison, Eur. J. Ope'r. Res.,4 (1980) 1-6. 18 M. D. Devine, M. J. Kumin and A. A. Aly, A survey of optimization and stochastic process techniques applied to solar energy systems, in B. Lev (ed.), Energy Models and Studies, North-Holland, New York, 1983. 19 G. F. Lameiro, Stochastic modcls of solar energy systems, Ph.D. Dissertation, Colorado State University, Fort Collins, CO, 1977. 20 G. F. Lameiro and W. S. Duff, A Markov model of solar energy space and hot water heating systems, Solar Energy, 22 (1979) 211. 21 J. B. Kemeny and J. L. Snell, Finite Markov Chains, Springer, New York, 1976. 22 J.-L. Scartezzini, F. Bottazzi and M. Nyg~rd-Ferguson, Application des mdthodes stochastiques: dimensionnement et rdgulation, rapport ddtailld (1987-1989), Laboratoire d'Energie Solaire et de Physique du B~timent et Chaire de Recherche Op6rationnelle, EPFL, Lausanne, 1989. 23 P. Billingsley, Statistical methods in Markov chains, Ann. Math. Stat., 35 (1961) 12. 24 M. Gondran and M. Minoux, Graphes et Algorithmes, Eyrolles, Paris, 1979. 25 R. E. Tarjan, Finding optimum branchings, Networks, 7 (1977) 2 5 - 3 5 .