Skip to content

Stochastic systems additions#714

Merged
murrayrm merged 13 commits intopython-control:masterfrom
murrayrm:stochsys
Apr 2, 2022
Merged

Stochastic systems additions#714
murrayrm merged 13 commits intopython-control:masterfrom
murrayrm:stochsys

Conversation

@murrayrm
Copy link
Member

@murrayrm murrayrm commented Mar 23, 2022 •
edited
Loading

This PR adds new functionality for supporting analysis and simulation of stochastic systems, based on a class that I taught at Caltech:

  • Adds two new functions supporting random signals: white_noise, which creates a white noise vector in continuous or discrete time, and correlation, which calculates the correlation function (or [cross-] correlation matrix), R(tau).

  • Adds a new function create_estimator_iosystem that matches the style of create_statefbk_iosystem (I/O system enhancements #710) and creates an I/O system implementing an estimator (including covariance update).

  • Adds the ability to specify initial conditions for input_output_response as a list of values, so that for estimators that keep track of covariance you can set the initial conditions as [X0, P0]. In addition, if you specify a fewer number of initial conditions than the number of states, the remaining states will be initialized to zero (with a warning if the last initial condition is not zero). This allows the initial conditions to be given as [X0, 0].

  • Adds the ability to specify inputs for input_output_response as a list of variables. Each element in the list will be treated as a portion of the input and broadcast (if necessary) to match the time vector. This allows input for a system with noise as [U, V] and inputs for a system with zero noise as [U, np.zero(n)] (where U is an input signal and np.zero(n) gets broadcast to match the time vector).

  • Added some new Jupyter notebooks demonstrate the use of these functions:

    • stochresp.ipynb: illustrates the implementation of random processes and stochastic response.

    • pvtol-outputfbk.ipynb: illustrates the implementation of an extended Kalman filter and the use of the estimated state for LQR feedback of a vectored thrust aircraft model.

    • kincar-fusion.ipynb: works through estimation of the state of a car changing lanes with two different sensors available: one with good longitudinal accuracy and the other with good lateral accuracy.

@murrayrm

This comment was marked as outdated.

Sign in to view
@murrayrm
Copy link
Member Author

murrayrm commented Apr 2, 2022

This PR adds new functionality and doesn't change any old functionality, so going ahead and merging (since it has been sitting for a while).

@murrayrm murrayrm merged commit 983726c into python-control:master Apr 2, 2022
@murrayrm murrayrm deleted the stochsys branch April 16, 2022 01:03
@murrayrm murrayrm added this to the 0.9.2 milestone Apr 17, 2022
@jbayless
Copy link

jbayless commented Jul 8, 2024 •
edited
Loading

Sorry to jump in a little late here...

Doesn't white_noise scale Q the wrong way for a continuous-time system? Right now it scales the standard deviation by 1/sqrt(dt). But it should scale by sqrt(dt) for the integral of the covariance to be Q over the simulation sample time (hence the noise integral from 0 to dt, and the std deviation, should decrease to 0 as dt-->0, instead of increasing to infinity).

Also, since this is generating multivariate noise, why not allow Q to be a covariance matrix with correlations between variables? This can be done by eg:

Q = np.array(Q)
if Q.ndim == 0:
    # scalar
    Q = np.ones((1,1), dtype = np.float64)*Q
    
elif Q.ndim == 1:
    # uncorrelated signals: diagonal covariance matrix
    Q = np.diag(Q)
    
elif Q.ndim == 2:
    # full covariance matrix
    if Q.shape[0] != Q.shape[1]:
        raise ValueError("Noise covariance matrix is not square!")
else:
    # error
    raise ValueError("Noise covariance matrix has >2 dimensions!")

means = np.zeros((Q.shape[0],), dtype = np.float64)
W = rng.multivariate_normal(means, Q*np.sqrt(dt), size = T.size, method = "cholesky")

@murrayrm
Copy link
Member Author

murrayrm commented Jul 8, 2024

The current code handles the multi-variate case, though in the case of a diagonal $Q$ you have to pass call np.diag yourself (which seems fine).

In terms of the scaling: we are essentially approximating the continuous white noise process by a piecewise constant (discrete) noise process. You have to scale the discrete-time signal by $1/\sqrt{dt}$ since the covariance is $R(t_1, t_2) = \mathbb{E}(W(t_1) W^\mathtt{T}(t_2)) = Q \delta(t_2 - t_1)$ (similar to the way that an impulse function can be approximated as a pulse of height $1/dt$ so that the integral is 1).

You can find an example of using these functions (and some insight into the scaling) here: https://github.com/python-control/python-control/blob/main/examples/stochresp.ipynb (Google Colab version).

@jbayless
Copy link

jbayless commented Jul 8, 2024 •
edited
Loading

OK, I see. You mean that you're scaling it such that after you integrate the output of white_noise, the result will be scaled correctly. I read the documentation a different way: that white_noise is scaled such that each sample of the output represents the integral of the noise over the previous dt. (Similar to the way that the integral of a pulse of height 1 has an magnitude of dt). But that was a misreading on my part.

(I also overlooked the last line of your function!).

Never mind my comment then!

@murrayrm murrayrm modified the milestones: 0.9.2, 0.10.1 Aug 8, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

No reviews

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

0.10.1

Development

Successfully merging this pull request may close these issues.

2 participants

@murrayrm @jbayless