Annals of Nuclear Energy 27 (2000) 1543±1575 www.elsevier.com/locate/anucene
A variational principle in probabilistic dynamics P.E. Labeau * Service de MeÂtrologie NucleÂaire (CP 165/84), Universite Libre de Bruxelles, 50, Av. F.D. Roosevelt, B-1050 Brussels, Belgium Received 31 January 2000; accepted 5 February 2000
Abstract Probabilistic dynamics oers a theoretical framework for an integrated treatment of accident sequences. But such realistic modelling in risk analysis calls for the development of ecient numerical schemes, in order to bene®t from the conceptual advantages it provides. Up to now, mostly simulation techniques have been envisioned to deal with this kind of problem. However, a suciently exhaustive search for scenarios and the low probability of the situations of interest limit the eciency of these methods. This work, therefore, explores an alternative approach, based on a variationally processed correction of a reliability characteristic, ®rst estimated with an approximation of the probability ®eld. This ®rst estimation can result either from the solution of a simpli®ed problem or from the output of another solution scheme. A variational principle for the evaluation of a functional of the probability distribution in this dynamic context is established. The possible achievements as well as the limitations of this new approach are discussed, starting from the results obtained for two test-cases. # 2000 Elsevier Science Ltd. All rights reserved.
1. Introduction The theory of continuous event trees, or probabilistic dynamics (Devooght and Smidts, 1992), has arisen from the need of improving classical risk analysis methodologies when applied to dynamic systems. Such systems (Siu, 1994) are characterized by a time-evolving response to an initial solicitation. This requires a timedependent treatment of the accident scenarios, integrating the dynamic evolution of the system in the course of a transient within the branching process in the event tree, corresponding to transitions between hardware states. * Postdoctoral Researcher, National Fund for Scienti®c Research, Belgium. Tel.: +32-2-650-2061; fax: +32-2-650-4534. E-mail address:
[email protected] 0306-4549/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved. PII: S0306-4549(00)00015-3
1544
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
For realistic applications, however, this problem turns out to be highly dimensional. The probability density of ®nding the system in a given hardware state, with given values of the process variables describing the system evolution, as a function of time, is de®ned in a space made of replicas of Rn , if is the number of hardware states and n that of process variables. Even though evolution equations can be written for these probability density functions, their solution with usual numerical schemes is hopeless for industrial plants. Speci®c methods are then required, and many authors have resorted to simulation techniques to deal with this issue. Discrete dynamic event trees (DDETs) consist in performing a systematic simulation of ``all'' sequences of the event tree, while restricting the possible branching times to a discrete set of epochs. The issue of exhaustivity in the search for scenarios is dealt with, in principle, but a combinatorial explosion of branches can limit the potentialities of such an analysis. On the other hand, Monte Carlo (MC) simulation has also been envisioned to solve dynamic reliability problems, thanks to its inherent ¯exibility allowing the inclusion of realistic assumptions in the modelling and to its ability to treat large-size applications. Nonetheless, the probability of rare events has to be assessed, and analogue games then turn out to be inecient. They must thereby be discarded in favour of biased algorithms, whose eciency is often problem-dependent, however. A comprehensive survey of these simulation techniques can be found in (Labeau et al., in press), as well as numerous references therein. A recent trend in dynamic reliability consists in coupling methods to improve the eciency of the calculations. A hybrid approach between DDETs and MC simulation was for instance propounded in Marchand et al. (1998). In this work, a simpli®ed problem where only transitions on demand (i.e. corresponding to the actuation of control/protection devices) are allowed is ®rst solved with a DDET technique, which provides in this case an exact result. When transitions in operation are again permitted, the solution of this ®rst problem corresponds to the underlying deterministic sub-structure (``skeleton'') of the whole continuous event tree. New sequences grafted on to this skeleton are accounted for via MC simulation. The same idea can be exploited in dynamic reliability calculations based on a system representation using event sequence diagrams (ESDs). This kind of description of the system logic leads to a subset of the whole event tree, which can be considered as another form of skeleton. It could also be complemented by a random sampling of new sequences added to this ®rst solution (Smidts and Swaminathan, 1999). In the present work, we explore the same kind of philosophy, in a somewhat different way. A variational principle of the Lagrange type is introduced in order to improve the accuracy of the estimation of a reliability characteristic. The latter is expressed as a functional of the probability ®eld, an approximation of which can be obtained by an (inherently imperfect) numerical scheme, or even as the solution of a simpli®ed problem. The Lagrangian is de®ned by adding an error-compensating term to the functional. This corrective term should vanish if the exact probability distribution were used in the estimation of the functional. It should also be stationary with respect to variations in the probability ®eld. Both requirements will be introduced in de®ning the variational principle. It is important to notice that this variational method was previously exploited in Markovian systems involving stochastic neutron
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1545
distributions (Lewins, 1978). It was recently applied in Markovian reliability (Lewins et al., 1998) and extended to the pseudo-Markov case (Lewins and Chang, 1999). This paper is structured in the following manner. Section 2 presents the mathematical formulation of probabilistic dynamics, and how it can be included in a variational principle. The meaning of this approach and its scope of applications in the frame of dynamic reliability are then discussed. The application of the method is consequently considered as a substitute for MC simulation in the coupling with DDETs mentioned above. The accuracy of the variational approach is assessed in Section 3, where it is applied to a test-case for which an analytical solution exists. Our method is then compared with the coupling DDET-MC simulation on a second application presented in Section 4. Some ®nal remarks and perspectives conclude the work. For the sake of readability of the paper, several mathematical developments are written in the appendices. 2. A variational principle for dynamic reliability calculations 2.1. Mathematical framework of probabilistic dynamics We ®rst provide an overview of the mathematical theory of dynamic reliability (Devooght and Smidts, 1992), in which we underline those aspects that are necessary for the variational principle that is introduced later. As stated in the introduction, dynamic reliability gives allowance to the mutual interaction between the dynamic evolution of the system during a transient and the sequence of transitions (i.e. branchings) de®ning a sequence. Then let x 2 Rn be the set of process variables characterizing the system dynamics. The evolution of this vector is in¯uenced by the hardware state the system lies in. Hence we label the system dynamics by an index i, current system state: dx fi
x ; i 1 . . . dt
1
where is the total number of distinct con®gurations of the system. Reciprocally, the value of x can aect the probability of a transition between two states. This statement is obvious when the state change is due to the successful actuation of a control device, triggered by the crossing of a setpoint in the process variables space. But the probability of a component failure can be increased when severe working conditions are met. For these reasons, we explicitly write a dependence on x for the transition rates li
x out of state i, and pji
x between states j and i. The condition li
x
X pij
x
2
j6i
must be satis®ed. An accidental transient thus develops in the following way: once an initiating event has occurred, the system leaves the steady-state conditions it was in, and enters
1546
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
a new state i in which the process variables evolve according to Eq. (1). The system trajectory is thereby deterministic as long as no transition out of i takes place. But when the system switches from i to another state j, a new dynamics, corresponding to the new state j, is obeyed up to the next transition, etc. This behaviour is characteristic of a piecewise-deterministic stochastic process (Davis, 1993). In order to quantify it, we introduce the conditional p.d.f.
x; i; tjx o ; k; to to ®nd the system in state i at time t, with process variables x, provided it was in state k at to , with variables x o . If the transitions between states are Markovian, this conditional distribution obeys a Chapman-Kolmogorov (CK) equation (Devooght and Smidts, 1992): X @ divx fi
x : li
x
x; i; tjx o ; k; to pji
x
x; j; tjx o ; k; to
3 @t j6i with the initial condition:
x; i; to jx o ; k; to ik
x ÿ x o
4
This is a forward form of the CK equation, as the evolution of the distribution is considered with respect to the ®nal coordinates. In contrast, a backward CK equation can also be deduced: ÿ
X @ ÿ fk
x o :r xo lk
x o
x; i; tjx o ; k; to pkj
x o
x; i; tjx o ; j; to @to j6k
5
which turns out to be the adjoint of Eq. (3). The failure criteria correspond to the system entering an unacceptable state or escaping a safety domain D de®ned in the process variables space. The failure probability will then be assessed by making all unacceptable states absorbing and by adding the following boundary condition:
x s ; i; tjx o ; k; to 0
8x s 2 ÿ
D
s:t:
n:fi
x s < 0
6
where ÿ is the border of D and n is the exterior normal to D. Let us introduce now the following matrix notations (Labeau and Lewins, 1999): . P; ; s:t:Pik
x; tjx o ; to
x; i; tjx o ; k; to . D; ; s:t:D
x diag div f1
x : ; div f2
x : ; . . . ; div f
x : . ÿ; ; s:t:
ÿii
x ÿij
x
ÿli
x ; pji
x ;
i 1... i 6 j 1 . . .
These matrices can be used to rewrite Eqs. (3) and (5) in a compact form, in order to simplify the writing of the Lagrangian (see Section 2.3). Indeed, we easily ®nd that
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
@ D
x ÿ ÿ
x P
x; tjx o ; to 0 @t
for the forward CK equation, while the backward CK equation becomes @ D
x o ÿ ÿT
x o PT
x; tjx o ; to 0 ÿ @to
1547
7
8
where ``T '' and `` '' stand for the transpose and the adjoint, respectively. 2.2. Expression of reliability characteristics The distribution
x; i; tjx o ; k; to provides a huge quantity of information to the risk analyst. However, there is no speci®c interest in knowing its detailed form deep inside the safety domain. Meaningful characteristics of the severe transients we want to analyse in dynamic reliability are e.g. the probability associated with each failure mode, the mean time to escape D or the distribution of the ``damage'' (whatever the physical reality beyond this term) as a function of time. Let c
t denote such a key characteristic, which is expressed as a functional of the probability ®eld: c
t h T
Rn
P
x; tjx o ; to
x o ; to dx o
9
where h is a vector operator and
x o ; to is a vector, the ith component of which represents the initial probability distribution of the process variables in state i. We assume for the sake of simplicity that
x; to
to
x ÿ x o
10
In the latter equation
to is the initial probability vector, when the system leaves its steady-state conditions, following the occurrence of the initiating event. If one is interested in the failure probability due to the crossing of the safety border, the components of operator h will all take the same form:
dx; i 1 . . .
11 hi Rn =D
2.3. De®nition of the Lagrangian As mentioned in the introduction, the solution of the system of PDEs (3) is a dif®cult numerical problem, due to its dimensionality. Any solution method will lead to an approximate probability ®eld P~ , and consequently to an approximation c~ of the characteristic of interest.
1548
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
The propounded variational procedure consists in de®ning a Lagrangian L
t, which adds to the estimation c~ of the functional c a supplementary term in order to compensate for the error on the matrix of conditional probabilities. L
t is written as (Lewins, 1978; Lewins et al., 1998; Labeau and Lewins, 1999)
t
Q
x; tjx 0 ; L
t c~
t h T n to R
12 @ ~ 0 0 0 P
x ; jx o ; to
to dx 0 d ÿ
x ÿ D
x ÿ @t Several important comments can be made on this expression. First, it can be observed that the corrective term in Eq. (12) vanishes if the exact solution is known, i.e. if P~ P, since the term between brackets is then identically zero. Secondly, the supplementary term in L
t is truly error-compensating, because Eq. (7) appears with a negative sign in the integral. Finally, a still to be determined weighting matrix Q
x; tjx 0 ; has been introduced in Eq. (12). This allows us to account for a supplementary condition, which is the stationarity of L
t to variations in P. Therefore, without writing explicitly the tilde symbol, we have:
t
Q
x; tjx 0 ; L
P h T P
x; tjx o ; to
to h T to R n
13 @ P
x 0 ; jx o ; to
to dx 0 d 0 ÿ
x 0 ÿ D
x 0 ÿ @ The last term in Eq. (13) can be integrated by parts. Furthermore, having P
x 0 ; to jx o ; to 0 because of Eq. (4) and setting Q
x; tjx 0 ; t I
x ÿ x 0 by assumption, Eq. (13) is satis®ed if the following condition on Q holds:
ÿT
x 0 ÿ D
x 0
@ QT
x; tjx 0 ; 0 @
14
Eq. (14) is nothing but the adjoint of Eq. (6), and consequently Q
x; tjx 0 ; P
x; tjx 0 ; . This de®nes a pseudo self-adjoint variational principle (Lewins et al., 1998) for the estimation of dynamic reliability functionals. 2.4. Applicability The object of the variational principle we have just presented is to improve the accuracy of the estimation of a reliability functional. Though it has not been demonstrated, numerical results in Markovian reliability support the idea that the standard deviation of the variationally-processed estimation evolves as N1 , instead of p1 in MC , if N is the number of histories (or trial functions) (Lewins and Chang, N
1999). This additional accuracy comes from the more complete use of the information contained in the conditional distribution.
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1549
However, this potential gain of the method should not be counterbalanced by the supplementary workload necessary to compute the error-compensating term. In the peculiar case of dynamic reliability, this situation is not obvious at all : Eq. (12) shows the convolution nature, both on time and process variables, of the term that needs to be postprocessed. This task appears even more numerically demanding when one considers the dierential operators (acting both on t and x) in the expression to be integrated. This issue can be partly dealt with in the following case. Consider that the matrix of transition rates can be split as ÿ
x ÿ~
x ÿ
x
15
Let the approximate probability ®eld be the exact solution of the problem de®ned by matrix ÿ~. In this case, the error-compensating term reduces to L
t ÿ c~
t h T
t
to
Rn
P~
x; tjx 0 ; ÿ
x 0 :P~
x 0 ; jx o ; to
to dx 0 d
16
thanks to Eqs. (14) and (15). This is a perturbation approach, which allows the consideration of the use of higher order theory (Lewins, in press). In Eq. (16), no dierential operator has to be considered in the postprocessing stage. This idea was ®rst exploited on a simple two-state one-dimensional test-case (Labeau and Lewins, 1999): the solution of the non-repairable case was used as approximate solution for the repairable case. Yet the possibility of repair is not clear in severe transients. Moreover, the application of this approach is still strongly limited when the dimension of x increases. Therefore, we consider in the following a more realistic partition of ÿ
x [see Eq. (15)] between the transitions on demand and those in operation. As discussed in the introduction, the problem based only on transitions on demand has a deterministic solution [which is the so-called skeleton (Marchand et al., 1998)]. The associated probability ®eld is then expressed as a linear combination of Dirac peaks, the integration of which in Eq. (16) is much simpler. Consequently, the variational principle can be envisioned as a substitute for MC in estimating the perturbation of the skeleton due to transitions in operation. Mathematically speaking, this decomposition can be understood as follows. Let P S KP
17
be a symbolic writing for the integral form of Eq. (7) [see (B.1)]. Its solution can be expanded into a Neumann series: P
I ÿ Kÿ1 S S KS K2 S . . .
18
If K KD KS , where KD is the integral operator associated with the problem with transitions on demand only, and KS the additional operator when allowance is given to the stochastic transitions, the approximate probability ®eld is
1550
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
P~
I ÿ KD ÿ1 S
19
Multiplying Eq. (17) by
I ÿ KD ÿ1 , we ®nd: or
I ÿ KD ÿ1
I ÿ KP P~
20
I ÿ
I ÿ KD ÿ1 KS P P~
21
Therefore, we can write: 2 P P~
I ÿ KD ÿ1 KS P~
I ÿ KD ÿ1 KS P~ . . .
22
The computation of the Lagrangian amounts to accounting for the second term in the Neumann series on the r.h.s. of Eq. (22). In the next section, we present a problem for which an analytical solution is available, so that the potential capabilities of the variational principle can be assessed. A second application will then compare our method to the hybrid DDET-MC approach. 3. Application 1: an analytical benchmark 3.1. Description of the problem Consider a system described by one process variable x, with its steady-state value xo ; at time to , a transient is initiated in state 1, and the evolution of x follows an exponential law. When the level x ` is reached, a protection device is solicited; it can either fail (with a probability po ), work correctly (prob. 1 ÿ po ÿ p1 , state 2) or be imperfectly triggered (prob. p1 , state 3). In the last two cases, x starts to decrease, and a safe shutdown is reached as soon as x < xo . But if x > L, the system fails. This situation could correspond to a simple modelling of the pressurisation of a tank, with a relief valve designed to open at the setpoint ` before the pressure of rupture L is reached. We add to this description the possibility of a supplementary component failure accelerating the transient (transition from state i to i 3, i 1 . . . 3). This more severe transient can still be mitigated in case of perfect working of the protection device, but it is only slowed down in case of partial triggering. This evolution of x is then as follows: dx dt
ai x; i 1; 4; 6 ;a > 0 ÿai x; i 2; 3; 5 i
8i
23
with a2 > a3 , a4 a1 , a5 a2 ÿ , a6 ÿa3 , and > 0. Some of the possible transients are represented in Fig. 1, while the state graph is shown in Fig. 2,
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1551
Fig. 1. Transient evolutions for application 1.
where solid edges are used for stochastic transitions and dashed edges for transitions on demand. The transition rate is assumed constant on [xo ; L], and zero outside this range. This makes both ®nal situations (failure or safe shutdown) absorbing. The estimation of the respective probabilities PL
t and Po
t of these events is the objective of the calculations. 3.2. Calculation of the Lagrangian When only transitions on demand are considered, the matrix of conditional p.d.f.'s has a block-diagonal form: 0 1 0 0 Pii 1 ~ 0 where P~ i @ P Pi1;i1 0 A; i 1; 4
24 P~ P i1;i 0 P~ 4 0 P P i2;i
i2;i2
It is easily checked that the elements of P~ have the following expression: ÿ P~ ii
x; tjx0 ; b
xjx0 : x ÿ x0 eai
tÿ ; i 1; 4
25
8 <1 with b
xjx0 po : 0
26
if x0 < x < ` or if x0 < `4x else
`4x0 < x
ÿ P~ ii
x; tjx0 ; x ÿ x0 eÿai
tÿ ; i 2; 3; 5
27
1552
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
Fig. 2. State graph for application 1.
ÿ P~ 66
x; tjx0 ; x ÿ x0 ea6
tÿ
28
P~ i1;i
x; tjx0 ; H
` ÿ x0 H
` ÿ x
1 ÿ po ÿ p1 ai1 =ai ! ` ÿai1
tÿ x ÿ `e ; i 1; 4 x0
29
P~ 31
x; tjx ; H
` ÿ x H
` ÿ xp1 : x ÿ `e
ÿa3
tÿ
P~ 64
x; tjx ; H
` ÿ x H
x ÿ `p1 : x ÿ `e
a6
tÿ
0
0
0
0
a3 =a1 ! ` 0 x
0 a6 =a4 ! x `
Matrix ÿ
x has also a block-diagonal structure: ÿÿs 0 :H
x ÿ xo :H
L ÿ x ÿ
x ÿs 0
30
31
32
with ÿs diag
l; l; l. Finally, we note that
to T
1; 0; 0; 0; 0; 0. Using results (24) and (32), one can directly ®nd that Eq. (16) becomes L
t c~
t h T :A
x; tjxo ; to c~
t h T
t to
d dx0 a
x; t; x0 ; ; xo ; to :H
x0 ÿ xo :H
L ÿ x0
33
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1553
with a
x; t; x0 ; ; xo ; to 3 2 ÿlP~ 11
x; tjx0 ; :P~ 11
x0 ; jxo ; to 6 ÿlP~
x; tjx0 ; :P~
x0 ; jx ; t ÿ lP~
x; tjx0 ; P~
x0 ; jx ; t 7 6 21 11 o o 22 21 o o 7 7 6 6 ÿlP~ 31
x; tjx0 ; :P~ 11
x0 ; jxo ; to ÿ lP~ 33
x; tjx0 ; P~ 31
x0 ; jxo ; to 7 7 6 6 7 7 6 lP~ 44
x; tjx0 ; :P~ 11
x0 ; jxo ; to 7 6 4 lP~ 54
x; tjx0 ; :P~ 11
x0 ; jxo ; to lP~ 55
x; tjx0 ; P~ 21
x0 ; jxo ; to 5 lP~ 64
x; tjx0 ; :P~ 11
x0 ; jxo ; to lP~ 66
x; tjx0 ; P~ 31
x0 ; jxo ; to
34
The explicit form of A
x; tjxo ; to , which is easily obtained thanks to the Dirac peaks in Eqs. (25)±(31), is given in Appendix A. We are interested in the estimation of two functionals: the failure probability PL
t and that of a safe shutdown Po
t. The vector operators associated to them are:
1 dx; i 1 . . . 6
35a hL;i L
ho;i
xo o
dx; i 1 . . . 6
35b
respectively. After some lengthy but trivial calculations, we obtain the expression of the Lagrangian associated with PL
t: ÿ LL
t po H xo ea1
tÿto ÿ L ÿ ÿ lpo L a4
t ÿ to ÿ `n :H xo ea4
tÿto ÿ L :H L ÿ xo ea1
tÿto a4 ÿ a1 xo " x a6 =a4 la4 p1 L a6 ` o a6
t ÿ to ÿ `n ÿ `n ÿL :H `ea6
tÿto a6
a4 ÿ a1 ` ` a4 xo # a6 =a1 a6 =a1 1 1 ` a6
tÿto xo a6
tÿto xo ÿ ÿL a6 :`n :H `e :H L ÿ `e ` ` a1 a4 xo " x a6 =a1 lp1 L a6 ` o : a6
t ÿ to ÿ `n ÿ `n ÿL :H `ea6
tÿto a3 a6 ` ` a1 xo 1 0 1 1 a6 a3 a6 ` B a6
tÿto xo a a 1 3 C `n H@L ÿ xo e A ` a3 xo 0
x a6 o B :H@xo ea6
tÿto `
1 1 1 # a1 a3 ÿLC A
36
1554
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
As for the Lagrangian associated with Po
t, we have: a a
tÿt o 1 2 l ` a1 a2 Lo
t
1 ÿ po ÿ p1 1 ÿ `n ÿ` :H xo e a2 xo a a
tÿt o 1 3 l ` l
1 ÿ po ÿ p1 p1 1 ÿ `n :H xo e a1 a3 ÿ ` a3 xo a4 ÿ a1 " a5 =a4 a4 a5 ` a5
tÿto xo `n ÿ` a4
t ÿ to ÿ :H xo e a5 ` xo # a5 =a1 a ÿ a a5 =a1 ` 4 1 a5
tÿto xo a5
tÿto xo `n :H xo e ÿ` :H ` ÿ xo e ` a1 ` xo " a a
tÿt o 1 2 l
1 ÿ po ÿ p1 1 1 ` a1 a2 ÿ` `n :H xo e a5
t ÿ to ÿ a5 a2 ÿ a5 a1 a2 xo # a5 =a1 a ÿ a a5 =a1 ` 2 5 a5
tÿto xo a5
tÿto xo :`n :H xo e ÿ` :H ` ÿ xo e ` a5 ` xo
37
3.3. Comparison with the analytical solution To assess the quality of the variationally-processed estimations (36) and (37), we shall now compare them with the analytical results which can be obtained for this problem (see Appendix B). Table 1 gathers the numerical values given to the parameters of the application: Table 1 Numerical values of the parameters xo ` L po p1
1. 3. 4. 0.02 0.04
a1 a2 a3 a4 a5 a6
0.2 0.25 0.1 0.35 0.1 0.05
The transition rate l has not yet been speci®ed. A characteristic time of the transient is 1 ` 1 `
38 tc `n `n a 1 xo a 2 x o which is the time elapsed until the accident is mitigated in the case of correct working of the protection device. Indeed, the ®rst term on the r.h.s. of (38) is the time to reach the control setpoint in state 1, while the second term corresponds to the time interval required for x to go below xo while evolving in state 2. Hence ltc is a measure of the probability of having a stochastic transition in the course of the transient. We will thereby analyse the numerical results with respect to an adimensional
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1555
parameter l=
a1 po , which gives the relative importance of the transitions in operation, compared with the transitions on demand. Figs. 3±7 present the evolution with time of the failure probability, obtained both analytically and variationally. These are compared with the result of the simpli®ed
Fig. 3. Application 1 Ð comparison of the exact failure probability with the variationally processed estimation
0:2.
Fig. 4. Application 1 Ð comparison of the exact failure probability with the variationally processed estimation
0:5.
1556
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
Fig. 5. Application 1 Ð comparison of the exact failure probability with the variationally processed estimation
1.
Fig. 6. Application 1 Ð comparison of the exact failure probability with the variationally processed estimation
2.
problem, referred to as the approximate solution. There is an obvious agreement between both results, even for rather large values of , i.e. for situations where the dierence between the exact and simpli®ed problems becomes quantitatively important. The variational method turns out to be inaccurate only for 10, but in
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1557
this case, the failure probability at the end of the transient is twice as large as that of the simpli®ed case without transitions in operation. However, in all cases, one can clearly observe the contribution to the failure risk due to an incorrect working of the protection device after a ®rst transition in operation.
Fig. 7. Application 1 Ð comparison of the exact failure probability with the variationally processed estimation
10.
Fig. 8. Application 1 Ð comparison of the exact shutdown probability with the variationally processed estimation
0:2.
1558
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
When considering the estimations of the shutdown probability (Figs. 8±10), we quickly face, for increasing values of , one of the limitations of our method. If the variational calculation of Po
t is still acceptable for limited values of , it can be seen that the accuracy of the results decreases much quicker with than in the previous case. Moreover, we observe for large that some estimations of Po
t exceed 1! This
Fig. 9. Application 1 Ð comparison of the exact shutdown probability with the variationally processed estimation
0:5.
Fig. 10. Application 1 Ð comparison of the exact shutdown probability with the variationally processed estimation
ÿ 10.
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1559
highlights the ``bias'' of the Lagrangian for large errors in the trial function, since L
t is second order accurate. Hence, the approximation of the functional has to be ``suciently good'' before applying the variational approach. However, this problem appears for the estimation of a large probability, which is somehow out of the scope of this work. 4. Application 2: a comparison with the hybrid DDET-MC approach 4.1. Description of the problem Let us bring some modi®cations to the case we considered in the previous section. As before, the transient evolution is ®rst exponential. When x `1 , a ®rst control device (c1 ) is solicitated. Again, it can either fail, work or be incorrectly actuated, with respective probabilities po , 1 ÿ po ÿ p1 and p1 . But this time, we assume in the latter case that the transient is only slowed down by the partial operation of c1 . For this reason, a second control unit c2 is devised to mitigate the transient, whatever the working mode of c1 . This second protection device comes into play with a probability 1 ÿ q when crossing the threshold x `2 . The system failure is again de®ned by x exceeding L. To this set of transitions on demand are grafted two possible transitions in operation: . the failure of an additional component Q, with a rate lq , making the transient more severe (term x in the dynamics); . the spurious actuation of c2 , with a rate l2 . Both transition rates vanish outside xo ; L. The eect of the failure of Q on the whole system is the following : none of the control units alone can mitigate the transient when Q fails, but their joint action, even in case of imperfect actuation of c1 , is sucient to bring the system back to a safe situation. Table 2 displays the various combinations of components states and the resulting dynamics, while Figs. 11 and 12 show some of the possible accident scenarios. The state graph of the system is presented in Fig. 13. The objective of the calculations is the estimation of PL
t, probability of a system failure. 4.2. Calculation of the Lagrangian In this problem again, the matrix of conditional p.d.f.'s P~ , when no transitions in operation are considered, has a block-diagonal form: P~ with
P~ 1 O
O P~ 7
39
1560
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
0
P11 BP B 21 B B P31 P~ 1 B BP B 41 B @ O
O P22
O O
O O O
P61
O
0
P77
B P87 B B B P97 7 ~ P B BP B 10;7 B @ P11;7 P12;7
O O
O O
P33
O
O
O O
P44 O
O P55
P63
O
O
1 O O C C C O C C and O C C C O A P66
O
O
O
O
P88
O
O
O
O O
P99 O
O P10;10
O O
P11;8 O
O P12;9
P11;10 P12;10
P11;11 O
O
40
1
O C C C O C C O C C C O A P12;12
For the sake of readability of the paper, the explicit expression of the components of matrix P~ is not provided. It is worth underlying the fact that P~ contains much more information than that pertaining to the simpli®ed problem alone. Indeed, when the stochastic transitions are not accounted for, only 5 states can physically be reached (see Table 2). This situation corresponds to the ®rst column of P~ 1 . However, in the application of the variational principle, an information relative to all system states in the exact problem must be available. The approximate solution P~ is thereby obtained by completing the solution of the simpli®ed problem with the control-based evolution from all starting conditions. In spite of this, the additional workload is not too important, Table 2 System states and associated dynamics Component states
System states
Dynamics
Q
c1
c2
i
fi
x
1 1 1 1 1 1 2 2 2 2 2 2
1 2 3 1 2 3 1 2 3 1 2 3
1 1 1 2 2 2 1 1 1 2 2 2
1 2 3 4 5a 6 7a 8a 9a 10a 11a 12a
1 x a 1 x
1 ÿ 2 x ÿa2 x
1 ÿ 3 x a3 x
1 ÿ 4 x ÿa4 x
1 ÿ 2 ÿ 4 x ÿa5 x
1 ÿ 2 ÿ 4 x ÿa6 x
1 x a7 x
1 ÿ 2 x a8 x
1 ÿ 3 x a9 x
1 ÿ 4 x a10 x
1 ÿ 2 ÿ 4 x ÿa11 x
1 ÿ 3 ÿ 4 x ÿa12 x
a
State not considered in the approximate solution, ai , j , > 0, i 1 . . . 12, j 1 . . . 4.
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1561
Fig. 11. Transient evolutions (a) for application 2.
Fig. 12. Transient evolutions (b) for application 2.
as the deduction of the diagonal terms in P~ is straightforward and because the matrix P~ is sparse. Consider now the matrix ÿ
x. Denoting by In the n n identity matrix, we have
1562
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
Fig. 13. State graph for application 2.
ÿ
x
ÿ1 lq I 6
o :H
x ÿ xo :H
L ÿ x ÿ7
41
with ÿ ÿ l2 lq I 3 ÿ l2 I 3 1
O ÿlq I3
and
7
ÿ
ÿl2 I3 l2 I 3
O O
42
Being interested in the Lagrangian LL
t associated with the failure probability, we need only to integrate the components of vector P~ :ÿ:P~ :
to [see Eq. (16)] corresponding to critical states, i.e. states i such that fi
L > 0. The expression of the Lagrangian is given in Appendix C. Though this ®nal result is rather cumbersome, LL
t has to be understood as the sum of the (approximate) contributions to the failure probability due to all additional (with respect to the simpli®ed problem) accident scenarios accounted for through the variational procedure [see Eq. (C.2) for the list of such sequences].
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1563
4.3. Results and comments The estimation of PL
t obtained with the Lagrangian (C.4) has been compared with the results given by the hybrid DDET-MC approach, and those without variational processing. In this method, the deterministic branches of the event tree, articulated on the control thresholds, are ®rst calculated. Then, situations departing from this skeleton are systematically considered, by forcing stochastic transitions to be sampled on each section of this deterministic event tree. For this simple test-case, this method provides a reference solution for the comparison. The numerical values given to the parameters can be found in Table 3. Figs.14 and 15 display the results obtained for two values of the adimensional parameter (see Section 3.3), which was de®ned for this problem with the transition rate lq . A satisfying agreement between both estimations can be observed, especially for the smallest value of , unsurprisingly. However, a growing gap between these two curves is noticed on the last part of the accident duration. What are the possible reasons for such a divergent behaviour? A ®rst idea that can be investigated is the lack of exhaustiveness in the accident sequences accounted for in the variational procedure. The Langrangian indeed gives Table 3 Numerical values of the parameters `1 `2 L po p1 q l2
3 3.5 4 0.02 0.04 0.01 lq =2
xo 1 2 3 4
Fig. 14. Application 2 Ð comparison between the variational and MC estimations
0:05.
1 0.2 0.3 0.1 0.35 0.2
1564
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
Fig. 15. Application 2 Ð comparison between the variational and MC estimations
0:5.
allowance to new scenarios corresponding to a ®rst-order improvement of the approximate solution. In this context, it means that transient evolutions combining two sections of the skeleton through one stochastic transition are included in the analysis. Yet, from the problem description, it is clear that a system failure can occur when both possible transitions in operation take place, together with the failure of c1 . However, the likelihood of these stochastic transitions is limited in this example. Moreover, the failure of Q and the spurious triggering of c2 drive the system in opposite directions, relatively to safety in the system operation. Indeed, the ®rst failure event makes the transient more severe, while the second one speeds up the transient mitigation. Therefore, the omission of the combined occurrence of these two events in the same sequence induces two small errors that are likely to compensate for each other. Another source of error in the application of the variational procedure lies in the contribution to LL
t from physically unrealistic sequences generated in calculating the Lagrangian. Consider in Eq. (C.2) the second scenario leading to the presence of the system in state 10 (second term in the 10th component of z): transitions on demand bring the system into state 4, what corresponds to the failure of c1 and the successful actuation of c2 ; while x is decreasing, a stochastic failure of component Q occurs, and the value of the process variable starts increasing again, without any possible means of mitigation. Actually, this description corresponds to the realistic development of the sequence, as analysed with the DDET-MC solution scheme. In the variational procedure, if the failure of Q occurs when x < `1 , P~ 10;10 will include the probability of c1 not responding to the solicitation at x `1 ; indeed, the de®nition of the approximate solution is done on the basis of the possible transitions on demand alone, independently of the constraints introduced by the scenario having developed so far. The Lagrangian hence under-estimates the probability of a system failure, because this sequence accounts for an additional unrealistic control failure.
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1565
If we now correct the expression of P~ 10;10 to forbid the second solicitation of c1 , the distance between the DDET-MC method and the variational approach is dramatically reduced. This can be observed in Figs. 14 and 15, where the results of this so-called modi®ed variational method are also displayed. With this modi®cation of the approximate solution of P~ , the agreement with the reference solution is obviously better. 5. Conclusions This paper has dealt with the application of a variational principle to the calculation of functionals in dynamic reliability. A Lagrangian has been associated with the characteristic of interest, the error-compensating term making use of the solution of the adjoint CK equation to improve the accuracy of the calculations. In practice however, this additional accuracy can be useless if it is obtained only at the cost of a large supplementary numerical workload. At ®rst, this seems to be the case in dynamic reliability, because a convolution product on the process variables and time has to be performed, after having applied dierential operators to the ®rst approximation of the probability ®eld. Nonetheless, this diculty can be, at least partly, overcome if the approximate p.d.f.'s are solution of a simpli®ed problem in which only transitions on demand are accounted for. In this case indeed, the expression to be integrated reduces to a matrix product, without dierential operators, and the p.d.f.'s are linear combinations of Dirac peaks, the integration of which is straightforward. This approach was applied to two test-cases, one of which having a fully analytical solution and the other one being also solved by a hybrid DDET-MC method. The variational scheme turned out to provide very good approximations of the results, as long as the eect of transitions in operation was kept limited in comparison with that of the transitions on demand. The biasedness of the Lagrangian was also highlighted in extreme situations. It is expected that problems of a larger size could be eciently treated by resorting to symbolic calculus. However, the capability of this variational method to compete with the hybrid DDET-MC approach is still questionable for realistic cases. In future work, the following approaches should be further investigated, in order to reduce the cost of postprocessing. First, as the most interesting functionals refer to an exit out of the safety domain (failure probability, damage...), one should account only for the contribution of the error-compensating term in critical states. These are the states i for which n:fi
x s > 0, if x s 2 @D, border of the safety domain, and n is the exterior normal to D. Secondly, when MC simulation is used to obtain the approximate probability ®eld, it could be interesting to estimate directly the error-compensating term as an objective of the MC game, instead of postprocessing the MC trial functions via a classical integration scheme. Finally, we would like to comment on a possible by-product of the variational developments presented in this paper, which could give a basis for improving simulation schemes in dynamic reliability. This results from the expression of the
1566
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
Lagrangian as a matrix product, and from the subsequent combination of scenarios that is encompassed in this formulation. Indeed, when estimating L
t, we ®rst compute @ ~ 0
43 P
x ; jx o ; to
to P~
x; tjx 0 ; ÿ
x 0 ÿ D
x 0 ÿ @ and then integrate it over x 0 and , before applying the vector operator h to the result [see Eqs. (12) and (14)]. The contribution in state i is linked to the ith component of , whose explicit form writes XXX XX @ ~ Pjk k
to P~ ij ÿjk P~ k` `
to ÿ P~ ij Djj
44 i @ j j ` k k if we omit most dependences on x and t for the sake of clarity. What is the eect of the variational procedure? It combines a scenario (identi®ed with the approximate solution method) leading the system from state ` to state k, with a transition between k and j, and ®nally with another possible scenario (highlighted by the approximate solution scheme) from j to i. The same kind of reasoning can be done for the second term on the r.h.s. of Eq. (44). The non-vanishing terms in the sums (44) provide thereby the new sequences reaching a state i, that are added by the variational procedure as a ®rst-order improvement of the approximate solution. The identi®cation of the additional sequences can be performed for all critical states, while skipping the time-consuming estimation of the Lagrangian. Instead, this information could be exploited to bias a MC algorithm towards these new sequences, and to estimate the quantitative impact of these scenarios on the total failure risk or damage. This would allow to use the variational approach in order to bene®t from the advantages of MC simulation in high-dimensional applications. Acknowledgements The author is indebted to Dr. Jeery Lewins, University of Cambridge (UK), for ongoing support and advice throughout this work. Mrs. JoseÂe Immers is also warmly thanked for her help in typing this paper. Appendix A We give below the explicit expression of vector A
x; tjxo ; to , introduced in Eqs. (33) and (34). The components of A correspond to the correction of the estimation of the functional associated to the simpli®ed problem, when allowance is given to the possible stochastic transitions via the Lagrangian. We have: ÿ 1 L A1
x; tjxo ; to ÿlb
xjxo x ÿ xo ea1
tÿto : min t ÿ to ; `n
A:1 a1 xo
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
ÿa2
tÿto
A2
x; tjxo ; to ÿl
1 ÿ po ÿ p1 x ÿ `e
a2 =a1 ! ` xo
ÿ 1 ` 1 ` :H x ÿ xo ea1
tÿto : min t ÿ to ÿ `n ; `n a1 xo a2 xo
1567
A:2
a3 =a1 ! " ` 1 ` min t ÿ to ; `n A3
x; tjxo ; to ÿlp1 x ÿ `e xo a1 xo # ÿ 1 ` 1 ` :H x ÿ xo ea1
tÿto :H
` ÿ x min t ÿ to ÿ `n ; `n a1 xo a3 xo
A:3
lb
xjxo
a4 ÿ a1 x x a4 =a1 ÿ ÿ o :H x ÿ xo ea1
tÿto :H xo ea4
tÿto ÿ x :H x ÿ Lea4
tÿto L
A:4
ÿa3
tÿto
A4
x; tjxo ; to
a5 =a1 ! l
1 ÿ po ÿ p1 ÿa5
tÿto ` ÿx :H
` ÿ x:H `e A5
x; tjxo ; to x xo " a5 =a4 ! a5 =a4 ! a a
tÿt a4 ` ÿ 1 5a o ÿa5
tÿto ` 4 :H x ÿ `e ÿx : :H `e a5
a4 ÿ a1 xo xo 0 1 # a2 =a1 ! a5 1 1 B C 1 ` ` a a ÿa2
tÿto ÿa
tÿt 1 2 5 o C :H x ÿ `e : :HB @x ÿ xo e A a2 ÿ a5 xo xo a6 =a4 la4 p1 a6
tÿto xo :H
x ÿ `:H `e ÿx A6
x; tjxo ; to a6
a4 ÿ a1 x ` x a6 =a1 a6 a1
tÿto x a6 =a4 o o H x ÿ `e a4 :H x ÿ `ea6
tÿto ` ` a3 =a1 ! lp1 ` ÿa3
tÿto :H x ÿ `e
a3 a6 x xo 1 0 1 1 x a6 a6 =a1 o B a6
tÿto xo a1 a3 C ÿx :H `e :H@x ÿ xo ea6
tÿto A ` `
A:5
A:6
Appendix B The analytical expression of the probability distributions corresponding to the problem presented in Section 3 can be obtained from the integral formulation of the CK equation (Devooght and Smidts, 1992):
1568
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
t
x; i; t
u; i; o
x ÿ g i
t; u : exp ÿ li
g i
s; u ds du 0
tÿ X
t pji
u
x ÿ g i
t ÿ ; u exp ÿ li
g i
s; u ds
u; j; dud
j6i
0
0
B:1 where g i
t; u is the explicit solution of Eq. (1), with the initial condition u g i
o; u . In our one-dimensional case with linear dynamics, this solution reduces to gi
t; u ueai t
B:2
where the sign of the argument in the exponential depends on the state i (see Fig. 1). Setting
x; i; o i1
x ÿ xo , we can deduce from Eq. (B.1) the distributions in the dierent states. State 1 Accounting for the probability po to remain in state i after the crossing of the control setpoint [see Eq. (26)], we have:
t
ÿ
x; 1; t du
uÿxo xÿuea1 t :b
xjxo : exp ÿ lH
uea1 s ÿxo :H
Lÿuea1 s ds 0 ÿ 1 L a1 t x ÿ xo e :b
xjxo exp ÿl: min t; `n a 1 xo
B:3 State 2 Having p12
u a1 `
1 ÿ po ÿ p1
u ÿ `, we can write, using Eq. (B.3):
t ÿ
x; 2; t a1 `
1 ÿ po ÿ p1
`; 1; : x ÿ `eÿa2
tÿ 0
t ÿa2 s ÿ xo :H
L ÿ eÿa2 s `ds d : exp ÿ lH
`e 0 a2 =a1 ! x l=a1 ÿ ` o ÿa2 t : x ÿ `e :H xo ea1 t ÿ `
1 ÿ po ÿ p1 : ` xo 1 ` 1 ` : exp ÿl: min t ÿ `n ; `n a 1 xo a 2 xo
B:4
State 3 The distribution in this state is easily obtained from Eq. (B.4) by symmetry considerations: a2 is to be replaced with a3 , and
1 ÿ po ÿ p1 with p1 .
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1569
State 4 This state can be reached only from state 1. The system will keep evolving in it from u to x with a probability b
xju. Therefore
t
ÿ
x; 4; t d dulH
u ÿ xo H
L ÿ u:
u; 1; : x ÿ uea4
tÿ b
xju 0
t
0
ÿ lH
L ÿ xo ea1 x ÿ xo ea1 a4
tÿ eÿl d:b
xjxo
B:5
The argument of the Dirac peak in Eq. (B.5) vanishes for
1 x a4 t ÿ `n a4 ÿ a1 xo
B:6
which must belong to the interval 0; t. Considering also the supplementary constraint induced by the stepfunction in Eq. (B.5), we ®nd: l ÿa4 t a ÿa ÿ 4 1 l xe : :H x ÿ xo ea1 t
x; 4; t
a4 ÿ a1 x xo
B:7 a4 =a1 ÿ a4 t a4 t xo :b
xjxo :H xo e ÿ x :H x ÿ Le L State 5 Two situations have to be considered here: a transition on demand from state 4, where the system reaches the setpoint at x `, and a stochastic transition from state 2. We thus write:
t ÿ
x; 5; t da4 `
1 ÿ po ÿ p1
`; 4; x ÿ `eÿa5
tÿ 0
t
ÿ d dulH
u ÿ xo H
L ÿ u
u; 2; x ÿ ueÿa5
tÿ
x; 5; t
x; 5; t 0
B:8 The ®rst integral is easily evaluated, thanks to the Dirac peak, whose argument vanishes for t ÿ a15 `n x` . Since we must have 2 0; t, we obtain, with the constraints introduced by the stepfunctions in Eq. (B.7): l " #a ÿa la4
1 ÿ po ÿ p1 `eÿa4 t ` a4 =a5 4 1
x; 5; t
a4 ÿ a1 a5 x xo x
B:9 a5 =a4 ! a5 =a1 ! ` ÿa5 t ` ÿa5 t :H
` ÿ x:H x ÿ `e ÿx :H `e xo xo
1570
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
The calculation of
x; 5; t is also straightforward, as a second Dirac peak is introduced in
u; 2; [see Eq. (B.4)]. Thereby l a2 =a1 ! l
1 ÿ po ÿ p1 xea5 t xo a2 =a1 a2 ÿa5 ` ÿa2 t
x; 5; t :H x ÿ `e :
a2 ÿ a5 x ` ` xo 0 1 a5 1 1 a5 =a1 ! B C ` a a ÿa t ÿa t 1 2 C:H `e 5 ` 5 ÿx :HB @x ÿ xo e A xo xo
B:10
State 6 The distribution
x; 6; t is again the sum of two contributions:
t ÿ
x; 6; t da4 `p1 :
`; 4; x ÿ `ea6
tÿ 0
t
ÿ d dulH
u ÿ xo H
L ÿ u
u; 3; x ÿ uea6
tÿ
x; 6; t
x; 6; t 0
B:11 By directly integrating the Dirac distributions, while paying careful attention to the constraints determining the support of the distributions, we ®nd after some calculations: ÿa4 t l la4 p1 `e x a4 =a6 a4 ÿa1
x; 6; t :H
x ÿ `
a4 ÿ a1 a6 x xo `
B:12 x a6 =a1 x a6 =a4 o o ÿx :H x ÿ `ea6 t :H `ea6 t ` `
and
ÿa6 t l a3 =a1 ! lp1 xe xo a3 =a1 a3 a6 ` ÿa3 t
x; 6; t :H x ÿ `e
a3 a6 x ` ` xo 1 0 1 1 x a6 x a6 =a1 o o B a1 a3 C ÿx :H@x ÿ xo ea6 t A:H `ea6 t ` `
B:13
Failure probability The probability of exceeding level L is given by:
1
1
1
x; 1; tdx
x; 4; tdx
x; 6; tdx PL
t L
L
L
From Eqs. (B.3), (B.7), (B.12) and (B.13), we obtain:
B:14
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1571
! l ÿa4 t a ÿa ÿ a1 t ÿ ÿ 4 1 Le PL
t po H xo e ÿ L po 1 ÿ :H L ÿ xo ea1 t H xo ea4 t ÿ L xo 8 2 3 > l > a4 =a6 !a4 ÿa < a6 =a4 1 ÿa4 t `e L 5:H `ea6 t xo ÿL p1 41 ÿ > xo ` ` > :
9 = x a6 =a1 x l=a1 x a6 =a1 o o o ÿL :H L ÿ `ea6 t 1ÿ :H `ea6 t ; ` ` ` 8 3 >2 l a6 =a1 !a3 a < x l=a1 > a6 =a1 6 ÿa6 t Le ` o 5:H `ea6 t xo 41 ÿ ÿL p1 > ` ` ` xo > : 0 0 1 19 = x a6 1 1 x l=a3 x a6 1 1 a1 a3 A a1 a3 o o o :H@L ÿ xo ea6 t :H@xo ea6 t 1ÿ ÿLA ; ` ` `
B:15
Safe shutdown probability The system will be brought down at a level lower than xo with a probability: P o
t
xo o
x; 2; tdx
xo o
x; 3; tdx
xo o
x; 5; tdx
B:16
Therefore, using results (B.4), (B.9) and (B.10), we have: ! 1 1 x l 1 1 a1 a2 o a1 a2 t Po
t
1 ÿ po ÿ p1 ÿ` :H xo e ` ! 1 1 x l 1 1 a1 a3 o a1 a3 t ÿ` p1 : :H xo e ` 82 3 l > 1a4 =a5 !a4 ÿa a5 =a4 ! < 1 ` ` 6 7 ÿa4 t ÿa5 t
1 ÿ po ÿ p1 41 ÿ e 5:H xo ÿ `e > x x o o : 9 ! a5 =a1 a5 =a1 !> = l=a 1 ` xo ` ÿxo 1 ÿ :H xo ÿ `eÿa5 t :H `eÿa5 t > ` xo xo ;
1572
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
3 (2 l x l=a1 x 1a5 =a1 a2 ÿa x l=a2 5 o o o 4 ea5 t 5
1 ÿ po ÿ p1 ÿ ` ` ` 0 1 ! a5 1 1 a5 =a1 a1 a2 C ` ` B a5 t ÿa5 t ÿxo :H@e ÿ A:H `e xo xo
B:17
a5 =a1 !) x l=a2 ` o ÿa5 t :H xo ÿ `e 1ÿ ` xo
Appendix C In the error-compensating term of the Lagrangian, the components of vector z, such that z
x; t; x0 ; ; xo ; to :H
x0 ÿ xo :H
L ÿ x0 P~
x; t; x0 ; ÿ
x0 :P~
x0 ; ; xo ; to
to ;
C:1
must be evaluated and integrated over x0 and . Using Eqs. (39) to (43), and assuming that the transient is initiated in state 1, we obtain: 3 2 ÿ P~ 11 :P~ 11 ÿ l 2 lq 7 6 ÿ 7 6 ÿ l2 lq P~ 21 :P~ 11 P~ 22 :P~ 21 7 6 7 6 ÿ 7 6 ÿ l2 lq P~ 31 :P~ 11 P~ 33 :P~ 31 7 6 7 6 ÿ 7 6 ÿ l2 lq P~ 41 :P~ 11 P~ 44 l2 P~ 11 ÿ lq P~ 41 7 6 7 6 7 6 l2 P~ 55 :P~ 21 7 6 ÿ 7 6 ÿ l l P~ :P~ P~ :P~ ~ ~ ~ 6 2 q 61 11 63 31 P 66 l2 P 31 ÿ lq P 63 7 7 6 z 6 l P~ :P~ 7: 7 6 q 77 11 7 6 7 6 lq P~ 87 :P~ 11 P~ 88 :P~ 21 7 6 7 6 7 6 lq P~ 97 :P~ 11 P~ 99 :P~ 31 7 6 7 6 7 6 l P~ :P~ P~ ~ 10;10 :P 41 7 6 q 10;7 11 7 6 7 6 l P~ :P~ P~ :P~ P~ ~ 7 6 q 11;7 11 11;8 21 11;10 :P 41 5 4 lq P~ 12;7 :P~ 11 P~ 12;9 :P~ 31 P~ 12;10 :P~ 41 P~ 12;12 :P~ 61
C:2
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1573
For the sake of clarity, we have not written in Eq. (C.2) the dependence on time and on the process variable of the P~ ij 's. However, in each term of this equation, the left factor depends on
x; t; x0 ; , and the right factor on
x0 ; ; xo ; to . A failure can only occur in states 1,3,7,8,9 and 10. Hence only the corresponding components of z need to be integrated, operator h [see Eq. (35a)] being applied to the result. With ! a3 =a1 ÿ a1
tÿto xo a3
tÿto ÿ L p1 q:H `1 e ÿL
C:3 c~L
t po q:H xo e `1 the solution of the simpli®ed problem, we ®nd for the Lagrangian the following expression: ÿ l2 L :H xo ea1
tÿto ÿ L LL po q 1 ÿ `n a 1 xo ! a3 =a1 ÿ 1 `1 1 L xo a3
tÿto p 1 q 1 ÿ l2 l q `n `n ÿL :H `1 e a1 xo a3 `1 `1 ÿ ÿ lq po q L a7
t ÿ to ÿ `n :H xo ea7
tÿto ÿ L :H L ÿ xo ea1
tÿto a7 ÿ a1 xo " ! a8 =a7 lq
1 ÿ po ÿ p1 q a7 L `1 x o a7
t ÿ to ÿ `n ÿ `n ÿL :H `1 ea8
tÿto a8 ` 1 xo `1 a7 ÿ a1 !# a8 =a1 ! a8 =a1 xo a8 ÿ a1 ` 1 xo a8
tÿto a8
tÿto `n :H `1 e ÿL :H L ÿ `1 e `1 a1 xo `1 " ! a8 =a1 lq
1 ÿ po ÿ p1 q a8 `1 L x o a8
t ÿ to ÿ `n ÿ `n ÿL :H `1 ea8
tÿto a2 a8 `1 a1 x o `1 0 1 a8 1 1 B xo a1 a2 C a8
tÿto C :HB @ L ÿ xo e A `1 13 a8 1 1 7 a2 a8 ` 1 B xo a1 a2 ÿLC C7 `n :HB xo ea8
tÿto @ A 5 a2 xo `1 0
! a9 =a7 a7 ` 1 `1 x o a7
t ÿ to ÿ `n ÿ `n ÿL :H `1 ea9
tÿto a9 xo xo `1 !# a9 =a1 ! a9 =a1 xo a7 ÿ a1 ` 1 xo a9
tÿto a9
tÿto :H L ÿ `1 e `n :H `1 e ÿL `1 a1 xo `1 lq p1 q a7 ÿ a1
"
1574
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
! a9 =a1 a9 ` 1 L xo a9
tÿto a9
t ÿ to ÿ `n ÿ `n ÿL :H `1 e `1 a1 xo `1 !# a3 =a1 ! a3 =a1 xo a9 ÿ a3 L xo a3
tÿto a3
tÿto `n :H `1 e ÿL :H L ÿ `1 e `1 a1 `1 `1 " ! a10 =a7 lq
1 ÿ q a7 L `2 x o a7
t ÿ to ÿ `n `n ÿL :H `2 ea10
tÿto a10 `2 xo `2 a7 ÿ a1 !# a10 =a1 ! a10 =a1 xo a7 ÿ a1 ` 2 xo a10
tÿto a10
tÿto `n H `2 e ÿL :H L ÿ `2 e `2 a1 xo `2 " ! a10 =a1 lq
1 ÿ qpo a10 `2 L x o a10
t ÿ to ÿ `n ÿ `n ÿL :H `2 ea10
tÿto a4 a10 a1 xo `2 `2 a10 =a4 a10 =a1 ! `1 xo a4 a10 `2 a10
tÿto `n :H L ÿ `1 e `2 `2 a4 `1 !# a10 =a4 a10 =a1 `1 xo :H `1 ea10
tÿto ÿL `2 `2 lq
1 ÿ qp2o a10 `2 a10 `2 L a10
t ÿ to ÿ `n ÿ `n ÿ `n a4 a10 a4 ` 1 a1 xo `1 0 1 ! a10 =a4 a10 =a1 a10 1 1 B `1 xo xo a1 a4 C C ÿL :HB L ÿ xo ea10
tÿto :H `1 ea10
tÿto @ A `2 `2 `2 lq p1 q a9 ÿ a 3
"
13 a10 1 1 7 a4 a10 `1 B xo a1 a4 ÿLC C7 `n :HB xo ea10
tÿto @ A 5 a4 xo `2 0
C:4
References Davis, M.H.A., 1993. Markov models and optimization. In: Monographs on Statistics and Applied Probability. No. 49, Chapman Hall. Devooght, J., Smidts, C., 1992. Probabilistic reactor dynamics. I. The theory of continuous event trees. Nucl. Sc. Eng. 111, 229. Labeau, P.E., Lewins, J.D., 1999. On the feasibility of variationally processed calculations in dynamic reliability. Proc. of Esrel'99 2, 987. Labeau, P.E., Smidts, C., Swaminathan, S., in press. Dynamic reliability: towards an integrated platform for probabilistic risk assessment, Rel. Engng. Syst. Safety. Lewins, J.D., 1978. Variational method in neutron statistics. Ann. Nucl. En. 5, 141. Lewins, J.D., Parks, G.T., Chang, M.C., 1998. A variational principle using Monte Carlo approximations for reliability studies. Proc. of ISSAT.
P.E. Labeau / Annals of Nuclear Energy 27 (2000) 1543±1575
1575
Lewins, J.D., Chang, M.C., 1999. Monte Carlo variational processing for semi-Markov systems. Proc. of Esrel'99, 2, 921. Lewins, J.D., in press. Classical perturbation theory for Monte Carlo studies of system reliability, to be published in the special issue of Nuclear Science and Engineering, in memory of G.Pomraging. Marchand, S., Tombuyses, B., Labeau, P.E., 1998. DDET and Monte Carlo simulation to solve some dynamic reliability problems. Proc. of PSAM-4 3, 2055. Siu, N., 1994. Risk assessment for dynamic systems: an overview. Rel. Engng. Syst. Safety 43, 43. Smidts, C., Swaminathan, S., 1999. An ESD-probabilistic dynamics hybrid approach for scenario generation. Proc. of M&C'99 2, 1349.