Add a New Comment (Rich Markup)

( Jump to the original submission )

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)

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

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

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.

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.

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)

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

>> [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.

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.

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

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

Attach Files: Comment:

Depends on the following items: None found

Items that depend on this one: None found

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.

Please enter the title of George Orwell's famous dystopian book (it's a date):

Follow 3 latest changes.

Copyright © 2023 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. Source Code

Powered by Savane 3.11