Hello,

during my work on atomic clock ensembling, which uses a lot of Control toolbox, I have encountered grave numerical errors during plotting of state-space (@lti/@ss) Bode plots, traced to _freqresp_.m.

In MATLAB it gave good plots out of the same data. I haven't looked into MATLAB and not into a SLICOT and tried a fix using high-level algebraic functions. Original _freqresp_ uses a textbook formula for on-the-fly transformation of ss LTI system into an operator transfer function, i.e. c* (e*s - a)^-1 *b, where e is identity for regular systems, and something for stiff systems.

Here, the matrix inverse yielded a lot of rcond warnings. Although I can agree, that my system is "ugly", and there might be lot of nearly cancelling modes, it turned out to be possible to fix freqresp by calling eigendecomposition first, and then to just fiddle with inversions on diagonal. Besides the more stable behaviour, it shall be much faster, where there are (typically) lots of frequency points to evaluate.

I am attaching my test data, together with new freqresp and test script bodetst.m.

What I have not done:

1. support for stiff systems: I can not imagine the function for case of non-identity e; however, stiff systems are rare enough, so in worst case a revert to old formula could be used;

2. it seems to me silly to call C*diag(..)*D, but I am unsure about internal optimizations, JIT etc., so I do not know, whether to use this, for cycles, or something else to get fastest result.


Besides the ss/freqresp function, it may be interesting to consider improvement of ss2tf, ss2zpk conversions. It seems, that these functions are outsourced to some SLICOT code (in what I do not dare to even look). Freqresp using these converted systems are even much worse, than the original.

It seems to me, that the conversion into tf & zpk respectively could be better made by the same eigendecomposition, which I have used for freqresp. But I am not so experienced LTI guru to imagine all possible scenarios, I mean, whether the SLICOT ones could be better in some situation (but I doubt it could).

Any advanced appreciated.


Greetings,
Marek