implement native @zpk class

[MitchellThompkins]

Description

This adds a first-class citizen zpk class which is intended to resolve #30.

The changes here are largely mechanical routing to move data around in z/p/k form that were already provided that way (as would be the case for a call to c2d with matched, you already have the system poles/zero/gain). The current zpk() implementation in 037fc5f is a shim that immediately calls poly() and returns a tf.

With these changes and a slightly modified test script from #30, I get a much improved accuracy from c2d.

Modified test script

if exist('OCTAVE_VERSION', 'builtin')
    pkg load control
    addpath (fullfile (fileparts (mfilename ('fullpath')), 'inst'));
    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\n', err);
end

~/workspace/pkg-control > octave --no-gui c2d_zpk_error.m
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: 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

Changes

• Adds a new @zpk class which mirrors the other system representation classes (I largely stole from @tf).
  • Most of the conversion shims that I added delegate directly to @tf. In some cases I think that's fine, in others it is probably less appropriate but I would consider that future work (basically audit them and make sure that for each shim it's doing the most numerically correct thing).
  • Adds several new tests which are mostly trivial validation (with the exception of the 25-pole test I concocted to test this new zpk c2d against)
• Applies a small patch to @zpk/__d2c__ that controls where wd was computed and compared in a while loop. As it was written it always compared against a stale wd.
• Also applied similar changes from https://github.com/gnu-octave/pkg-control/pull/29/changes to the __d2c__ path here because that patch only fixed the __c2d__ path.
• Removes TODOs from @tf and uses matched-z discretization on tf when passed a tf system. The matched-z is properly applied to the zpk system as well.
• Ignores temporary vim files

[MitchellThompkins]

@ttl-octave I'm not sure how appropriate or necessary it is to add my own name here. I was just trying to follow suit.

[ttl-octave]

This depends on the amount of changes or extension that you have contributed to the file. Im this case (and maybe in some other of this commit) I think there is too little originality. Moreover, your contribution is recorded in the git history anyways.

[MitchellThompkins]

@ttl-octave The formatting here is different than what matlab produces, largely because I think getting their formatting correct is hard and I've included the relevant information someone would want when printing this thing. For example for my 25-pole test system it produces this:

Zero/pole/gain model 'sys':

y1 <- u1:
gain:  1
zeros: (none)
poles:
  -0.0251327+25.1327i
  -0.0251327+50.2655i
  -0.0251327+75.3982i
  -0.0251327+100.531i
  -0.0251327+125.664i
  -0.0251327+150.796i
  -0.0251327+175.929i
  -0.0251327+201.062i
  -0.0251327+226.195i
  -0.0251327+251.327i
  -0.0251327+276.46i
  -0.0251327+301.593i
  -0.0251327+326.726i
  -0.0251327+351.858i
  -0.0251327+376.991i
  -0.0251327+402.124i
  -0.0251327+427.257i
  -0.0251327+452.389i
  -0.0251327+477.522i
  -0.0251327+502.655i
  -0.0251327+527.788i
  -0.0251327+552.92i
  -0.0251327+578.053i
  -0.0251327+603.186i
  -0.0251327+628.319i

Where matlab produces:

sys =
 
                                                                                           
                                                                                          1
                                                                                           
  ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
                                                                                                                                                                                    
  (s+(0.02513-25.13i)) (s+(0.02513-50.27i)) (s+(0.02513-75.4i)) (s+(0.02513-100.5i)) (s+(0.02513-125.7i)) (s+(0.02513-150.8i)) (s+(0.02513-175.9i)) (s+(0.02513-201.1i)) (s+        
                                                                                                                                                                                    
          (0.02513-226.2i)) (s+(0.02513-251.3i)) (s+(0.02513-276.5i)) (s+(0.02513-301.6i)) (s+(0.02513-326.7i)) (s+(0.02513-351.9i)) (s+(0.02513-377i)) (s+(0.02513-402.1i))        
                                                                                                                                                                                    
          (s+(0.02513-427.3i)) (s+(0.02513-452.4i)) (s+(0.02513-477.5i)) (s+(0.02513-502.7i)) (s+(0.02513-527.8i)) (s+(0.02513-552.9i)) (s+(0.02513-578.1i)) (s+(0.02513-603.2i))   
                                                                                                                                                                                    
                                                                                                                                                                (s+(0.02513-628.3i))

[MitchellThompkins]

@ttl-octave I'm pretty sure that

pkg-control/inst/@tf/__d2c__.m

Line 52 in 037fc5f

w_d = exp (w_c * tsam);

meant to walk wd in the while loop, so I changed it.

[MitchellThompkins]

@ttl-octave I applied the same fix from #29; let me know if you think that's incorrect/inappropriate.

[MitchellThompkins]

@ttl-octave I applied the same fix from #29; let me know if you think that's incorrect/inappropriate.

[ttl-octave]

Impressive. Thank you very much for the contribution.

I am going to make a review in the next couple of days. None of the code is based on any other software (see Contribution.md), right?

Since this is a quite big extension, could you please change the pull request such that it compares to the dev branch? I have just merged main into dev and I therefore don't expect any problems.

[MitchellThompkins]

Impressive. Thank you very much for the contribution.

I am going to make a review in the next couple of days. None of the code is based on any other software (see Contribution.md), right?

Since this is a quite big extension, could you please change the pull request such that it compares to the dev branch? I have just merged main into dev and I therefore don't expect any problems.

Thanks! I did read the guidelines and am familiar with GPL licensing in general. I didn't use any code reference other than the implementations already provided in this library (like I said it's largely just mechanical
routing of the data).

[MitchellThompkins]

@ttl-octave one other note here, I took the notice from https://github.com/MitchellThompkins/pkg-control/blob/037fc5f50fb9666518667de018fe84865bbbac29/inst/%40tf/__c2d__.m#L1-L23 but that doesn't actually match https://github.com/gnu-octave/pkg-control/blob/eaba61e11b22de7e9b48871f939423e618d12f0e/CONTRIBUTING.md#license-and-documentation. I'm not sure which style is right. (This applies to both the GPL notice and the textinfo section info.

[ttl-octave]

Regarding the texinfo section, the help text for internal function __xyz__ does not need to follow the general formatting starting with the @deftypefn statement.

The GPL notice of existing files still contains the hint to the name "LTI Syncope", which was introduced by Lukas Reichlin as far as I know. As this name was and is not used officially for the control package it is not referred to in the example comment in the guidelines for contributing (where I just saw the copy & paste error, "statistics package").