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

 
 

bug #53881: [octave forge] (control) impulse response differs greatly on i386 vs x86_64

Submitter:  Mike Miller <mtmiller>
Submitted:  Fri 11 May 2018 08:49:22 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  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

Sun 21 Jan 2024 02:29:58 PM UTC, comment #9: 

I have changed the algorithm of the control package's imp_invar for transfer functions and have also added the impulse invariant discretization for state space systems (dev branch). These changes will be included in the next minor release 4.1.0.

I am closing the report for now. It can be reopened if any new numerical issues show up in the new version.

Torsten Lilge <ttl>
Group Member
Sat 03 Jul 2021 02:14:00 PM UTC, comment #8: 

I have realized some numerical problems in imp_invar when it comes o systems with orders higher than 3 or 4, especially when the system has repeated poles. I have found a quite simple way to compute the impulse invariant discretization by firstly using c2d with the 'zoh' method for s*G(s) and multiplying afterwards with z/z-1 (see comment in help test).

I have attached the new inp_invar.m for testing.

A remark: Matlab and obviously also the impinvar routine in the signal package assume an impulse of height 1/T for the impulse response. imp_invar is assuming the unit impulse and therefore returns a discrete-time description multiplied by 1/T compared to impinvar. This was already the case before my change and, to my impression, is correct from a control perspective. Discrete-time systems are not necessarily connected to a sampling time and a test input 1/T makes no sense for system being discrete-time "by nature".


(file #51643)

Torsten Lilge <ttl>
Group Member
Sat 12 May 2018 06:24:54 PM UTC, comment #7: 

OK so matlab does not do the algebraic simplification before doing the impulse invariant math.

Doug Stewart <dastew>
Sat 12 May 2018 06:21:08 PM UTC, comment #6: 

Yes. and I understand why.
I want to see what matlab does.

Doug Stewart <dastew>
Sat 12 May 2018 06:20:33 PM UTC, comment #5: 

The Matlab impinvar function returns the following


bz = 1×2
10⁻³ ×

     1.000000000000000  -0.999000499833375

az = 1×3

     1.000000000000000  -1.998000999666750   0.998001998667333



Identical to the Octave signal impinvar function results.

Mike Miller <mtmiller>
Group Member
Sat 12 May 2018 06:04:14 PM UTC, comment #4: 

For comparison, here is Octave signal package impinvar result, which I hadn't listed before:


bz =

   1.000000000000000e-03  -9.990004998333751e-04

az =

   1.000000000000000e+00  -1.998000999666750e+00   9.980019986673331e-01


which is different from the control package imp_invar result.

Mike Miller <mtmiller>
Group Member
Sat 12 May 2018 12:08:32 PM UTC, comment #3: 

I Think we should find out what Matlab does with poles and zeroes that can be canceled.

what does matlab do with:
 [bz, az] = impinvar ([1 1], [1 2 1], 1000)


Doug Stewart <dastew>
Sat 12 May 2018 02:27:12 AM UTC, comment #2: 

As I discussed with Doug on IRC, the system tf([1, 1], [1, 2, 1]) is a particularly problematic system because it contains a pole and zero that cancel each other out at -1.

Should this be fixed or should it be accepted as reasonable numerical error in this particular system?

If the system is simplified to the expression tf(1, [1, 1]), then there is no error, both i386 and x86_64 give the same result.

The source of the error between i386 and x86_64 is the impulse invariant c2d transform, as shown here.

On i386


>> [A, B, C, D, dt] = ssdata (c2d (tf ([1 1], [1 2 1]), 1/1000, 'impulse'))
A =    9.989991670891507e-01
B =    9.989991670891507e-01
C =  1
D =  1
dt =    1.000000000000000e-03


On x86_64


>> [A, B, C, D, dt] = ssdata (c2d (tf ([1 1], [1 2 1]), 1/1000, 'impulse'))
A =    9.990076217466157e-01
B =    9.990076217466157e-01
C =  1
D =  1
dt =    1.000000000000000e-03


The error between these computed values A and B on the two systems (8.4e-6 in this case) is solely responsible for the compounding relative error between the impulse responses.

And this comes down to the imp_invar.m function in the control package returning different results on i386 vs x86_64.

On i386


>> [bz, az] = imp_invar ([1 1], [1 2 1], 1000)
bz =

   1   0

az =

   1.000000000000000e+00  -9.989991670891507e-01


On x86_64


>> [bz, az] = imp_invar ([1 1], [1 2 1], 1000)
bz =

   1   0

az =

   1.000000000000000e+00  -9.990076217466157e-01


The irony is that I am only looking into this because of a test failure in the signal package function impinvar, where the test compares the discrete filter obtained from impinvar against the impulse response calculated by the control package. The signal package impinvar results do not vary between i386 and x86_64, but the control package results do, and the error was more than the test tolerance was written for.

So if we leave this as reasonable numerical error, then users of impulse on continuous systems have to be aware that small numerical errors are possible, which compound when the impulse response is calculated over a large time series.

Mike Miller <mtmiller>
Group Member
Fri 11 May 2018 08:55:04 PM UTC, comment #1: 

Sorry, I was playing with different time steps. The scripts that I pasted should be


pkg load control
sys = tf ([1, 1], [1, 2, 1]);
y = impulse (sys, [0:999] ./ 1000).';
save results32.txt


and


pkg load control
sys = tf ([1, 1], [1, 2, 1]);
y = impulse (sys, [0:999] ./ 1000).';
results32 = load ('results32.txt');
fprintf (stdout, '%4u  %-17.15g  %-17.15g\n', [0:999; results32.y; y]);


to get the data arrays and divergence that I showed in the original post.

Mike Miller <mtmiller>
Group Member
Fri 11 May 2018 08:49:22 PM UTC, original submission:  

The impulse response of a simple system is very inaccurate when comparing the results on X86-32 (i386) to the results on X86-64.

I run the following script on a 32-bit Octave 4.4 with control 3.1.0


pkg load control
sys = tf ([1, 1], [1, 2, 1]);
y = impulse (sys, [0:99] ./ 100).';
save results32.txt


and the following script on a 64-bit Octave 4.4 with control 3.1.0


pkg load control
sys = tf ([1, 1], [1, 2, 1]);
y = impulse (sys, [0:99] ./ 100).';
results32 = load ('results32.txt');
fprintf (stdout, '%4u  %-17.15g  %-17.15g\n', [0:99; results32.y; y]);


and I get the following output. Even at the start of the impulse response the vectors are accurate to only three decimal places. By the end the accuracy is down to one decimal place.


   0  1                  1
   1  0.998999167089151  0.999007621746616
   2  0.997999335844817  0.998016228307829
   3  0.997000505264498  0.997025818706332
   4  0.996002674346696  0.996036391965785

 250  0.778541081031692  0.780190040994487
 251  0.777761891495348  0.779415797364297
 252  0.776983481797535  0.778642322076649
 253  0.776205851157766  0.777869614369055
 254  0.775428998796333  0.777097673479787

 500  0.606126214853996  0.60869650006698
 501  0.605519583790042  0.608092442897402
 502  0.604913559862421  0.607488985181024
 503  0.604308142463491  0.60688612632296
 504  0.60370333098622   0.606283865728916

 995  0.369232996221608  0.372355352593404
 996  0.368863455687218  0.371985835238959
 997  0.368494285001156  0.371616684585501
 998  0.368125483793267  0.371247900269124
 999  0.367757051693765  0.370879481926282


Mike Miller <mtmiller>
Group Member

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #51643:  imp_invar.m added by ttl (5KiB - text/x-objcsrc)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by ttl (Updated the item)
  • -email is unavailable- added by dastew (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
    2024-01-21 ttl StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2021-07-03 ttl Attached File- Added imp_invar.m, #51643
    2019-02-26 mtmiller StatusNone Confirmed
    2019-02-26 mtmiller Carbon-CopyRemoved 80942 -

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code