Maximum Entropy Multivariate Analysis of Uncertain Dynamical ...
Maximum Entropy Multivariate Analysis of Uncertain Dynamical Systems Based on the Wiener-Askey Polynomial Chaos Gabriele D’Antona1, Antonello Monti2, Ferdinanda Ponci2, Luca Rocca1 1
Dipartimento di Elettrotecnica , Politecnico di Milano P.zza Leonardo da Vinci 32, 20133 Milano, Italy Phone: +39-02-2399-3706, Fax: +39-02-2399-3703, Email: [email protected] 2 Department of Electrical Engineering , University of South Carolina 301 S. Main st., Columbia SC 29208, USA Phone: +1-803-7777365, Fax: +1-7778045, Email: [email protected]
AMUEM 6023
Abstract – Many measurements models are formalized in terms of a stochastic Ordinary Differential Equations (ODE) relating its solution to some given observables. The expression of the measurement uncertainty for the solution requires the determination of its (joint)Probability Density Function (pdf), evaluated in an assigned time window. Recently, Polynomial Chaos Theory (PCT) has been widely recognized as a promising technique in order to address the problem. The uncertainty estimation via PCT requires the use of a Monte Carlo integration sampling strategy. In this paper a novel approach will be presented in order to achieve the PCT uncertainty estimation on the basis of an analytical methodology, requiring only an optimization calculus.
Keywords – Measurement uncertainty, maximum entropy, uncertain systems, density function, stochastic circuits. I. INTRODUCTION Many measurements models are formalized in terms of dynamic systems, relating some given observables to an unknown quantity which has to be estimated. In this work we will restrict our attention to systems governed by linear Ordinary Differential Equations (ODE) of the form
x (t ) a x(t ) b ut
x(0) 0
(1)
Although it has been common practice in engineering to analyze systems, which are based on deterministic mathematical models of the type (1), in which all the quantities involved are deterministic, such ideal situations are rarely encountered in practice. As a matter of fact, the solution x(t) of the system (1) is in practice affected by uncertainty. This is because the knowledge of either the model parameters a and/or b or the input u(t) is uncertain. Furthermore, the initial condition x(0) of the ODE (1) can be uncertain as well. All these sources of uncertainty can be present one at time or even simultaneously. Under the more general case that all the sources of uncertainty are simultaneously present, the ordinary differential equation (1) has to be replaced by the following linear stochastic equation
x (t ) a x(t ) b u(t )
x(0) x0
(2)
where a and b are random variables and x(t) and u(t) are stochastic processes. Although, the solution x(t) of (2) is viable on the basis of a random sampling strategy based on Monte Carlo (MC) integration, it is however notoriously computationally intensive. In particular, in order to fully characterize the stochastic process x(t) (the solution of (2)), it is necessary to run a sequence of N MC simulations, relative to different time instants ti (i=1,…,N), starting from the prior knowledge of the stochastic process u(t) and the probability distributions of the random variables a and b. In order to overcome this computational limit, recently generalized Polynomial Chaos Theory (PCT) or Askey Chaos, has been proposed [1], and its application has been suggested by many authors coming from different branches of science and technology. For example, in [2], [3] applications to solve stochastic electrical networks have been suggested. According to PCT it is possible, under appropriate assumptions, to approximate the solution x(t) of (2) by means of the following analytical expression:
P
x(t ) ak (t ) f k (ξ)
(3)
k 0
In (3) fk(ξ) are assigned polynomials, which are considered function of an appropriate multidimensional random vector ξ=(ξ1, ξ2,…, ξM) and ak(t) are known (deterministic) coefficients. The main advantage of the PCT respect to a conventional MC method is given by the fact that in the PCT the stochastic feature of the solution x(t) is fully characterized by the random variables fk(ξ), while the time dependence of x(t) is fully described by the deterministic coefficients ak(t). Thus, the stochastic process x(t) can be fully characterized running a unique MC simulation, starting from the prior knowledge of the pdf relative to the random vector ξ. In this work we restrict specifically our attention to the case of a scalar random variable ξ only, although the reasoning we will propose can be further extended to take into account random vectors ξ as well. In particular, following the basic steps that have already been presented in [4], we will suggest a novel approach to express the uncertainty relative to the stochastic process x(t), based on an analytical estimation of its joint pdf, evaluated in correspondence of prescribed time instants t1,…,tN. The availability of an analytical expression for the PCT solution’s pdf has many advantages. First of all it is a natural final step of the PCT considering that the solution is expressed as a polynomial combination of a variable with known distribution. This is a significant advantage with respect to the MC analysis that provides only a set of samples, but the possibility of transforming this result in an analytical expression of the pdf dramatically enriches the value of the PCT. In effect, it allows for more sophisticated post-processing of the simulation results. The PCT approach can play a significant role in defining level of acceptance of a design affected by uncertainty. Thanks to analytical expression it is possible to quickly check many characteristics of the output pdf for different instants of time.
Last but not least the algorithm can be used to create realistic pictures of the evolution of the uncertain system representing the analytical curve as a function of time. In particular, the analytical expression of the pdf will result to be parameterized in terms of an appropriate set of coefficients, acting on which it will be possible to interpolate the whole shape of the pdf at prescribed times. This constitutes a significant improvement respect to the conventional MC analysis, as the number of parameters are far smaller that the number of MC samples. Furthermore, the MC analysis usually gives a noisy shape and that can make quite impossible to appreciate the evolution in time. The suggested analytical approach instead, is based on a regularized estimation of the joint pdf for the state x(t), which allows to define a well defined shaped pdf. The novel methodology is based on a higher order statistical moments analysis. As matter of fact, if we know the pdf of the random vector ξ, than we have an a priori knowledge of its moments of any order. Based on (3), it is possible to propagate the higher order moments of ξ to obtain those of x at any instants t1,…,tN. Once evaluated some of the higher order moments of x at t1,…,tN, the joint pdf can be estimated analytically on the basis of the maximum entropy (ME) approach [5]. According to the ME, the pdf of x is the one satisfying the prescribed set of moments constraints and maximizing the Shannon entropy, whose definition generalizes the well known one coming from statistical mechanics. Technically, the entropy maximization palys the role of regularizing the solution, thus enforcing a suitable well behaved shape for the estimated pdf. In order to prove the goodness of the suggested methodology, an application in the field of electrical science will be analyzed. The objective of the application is the estimation of a voltage transient in a first order electrical network, under the hypothesis that the values of some circuit parameters are uncertain. A comparison between the suggested approach and a conventional Monte Carlo sampling strategy will be given, on the basis of their relative statistical performances.
II. THE NOVEL APPROACH Focusing our attention to the case of a single random variable ξ, we first rewrite the Polynomial Chaos expansion (3) in the following compact matrix notation:
x A f (ξ)
(4)
where x and f(ξ) are random vectors defined as below:
x [x(t1 ),..., x(t N )]T
(5)
f (ξ) [ f 0 (ξ),..., f P (ξ)]T
(6)
and A is the following (fixed) matrix of Polynomial Chaos coefficients
a0 (t1 ) a p (t1 ) A a0 (t N ) a p (t N )
(7)
In order to estimate the nth order moment of the random vector x in (4), the first step is the propagation of the nth order moment mn,ξ of the random variable ξ through the vector function f(ξ). The nth order moment mn,f of the random variable f(ξ) is expressed, by definition, as:
mn, f E f (ξ)( n )
where the symbol (·)(n) denotes the n-fold Kronecker product [6]:
(8)
f ξ
n
f ξ f ξ ... f ξ
(9)
n times
The vector of polynomials defined in (6), can be written as reported below
P
f c k k
(10)
k 0
where in (10) ck are vectors of polynomials’ coefficients. Applying the definition (8), it can be verified that the nth order moment mn,f can be expressed as
mn , f
1 c k n mk1 ... kn , k1 0 k n 0 k1!...k n !
(11)
The relationship (11) allows the determination of the nth order moment mn,f of the random variable f(ξ), starting from the (prior) knowledge of the (k1+,…,kn)th order moment of the random variable ξ. Once the nth order moment mn,f of f(ξ) has been computed, the next step is its propagation through matrix transformation defined in (4), in order to determine the nth order moment mn,x of the random variable x. The first order moment can be computed directly from the linear system (4), employing elementary laws of statistics, as:
m1, x A m1, f
(12)
The estimation of higher order moments (n2) instead, requires a preliminary scaling of both members in the linear system (4) by the factor -m1,x:
x A[ f (ξ)]
(13)
where
x x m1, x
f (ξ) f (ξ) m1, f
and
(14)
It is trivial to verify that Δx is a zero mean random variable such that:
mn , x mn , x
n 2
(15)
The nth order moment mn,Δx can be computed starting from the knowledge of the nth order moment mn,Δf of the random variable Δf(ξ), employing the following higher order moments propagation formula [7]:
mn, x A( n ) mn, f
(16)
In (16), the nth order moment mn,Δf can be computed starting from the knowledge of the moments mk,f , with k = 0,…, n, determined on the basis of (11), making use of the following binomial power formula [8]:
n
[ f ( ) m1, f ]( n ) M k [ f k 0
n
(k )
( ) (m1, f )( n k ) ]
(17)
where {M0n,…, Mnn} are the matrix coefficients of the binomial power formula. Details of their computation are given in [8]. Applying the expected value operator to both members in (17), it is trivial to obtain the following relationship between mn,Δf and mk,f:
n
mn, f M k mk , f [m1, f ]( n k ) n
(18)
k 0
It is then possible to apply the propagation formula (16). At last, joining (15) and (12), we can estimate all moments up to an assigned order n of the random variable x. The pdf p(x) of x can finally be determined on account of the Maximum Entropy (ME) approach [5]. Derived from the ME, the joint pdf p(x) can be approximated by means of the following exponential function
n (i ) ~ p ( x) exp i x k 0
(19)
where λi is a row vector of unknown parameters, which are determined by finding a least squares solution for the following system of L+1 moments constraints
x
(n)
p( x)d x m n,x
n 0,..., L
(the zero order moment m0,x is assumed equal to 1).
(20)
III. THE APPLICATION In order to prove the goodness of the proposed approach, it has been applied to the study of a voltage transient in the RC electrical network shown in figure 1, containing a resistor, a DC voltage source, and a capacitor. The state of the system is represented by the voltage across the capacitance, which, under the hypothesis that all the circuit parameters are not affected by uncertainty, is governed by the following ordinary differential equation:
Vin vc (t ) RC
dvc (t ) dt
(21)
where Vin is the voltage source, vc(t) is the voltage across the capacitance, R is the resistance and c the capacitance. Under the hypothesis instead that some circuit parameters are affected by uncertainty, the ODE (21) has to be replaced by a suitable stochastic differential equation. In the following two distinct cases will be analyzed, in which the only uncertain variable in the system will be in turn the capacitance (case 1)[4], and the voltage source (case 2).
A. Case 1 The only uncertain variable in the system is the capacitance of the capacitor, modeled in terms of the random variable c. Since the value of the capacitance is uncertain, the ODE (21) has to be replaced by the following stochastic equation:
Vin v c (t ) Rc
dv c (t ) dt
(22)
where the symbol vc(t) is a stochastic process modeling the dynamics of the capacitance voltage. The stochastic differential equation (22) has been solved employing the Wiener-Askey polynomial
chaos expansion, leading to a solution for the capacitance voltage of the form (3), particularized as below for the case of our electrical network application:
P
v c (t ) vk ,C (t ) f k (ξ C )
(23)
k 1
In (23) ξc is a random variable such that, according to the Wiener-Askey theory, its pdf is a shifted and rescaled version of that relative to the capacitance c. Thus, assuming the pdf of c is a priori known, the moments of any order relative the random variable ξc result to be a priori known as well. Starting from a (finite) set of higher order moments mn,ξc of ξc, it is than possible to estimate some higher order moments mn,v relative to vc=[vc(t1),…, vc(tN)], on the basis of the propagation moments rules described in the previous paragraph. Once a set of higher order moments for vc has been estimated, it will be possible to estimate its pdf, on the basis of the analytical maximum entropy approximation (19). In our case study we have assumed that the parameters of the circuit are as follows:
R=1 Ω
C = 1 mF
Vin=10 V
with C being the mean of the capacitance. We have further assumed that the capacitance c is distributed according to a uniform pdf between [-E, +E], with E being the 30% of the central value
C . Thus, based on the Wiener-Askey theory [1], the polynomials fk(ξc) in (23) have been assumed to be Legendre polynomials. Figure 2 shows the first order moment of the voltage transient from 0 to 10 ms. together with the standard deviation relative to 20 uniformly distributed time instants. It can be noticed that as the
voltage vc reaches its steady state value of 10 V, the uncertainty is progressively reduced. This behavior can be well understood on the basis of elementary considerations relative to the physics of the RC circuit shown in figure 1. As a matter of fact, having assumed the circuit supplied with a DC voltage, under steady state conditions the voltage across the capacitance is not influenced by the uncertain parameter (the capacitance) and, thus, no uncertainty is associated with the knowledge of the voltage vc. Figures 3 and 4 show the voltage marginal pdf’s at the times t=1 ms and t=3 ms (see figure 2) estimated on the basis of the suggested approach, using the first 5 order moments m0,v, m1,v, m2,v, m3,v , m4,v, together with a comparison with the correspondent pdf’s recovered by a conventional PCT, running ten thousands of Monte Carlo trials of the random variable ξc in (23) [4]. Figures 5 and 6 provides the multivariate analysis. In particular, a couple of time instants, t1 and t2, is taken into account and the random voltages vc(t1) and vc(t2) are evaluated in correspondence. Since in the more general case vc(t1) and vc(t2) are not independent (with each other), in order to analyze possible correlations between them it is necessary to estimate the bivariate pdf of the full random vector [vc(t1) vc(t2)]T. To better compare the results, in figure 5 and 6 the bivariate pdf’s have been plotted with a scaled common mean (the origin). This is equivalent to scaling the random
vector [vc(t1) vc(t2)]T by its mean v c t1 v c t 2 in advance: T
v t v c t1 vc c 1 v c t 2 v c t 2
(24)
In order to help understand the contours and the extent of the estimated pdf’s, in both figures 5 and 6 the plots are color-coded projection of the correspondent 3D plots. In particular, in figure 5 it is assumed that the time instants t1 and t2 are “close” to each other: t1 = 0.9 ms, t2 = 1 ms (see also figure 2 in order to visualize the couple of time instants). In figure 6 instead the same analysis of figure 5 is repeated considering t1= 1 ms and t2 = 3 ms, which have a
wider distance in time (figure 2). In both cases, the correlation coefficients relative to vc(t1) and vc(t2) have been computed and have resulted to be equal to 1. Based on the multivariate analysis of figure 5 and 6, it is thus possible to conclude that as time elapses (from t2 = 1 ms to t2 = 3 ms) the variances of the random capacitance voltage tends to becomes smaller and smaller (compare the error voltages at t2 of figures 5 and 6). However, the correlation between voltages (at t1 and t2) is unchanged.
B. Case 2 In case 2, the only uncertain variable is the voltage source, modeled in terms of the random variable Vin. Thus the ODE (21) has to be replaced by the following stochastic equation
Vin v c (t ) RC
dv c (t ) dt
(25)
The solution of (25) obtained on the basis of the Wiener-Askey polynomial chaos expansion is given by the following expression
P
v c (t ) vk , Vin (t ) f k (ξ Vin ) k 1
where vk,Vin(t) and ξVin have an analogous meaning to that provided for vk,C(t) and ξC in (23). The circuit parameters are as follows:
V in =10 V
R=1Ω
(26)
C=1 mF
with V in being the mean of the voltage source. The voltage V in is distributed according to a uniform pdf between [-E, +E], with E being equal to the 20% of the central value V in . Figure 7 shows the first order moment of the capacitance voltage transient from 0 to 10 ms. together with the standard deviation relative to 20 uniformly distributed time instants. Figures 8 and 9 show the capacitance voltage marginal pdf’s at the times t=1 ms. and t=5 ms. estimated on the basis of the suggested approach, using the first 5 order moments m0,v, m1,v, m2,v, m3,v , m4,v, together with a comparison with the correspondent pdf’s recovered by a conventional PCT, running ten thousands of Monte Carlo trials of the random variable ξVin in (26). Analogously to what previously illustrated in figure 5 and 6, figure 10 and 11 provides the multivariate analysis. In particular figure 10 gives the joint pdf relative to the voltages at t1=0.9 ms. and t2=1 ms. (see also figure 7 to visualize the time instants). Figure 11 repeats the same analysis at t1=1 ms. and t2=5 ms. In a similar fashion to case 1, the correlation coefficients are equal to 1.
IV. CONCLUSION This paper has focused on the analysis of measurements models expressed in terms of a stochastic ode. In particular, a novel approach to express the uncertainty for its solution has been suggested, based on the Wiener-Askey Polynomial Chaos theory. While the PC allows the calculation of the solution of the stochastic differential equation very quickly and faster than using a traditional Monte Carlo analysis, the results are anyhow available in terms of a polynomial expansion. This polynomial expansion allows a rapid evaluation of certain properties of the resulting pdf as the mean value or an estimation of the worst case, but it does not support a complete and simple reconstruction of the final pdf.
A simple approach to the problem is the evaluation of the solution by means of a suitable sampling of the random quantity ξ. This process is very computationally efficient but it looses the possibility to have an analytical solution which is one of the significant advantages of PC. The method here proposed, is less computationally efficient (on a regular PC it may take about 30 times more than the sampling process), but it preserve the analytical aspect of the solution. In conclusion, in the authors’ opinion, both approaches can be valuable to post-process the PC solution. The random sampling is very attractive for quick display of the results while the method proposed by this paper is very important for any situation where the analytical form is necessary. For example the analytical solution better support the possibility to calculate different values of the coverage factors. Another interesting example is given by the need of using the results of the analysis of the measurement process for a second analysis step where the measurement uncertainty has to be used as input. A concrete example of such a situation is a control design performed in the PC space that includes the measurement as source of uncertainty. Future work will focus on the application of the suggested methodology to the case in which more than one uncertain parameter is present in the system. It will be thus be assumed that the Polynomial Chaos solution of each equation is obtained starting from (3), which, according to a conventional approach, requires running a Monte Carlo of the whole random vector ξ. AKNOWLEDGEMENT This work has been partially supported by the US Office of Naval Research under the grant N00014-02-1-0623. REFERENCES [1] D. Xiu & Karniadakis, G. E. 2002b., “The Wiener-Askey polynomial chaos for stochastic differential equations”, SIAM J. Sci. Comput., vol. 24, pp. 619-644, 2002.
[2] Monti, A., Ponci, F., Lovet T., “A polynomial chaos theory approach to the control design of a power converter”, IEEE Annual Power Electronics Specialists Conference, Vol. 6, pp.48094813, June 2004. [3] Monti, A., Ponci, F., Lovet T., Smith A., Dougal R., “Modeling of uncertainty and applications in monitoring and control of power electronics”, Proc. of the American Control Conference, Vol3, pp.2011-2016, 8-10 June 2005. [4] D’Antona G., Monti A., Ponci F., Rocca L, Maximum Entropy Analytical Solution for Stochastic Differential Equations Based on the Wiener Askey Polynomial Chaos, Proc. 21st IEEE Advanced Methods for Uncertainty Estimation in Measurement, Vol.1 , 20-21 April 2006, pp. 62 - 66. [5] D’Antona, G., “Expanded uncertainty and coverage factor computation by higher order moments analysis”, Proc. 21st IEEE Instr. and Meas. Technology Conf., Vol.1 , 18-20 May 2004, p.234 - 238. [6] Brewer, J.W., “Kronecker product and matrix calculus in system theory”, IEEE Trans. on Circuits and Systems, vol. cas-25, No. 9, September 1978. [7] André,T.F., Nowak, R.D., Van Veen, B. D., “Low-Rank Estimation of Higher Order Statistics”, IEEE Trans. on Signal Processing, Vol.45, No.3, March 1997. [8] Carravetta F., Germani A., Raimondi M., “Polynomial Filtering for Linear Discrete Time non Gaussian Systems”, SIAM J. Control and Optimization, Vol. 34, No. 5, pp. 1666-1690, September 1996.
R Vin
C
Fig. 1. The RC circuit.
12
10
voltage [V]
8
6
4
2
0
0
0.002
0.004 0.006 time [s]
0.008
0.01
Fig. 2. First and second order moments for the capacitance voltage.
0.8 0.7 0.6
pdf
0.5 0.4 0.3 0.2 0.1 0
0
2
4
6 voltage [V]
8
10
Fig. 3. Capacitance voltage pdf at t=1 ms.
2.5
2
pdf
1.5
1
0.5
0
0
2
4
6 voltage [V]
8
10
Fig. 4. Capacitance voltage at t=3 ms.
2 1.5
error volatge at t2
1 0.5 0 -0.5 -1 -1.5 -2 -2
-1.5
-1
-0.5 0 0.5 error voltage at t1
1
1.5
2
Fig. 5 The contour plot for the joint pdf at t1 = 0.9 ms and t2 = 1 ms.
2 1.5
error volatge at t2
1 0.5 0 -0.5 -1 -1.5 -2 -2
-1.5
-1
-0.5 0 0.5 error voltage at t1
1
1.5
2
Fig. 6 The contour plot for the joint pdf shown at t1 = 1 ms and t2 = 3 ms.
12
10
voltage [V]
8
6
4
2
0
0
0.002
0.004 0.006 time [s]
0.008
0.01
Fig. 7. First and second order moments for the capacitance voltage .
0.6
0.5
pdf
0.4
0.3
0.2
0.1
0
0
2
4
6 8 voltage [V]
10
12
14
Fig. 8. Capacitance voltage pdf at t=1 ms.
0.3
0.25
pdf
0.2
0.15
0.1
0.05
0
0
2
4
6 8 voltage [V]
10
12
14
Fig. 9. Capacitance voltage pdf at t=5 ms.
2 1.5
error volatge at t2
1 0.5 0 -0.5 -1 -1.5 -2 -2
-1.5
-1
-0.5 0 0.5 error voltage at t1
1
1.5
2
Fig. 10 The contour plot for the joint pdf at t1 = 0.9 ms and t2 = 1 ms.
2.5 2 1.5
error volatge at t2
1 0.5 0 -0.5 -1 -1.5 -2 -2.5
-2
-1
0 error voltage at t1
1
2
Fig. 11 The contour plot for the joint pdf at t1 = 1 ms and t2 = 5 ms.
Dipartimento di Elettrotecnica , Politecnico di Milano P.zza Leonardo da Vinci 32, 20133 Milano, Italy Phone: +39-02-2399-3706, Fax: +39-02-2399-3703, Email: [email protected] 2 Department of Electrical Engineering , University of South Carolina 301 S. Main st., Columbia SC 29208, USA Phone: +1-803-7777365, Fax: +1-7778045, Email: [email protected]
AMUEM 6023
Abstract – Many measurements models are formalized in terms of a stochastic Ordinary Differential Equations (ODE) relating its solution to some given observables. The expression of the measurement uncertainty for the solution requires the determination of its (joint)Probability Density Function (pdf), evaluated in an assigned time window. Recently, Polynomial Chaos Theory (PCT) has been widely recognized as a promising technique in order to address the problem. The uncertainty estimation via PCT requires the use of a Monte Carlo integration sampling strategy. In this paper a novel approach will be presented in order to achieve the PCT uncertainty estimation on the basis of an analytical methodology, requiring only an optimization calculus.
Keywords – Measurement uncertainty, maximum entropy, uncertain systems, density function, stochastic circuits. I. INTRODUCTION Many measurements models are formalized in terms of dynamic systems, relating some given observables to an unknown quantity which has to be estimated. In this work we will restrict our attention to systems governed by linear Ordinary Differential Equations (ODE) of the form
x (t ) a x(t ) b ut
x(0) 0
(1)
Although it has been common practice in engineering to analyze systems, which are based on deterministic mathematical models of the type (1), in which all the quantities involved are deterministic, such ideal situations are rarely encountered in practice. As a matter of fact, the solution x(t) of the system (1) is in practice affected by uncertainty. This is because the knowledge of either the model parameters a and/or b or the input u(t) is uncertain. Furthermore, the initial condition x(0) of the ODE (1) can be uncertain as well. All these sources of uncertainty can be present one at time or even simultaneously. Under the more general case that all the sources of uncertainty are simultaneously present, the ordinary differential equation (1) has to be replaced by the following linear stochastic equation
x (t ) a x(t ) b u(t )
x(0) x0
(2)
where a and b are random variables and x(t) and u(t) are stochastic processes. Although, the solution x(t) of (2) is viable on the basis of a random sampling strategy based on Monte Carlo (MC) integration, it is however notoriously computationally intensive. In particular, in order to fully characterize the stochastic process x(t) (the solution of (2)), it is necessary to run a sequence of N MC simulations, relative to different time instants ti (i=1,…,N), starting from the prior knowledge of the stochastic process u(t) and the probability distributions of the random variables a and b. In order to overcome this computational limit, recently generalized Polynomial Chaos Theory (PCT) or Askey Chaos, has been proposed [1], and its application has been suggested by many authors coming from different branches of science and technology. For example, in [2], [3] applications to solve stochastic electrical networks have been suggested. According to PCT it is possible, under appropriate assumptions, to approximate the solution x(t) of (2) by means of the following analytical expression:
P
x(t ) ak (t ) f k (ξ)
(3)
k 0
In (3) fk(ξ) are assigned polynomials, which are considered function of an appropriate multidimensional random vector ξ=(ξ1, ξ2,…, ξM) and ak(t) are known (deterministic) coefficients. The main advantage of the PCT respect to a conventional MC method is given by the fact that in the PCT the stochastic feature of the solution x(t) is fully characterized by the random variables fk(ξ), while the time dependence of x(t) is fully described by the deterministic coefficients ak(t). Thus, the stochastic process x(t) can be fully characterized running a unique MC simulation, starting from the prior knowledge of the pdf relative to the random vector ξ. In this work we restrict specifically our attention to the case of a scalar random variable ξ only, although the reasoning we will propose can be further extended to take into account random vectors ξ as well. In particular, following the basic steps that have already been presented in [4], we will suggest a novel approach to express the uncertainty relative to the stochastic process x(t), based on an analytical estimation of its joint pdf, evaluated in correspondence of prescribed time instants t1,…,tN. The availability of an analytical expression for the PCT solution’s pdf has many advantages. First of all it is a natural final step of the PCT considering that the solution is expressed as a polynomial combination of a variable with known distribution. This is a significant advantage with respect to the MC analysis that provides only a set of samples, but the possibility of transforming this result in an analytical expression of the pdf dramatically enriches the value of the PCT. In effect, it allows for more sophisticated post-processing of the simulation results. The PCT approach can play a significant role in defining level of acceptance of a design affected by uncertainty. Thanks to analytical expression it is possible to quickly check many characteristics of the output pdf for different instants of time.
Last but not least the algorithm can be used to create realistic pictures of the evolution of the uncertain system representing the analytical curve as a function of time. In particular, the analytical expression of the pdf will result to be parameterized in terms of an appropriate set of coefficients, acting on which it will be possible to interpolate the whole shape of the pdf at prescribed times. This constitutes a significant improvement respect to the conventional MC analysis, as the number of parameters are far smaller that the number of MC samples. Furthermore, the MC analysis usually gives a noisy shape and that can make quite impossible to appreciate the evolution in time. The suggested analytical approach instead, is based on a regularized estimation of the joint pdf for the state x(t), which allows to define a well defined shaped pdf. The novel methodology is based on a higher order statistical moments analysis. As matter of fact, if we know the pdf of the random vector ξ, than we have an a priori knowledge of its moments of any order. Based on (3), it is possible to propagate the higher order moments of ξ to obtain those of x at any instants t1,…,tN. Once evaluated some of the higher order moments of x at t1,…,tN, the joint pdf can be estimated analytically on the basis of the maximum entropy (ME) approach [5]. According to the ME, the pdf of x is the one satisfying the prescribed set of moments constraints and maximizing the Shannon entropy, whose definition generalizes the well known one coming from statistical mechanics. Technically, the entropy maximization palys the role of regularizing the solution, thus enforcing a suitable well behaved shape for the estimated pdf. In order to prove the goodness of the suggested methodology, an application in the field of electrical science will be analyzed. The objective of the application is the estimation of a voltage transient in a first order electrical network, under the hypothesis that the values of some circuit parameters are uncertain. A comparison between the suggested approach and a conventional Monte Carlo sampling strategy will be given, on the basis of their relative statistical performances.
II. THE NOVEL APPROACH Focusing our attention to the case of a single random variable ξ, we first rewrite the Polynomial Chaos expansion (3) in the following compact matrix notation:
x A f (ξ)
(4)
where x and f(ξ) are random vectors defined as below:
x [x(t1 ),..., x(t N )]T
(5)
f (ξ) [ f 0 (ξ),..., f P (ξ)]T
(6)
and A is the following (fixed) matrix of Polynomial Chaos coefficients
a0 (t1 ) a p (t1 ) A a0 (t N ) a p (t N )
(7)
In order to estimate the nth order moment of the random vector x in (4), the first step is the propagation of the nth order moment mn,ξ of the random variable ξ through the vector function f(ξ). The nth order moment mn,f of the random variable f(ξ) is expressed, by definition, as:
mn, f E f (ξ)( n )
where the symbol (·)(n) denotes the n-fold Kronecker product [6]:
(8)
f ξ
n
f ξ f ξ ... f ξ
(9)
n times
The vector of polynomials defined in (6), can be written as reported below
P
f c k k
(10)
k 0
where in (10) ck are vectors of polynomials’ coefficients. Applying the definition (8), it can be verified that the nth order moment mn,f can be expressed as
mn , f
1 c k n mk1 ... kn , k1 0 k n 0 k1!...k n !
(11)
The relationship (11) allows the determination of the nth order moment mn,f of the random variable f(ξ), starting from the (prior) knowledge of the (k1+,…,kn)th order moment of the random variable ξ. Once the nth order moment mn,f of f(ξ) has been computed, the next step is its propagation through matrix transformation defined in (4), in order to determine the nth order moment mn,x of the random variable x. The first order moment can be computed directly from the linear system (4), employing elementary laws of statistics, as:
m1, x A m1, f
(12)
The estimation of higher order moments (n2) instead, requires a preliminary scaling of both members in the linear system (4) by the factor -m1,x:
x A[ f (ξ)]
(13)
where
x x m1, x
f (ξ) f (ξ) m1, f
and
(14)
It is trivial to verify that Δx is a zero mean random variable such that:
mn , x mn , x
n 2
(15)
The nth order moment mn,Δx can be computed starting from the knowledge of the nth order moment mn,Δf of the random variable Δf(ξ), employing the following higher order moments propagation formula [7]:
mn, x A( n ) mn, f
(16)
In (16), the nth order moment mn,Δf can be computed starting from the knowledge of the moments mk,f , with k = 0,…, n, determined on the basis of (11), making use of the following binomial power formula [8]:
n
[ f ( ) m1, f ]( n ) M k [ f k 0
n
(k )
( ) (m1, f )( n k ) ]
(17)
where {M0n,…, Mnn} are the matrix coefficients of the binomial power formula. Details of their computation are given in [8]. Applying the expected value operator to both members in (17), it is trivial to obtain the following relationship between mn,Δf and mk,f:
n
mn, f M k mk , f [m1, f ]( n k ) n
(18)
k 0
It is then possible to apply the propagation formula (16). At last, joining (15) and (12), we can estimate all moments up to an assigned order n of the random variable x. The pdf p(x) of x can finally be determined on account of the Maximum Entropy (ME) approach [5]. Derived from the ME, the joint pdf p(x) can be approximated by means of the following exponential function
n (i ) ~ p ( x) exp i x k 0
(19)
where λi is a row vector of unknown parameters, which are determined by finding a least squares solution for the following system of L+1 moments constraints
x
(n)
p( x)d x m n,x
n 0,..., L
(the zero order moment m0,x is assumed equal to 1).
(20)
III. THE APPLICATION In order to prove the goodness of the proposed approach, it has been applied to the study of a voltage transient in the RC electrical network shown in figure 1, containing a resistor, a DC voltage source, and a capacitor. The state of the system is represented by the voltage across the capacitance, which, under the hypothesis that all the circuit parameters are not affected by uncertainty, is governed by the following ordinary differential equation:
Vin vc (t ) RC
dvc (t ) dt
(21)
where Vin is the voltage source, vc(t) is the voltage across the capacitance, R is the resistance and c the capacitance. Under the hypothesis instead that some circuit parameters are affected by uncertainty, the ODE (21) has to be replaced by a suitable stochastic differential equation. In the following two distinct cases will be analyzed, in which the only uncertain variable in the system will be in turn the capacitance (case 1)[4], and the voltage source (case 2).
A. Case 1 The only uncertain variable in the system is the capacitance of the capacitor, modeled in terms of the random variable c. Since the value of the capacitance is uncertain, the ODE (21) has to be replaced by the following stochastic equation:
Vin v c (t ) Rc
dv c (t ) dt
(22)
where the symbol vc(t) is a stochastic process modeling the dynamics of the capacitance voltage. The stochastic differential equation (22) has been solved employing the Wiener-Askey polynomial
chaos expansion, leading to a solution for the capacitance voltage of the form (3), particularized as below for the case of our electrical network application:
P
v c (t ) vk ,C (t ) f k (ξ C )
(23)
k 1
In (23) ξc is a random variable such that, according to the Wiener-Askey theory, its pdf is a shifted and rescaled version of that relative to the capacitance c. Thus, assuming the pdf of c is a priori known, the moments of any order relative the random variable ξc result to be a priori known as well. Starting from a (finite) set of higher order moments mn,ξc of ξc, it is than possible to estimate some higher order moments mn,v relative to vc=[vc(t1),…, vc(tN)], on the basis of the propagation moments rules described in the previous paragraph. Once a set of higher order moments for vc has been estimated, it will be possible to estimate its pdf, on the basis of the analytical maximum entropy approximation (19). In our case study we have assumed that the parameters of the circuit are as follows:
R=1 Ω
C = 1 mF
Vin=10 V
with C being the mean of the capacitance. We have further assumed that the capacitance c is distributed according to a uniform pdf between [-E, +E], with E being the 30% of the central value
C . Thus, based on the Wiener-Askey theory [1], the polynomials fk(ξc) in (23) have been assumed to be Legendre polynomials. Figure 2 shows the first order moment of the voltage transient from 0 to 10 ms. together with the standard deviation relative to 20 uniformly distributed time instants. It can be noticed that as the
voltage vc reaches its steady state value of 10 V, the uncertainty is progressively reduced. This behavior can be well understood on the basis of elementary considerations relative to the physics of the RC circuit shown in figure 1. As a matter of fact, having assumed the circuit supplied with a DC voltage, under steady state conditions the voltage across the capacitance is not influenced by the uncertain parameter (the capacitance) and, thus, no uncertainty is associated with the knowledge of the voltage vc. Figures 3 and 4 show the voltage marginal pdf’s at the times t=1 ms and t=3 ms (see figure 2) estimated on the basis of the suggested approach, using the first 5 order moments m0,v, m1,v, m2,v, m3,v , m4,v, together with a comparison with the correspondent pdf’s recovered by a conventional PCT, running ten thousands of Monte Carlo trials of the random variable ξc in (23) [4]. Figures 5 and 6 provides the multivariate analysis. In particular, a couple of time instants, t1 and t2, is taken into account and the random voltages vc(t1) and vc(t2) are evaluated in correspondence. Since in the more general case vc(t1) and vc(t2) are not independent (with each other), in order to analyze possible correlations between them it is necessary to estimate the bivariate pdf of the full random vector [vc(t1) vc(t2)]T. To better compare the results, in figure 5 and 6 the bivariate pdf’s have been plotted with a scaled common mean (the origin). This is equivalent to scaling the random
vector [vc(t1) vc(t2)]T by its mean v c t1 v c t 2 in advance: T
v t v c t1 vc c 1 v c t 2 v c t 2
(24)
In order to help understand the contours and the extent of the estimated pdf’s, in both figures 5 and 6 the plots are color-coded projection of the correspondent 3D plots. In particular, in figure 5 it is assumed that the time instants t1 and t2 are “close” to each other: t1 = 0.9 ms, t2 = 1 ms (see also figure 2 in order to visualize the couple of time instants). In figure 6 instead the same analysis of figure 5 is repeated considering t1= 1 ms and t2 = 3 ms, which have a
wider distance in time (figure 2). In both cases, the correlation coefficients relative to vc(t1) and vc(t2) have been computed and have resulted to be equal to 1. Based on the multivariate analysis of figure 5 and 6, it is thus possible to conclude that as time elapses (from t2 = 1 ms to t2 = 3 ms) the variances of the random capacitance voltage tends to becomes smaller and smaller (compare the error voltages at t2 of figures 5 and 6). However, the correlation between voltages (at t1 and t2) is unchanged.
B. Case 2 In case 2, the only uncertain variable is the voltage source, modeled in terms of the random variable Vin. Thus the ODE (21) has to be replaced by the following stochastic equation
Vin v c (t ) RC
dv c (t ) dt
(25)
The solution of (25) obtained on the basis of the Wiener-Askey polynomial chaos expansion is given by the following expression
P
v c (t ) vk , Vin (t ) f k (ξ Vin ) k 1
where vk,Vin(t) and ξVin have an analogous meaning to that provided for vk,C(t) and ξC in (23). The circuit parameters are as follows:
V in =10 V
R=1Ω
(26)
C=1 mF
with V in being the mean of the voltage source. The voltage V in is distributed according to a uniform pdf between [-E, +E], with E being equal to the 20% of the central value V in . Figure 7 shows the first order moment of the capacitance voltage transient from 0 to 10 ms. together with the standard deviation relative to 20 uniformly distributed time instants. Figures 8 and 9 show the capacitance voltage marginal pdf’s at the times t=1 ms. and t=5 ms. estimated on the basis of the suggested approach, using the first 5 order moments m0,v, m1,v, m2,v, m3,v , m4,v, together with a comparison with the correspondent pdf’s recovered by a conventional PCT, running ten thousands of Monte Carlo trials of the random variable ξVin in (26). Analogously to what previously illustrated in figure 5 and 6, figure 10 and 11 provides the multivariate analysis. In particular figure 10 gives the joint pdf relative to the voltages at t1=0.9 ms. and t2=1 ms. (see also figure 7 to visualize the time instants). Figure 11 repeats the same analysis at t1=1 ms. and t2=5 ms. In a similar fashion to case 1, the correlation coefficients are equal to 1.
IV. CONCLUSION This paper has focused on the analysis of measurements models expressed in terms of a stochastic ode. In particular, a novel approach to express the uncertainty for its solution has been suggested, based on the Wiener-Askey Polynomial Chaos theory. While the PC allows the calculation of the solution of the stochastic differential equation very quickly and faster than using a traditional Monte Carlo analysis, the results are anyhow available in terms of a polynomial expansion. This polynomial expansion allows a rapid evaluation of certain properties of the resulting pdf as the mean value or an estimation of the worst case, but it does not support a complete and simple reconstruction of the final pdf.
A simple approach to the problem is the evaluation of the solution by means of a suitable sampling of the random quantity ξ. This process is very computationally efficient but it looses the possibility to have an analytical solution which is one of the significant advantages of PC. The method here proposed, is less computationally efficient (on a regular PC it may take about 30 times more than the sampling process), but it preserve the analytical aspect of the solution. In conclusion, in the authors’ opinion, both approaches can be valuable to post-process the PC solution. The random sampling is very attractive for quick display of the results while the method proposed by this paper is very important for any situation where the analytical form is necessary. For example the analytical solution better support the possibility to calculate different values of the coverage factors. Another interesting example is given by the need of using the results of the analysis of the measurement process for a second analysis step where the measurement uncertainty has to be used as input. A concrete example of such a situation is a control design performed in the PC space that includes the measurement as source of uncertainty. Future work will focus on the application of the suggested methodology to the case in which more than one uncertain parameter is present in the system. It will be thus be assumed that the Polynomial Chaos solution of each equation is obtained starting from (3), which, according to a conventional approach, requires running a Monte Carlo of the whole random vector ξ. AKNOWLEDGEMENT This work has been partially supported by the US Office of Naval Research under the grant N00014-02-1-0623. REFERENCES [1] D. Xiu & Karniadakis, G. E. 2002b., “The Wiener-Askey polynomial chaos for stochastic differential equations”, SIAM J. Sci. Comput., vol. 24, pp. 619-644, 2002.
[2] Monti, A., Ponci, F., Lovet T., “A polynomial chaos theory approach to the control design of a power converter”, IEEE Annual Power Electronics Specialists Conference, Vol. 6, pp.48094813, June 2004. [3] Monti, A., Ponci, F., Lovet T., Smith A., Dougal R., “Modeling of uncertainty and applications in monitoring and control of power electronics”, Proc. of the American Control Conference, Vol3, pp.2011-2016, 8-10 June 2005. [4] D’Antona G., Monti A., Ponci F., Rocca L, Maximum Entropy Analytical Solution for Stochastic Differential Equations Based on the Wiener Askey Polynomial Chaos, Proc. 21st IEEE Advanced Methods for Uncertainty Estimation in Measurement, Vol.1 , 20-21 April 2006, pp. 62 - 66. [5] D’Antona, G., “Expanded uncertainty and coverage factor computation by higher order moments analysis”, Proc. 21st IEEE Instr. and Meas. Technology Conf., Vol.1 , 18-20 May 2004, p.234 - 238. [6] Brewer, J.W., “Kronecker product and matrix calculus in system theory”, IEEE Trans. on Circuits and Systems, vol. cas-25, No. 9, September 1978. [7] André,T.F., Nowak, R.D., Van Veen, B. D., “Low-Rank Estimation of Higher Order Statistics”, IEEE Trans. on Signal Processing, Vol.45, No.3, March 1997. [8] Carravetta F., Germani A., Raimondi M., “Polynomial Filtering for Linear Discrete Time non Gaussian Systems”, SIAM J. Control and Optimization, Vol. 34, No. 5, pp. 1666-1690, September 1996.
R Vin
C
Fig. 1. The RC circuit.
12
10
voltage [V]
8
6
4
2
0
0
0.002
0.004 0.006 time [s]
0.008
0.01
Fig. 2. First and second order moments for the capacitance voltage.
0.8 0.7 0.6
0.5 0.4 0.3 0.2 0.1 0
0
2
4
6 voltage [V]
8
10
Fig. 3. Capacitance voltage pdf at t=1 ms.
2.5
2
1.5
1
0.5
0
0
2
4
6 voltage [V]
8
10
Fig. 4. Capacitance voltage at t=3 ms.
2 1.5
error volatge at t2
1 0.5 0 -0.5 -1 -1.5 -2 -2
-1.5
-1
-0.5 0 0.5 error voltage at t1
1
1.5
2
Fig. 5 The contour plot for the joint pdf at t1 = 0.9 ms and t2 = 1 ms.
2 1.5
error volatge at t2
1 0.5 0 -0.5 -1 -1.5 -2 -2
-1.5
-1
-0.5 0 0.5 error voltage at t1
1
1.5
2
Fig. 6 The contour plot for the joint pdf shown at t1 = 1 ms and t2 = 3 ms.
12
10
voltage [V]
8
6
4
2
0
0
0.002
0.004 0.006 time [s]
0.008
0.01
Fig. 7. First and second order moments for the capacitance voltage .
0.6
0.5
0.4
0.3
0.2
0.1
0
0
2
4
6 8 voltage [V]
10
12
14
Fig. 8. Capacitance voltage pdf at t=1 ms.
0.3
0.25
0.2
0.15
0.1
0.05
0
0
2
4
6 8 voltage [V]
10
12
14
Fig. 9. Capacitance voltage pdf at t=5 ms.
2 1.5
error volatge at t2
1 0.5 0 -0.5 -1 -1.5 -2 -2
-1.5
-1
-0.5 0 0.5 error voltage at t1
1
1.5
2
Fig. 10 The contour plot for the joint pdf at t1 = 0.9 ms and t2 = 1 ms.
2.5 2 1.5
error volatge at t2
1 0.5 0 -0.5 -1 -1.5 -2 -2.5
-2
-1
0 error voltage at t1
1
2
Fig. 11 The contour plot for the joint pdf at t1 = 1 ms and t2 = 5 ms.
