[Octave-patch-tracker] [patch #8778] Control pkg: improvement of lti/ss/freqresp numerics, possibly also a way to better ss2tf

[anonymous]

URL:
  <http://savannah.gnu.org/patch/?8778>

                 Summary: Control pkg: improvement of lti/ss/freqresp
numerics, possibly also a way to better ss2tf
                 Project: GNU Octave
            Submitted by: None
            Submitted on: Thu 22 Oct 2015 05:08:23 PM UTC
                Category: None
                Priority: 5 - Normal
                  Status: None
                 Privacy: Public
             Assigned to: None
        Originator Email: Marek Peca <address@hidden>
             Open/Closed: Open
         Discussion Lock: Any

    _______________________________________________________

Details:

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




    _______________________________________________________

File Attachments:


-------------------------------------------------------
Date: Thu 22 Oct 2015 05:08:23 PM UTC  Name: lti_ss_freqresp_improve.png 
Size: 88kB   By: None

<http://savannah.gnu.org/patch/download.php?file_id=35263>
-------------------------------------------------------
Date: Thu 22 Oct 2015 05:08:23 PM UTC  Name: lti_ss_xfer_improve.tar.gz  Size:
12kB   By: None

<http://savannah.gnu.org/patch/download.php?file_id=35264>

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?8778>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/