jonathf / chaospy

Chaospy - Toolbox for performing uncertainty quantification.
https://chaospy.readthedocs.io/
MIT License
447 stars 87 forks source link

Continuous evaluation of the uncertainty using PCE #378

Open 2529481293 opened 2 years ago

2529481293 commented 2 years ago

Hi Jonathan,

Hope you are doing great! I have a question recently about implementing PCE about whether it is possible to achieve the continuous evaluation of the uncertainty over a time period for a dynamic system.

As shown in the attached fig, previously I have used the PCE method (Eq. 1) to evaluate the uncertainties of a dynamic system over discrete-time T = {t0, t1, t2,..., tn} through Chaospy, the results is very satisfactory. The \vec {\theta} is the uncertain parameters for the system that follows normal distribution; the \vec {\Phi_i} is the Hermite polynomials of i_th order; t is the time. for each time t in T, we conducted the PCE to obtain the coefficients \alpha_i(t) and uncertainties of f(t, \vec {\theta}) corresponding to time t.

I was wondering if it is possible to implement something like Eq. 2 to get the coefficients independent of the time t, so I can simply do interpolation of t to get the uncertainty for t at any time within the time period without evaluating PCE over and over again? I am not really sure if this is possible, though.

let me know your thoughts!

Best!

Formula

jonathf commented 2 years ago

Yeah, that should be possible. at least if you use point collocation method. But it isn't supported out of the box, so you need a little bit of tinkering.

Constructing the \Phi functions sould be straight forward, just make the distribution one dimentional higher to account for \Phi_j(t). Classical PCE already assumes that each Fourier term is seperable in each dimension. Just remember to choose the distribution where all values in T are represented. For simplicity of notation we place the dimension for t at the end.

Estimating the a_{ij} is a bit more involved. You need to:

  1. Sample from joint distribution [D+1, N] samples (one extra dimension for time).
  2. Use np.tile to repeat the samples T times (number of time steps) giving [D+1, N, T] samples.
  3. Fill inn samples[-1, :, :] with the time steps (as a mesh between t and \theta)
  4. Reshape to [D+1, N*T].

At the same time, you need to create evalutations.

  1. take the samples from (1).
  2. Remove the last dimension resulting in shape [D, N]
  3. Evaluate f N time resulting in output shape [N, T].
  4. Reshape to match (4): [N*T]

If everything is done right, you should be able to use cp.fit_regression to but it all together. Do a retall=True to get the Fourier coefficients a_{ij}. Just remember to reshape [N*T] back into [N, T] to get the ij indexing right.

This is a quick skimover on my part, so I appologize if I got the shapes incorrect. Hopefully you get the gist of it.

2529481293 commented 2 years ago

Hi Jonathan,

Thank you very much for the detailed reply! I tried ur suggestions and it went well until the last step to get the coefficients through cp.fit_regression, where I would need to specify the orthogonal polynomials generated through cp.orth_ttr(M,dist). M is the polynomial order and dist is the joint distribution of t and \theta. Based on the number of the coefficients, I think the problem become the following formula Eq. (3) instead of Eq. (2), right? (unless we can somehow specify the polynomials differently). I tried Eq. (3) with different orders, but it seems using Eq. (3) can't well approximate the function f though.

Let me know if I didn't understand your reply properly!

Thanks and regards!

image

jonathf commented 2 years ago

Well, yes that is the case, but in our example (2) and (3) are just semantic reformulations of eachother. Doing the reshaping between [N, T] and [N*T] is equivalent to going between (2) and (3), if you account for the new indexing and let \Phi_j(t) \Phi_i(\theta) = \Phi_n(t, \theta) when (i, j) and n refer to the same element.

2529481293 commented 2 years ago

Hi Jonathan,

The other thing I am not very sure about is whether it is correct to use the product of \Phi_j (t) and \Phi_i (\theta) in Eq. (2). Comparing the Eq. (2) and Eq. (3), the difference I can notice apart from what you've pointed out in your last reply is that the polynomials generated through Eq. (2) and Eq. (3) are quite different. I am not sure how this will influence the results, but according to the result Eq. (3) does a better job than Eq. (2) compared to the standard MC result. I think I will probably go with Eq. (3) instead of Eq. (2).

Thanks and Regards!

jonathf commented 2 years ago

In general (2) is a subset of (3), so yes, more complicated forms can be formed using the latter. But considering the original question, I have assumed (2) is what you wanted.

If you want a more general solution that covers (3) outside the space that (2) provides, things get complicated quite fast. Like e.g. canonical PCE itself makes an assumption that \Phi(\theta_0, \theta_1, \theta_2) can be decomposed into the dimensional parts \Phi(\theta_0) \Phi(\theta_1) \Phi(\theta_2). It can not be used to create anything outside of (2).