c2d with matched-z introduces numerical instability from ZPK-to-TF polynomial round-trip

[MitchellThompkins]

Version of the Control Package, version of Octave

~/workspace/pkg-control > octave --no-gui --quiet --eval "disp(version()); pkg load control; pkg list control"
10.3.0
Package Name  | Version | Installation directory
--------------+---------+-----------------------
     control *|   3.6.1 | .../share/octave/octave_packages/control-3.6.1

What happens?

I have been working on my own constexpr C++ filter library constfilt which uses octave to generate references against which the library is tested. My library supports ZOH, matched-z, and Tustin discretization.

I started to implement matched-z biquad filters as I drove to higher orders and started to generate references, and I started to see the behavior in octave degrade.

I'm pretty sure the issue is that @zpk today is just a wrapper function that converts poles/zeros to @tf. So by the time c2d sees the system, it's already a TF object with polynomial coefficients. The matched-z impl in @tf/__c2d__.m has no choice but to call zpkdata() which root-finds that polynomial to recover teh poles.

[z_c, p_c, k_c] = zpkdata (sys, "vector");

There's actually a TODO surrounding this today in the implementation. It looks liek this has always been planned but hasn't happened yet. By not having a zpk container though whose shape matches data born as poles, zeros, and gain, the impl today stuffs that data into a stand-in container which ends up storing the data in the wrong shape (tf), and the reshaping is destructive at high orders b/c it ends up performing root-finding on something ill-conditioned.

What is the expected result?

c2d with the matched-z method should map poles directly via exp(p*Ts) without converting to a polynomial representation. In MATLAB, zpk() returns a native ZPK object and c2d uses it directly, so this problem does not occur there.

Steps to reproduce

I threw together this script to demonstrate this:

if exist('OCTAVE_VERSION', 'builtin')
    pkg load control
    fprintf ('Running with octave\n');
else
  fprintf ('Running with matlab\n');
end

N  = 25; % fitler order
Ts = 1/1000;

% N poles in the left half-plane which are all close to the imaginary axis (i.e.
% stable but their real parts are small relative to their imaginary parts,
% making poly(ps) ill-conditioned for root finding and placing the digital
% poles close to the unit circle).
ps = (-0.001 + 1i * linspace(1, 25, N)).' * 2*pi*4;

%% analytic matched-z transform per definition
pd_matchedz = exp(ps * Ts);

fprintf('\nAnalog poles start stable: all Real(p) < 0 = %d\n', all(real(ps) < 0));
fprintf('Digital poles are stable: all |pd| < 1 = %d\n', all(abs(pd_matchedz) < 1));
fprintf('Digital poles are all very close to unit circle: max|pd| = %.10f\n\n', max(abs(pd_matchedz)));

%% what does zpk + c2d produce
sys = zpk([], ps, 1);
sys_d = c2d(sys, Ts, 'matched');
[~, pd, ~] = zpkdata(sys_d, 'v');

fprintf('zpk class: %s\n', class(sys));
fprintf('c2d digital poles are stable: all |pd| < 1 = %d\n', all(abs(pd) < 1));
fprintf('c2d digital poles max magnitude: max|pd| = %.10f\n\n', max(abs(pd)));

%% results
if any(abs(pd) > 1)
    fprintf('incorrect: c2d produced unstable poles from a stable system.\n');
else
    fprintf('correct: c2d produced all stable poles from a stable system.\n');
    [~, ia] = sort(abs(pd));
    [~, ib] = sort(abs(pd_matchedz));
    err = max(abs(pd(ia) - pd_matchedz(ib)));
    fprintf('max error between analytically calculated matched-z discrete poles and zpkdata calculated %.16e', err);
end

Matlab

>> c2d_zpk_error
Running with matlab

Analog poles start stable: all Real(p) < 0 = 1
Digital poles are stable: all |pd| < 1 = 1
Digital poles are all very close to unit circle: max|pd| = 0.9999748676

Warning: This zpk model has a complex gain or some complex zeros or poles that do not come in conjugate pairs.
> In zpk (line 388)
In zpk_test (line 39)
zpk class: zpk
c2d digital poles are stable: all |pd| < 1 = 1
c2d digital poles max magnitude: max|pd| = 0.9999748676

correct: c2d produced all stable poles from a stable system.
max error between analytically calculated matched-z discrete poles and zpkdata calculated 0.0000000000000000e+00>>

Octave

~/workspace/pkg-control > octave --no-gui c2d_zpk_error.m
octave: X11 DISPLAY environment variable not set
octave: disabling GUI features
Running with octave

Analog poles start stable: all Real(p) < 0 = 1
Digital poles are stable: all |pd| < 1 = 1
Digital poles are all very close to unit circle: max|pd| = 0.9999748676

zpk class: tf
c2d digital poles are stable: all |pd| < 1 = 0
c2d digital poles max magnitude: max|pd| = 178.9002929728

incorrect: c2d produced unstable poles from a stable system.

[ttl-octave]

Your test produces the correct result after merging pull request #31. Can this issue be closed?

[MitchellThompkins]

Your test produces the correct result after merging pull request #31. Can this issue be closed?

I'm good with closing if you are. My changes address the specific problem I identified here (although there's some lite cleanup I haven't had the chance to sit down and do yet, but that's unrelated to this issue). The changes landed in dev though so I'm not sure if that matters (or when dev it gets merged to main).

[ttl-octave]

Very good. The plan is to publish a 4.2.3 release based on the current changes in the main branch soon and then go for a 4.3.0 release with the dev branch merged to main.