8000 Stochastic systems additions by murrayrm · Pull Request #714 · python-control/python-control · GitHub
[go: up one dir, main page]
Skip to content

Stochastic systems additions #714

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom
Merged
merged 13 commits into from
Apr 2, 2022
Merged

Stochastic systems additions #714

merged 13 commits into from
Apr 2, 2022

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 con 10000 ditions 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
0