bugGNU Octave - Bugs: bug #66651, control package, MIMO systems...

 
 

bug #66651: control package, MIMO systems wrong created C matrix for state space representation using ss

Submitter:  None
Submitted:  Fri 10 Jan 2025 12:07:22 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  None
Status:  Fixed Assigned to:  ttl
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 7.1.0
Operating System:  * GNU/Linux Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Sun 19 Jan 2025 12:39:09 PM UTC, comment #5: 

I have revised the algorithm for computing a state space representation from a transfer function with this commit: https://github.com/gnu-octave/pkg-control/commit/9b4481b280829d6dbfd98cf312dd688b176015bb . The new algorithm uses an automated "manual" composition and looks for a minimal realization afterwards. The results of your example and other systems are as expected.

The new method will be included in the upcoming release 4.1.1. of the control package.

Thanks again for reporting the issue and providing the example.

 Torsten Lilge <ttl>
Group Member
Sat 11 Jan 2025 12:10:15 PM UTC, comment #4: 

Thanks for the code, now I see the problem. From your original post, where you were referring to a wrong C Matrix, I supposed that only the structure of the state space representations were different.

There are obviously some numerical issues for systems of this order. I will think about an improved method for ss().

 Torsten Lilge <ttl>
Group Member
Fri 10 Jan 2025 09:55:24 PM UTC, comment #3: 

i tried the provided code and it is still not working, i provide code with the correct step response using manual transformation to MIMO so you can test and make ss more robust





pkg load control
pkg load signal
clear

t=u1y1 =        tf([  4.74972290301299  0.0962858385264363  0.000679935548113776  3.72643701613439e-07  ],[  1  0.0207380329272262  0.000358651716393651  1.53426197558325e-06  1.60694141834117e-09  ]);
numu1y1_m1=t.num{1};
denu1y1_m1=t.den{1};


t=u1y2 =               tf([  -2.11493904790791e-05  -3.98828140784831e-07  -1.66475746406136e-08  -1.83748125520122e-10  -3.17372735779801e-12  -1.51387201191368e-14  -7.95515361277902e-18  ],[  1  0.0237440013869291  0.00093839286173024  1.36557230877558e-05  2.17818987171716e-07  1.82722920274503e-09  5.29889671308072e-12  4.97347909522557e-15  ]);
numu1y2_m1=t.num{1};
denu1y2_m1=t.den{1};


t=     u2y1=             tf([  0.0337260175091407  0.414338784602062  0.00343501567979344  0.00132279342826428  8.05249078940347e-06  1.39424444953855e-06  5.51151054539321e-09  5.10115766190411e-10  9.91853470677749e-13  2.88997553994485e-14  ],[  1  0.153889837566522  0.00677224375967168  0.000548671957421157  1.41722799246167e-05  6.85514842206994e-07  1.19614786183711e-08  3.46611370329274e-10  3.77280525883548e-12  5.76999094807308e-14  2.00325109640768e-16  ]);
numu2y1_m1=t.num{1};
denu2y1_m1=t.den{1};


t=   u2y2 =             tf([  -1.48787383837688e-06  -3.24705706205277e-08  -6.54751170609613e-10  -7.46171917997884e-12  ],[  1  0.0134226252067825  0.000316273404177729  2.74708154815059e-06  5.9854886873292e-09  ]);
numu2y2_m1=t.num{1};
denu2y2_m1=t.den{1};


t=u3y1 =        tf([  -8318.14771695403  -1005.17247836674  -69.388714502427  -2.13567218084595  -0.057214617650067  -0.00051673284923126  -5.21707731497307e-06  -3.10877494268659e-09  ],[  1  0.127298409756932  0.0105061762139007  0.000356618210494435  1.16880708107467e-05  1.51741412632152e-07  2.05485336227224e-09  2.72630351738011e-12  ]);
numu3y1_m1=t.num{1};
denu3y1_m1=t.den{1};


t=u3y2 =       tf([  6.58550303958022e-07  1.12032568300048e-06  4.84240826234966e-08  2.06589443851401e-09  3.25197618527417e-11  2.6086983962496e-14  ],[  1  0.0340866722611228  0.0019478738998518  3.09342222656595e-05  4.969573841421e-07  2.19170977290605e-09  2.86190017229284e-12  ]);
numu3y2_m1=t.num{1};
denu3y2_m1=t.den{1};

%tf matrix
sys=[u1y1 u2y1 u3y1 ;
u1y2 u2y2 u3y2];

sys_ss=ss(sys);%non working ss
A=sys_ss.a;
B=sys_ss.b;
C=sys_ss.c;
D=sys_ss.d;

%SISO ss
[A_u1y1,B_u1y1,C_u1y1,D_u1y1]=tf2ss(u1y1)
[A_u2y1,B_u2y1,C_u2y1,D_u2y1]=tf2ss(u2y1)
[A_u3y1,B_u3y1,C_u3y1,D_u3y1]=tf2ss(u3y1)
[A_u1y2,B_u1y2,C_u1y2,D_u1y2]=tf2ss(u1y2)
[A_u2y2,B_u2y2,C_u2y2,D_u2y2]=tf2ss(u2y2)
[A_u3y2,B_u3y2,C_u3y2,D_u3y2]=tf2ss(u3y2)

%manual transform to MIMO ss
A_ss = blkdiag(A_u1y1, A_u2y1, A_u3y1, A_u1y2, A_u2y2, A_u3y2)
B_ss = [ blkdiag(B_u1y1, B_u2y1, B_u3y1) ; blkdiag(B_u1y2, B_u2y2, B_u3y2) ]
C_ss = blkdiag([C_u1y1, C_u2y1, C_u3y1] , [C_u1y2, C_u2y2, C_u3y2])
D_ss = [D_u1y1, D_u2y1, D_u3y1 ; D_u1y2, D_u2y2, D_u3y2 ]

close all
s=ss(A_ss,B_ss,C_ss,D_ss)
step(s)%correct step response


%provided code

# Similarity trnasformation
[T, lambda] = eig (sys_ss.a);
A = inv(T)*sys_ss.a*T;
B = inv(T)*sys_ss.b;
C = sys_ss.c*T;
# Set very small values to 0 (finite
# numerical accuracy) for a nicer output
A(abs(A)<1e-12)= 0
B(abs(B)<1e-12)= 0
C(abs(C)<1e-12)= 0

figure
s=ss(A,B,C,D)
step(s)%incorrect step response

 Anonymous
Fri 10 Jan 2025 08:32:20 PM UTC, comment #2: 

There exists an indefinite number of state space representations with identical transfer function. You are expecting a diagonal form. You can try to transform the result of ss() into a diagonal form by

# Similarity trnasformation
[T, lambda] = eig (sys_ss.a);
A = inv(T)*sys_ss.a*T;
B = inv(T)*sys_ss.b;
C = sys_ss.c*T;
# Set very small values to 0 (finite
# numerical accuracy) for a nicer output
A(abs(A)<1e-12)= 0
B(abs(B)<1e-12)= 0
C(abs(C)<1e-12)= 0

However, this method is only possible for systems with single eigenvalues.

Could you please try this and report back whether you get the expected results? The order of non-zero entries in B and C might differ because the order of the subsystems in A follows the magnitude of their eigenvalues.

 Torsten Lilge <ttl>
Group Member
Fri 10 Jan 2025 06:15:42 PM UTC, comment #1: 


> ... i dont have access to your source code ...

Au contraire, you do. Octave is an OSS project, remember?  :-)
You can type 'which ss' in the octave command window; it'll tell you where the file ss.m can be found, and if you do 'edit ss' you'll have the source code in the editor.

Starting from https://packages.octave.org, the source code of 'ss' with all its helper functions is easily found:

control package repository:
https://github.com/gnu-octave/pkg-control

'ss':
https://github.com/gnu-octave/pkg-control/tree/main/inst/%40ss

Turns out to be an 'old-style' class, but ss.m is here:
https://github.com/gnu-octave/pkg-control/blob/main/inst/%40ss/ss.m

Philip Nienhuis <philipnienhuis>
Group Member
Fri 10 Jan 2025 12:07:22 PM UTC, original submission:  

thanks for the control package


i have a MIMO system like this, where u1y1 u2y1 etc are transfer functions:

sys=[u1y1 u2y1 u3y1 ;
u1y2 u2y2 u3y2];

sys_ss=ss(sys);

this results in the wrong state space representation of the system, i dont have access to your source code but from the result of your function it seems the C matrix is created wrongly


it should be created like this:

[A_u1y1,B_u1y1,C_u1y1,D_u1y1]=tf2ss(u1y1)
[A_u2y1,B_u2y1,C_u2y1,D_u2y1]=tf2ss(u2y1)
[A_u3y1,B_u3y1,C_u3y1,D_u3y1]=tf2ss(u3y1)
[A_u1y2,B_u1y2,C_u1y2,D_u1y2]=tf2ss(u1y2)
[A_u2y2,B_u2y2,C_u2y2,D_u2y2]=tf2ss(u2y2)
[A_u3y2,B_u3y2,C_u3y2,D_u3y2]=tf2ss(u3y2)


A_ss = blkdiag(A_u1y1, A_u2y1, A_u3y1, A_u1y2, A_u2y2, A_u3y2)
B_ss = [ blkdiag(B_u1y1, B_u2y1, B_u3y1) ; blkdiag(B_u1y2, B_u2y2, B_u3y2) ]
C_ss = blkdiag([C_u1y1, C_u2y1, C_u3y1] , [C_u1y2, C_u2y2, C_u3y2])
D_ss = [D_u1y1, D_u2y1, D_u3y1 ; D_u1y2, D_u2y2, D_u3y2 ]




Anonymous

 

(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)

    •  

    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 5 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2025-01-19 ttl StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2025-01-11 ttl Assigned toNone ttl
    2025-01-11 ttl StatusNeed Info Confirmed
    2025-01-10 ttl StatusNone Need Info

    Back to the top

    Copyright © 2025 Free Software Foundation, Inc.
    Verbatim copying and distribution of this entire article is permitted in any medium, provided this notice is preserved.
    The Levitating, Meditating, Flute-playing Gnu logo is a GNU GPL'ed image provided by the Nevrax Design Team.
    Page source code

    Powered by Savane 3.14-708e.
    Corresponding source code