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 Forge Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Need Info Assigned to:  None
Originator Name:  Luiz Antonio Maccari Jr. Open/Closed:  * Open
Release:  * 7.2.0 Operating System:  * Microsoft Windows
* Mandatory Fields

Add a New Comment (Rich Markup)
   

Comment Type & Canned Response:
       

 

( Jump to the original submission )

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>
Project 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 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 project members can vote.

     

     

     

     

    Follow 2 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    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.9