bugGNU Octave - Bugs: bug #63052, [octave forge] (control) Inverse...

 
 

bug #63052: [octave forge] (control) Inverse operation with transfer function variables resulting in wrong results

Submitter:  Luiz Antonio Maccari Junior <luiz>
Submitted:  Mon 12 Sep 2022 03:51:36 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Ready For Test Assigned to:  None
Originator Name:  Luiz Antonio Maccari Jr. Open/Closed:  * Open
Release:  * other Operating System:  * Microsoft Windows
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 26 Jan 2024 09:43:34 PM UTC, comment #13: 

The current dev branch of the control repo at github contains some hacks in @tf/__sys2tf__ similar to the one proposed by Doug in comment #6. With this, the example produces the same result for FT compared to FT2.

Torsten Lilge <ttl>
Group Member
Sat 15 Oct 2022 09:52:40 PM UTC, comment #12: 

This is definitely A control pkg issue, not a version of octave.


The Real issue is:
1)  some round off errors in a FORTRAN library
     which I can apply my eps hack to, to lessen the effect, of
     the round off errors.
2) the method that the OP used looks correct in text books,
    but is very susceptible to round off errors.

The normal way to do this is using:
[num,den] =ss2tf (A,B,C,D);
FT2 = tf (num,den)

always gives the correct answer. So I am hesitant to make any changes.

I tried many other examples and they worked ok. The OP just tried one that has round off errors

Any way the changes would only be in the control pkg

Doug Stewart <dastew>
Sat 15 Oct 2022 06:52:13 PM UTC, comment #11: 

it seems that this is an issue constrained to a package issue rather than a particular octave version.  I'm going to change the Release tag.
@Doug, is the change to _zero_.m ready to make a patch from or are we still seeking feedback on whether it works?

Nicholas Jankowski <nrjank>
Group Member
Wed 14 Sep 2022 07:46:53 PM UTC, comment #10: 

It seems apparent that this method

FT = C*inv(s*eye(3)-A)*B+D
FT= minreal (FT)

is sensitive to round off errors.


but this one is not;

[num,den] =ss2tf (A,B,C,D);
FT2 = tf (num,den)

Most users will use the second method, so I don't think that this is a high priority. but it should work.


Your question about getting different answers on different systems must be because of the different compilers for the Fortran library.

Doug Stewart <dastew>
Wed 14 Sep 2022 07:26:51 PM UTC, comment #9: 

Notice that when I said that FT was equal FT2 I was wrong.

Before I had

Transfer function 'FT' from input 'u1' to output ...

      2 s^3 + 17 s^2 + 332.4 s + 1222
 y1:  -------------------------------
         s^3 + 8 s^2 + 157 s + 580

Continuous-time model.

Transfer function 'FT2' from input 'u1' to output ...

      2 s^3 + 17 s^2 + 331 s + 1180
 y1:  -----------------------------
        s^3 + 8 s^2 + 157 s + 580

They were different.
 
Now with your modification I've got


Transfer function 'FT' from input 'u1' to output ...

      2 s^3 + 17 s^2 + 331.2 s + 1179
 y1:  -------------------------------
         s^3 + 8 s^2 + 157 s + 580

Continuous-time model.

Transfer function 'FT2' from input 'u1' to output ...

      2 s^3 + 17 s^2 + 331 s + 1180
 y1:  -----------------------------
        s^3 + 8 s^2 + 157 s + 580

Now the results are closer. I have tested on Control pkg 3.40 on linux. Could someone test the modification on windows?

Other question: How could the result being different using the same packages on different OS?

Luiz Antonio Maccari Junior <luiz>
Wed 14 Sep 2022 07:00:32 PM UTC, comment #8: 

on my comp.:

which @ss/__zero__
'@ss/__zero__' is a function from the file /home/doug/.local/share/octave/api-v57+/packages/control-3.4.0/@ss/__zero__.m

>>

Doug Stewart <dastew>
Wed 14 Sep 2022 06:45:28 PM UTC, comment #7: 

Could you give the location of the file zero.m ? I am looking for it, but i'm not finding. Is it part of the control package?

Luiz Antonio Maccari Junior <luiz>
Wed 14 Sep 2022 12:36:11 AM UTC, comment #6: 

I found that ,in _zero_.m
this line:
[zer, rank, infz, kronr, kronl] = _sl_ag08bd_ (sys.a, sys.e, sys.b, sys.c, sys.d, sys.scaled);

is returning zer close to eps, when it should return 0;

I added:

if(abs(zer)<2*eps)
     zer=0;
   endif

This definitely reduced (removed) the spurious wrong numbers, that I was seeing.

So the FORTRAN function _sl_ag08bd_  seems to have some round off problems.


Doug Stewart <dastew>
Tue 13 Sep 2022 01:01:16 AM UTC, comment #5: 

Wow that Matlab output looks totally wrong.

I think there is a problem in the control pkg.


clear
s=tf('s')

A=[ 0 0 0
    1 2 3
    0 0 0]

sii=[s 0 0
     0 s 0
     0 0 s]

sia=sii-A


This should be easy to do, but I get things like this

Transfer function 'sys' from input 'u9' to output ...

 y1:  0

 y2:  0

 y3:  1 s - 1.57e-16

I don't know where the - 1.57e-16 comes from --- it should be :

y3:   1 s

I am still working on it but help would be accepted :-)
To see the round off error show up set a break point in
mtimes at line 58 and look at sys before and after that line.

Doug Stewart <dastew>
Mon 12 Sep 2022 07:51:48 PM UTC, comment #4: 

what control and signal package versions did you have installed on each system?

here's what I get with Matlab 2022a:


FT =

  2 s^9 + 49 s^8 + 1359 s^7 + 2.025e04 s^6 + 274280 s^5 + 2.603e06 s^4 + 1.933e07 s^3 + 1.06e08 s^2 + 3.263e08 s + 3.97e08
  ------------------------------------------------------------------------------------------------------------------------
     s^9 + 24 s^8 + 663 s^7 + 9788 s^6 + 131931 s^5 + 1.249e06 s^4 + 9.25e06 s^3 + 5.096e07 s^2 + 1.584e08 s + 1.951e08

Continuous-time transfer function.


FT =

  2 s^9 + 49 s^8 + 1359 s^7 + 2.025e04 s^6 + 2.743e05 s^5 + 2.603e06 s^4 + 1.933e07 s^3 + 1.06e08 s^2 + 3.262e08 s + 3.97e08
  --------------------------------------------------------------------------------------------------------------------------
     s^9 + 24 s^8 + 663 s^7 + 9788 s^6 + 1.319e05 s^5 + 1.249e06 s^4 + 9.25e06 s^3 + 5.096e07 s^2 + 1.584e08 s + 1.951e08

Continuous-time transfer function.


FT2 =

  2 s^3 + 17 s^2 + 331 s + 1180
  -----------------------------
    s^3 + 8 s^2 + 157 s + 580

Continuous-time transfer function.


and on my Windows 10, Octave 7.2.0, control 3.4.0, signal 1.4.2:

>> FT = C*inv(s*eye(3)-A)*B+D
warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    __sys_connect__ at line 103 column 5
    mtimes at line 58 column 7

warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    __sys_connect__ at line 103 column 5
    mtimes at line 58 column 7
    minus at line 45 column 7

warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    __sys_connect__ at line 103 column 5
    mtimes at line 58 column 7
    minus at line 45 column 7

warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    __sys_connect__ at line 103 column 5
    mtimes at line 58 column 7

warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    __sys_connect__ at line 103 column 5
    mtimes at line 58 column 7


Transfer function 'FT' from input 'u1' to output ...

      2 s^2 + 80.76 s + 291
 y1:  ---------------------
      s^2 + 38.66 s + 142.5

Continuous-time model.
>> FT= minreal (FT)

Transfer function 'FT' from input 'u1' to output ...

      2 s^2 + 80.76 s + 291
 y1:  ---------------------
      s^2 + 38.66 s + 142.5

Continuous-time model.
>> [num,den] =ss2tf (A,B,C,D);
>> FT2 = tf (num,den)

Transfer function 'FT2' from input 'u1' to output ...

      2 s^3 + 17 s^2 + 331 s + 1180
 y1:  -----------------------------
        s^3 + 8 s^2 + 157 s + 580

Continuous-time model.
-verbatm-

Nicholas Jankowski <nrjank>
Group Member
Mon 12 Sep 2022 07:34:04 PM UTC, comment #3: 

Can you show me Matlabs answer?
I tried some examples from the internet and Octave gave the correct answers.

Doug Stewart <dastew>
Mon 12 Sep 2022 07:15:37 PM UTC, comment #2: 

This is interesting. Could the problem have been transported to the beta version (8.0.0)?

OBS: I have tried the same code on MATLAB and the result is different but also wrong.

Luiz Antonio Maccari Junior <luiz>
Mon 12 Sep 2022 04:11:50 PM UTC, comment #1: 


I just ran your code on Linux Ubuntu 22.04
Octave 8.0.0   id  f5755dbacd8d

I got the same as your widows results : different answers for
ft and ft2

Doug Stewart <dastew>
Mon 12 Sep 2022 03:51:36 PM UTC, original submission:  

Hello,

I was trying to obtain the transfer function from a third order state-space system by means of the equation C*inv(s*eye(3)-A)*B+D. However, the results in linux (GNU Octave 7.2) are correct and in windows 10 are wrong (also GNU Octave 7.2).

The test code is:

%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
clc;

pkg load control;
pkg load signal;

A = [-1 0 10;
         2 -4 6;
        -9 -8 -3];

B= [1;1;1];
C = [1 0 0];
D=2;

sis = ss(A,B,C,D);

s = tf([1 0],[1]);

FT = C*inv(s*eye(3)-A)*B+D
FT= minreal (FT)

[num,den] =ss2tf (A,B,C,D);
FT2 = tf (num,den)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In linux the result is:

warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    _sys_connect_ at line 103 column 5
    mtimes at line 58 column 7
    asdasd at line 20 column 4

warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    _sys_connect_ at line 103 column 5
    mtimes at line 58 column 7
    minus at line 45 column 7
    asdasd at line 20 column 4

warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    _sys_connect_ at line 103 column 5
    mtimes at line 58 column 7
    minus at line 45 column 7
    asdasd at line 20 column 4

warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    _sys_connect_ at line 103 column 5
    mtimes at line 58 column 7
    asdasd at line 20 column 4

warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    _sys_connect_ at line 103 column 5
    mtimes at line 58 column 7
    asdasd at line 20 column 4


Transfer function 'FT' from input 'u1' to output ...

      2 s^3 + 17 s^2 + 332.4 s + 1222
 y1:  -------------------------------
         s^3 + 8 s^2 + 157 s + 580

Continuous-time model.

Transfer function 'FT' from input 'u1' to output ...

      2 s^3 + 17 s^2 + 332.4 s + 1222
 y1:  -------------------------------
         s^3 + 8 s^2 + 157 s + 580

Continuous-time model.

Transfer function 'FT2' from input 'u1' to output ...

      2 s^3 + 17 s^2 + 331 s + 1180
 y1:  -----------------------------
        s^3 + 8 s^2 + 157 s + 580

Continuous-time model.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This result above is correct, FT = FT2.

In windows 10 the result is:

Continuous-time model.
warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    _sys_connect_ at line 103 column 5
    mtimes at line 58 column 7
    teste at line 23 column 4


warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    _sys_connect_ at line 103 column 5
    mtimes at line 58 column 7
    minus at line 45 column 7
    teste at line 23 column 4


warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    _sys_connect_ at line 103 column 5
    mtimes at line 58 column 7
    minus at line 45 column 7
    teste at line 23 column 4


warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    _sys_connect_ at line 103 column 5
    mtimes at line 58 column 7
    teste at line 23 column 4


warning: tf: converting to minimal state-space for MIMO TF interconnections
warning: called from
    _sys_connect_ at line 103 column 5
    mtimes at line 58 column 7
    teste at line 23 column 4


Transfer function 'FT' from input 'u1' to output ...


      2 s^2 + 80.76 s + 291
 y1:  ---------------------
      s^2 + 38.66 s + 142.5


Continuous-time model.


Transfer function 'FT' from input 'u1' to output ...


      2 s^2 + 80.76 s + 291
 y1:  ---------------------
      s^2 + 38.66 s + 142.5


Continuous-time model.


Transfer function 'FT2' from input 'u1' to output ...


      2 s^3 + 17 s^2 + 331 s + 1180
 y1:  -----------------------------
        s^3 + 8 s^2 + 157 s + 580


Continuous-time model.

%%%%%%%%%%%%%%%%%%

As one can see FT is different from FT2. It was expected the same result as the one obtained on linux.

Luiz Antonio Maccari Junior <luiz>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by ttl (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by dastew (Posted a comment)
  • -email is unavailable- added by luiz (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2024-01-26 ttl StatusNeed Info Ready For Test
    2022-10-15 nrjank Release7.2.0 other
    2022-09-12 nrjank StatusNone Need Info
        Summary(Control) Inverse operation with transfer function variables resulting in wrong results [octave forge] (control) Inverse operation with transfer function variables resulting in wrong results

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code