Repository navigation
Improve numerical stability in computation of Cholesky-factored controllability gramian #1123
Description
Activity
Thanks @mccullerlp. Is it possible to give a simple example that demonstrates this issue, so that we can include it in a unit test?
It may be, but will take me a while. I work on this sporadically and am principally working out of my developmental code https://github.com/wieldphysics/wield-control and my student's code "buzz" https://git.ligo.org/ian.macmillan/buzz/-/tree/dev?ref_type=heads.
(I'm not using python-control for a handful or reasons, mostly that I want support for descriptor state space due to atypical and diabolical numerical stability needs and have various philosophies on indexing).
Numerical issues can be tough to make unit tests for. In our case we see the jitter as we compute the Pareto front for an optimal controller that balances low and high frequency noise injections. I can export the systems of that Pareto front and we could compute the H2 norms using the gramians. I put this note here since I reviewed this implementation to speed up my use of the slycot routine and noticed the improvement in my implementation (where a unit test is much easier due to familiarity). I just wanted to post to upstream the fix, but can try to get you data from where we saw improvement.
Just an observation from utilizing sb030d in some models with extreme dynamic range (LIGO systems). There I had to resort to using Cholesky factored Grammians due to range issues in the output covariance matrix. In these models my A matrix is extremely non-normal and its Schur form has large off-diagonal elements even after balancing.
In this case, general orthogonal transformations are dangerous and the Schur form needs to be carefully preserved. So I noticed that the
python-control/control/statefbk.py
Line 1212 in f73e893
should actually use tra='T' for type='cf' and not transpose A. If you transpose A, then internally it becomes lower-triangular and the routine re-runs QR-algo to make Schur, which adds epsilon-errors all over and can smear out the response on sensitive systems. You can test this by keeping fact='F' and getting an error since it checks the triangularity.