[doc] Add the PBH test in `obsv` and `ctrb` docstring

[franckgaga]

I suggest to add the code for the PBH tests in the documentation to help the users.

While writing my package, I struggle to find a simple method to assess the observability of a linear system that is both numerically robust and universal, that is, also works on unstable systems (the gramian is only defined for stable A). There is many complex methods from what I found in the literature. A simple method that I just learned the existence is the PBH Test. That's what MATLAB uses internally to quickly check observability.

Keep up the good work!

Francis

[mzaffalon]

I (strongly) dislike the use of rank to determine the observability/controllability although every textbook and course I attended use it. rank computes the SVD and counts how many singular values are above a certain threshold, which is usually proportional to the type's machine precision, 1e-15 for Float64, hardly a relevant precision for engineering measurements! So perhaps adding a threshold for the rank computation would be a good idea.

Just to be clear, my comment has nothing to do with the quality of your contribution.

[franckgaga]

I agree that using rank is still not optimal, but, theoretically, it is a better method than checking the rank of the observability matrix, since it does not use powers of A. It's a great idea to add an explicit tolerance parameter.

[baggepinnen]

Thanks for the PR :) Would it make sense to perform balance_statespace or some similar numerical balancing before computing the rank here? If it can be made sensible and robust, we could perhaps include a function for it. Maybe one of the methods from this PR would be useful? https://github.com/andreasvarga/MatrixPencils.jl/commit/d86bdd5109e2c352da7706f5e184b0fa4426ff4a#diff-b2153cfad28eeac90ac98310abb0817e11261146ce234fd6c5bcb554d94dd3cdR1032

[baggepinnen]

Internally, minreal does perform a check for observability and/or controlability before deciding which modes to remove, what condition are they using? @andreasvarga do you have any opinion on what the numerically best method to check for observability/controlability is, that is applicable also for unstable systems?

[franckgaga]

Thanks for the PR :) Would it make sense to perform balance_statespace or some similar numerical balancing before computing the rank here? If it can be made sensible and robust, we could perhaps include a function for it. Maybe one of the methods from this PR would be useful? andreasvarga/MatrixPencils.jl@d86bdd5#diff-b2153cfad28eeac90ac98310abb0817e11261146ce234fd6c5bcb554d94dd3cdR1032

I'm not sure it's a good idea to transform the state-space equations beforehand. We want to test the controllability/observability of the original state-space representation, not the balanced one. Or maybe it would works? Since it both modifies A and C matrices in a way that keeps the I/O response identical ? @baggepinnen

[andreasvarga]

Which is the best method to check controllability? The short answer is: it depends!

For a yes/no check, probably the best method is the PBH test. This test is numerically very reliable. The reason is that the eigenvalues of A used in the test correspond to a slightly perturbed A, i.e., to an A+E with E of size eps()A and the SVD based rank test is thus very reliable. Moreover, the size of smallest singular value of [A-eI B] can be interpreted as a degree of controllability of the eigenvalue e. But, this test is not cheap! It has complexity O(n^4), where n is the order of A. Alternative methods have complexity O(n^3), the usual complexity for linear algebra computations.

For well-scalled models, the orthogonal reduction to a Kalman controllability form is an efficient method to check controllabiliy and to separate the uncontrollable part. The staircase reduction algorithms implemented in the MatrixPencils package are well-suited for this purpose.While numerically stable, this method employs a unique tolerance to determine the ranks of various submatrices resulting during the reduction. Using a single tolerance for potentially high number of rank determinations can be an issue, because different tolerances can provide different ranks, and at worse, there could not exist a single tolerance which ensures that all rank decisions are correct, excepting the single input case.

For stable systems, the best method is to evaluate the rank of the controllabiliy gramian, or even better, the rank of its Cholesky factor, which can be computed by methods available in the MatrixEquations package. Once again, the size of least singular value provides a quantitive measure of the degree of controllability. And, for yes/no checks, this method can be applied for unstable systems too, by simply shifting the eigenvalues of A to A-eI in the continuous time case or using a suitable scaling Ak for discrete time systems.

[franckgaga]

Now that's an interesting discussion! My conclusion to that would be: It would make sense to create two functions, iscontrollable and isobservable. They both receive the system, a tolerance and a method arguments. The first implemented method could be the PBH test since it's a one-liner. An we could add other methods afterward.

What's your opinion ?

[baggepinnen]

Andreas, thanks for the very detailed response!

It does indeed seem like it would be worthwhile to have a function for the PHB test.

If we are to turn it into a function. perhaps this function could return not only true / false, but also some additional info that was computed, seeing how this information can be useful and is very expensive to recompute for large systems (PDE discretizations etc.). The poles and rank of $[A - \lambda I \quad B]$ for each pole appears to be useful. The nullspace associated with the rank loss tells you what actuator directions you must span to make the system controllable, which could perhaps also be returned.

I assume that it's sufficient to compute the rank for only one of the poles in a complex conjugated pole pair?

[andreasvarga]

For real data, it is sufficient to check only for a complex eigenvalue in each complex conjugated pair.

[franckgaga]

well for now we could just merge the tips in the documentation, that's a start!

[franckgaga]

on another subject, is anyone here is going to adchem (ifac) 2024 toronto ?

[franckgaga]

good job @baggepinnen! I'll use them in my package :)