bugGNU Octave - Bugs: bug #58538, [octave forge] (control) step and...

 
 

bug #58538: [octave forge] (control) step and lsim fail on high order transfer functions

Submitter:  None
Submitted:  Wed 10 Jun 2020 05:16:56 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  3 - Low Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  gs Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * dev
Operating System:  * GNU/Linux Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sat 03 Jul 2021 07:57:21 PM UTC, comment #9: 

The problem for the impulse response would be much harder to solve. The impulse response requires the discretization of the system in tf-form. Although the discretizaiton itself is done in state space representation, computing the tf from the discretized system can not be done with sufficient accuracy for very high order systems. This effect can also be seen in other tools, see comment #5.

Closing this report since the original issues with lsim() and step() are fixed.

Torsten Lilge <ttl>
Group Member
Fri 19 Jun 2020 07:49:27 PM UTC, comment #8: 

Changeset 1618:9bfa416611fd solves the issue for lsim ().

Torsten Lilge <ttl>
Group Member
Fri 19 Jun 2020 07:19:10 PM UTC, comment #7: 

I have pushed changeset 1617:104e9c0fc974 which solves the problem for step () by converting a continuous-time system into a state space representation before using c2d for discretization.

The issue is still present in lsim() and impulse(). In the latter case, the discretization is done using the "impulse" option of c2d, which obviously leads to an inaccurate discrete-time model even for a state space representation.

Torsten Lilge <ttl>
Group Member
Wed 10 Jun 2020 08:52:04 PM UTC, comment #6: 

Nicholas, thanks for posting.
Just to avoid confusions: I changed the order of the two cases before posting here. So the results posted by in comment #5 are in different order compared to my comment #2.

Torsten Lilge <ttl>
Group Member
Wed 10 Jun 2020 08:31:20 PM UTC, comment #5: 

(adding details from the mailing list)
running step and lsim commands in comment #0  result in the following warning being displayed repeatedly (for step, it appears every time the figure is given focus until a ctrl-c)

warning: opengl_renderer: data values greater than float capacity.  (1) Scale data, or (2) Use gnuplot

switching to gnuplot results in the output attached as file OctaveControlGnuplot.png

looking at the output from Torsten's tests:

output on Octave:

>> pkg load control
>> q=tf([  -1.46668928632582  38.5234944088261  -1023.42433722006  18055.8729050165  -275471.128598611  3680685.32154243  -40413629.6869944  433109318.542883  -3699840188.09672  32827566143.7738  -225798332168.488  1689921426466.71  -9494747382069.7  60613579740811.5  -278626623858487  1.52747193482069e+15  -5.6883783399098e+15  2.69067309068308e+16  -7.92226372994905e+16  3.25360727881732e+17  -7.2059792475876e+17  2.61175754638012e+18  -3.8909506684981e+18  1.31637705830625e+19  -9.12640045084024e+18  3.78290213329238e+19  1.41630594163597e+19  5.02932079881115e+19  1.29547036349465e+20  1.05983336007477e+19  2.43633391646665e+20  -1.31208497581345e+19  1.04713601706961e+20  7.83465823853832e+18  ],[  1  17.1698120499521  542.934121038049  7294.98595152311  126155.971925398  1377306.37948524  16829405.6386297  152789570.438515  1448656865.6342  11095921096.5808  85365123037.7488  556372911296.259  3555399791580.02  19799716986834.5  106353561602035  506292672690482  2.29624277698132e+15  9.31618125768562e+15  3.56349262092232e+16  1.22398996473335e+17  3.92373264437545e+17  1.12860046096021e+18  2.99632129507514e+18  7.09934073973815e+18  1.53176948230893e+19  2.91760427412351e+19  4.97010301791734e+19  7.33491370510952e+19  9.42777562297368e+19  1.01604047339533e+20  9.12510568398965e+19  6.44576522525933e+19  3.46367980244933e+19  1.23309424846899e+19  2.4145011614912e+18  ]);
>> sys_ss_ct = ss (q);
>> sys_dt2   = c2d (sys_ss_ct, 0.01);
>> eig (sys_dt2.a)
ans =

   0.99422 + 0.08037i
   0.99422 - 0.08037i
   0.99346 + 0.07615i
   0.99346 - 0.07615i
   0.99333 + 0.07155i
   0.99333 - 0.07155i
   0.99368 + 0.06688i
   0.99368 - 0.06688i
   0.99488 + 0.00426i
   0.99488 - 0.00426i
   0.99458 + 0.00842i
   0.99458 - 0.00842i
   0.99435 + 0.01262i
   0.99435 - 0.01262i
   0.99443 + 0.01718i
   0.99443 - 0.01718i
   0.99379 + 0.02253i
   0.99379 - 0.02253i
   0.99340 + 0.02768i
   0.99340 - 0.02768i
   0.99335 + 0.03249i
   0.99335 - 0.03249i
   0.99329 + 0.03711i
   0.99329 - 0.03711i
   0.99324 + 0.04174i
   0.99324 - 0.04174i
   0.99327 + 0.04646i
   0.99327 - 0.04646i
   0.99345 + 0.05142i
   0.99345 - 0.05142i
   0.99368 + 0.05637i
   0.99368 - 0.05637i
   0.99406 + 0.06152i
   0.99406 - 0.06152i

>> sys_tf_dt = c2d (q, 0.01);
>> sys_dt1   = ss (sys_tf_dt);
>> eig (sys_dt1.a)
ans =

   1.42735 + 0.92459i
   1.42735 - 0.92459i
   1.17401 + 0.91528i
   1.17401 - 0.91528i
   0.97044 + 0.84467i
   0.97044 - 0.84467i
   0.81966 + 0.74285i
   0.81966 - 0.74285i
   0.70737 + 0.63075i
   0.70737 - 0.63075i
   0.61757 + 0.53082i
   0.61757 - 0.53082i
   0.55208 + 0.43837i
   0.55208 - 0.43837i
   0.50582 + 0.34596i
   0.50582 - 0.34596i
   0.47599 + 0.25501i
   0.47599 - 0.25501i
   0.45893 + 0.16733i
   0.45893 - 0.16733i
   0.45075 + 0.08282i
   0.45075 - 0.08282i
   0.44843 + 0.00000i
   0.79468 + 0.52856i
   0.79468 - 0.52856i
   0.64503 + 0.00000i
   0.93407 + 0.00000i
   1.70651 + 0.84064i
   1.70651 - 0.84064i
   1.96533 + 0.64519i
   1.96533 - 0.64519i
   2.14777 + 0.35049i
   2.14777 - 0.35049i
   2.21291 + 0.00000i

>> diary off



output on Matlab 2020a:

q=tf([  -1.46668928632582  38.5234944088261  -1023.42433722006  18055.8729050165  -275471.128598611  3680685.32154243  -40413629.6869944  433109318.542883  -3699840188.09672  32827566143.7738  -225798332168.488  1689921426466.71  -9494747382069.7  60613579740811.5  -278626623858487  1.52747193482069e+15  -5.6883783399098e+15  2.69067309068308e+16  -7.92226372994905e+16  3.25360727881732e+17  -7.2059792475876e+17  2.61175754638012e+18  -3.8909506684981e+18  1.31637705830625e+19  -9.12640045084024e+18  3.78290213329238e+19  1.41630594163597e+19  5.02932079881115e+19  1.29547036349465e+20  1.05983336007477e+19  2.43633391646665e+20  -1.31208497581345e+19  1.04713601706961e+20  7.83465823853832e+18  ],[  1  17.1698120499521  542.934121038049  7294.98595152311  126155.971925398  1377306.37948524  16829405.6386297  152789570.438515  1448656865.6342  11095921096.5808  85365123037.7488  556372911296.259  3555399791580.02  19799716986834.5  106353561602035  506292672690482  2.29624277698132e+15  9.31618125768562e+15  3.56349262092232e+16  1.22398996473335e+17  3.92373264437545e+17  1.12860046096021e+18  2.99632129507514e+18  7.09934073973815e+18  1.53176948230893e+19  2.91760427412351e+19  4.97010301791734e+19  7.33491370510952e+19  9.42777562297368e+19  1.01604047339533e+20  9.12510568398965e+19  6.44576522525933e+19  3.46367980244933e+19  1.23309424846899e+19  2.4145011614912e+18  ]);
sys_ss_ct = ss (q);
sys_dt2   = c2d (sys_ss_ct, 0.01);
eig (sys_dt2.a)

ans =

   0.9942 + 0.0804i
   0.9942 - 0.0804i
   0.9935 + 0.0762i
   0.9935 - 0.0762i
   0.9933 + 0.0716i
   0.9933 - 0.0716i
   0.9937 + 0.0669i
   0.9937 - 0.0669i
   0.9941 + 0.0615i
   0.9941 - 0.0615i
   0.9937 + 0.0564i
   0.9937 - 0.0564i
   0.9935 + 0.0514i
   0.9935 - 0.0514i
   0.9933 + 0.0465i
   0.9933 - 0.0465i
   0.9932 + 0.0417i
   0.9932 - 0.0417i
   0.9933 + 0.0371i
   0.9933 - 0.0371i
   0.9934 + 0.0325i
   0.9934 - 0.0325i
   0.9934 + 0.0277i
   0.9934 - 0.0277i
   0.9938 + 0.0225i
   0.9938 - 0.0225i
   0.9944 + 0.0172i
   0.9944 - 0.0172i
   0.9943 + 0.0126i
   0.9943 - 0.0126i
   0.9946 + 0.0084i
   0.9946 - 0.0084i
   0.9949 + 0.0043i
   0.9949 - 0.0043i

sys_tf_dt = c2d (q, 0.01);
sys_dt1   = ss (sys_tf_dt);
eig (sys_dt1.a)

ans =

   2.0691 + 0.0000i
   2.0125 + 0.2961i
   2.0125 - 0.2961i
   1.8581 + 0.5393i
   1.8581 - 0.5393i
   1.6475 + 0.6911i
   1.6475 - 0.6911i
   1.4552 + 0.7483i
   1.4552 - 0.7483i
   1.2923 + 0.8011i
   1.2923 - 0.8011i
   1.1007 + 0.8116i
   1.1007 - 0.8116i
   0.9311 + 0.7639i
   0.9311 - 0.7639i
   0.8008 + 0.6860i
   0.8008 - 0.6860i
   0.7056 + 0.5990i
   0.7056 - 0.5990i
   0.6352 + 0.5112i
   0.6352 - 0.5112i
   0.5848 + 0.4239i
   0.5848 - 0.4239i
   0.5529 + 0.3406i
   0.5529 - 0.3406i
   0.5346 + 0.2680i
   0.5346 - 0.2680i
   0.5163 + 0.2043i
   0.5163 - 0.2043i
   0.4994 + 0.1373i
   0.4994 - 0.1373i
   0.4897 + 0.0687i
   0.4897 - 0.0687i
   0.4866 + 0.0000i

diary off



Nicholas Jankowski <nrjank>
Group Member
Wed 10 Jun 2020 08:12:55 PM UTC, comment #4: 

comment #1:

> So, to clarify because it's not entirely clear from your bug description, the bug you are reporting is that 'step(sys)' tries to plot too many points and overloads the Octave OpenGL plotting system, while 'step(ss(sys))' works as a workaround and does plot the correct figure. Is that accurate?



step() and lsim() create simulations that are completely wrong, the simulations seem to be unstable which create very large values and the plot window errors.

so 2 problems:
1: simulation is incorrect, gives incorrect values, not even close
2. plot window can't handle very large (probably inf) values

Anonymous
Wed 10 Jun 2020 07:50:44 PM UTC, comment #3: 

Changing from "Matlab Compatibility" to "Incorrect Result" and raising severity, since this is not just a problem in plotting the result.

Torsten Lilge <ttl>
Group Member
Wed 10 Jun 2020 07:47:25 PM UTC, comment #2: 

No, there is a problem with the internal system representation used by step (). As posted in the help list, this seems to be related to a problem when discretizing a continuous-time transfer function. To test:


# Get q from the original post

# discretization of the transfer function
# unstable eigenvalues (|lambda| > 1)
sys_tf_dt = c2d (q, 0.01);
sys_dt1   = ss (sys_tf_dt);
lambda = eig (sys_dt1.a)
stable = isstable (sys_dt1)

# discretization of the state space model
# => stable eigenvalues
sys_ss_ct = ss (q);
sys_dt2   = c2d (sys_ss_ct, 0.01);
lambda = eig (sys_dt2.a)
stable = isstable (sys_dt2)


I have attached a script that reads the transfer function q and then does the test given above.

(file #49237)

Torsten Lilge <ttl>
Group Member
Wed 10 Jun 2020 07:25:29 PM UTC, comment #1: 

So, to clarify because it's not entirely clear from your bug description, the bug you are reporting is that 'step(sys)' tries to plot too many points and overloads the Octave OpenGL plotting system, while 'step(ss(sys))' works as a workaround and does plot the correct figure. Is that accurate?

Mike Miller <mtmiller>
Group Member
Wed 10 Jun 2020 05:16:56 PM UTC, original submission:  

%Octave 5.2.1 and Control 3.2.0


pkg load control
q=tf([  -1.46668928632582  38.5234944088261  -1023.42433722006  18055.8729050165  -275471.128598611  3680685.32154243  -40413629.6869944  433109318.542883  -3699840188.09672  32827566143.7738  -225798332168.488  1689921426466.71  -9494747382069.7  60613579740811.5  -278626623858487  1.52747193482069e+15  -5.6883783399098e+15  2.69067309068308e+16  -7.92226372994905e+16  3.25360727881732e+17  -7.2059792475876e+17  2.61175754638012e+18  -3.8909506684981e+18  1.31637705830625e+19  -9.12640045084024e+18  3.78290213329238e+19  1.41630594163597e+19  5.02932079881115e+19  1.29547036349465e+20  1.05983336007477e+19  2.43633391646665e+20  -1.31208497581345e+19  1.04713601706961e+20  7.83465823853832e+18  ],[  1  17.1698120499521  542.934121038049  7294.98595152311  126155.971925398  1377306.37948524  16829405.6386297  152789570.438515  1448656865.6342  11095921096.5808  85365123037.7488  556372911296.259  3555399791580.02  19799716986834.5  106353561602035  506292672690482  2.29624277698132e+15  9.31618125768562e+15  3.56349262092232e+16  1.22398996473335e+17  3.92373264437545e+17  1.12860046096021e+18  2.99632129507514e+18  7.09934073973815e+18  1.53176948230893e+19  2.91760427412351e+19  4.97010301791734e+19  7.33491370510952e+19  9.42777562297368e+19  1.01604047339533e+20  9.12510568398965e+19  6.44576522525933e+19  3.46367980244933e+19  1.23309424846899e+19  2.4145011614912e+18  ]);
 
step(q)
%lsim(q,ones(1,1001),[0:.01:10])  %another simulator, same result

***************************************

these workarounds work which indicate that the internal tf2ss step needs correction:

1:
pkg load signal
[a,b,c,d] = tf2ss(q);
step(ss(a,b,c,d))

2:
only with the control pkg
step(ss(q))

************************************

picture show result from Matlab 2020a, matlab can handle the tf's. Was only possible to upload one picture.

see this thread for mote info:
https://octave.1599824.n4.nabble.com/control-package-can-t-simulate-complex-systems-td4697481.html




Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #49238:  OctaveControlGnuplot.png added by nrjank (110KiB - image/png)
file #49237:  bug58538_test.m added by ttl (2KiB - text/x-objcsrc)
file #49236:  MatlabControlOutput.png added by None (103KiB - image/png)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by nrjank (Updated the item)
  • -email is unavailable- added by ttl (Updated 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 14 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-07-03 ttl StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2020-06-19 ttl StatusConfirmed In Progress
    2020-06-10 nrjank Attached File- Added OctaveControlGnuplot.png, #49238
    2020-06-10 ttl Severity2 - Minor 3 - Normal
        Item GroupMatlab Compatibility Incorrect Result
    2020-06-10 ttl Attached File- Added bug58538_test.m, #49237
    2020-06-10 mtmiller CategoryNone Octave Package
        Severity3 - Normal 2 - Minor
        Priority5 - Normal 3 - Low
        StatusNone Confirmed
        Release5.2.0 dev
        Summarycontrol package, step() and lsim() can't simulate complex transfer functions [octave forge] (control) step and lsim fail on high order transfer functions
    2020-06-10 None Attached File- Added MatlabControlOutput.png, #49236

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code