User:Egm6936.f10/Stochastic systems

The content of this section is mainly from Xiu 2010. For more details, one can refer to Xiu 2010, Chapter 4, p.44

After establishing a deterministic model for a physical system, prior to any simulation being taken, one must set up a proper stochastic model to study the effect of uncertainties in the inputs to the system. This becomes a key step in the stochastic system formulation. More work is devoted here to study some general aspects of this topic.

=Input parameterization=

How to parameterization the input random variables and reduce the dimension of these variables are two important aspects in stochastic processes. we will focus on these two aspects in this section.

Random parameters
When the random inputs to a system are the system parameters, the more important issue is to identify the independent parameters in the set. Such problem can be presented as:
 * $$\displaystyle X = (X_1,X_2,...,X_n), n > 1$$ is the system parameters with a prescribed distribution function $$\displaystyle F_X(x) = P(X \leq x), x \in \mathbb R^n$$, find a set of mutually independent random variables $$\displaystyle Y = (Y_1,Y_2,...,Y_k)\in \mathbb R^k, 1 \leq k \leq n $$, such that $$\displaystyle X = T(Y)$$ for a suitable transformation function $$\displaystyle T$$.

We use a simply example to demonstrate above idea. Suppose the inputs to a system are $$\displaystyle x$$ and $$\displaystyle y$$, and the output of the system is $$\displaystyle U(x,y)$$. Thus, the input random variables are $$\displaystyle X(\omega) = (x, y)$$. Surely, the output $$\displaystyle U(x,y)$$ depends on the input variables $$\displaystyle x$$ and $$\displaystyle y$$.
 * {| style="width:100%" border="0"

$$\displaystyle U(x,y): \Omega \rightarrow \mathbb R $$ (ab1) If $$\displaystyle x$$ and $$\displaystyle y$$ are mutually independent, then we simply let $$\displaystyle Y(\omega) = X(\omega)$$. The output field can be presented as
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$\displaystyle U(x,y): \mathbb R \times \mathbb R \rightarrow \mathbb R $$ (ab2) If $$\displaystyle x$$ and $$\displaystyle y$$ are not independent of each other, then there exist a function such that $$\displaystyle f(x,y) = 0$$. It is possible to find a random variable $$\displaystyle Y(\omega)$$ to parameterize the relation such that
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * $$\displaystyle x(\omega) = a(Y(\omega))$$, $$\displaystyle y(\omega) = b(Y(\omega))$$ and $$\displaystyle f(a,b) = 0$$.

Then, the random inputs via $$\displaystyle x$$ and $$\displaystyle y$$ can now be expressed by a single random parameter $$\displaystyle Y$$. The output field becomes
 * {| style="width:100%" border="0"

$$\displaystyle U(x,y): \mathbb R \rightarrow \mathbb R $$ (ab3)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

In practice, there are many input parameters in a system, finding the exact form of the functional dependence among all of them can be a challenging (and unnecessary) task. It is obvious when the only available information for the parameters is the joint probability distributions. Next we will focus on transforming the parameters to a set of independent random parameters by using their distribution functions

Gaussian parameters
When the system parameter has a Gaussian distribution, the problem becomes easy by Theorem 1. Let $$\displaystyle \mathbf X = (X_1,X_2,...,X_n)$$ be a random vector, and $$\displaystyle \mathbf X \doteqdot N(0, \mathbf C)$$, where $$\displaystyle \mathbf C \in \mathbb R^{n \times n}$$ is the covariance matrix and the expectation is assumed to be zero. Let $$\displaystyle \mathbf Z \doteqdot N(0, \mathbf I)$$, where $$\displaystyle \mathbf I $$ is the $$\displaystyle n \times n$$ identity matrix, be a uncorrelated Gaussian vector of size $$\displaystyle n $$. Let $$\displaystyle \mathbf A$$ be an $$\displaystyle n \times n$$ matrix; then according to Theorem 1, $$\displaystyle \mathbf {AZ} \doteqdot N(0, \mathbf {AA}^T)$$. Therefore, if one finds a matrix $$\displaystyle \mathbf A$$ such that $$\displaystyle \mathbf {AA}^T = \mathbf C$$, then $$\displaystyle \mathbf X = \mathbf {AZ}$$ will have the given distribution $$\displaystyle N(0,\mathbf C)$$. Since $$\displaystyle \mathbf C$$ is real and symmetric, solving the problem $$\displaystyle \mathbf {AA}^T = \mathbf C$$ can be readily done via Cholesky's decomposition.

Non-Gaussian parameters
Mathematically, for non-Gaussian parameters, we can utilize the Rosenblatt transformation to accomplish our goal. It's theoretically simply and powerful, but not easy to carry out in practice. Let $$\displaystyle \mathbf X = (X_1,X_2,...,X_n)$$ be a random vector with a (non-Gaussian) pdf $$\displaystyle F_X(x)=P(X \leq x)$$ and let $$\displaystyle \mathbf z = (z_1,z_2,...,z_n)= Tx = T(x_1,x_2,...,x_n)$$ be a transformation defined as
 * $$\displaystyle z_1 = P(X_1 \leq x_1) = F_1(x_1)$$,
 * $$\displaystyle z_2 = P(X_2 \leq x_2 \mid X_1 = x_1) = F_2(x_2 \mid x_1)$$,
 * $$\displaystyle z_n = P(X_n \leq x_n \mid X_{n-1} = x_{n-1},...,X_1 = x_1) = F_n(x_n \mid x_{n-1},...,x_1)$$,
 * $$\displaystyle z_n = P(X_n \leq x_n \mid X_{n-1} = x_{n-1},...,X_1 = x_1) = F_n(x_n \mid x_{n-1},...,x_1)$$,

Then, It can be shown that
 * $$\displaystyle P(Z_i \leq z_i; i = 1,2,...,n)$$
 * $$\displaystyle = \int_...\int d_{x_n}F_n(x_n \mid x_{n-1},...,x_1)...d_{y_1}F_1(x_1)$$
 * $$\displaystyle \int_{0}^{z_n}...\int_{0}^{z_1}d_{z_1}...d_{z_n} = \prod_{i=1}^{n}z_{i}$$,

where $$\displaystyle 0 \leq z_i \leq 1, i = 1,2,...,n)$$. Hence, $$\displaystyle \mathbf Z = (Z_1, Z_2,...,Z_n)$$ are independent and identically distributed random variables with uniform distribution in $$\displaystyle [0,1]^n$$.

In practice, numerical approximations of Rosenblatt transformation are required. And this topic is still understudied.

Dimension reduction
In many cases, the random inputs are stochastic processes. The inputs could be a time-dependent or space-dependent random variables, such as velocity and conductivity. We restate the inputs parameterization problem in stochastic process as follows.
 * Let $$\displaystyle (X_t, t\in \Omega)$$ be a stochastic process that models the random inputs. Find a appropriate transformation function $$\displaystyle T $$ such that $$\displaystyle X_t = T(Z)$$, where $$\displaystyle Z = (Z_1, Z_2,...,Z_d), 1 \leq d$$, are mutually independent.

Note that the index set $$\displaystyle \Omega$$ can be in either the time space or the space domain, and usually is an infinite-dimensional object. Since we require $$\displaystyle d $$ to be a finite integer, the transformation cannot be exact. $$\displaystyle X_t \approx T(Z)$$ in a proper norm or metric, and the accuracy of the approximation will be problem-dependent. A straightforward approach is discretize the index domain $$\displaystyle \Omega$$ into a set of indices and then study the process
 * $$\displaystyle (X_{t_1}, X_{t_2},..., X_{t_n}), t_1, t_2,...,t_n \in \Omega$$.

Above discretization is obviously an approximation. The finer this discretization, the better the approximation. However, this is not desired because a finer discretization leads to a larger dimension of $$\displaystyle n$$ and can significantly increase the computational burden. Some kinds of dimension reduction techniques are required to make the dimension as low as possible while maintaining a satisfactory accuracy.

Karhunen-Loeve expansion
The Karhunen-Loeve is one of the most widely used techniques for reducing dimension in representing random processes.

Let $$\displaystyle \mu_X(t) $$ be the mean of the input process $$\displaystyle X_t $$ and $$\displaystyle C(t,s) = cov(X_t,X_s)$$ be its covariance function. The Karhunen-Loeve expansion of $$\displaystyle X_t$$ is


 * {| style="width:100%" border="0"

$$\displaystyle X_t(\omega) = \mu_X(t)+\sum_{i=1}^{\infty}\sqrt{\lambda_i}\psi_i(t)X_i(\omega) $$ (ab4) where $$\displaystyle \psi_i $$ are the orthogonal eigenfunctions and $$\displaystyle \lambda_i$$ are the corresponding eigenvalues of the eigenvalue problem
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$\displaystyle \int_{T}C(t,s)\psi_i(s)ds = \lambda_i\psi_i(t), t\in T, $$ (ab5) and $$\displaystyle {X_i(\omega)}$$ are mutually uncorrelated random variables satisfying
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$\displaystyle \mathbb E[X_i] = 0, \;\; \mathbb E[X_iX_j] = \delta_{ij} $$ (ab6) and defined by
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$\displaystyle X_i(\omega) = \frac{1}{\sqrt{\lambda_i}}\int_{T}(X_t(\omega)-\mu_X(t))\psi_i(t)dt, \forall i $$ (ab7) In practice, one adopts a finite series expansion, e.g.,
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$\displaystyle X_t(\omega) = \mu_X(t)+\sum_{i=1}^{d}\sqrt{\lambda_i}\psi_i(t)X_i(\omega), 1\leq d $$ (ab8)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

When to truncate the series is problem-dependent. One important property of KL expansion is Eigenvalues $$\displaystyle \lambda_i$$ decay as index $$\displaystyle i$$ increase. For a given covariance function, the decay rate of the eigenvalues depends inversely on the correlation length. Long correlation length implies that the process is strongly correlated and results in a fast decay of the eigenvalues. The limit of this, infinitely long correlation length, is the fully correlated case where the eignevalues decay to zero immediatedly. Contrary, a weakly correlated process has short correlation length and results in a slow decay of the eignevalues. The limit of this, the uncorrelated process with zero correlation length, has no eigenvalue decay.

$$\displaystyle $$

Gaussian processes
For Gaussian random variables, they are uncorrelated with each other. Hence, we can directly use KL expansion to parameterize a Gaussian process by a finite number of independent Gaussian random variables.

Non-Gaussian processes
When the input random processes are non-Gaussian, their parameterization and dimension reduction are significantly more challenging. The main problem is that, for non-Gaussian distribution, uncorrelation of the random variables $$\displaystyle X_i$$ in (ab4) doesn't imply independence. Hence the KL expansion doesn't provide a way of parameterization with independent variables. Currently, there are not many practical methods for achieving the parameterization procedure. In many practical computations, one often still uses the KL expansion for an input process and then further assumes that the $$\displaystyle X_i$$ are independent. This field is still an active research topic.

 I'm a little confused about the KL expansion and input parameterization. As far as I understand, the KL expansion is to reduce the dimension of the stochastic processes, while the input parameterization is to parameterize the random variables such that they are independent with each other. A general procedure is that we adopt input paramerization techniques, such as Gaussian parameter and Rosenblatt transformation method, to make all the variables independent or uncorrelated. Then we use KL expansion to make the infinite stochastic process domain finite at certain accuracy level. Finally, solve the deterministic PDEs or ODEs. The KL expansion and parameterization are just two different steps. Am i right? Hylon.Chen 21:33, 1 February 2012 (UTC)

yes, these are different steps. you want to do the input parameterization first to define your model. then you use the KL expansion, with a truncation, to obtain a finite-dimensional system. solving the eigenvalue problem (ab5) would result in an infinite number of eigenpairs; only a few selected ones would be used in the finite series of the KL expansion in (ab8). Egm6322.s12 13:30, 2 February 2012 (UTC)

=Formulation of stochastic systems=

We briefly outline the main steps in formulating a stochastic system by taking into account random inputs to a well-established deterministic system. The general procedure is as follow:
 * 1. Parameterize the random input variables;
 * This step is done via some transformation techniques to make all inputs independent with each other.
 * 2. Reduce the dimension of the stochastic process at certain accuracy level;
 * KL expansion can be adopted in this step for uncorrelated random variables. As for correlated random variables, this topic is still active in recent research.
 * 3. Solve the deterministic equations system using traditional methods.

We introduce one partial differential equations (PDEs) as a basic example to illustrate above procedures for formulating a stochastic system. Let's consider a system of PDEs defined in a spatial domain $$\displaystyle D \subset \mathbb R^k,\;\; k = 1,2,3$$, and a time domain $$\displaystyle [0, T], \;\; T > 0 $$,
 * {| style="width:100%" border="0"

$$\displaystyle u_t(x,t,\omega) = \mathcal L (u), \;\; D\times (0,T]\times \Omega$$, $$\displaystyle \mathcal B (u) = 0, \;\; \partial D \times [0, T]\times \Omega$$, $$\displaystyle u = u_0, \;\; D \times {t = 0}\times \Omega, $$ (ab9)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Where $$\displaystyle \mathcal L $$ is a nonlinear differential operator, $$\displaystyle \mathcal B$$ is the boundary condition operator, $$\displaystyle u_0$$ is the initial condition, and $$\displaystyle \omega \in \Omega$$ denotes the random inputs of the system in a properly defined probability space $$\displaystyle (\Omega, \mathcal F, P) $$. The solution is therefore a random quantity,
 * {| style="width:100%" border="0"

$$\displaystyle u(x,t,\omega): \bar{D}\times[0, T]\times \Omega \rightarrow \mathbb R^{n_u} $$ (ab10) where $$\displaystyle n_u \geq 1$$ is the dimension of $$\displaystyle u$$. The random input to (ab9) take the form of random parameters and random processes. Let us assume that they can all be properly parameterized by a set of independent random variables using the techniques mentioned in the previous sections. Let $$\displaystyle Z = (Z_1, Z_2,...,Z_n) \in \mathbb R^d, d\geq 1$$, be the set of independent random variables characterizing the random inputs. The system (ab9) can be rewritten as
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$\displaystyle u_t(x,t,Z) = \mathcal L (u), \;\; D\times (0,T]\times \mathbb R^d $$, $$\displaystyle \mathcal B(u) = 0, \;\; \partial D \times [0, T]\times \mathbb R^d $$, $$\displaystyle u = u_0,\;\; D \times {t = 0}\times \mathbb R^d, $$ (ab11)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

The solution is now
 * {| style="width:100%" border="0"

$$\displaystyle u(x,t,Z): \bar{D}\times[0, T]\times \mathbb R^d \rightarrow \mathbb R^{n_u} $$ (ab12) Up to this point, the system (ab12) can be deterministically solved.
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

$$\displaystyle $$

=Traditional numerical methods=

We briefly review the traditional methods for solving practical systems with random inputs. Example used in above section will reused in this section for the purpose of illustration.

Monte Carlo Sampling
Monte Carlo sampling(MCS) is a statistical sampling method. The general procedure is as follows.
 * 1. Generate identically and independently distributed random numbers $$\displaystyle Z_i, \;\; i=1, 2,..., N $$, according to the pdf;


 * 2. For each $$\displaystyle i=1, 2,..., N $$, solve the governing equation;


 * 3. Estimate the required solution statistics.

It is obvious that the step 1 and step 3 are preprocessing and postprocessing steps, respectively. Only step 2 requires solution of the original problem, and it involves repetitive simulations of the deterministic counterpart of the problem.

Error estimate of MCS follows immediately from the central limit theorem. The widely adopted concept that the error convergence rate of MCS is inversely proportional to the square root of the number of realizations, which is relatively slow. One advantage of MCS is the convergence rate is independent of the dimension of the random space.

For more information about Monte Carlo Sampling, one can refer to wiki article Monte Carlo Sampling.

Moment equation approach
The objective of the moment equation approach is to compute the moments of the solution directly because in many cases the moments of the solution are what one needs. Consider a ordinary differential equation as
 * {| style="width:100%" border="0"

$$\displaystyle \frac{du}{dt}(t,\omega) = -\alpha(\omega)u,\;\;\; u(0,\omega) = \beta(\omega) $$ (ab13) Taking the expectation of both sides of (ab13) we obtain
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$\displaystyle \frac{d\mathbb E[u]}{dt}(t) = -\mathbb E[\alpha u],\;\;\; \mu(0) = \mathbb E[\beta] $$ (ab14) This requires a knowledge of $$\displaystyle \mathbb E[\alpha u]$$, which is unknown. We then attempt to derive an equation for this new quantity. Multiplying (ab13) by $$\displaystyle \alpha$$ and taking the expectation,
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$\displaystyle \frac{d}{dt}\mathbb E[\alpha u](t) = -\mathbb E[\alpha^2 u],\;\;\; \mathbb E[\alpha u](0) = \mathbb E[\alpha\beta] $$ (ab15) A new quantity $$\displaystyle \mathbb E[\alpha^2 u]$$ appears. In general, if we keep multiplying the original system by $$\displaystyle \alpha^k,\;\;\;k=0,1,2,... $$, then we get
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$\displaystyle \frac{}{dt}\mathbb E[\alpha^k u](t) = -\mathbb E[\alpha^{k+1} u],\;\;\; \mathbb E[\alpha^k u](0) = \mathbb E[\alpha^k\beta] $$ (ab16) The system of equations thus cannot be closed, since new variables keep being introduced into the exist system. This is the well-known closure problem. The typical approach to remedy the difficulty is to assume, for a fixed $$\displaystyle k$$, that the higher-order term is a function of the lower-order ones, that is,
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">
 * }
 * {| style="width:100%" border="0"

$$\displaystyle \mathbb E[\alpha^{k+1} u]= g(\mathbb E[u](0), \mathbb E[\alpha u](0), ..., \mathbb E[\alpha^k u](0)) $$ (ab17) where the form of the function $$\displaystyle g$$ is determined on a problem-dependent basis with sufficient justification. There exists, to date, no effective general strategy for solving the closure problem. Also, the error caused by most, if not all, closure techniques is not easy to quantify.
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">
 * }

Perturbation method
In the perturbation method, the random part of the solution is treated as a perturbation. The fundamental assumption is that such a perturbation is small. And this is typically achieved by requiring the standard deviation of the solution to be small.

One can refer to Perturbation method for more details about this method.

$$\displaystyle $$

=References=