[MATLAB Grader] Comparing Transfer Functions

Just wondering if there is a good way to compare transfer functions in MATLAB Grader?

I am testing a pseudo student response against the reference solution, purposely adding some slight deviation to check the test function.

Both the pseudo and reference solution output the same transfer function (at least to 4 dp), however from my testing, I believe the the condition:

abs(expected-actual)

cant be met, as it is the incorrect argument data type. Further to this,

(expected-actual)

creates an error function with higher order terms, I believe

minreal(expected-actual)

would make more sense.

Any help appreciated.

[Paul]

Further discussion here that raises some issues that might be of interest.

Perhaps @Bill Tubbs has found a battle tested solution he's willing to share in addition to that shown in the linked Question.

[Bill Tubbs]

Im afraid i didnt. I asked my supervisors about this and they suggested to compare the step responses. The Solution by Matt Rich below looks better.

[Matt Rich]

@Paul The norm of the difference of the two transfer functions will return exactly zero whenever the accepted answer in the linked question would work. A nice advantage is that bounding it by a small tolerance will also account for the computational numerical precision issues, as well as the ss/zpk/tf or other class conversion precision issues mentioned there. Even better, it will work directly on ss or zpk or even a mixture of LTI system classes, it will work for systems with zeros, it will work for MIMO not just SISO, and will handle unstable systems. The only thing I can think of right now that would cause any problems is systems with internal delays.

@Bill Tubbs For stable SISO systems comparing the impulse or step response is going to be good too. I think it is probably the best way to go if you have delays.

[Paul]

@Matt Rich

I do like the L-inf approach, but it all really depends on what @BW considers equivalent.

For example:

sys1 = zpk([],-1,1)

sys1 =
 
    1
  -----
  (s+1)
 
Continuous-time zero/pole/gain model.

sys2 = zpk(1,[-1 1],1)

sys2 =
 
     (s-1)
  -----------
  (s+1) (s-1)
 
Continuous-time zero/pole/gain model.

Of course the difference from an L-inf norm perspective is going to be very small

norm(sys1-sys2,Inf)

ans = 3.2149e-16

But from a grading perspective, whether or not sys1 is equivalent to sys2 depends on whether or not the problem requires the student to perform the pole zero cancellation.

Perhaps a more concrete example would be

sys1 = tf(1,[1 1]);
sys2 = tf([1 2],[1 3]);
sysfb1 = feedback(sys1,sys2)

sysfb1 =
 
      s + 3
  -------------
  s^2 + 5 s + 5
 
Continuous-time transfer function.

sysfb2 = sys1/(1 + sys1*sys2)

sysfb2 =
 
      s^2 + 4 s + 3
  ----------------------
  s^3 + 6 s^2 + 10 s + 5
 
Continuous-time transfer function.

Here, sysfb1 and sysfb2 are identical by the L-inf test

norm(sysfb1 - sysfb2,Inf)

ans = 0

But, if the point of exercise is learn how to use feedback() instead of doing transfer function algebra in Matlab, then perhaps sysfb2 is to be considered incorrect.

And then of course we have to make sure the tolerance used in the L-inf test is smaller than the gain of the systems in question. Are sys1 and sys2 suppsed to be equal.

sys1 = 1e-10*tf(1,[1 1]);
sys2 = 1.1e-10*tf(1,[1 1]);
norm(sys1-sys2,inf)

ans = 1.0000e-11

Finally, as to your excellent point about delays, I think it really extends even to systems with input/ouput delays, because even those can become internal delays after the subtraction

sys1 = tf(1,[1 1],'InputDelay',1);
sys2 = sys1;
norm(parallel(sys1,-sys2),Inf)  % works ok

ans = 0

sys2 = tf(1,[1 1],'InputDelay',2);  % doesn't work ok
norm(parallel(sys1,-sys2),Inf)

Error using DynamicSystem/norm
The "norm" command cannot be used for continuous-time models with internal delays. Use the "pade" command to approximate such delays.

[BW]

@Matt Rich, out of interest did you test or change anything for unstable zpk systems using the approach?

For example - the TF (say )

-7.0554 (s+2.919) (s^2 + 0.04626s + 0.3162)
-----------------------------------------------
s (s+2.924) (s+0.01636) (s^2 + 0.5739s + 6.944)

will have an infinite norm. Given a second TF (say ), where the difference has a non-zero gain:

   -0.035277 (s+2.919) (s^2 + 0.04626s + 0.3162)
  -----------------------------------------------
  s (s+2.924) (s+0.01636) (s^2 + 0.5739s + 6.944)

The norm of this difference will still be infinite, so the only way to assess it as being correct is:

  assert(norm( dr2psi_SYSzpk_SA - dr2psi.SYSzpk , 'inf' ) <= inf , "Incorrect Transfer Function")

Obviously this is not a representative test condition as any submitted answer would pass.

[Paul]

Shouldn't G1 and G2 be assessed as incorrect ? They don't look the same to me because -7.0554 ~= -0.035277.

But this brings up an excellent point about how to compare systems with poles on the imaginary axis using the L-inf norm.

[BW]

Hi Paul,

The question was maybe ill-posed,

In the above, G2 was set to 1.005*G1 to test a pseudo student response with 0.5% error against the reference solution, therefore the gain -0.035277 came from the difference between them (G2-G1 = -7.090677--7.0554).

[Paul]

Ah, I didn't realize that the second tf was the difference G1 - G2.

[Matt Rich]

Yes, @BW @Paul excellent point about systems with poles on the imaginary axis! In those cases, one potential modification to the approach I can think of right now is to evaluate for some nonzero, conveniently chosen U which cancels the poles of the error system that are on the imaginary axis.

For example, form two "nearly identical" systems, each with 3 poles on the imaginary axis:

s = tf('s'); 
num = -7.0554*(s+2.919)*(s^2 + 0.04626*s + 0.3162);
den = s*(s+2.924)*(s+0.01636)*(s^2 + 0.5739*s + 6.944)*(s^2+1);
G1 = num/den

G1 =
 
                   -7.055 s^3 - 20.92 s^2 - 3.184 s - 6.512
  --------------------------------------------------------------------------
  s^7 + 3.514 s^6 + 9.679 s^5 + 23.96 s^4 + 9.011 s^3 + 20.45 s^2 + 0.3322 s
 
Continuous-time transfer function.

% duplicate the transfer function above introducing small errors 
err = 100*eps; 
num2 = (-7.0554 + err )*(s+2.919) *(s^2 + 0.04626*s + 0.3162);
den2 =   s *(s+2.924) *(s+0.01636 + err) *(s^2 + 0.5739*s + 6.944)*(s^2+1+err);
G2 = num2/den2

G2 =
 
                   -7.055 s^3 - 20.92 s^2 - 3.184 s - 6.512
  --------------------------------------------------------------------------
  s^7 + 3.514 s^6 + 9.679 s^5 + 23.96 s^4 + 9.011 s^3 + 20.45 s^2 + 0.3322 s
 
Continuous-time transfer function.

Check (expected to be ∞ in this case since both systems have a pole on the imaginary axis)

norm(G1-G2,inf)

ans = Inf

Now, calculate the poles of each system, and extract those on the imaginary axis (allowing for some numerical error):

p1 = pole(G1);
p1imag = p1( abs(real(p1)) < 10*eps )

p1imag = 
   0.0000 + 0.0000i
   0.0000 + 1.0000i
   0.0000 - 1.0000i

p2 = pole(G2);
p2imag = p2( abs(real(p2)) < 10*eps )

p2imag = 
   0.0000 + 0.0000i
   0.0000 + 1.0000i
   0.0000 - 1.0000i

Take the union of the two sets of poles:

pImag = union(p1imag,p2imag)

pImag = 
   0.0000 + 0.0000i
   0.0000 - 1.0000i
   0.0000 + 1.0000i
   0.0000 - 1.0000i
   0.0000 + 1.0000i

Use this set of poles to create a "convenient" (in this case such that ) nonzero U:

U = zpk(pImag,-1*ones(size(pImag)) ,1)

U =
 
  s (s^2 + 1)^2
  -------------
     (s+1)^5
 
Continuous-time zero/pole/gain model.

U = tf(U)

U =
 
  s^5 - 1.11e-15 s^4 + 2 s^3 - 1.11e-15 s^2 + s
  ---------------------------------------------
     s^5 + 5 s^4 + 10 s^3 + 10 s^2 + 5 s + 1
 
Continuous-time transfer function.

Finally, check

norm((G1-G2)*U,inf)

ans = 2.7102e-11

I have not tested this too thoroughly yet with different amounts of error, MIMO systems, and need to think if there's a nice way to make it work for discrete without adding much code for logical statements but it seems promising.

[BW]

Thanks you once again,

Another option which worked in this case was the ratio of transfer functions, which should be close to unity (using minreal)

Answer 1

[Matt Rich]

The most robust approach to assess LTI systems in MATLAB Grader is probably to use the norm of the difference between the reference and learner solution: for some small tolerance.

% check existence

assert( exist('G', 'var') , "The submission must contain a variable named G.")

% check class is a type of LTI system object

assert( isa(G,'tf') | isa(G,'ss') | isa(G,'zpk') , "The variable G must be an LTI system object (tf, ss, or zpk).")

% check IO size

assert( isequal( size(G), size(referenceVariables.G) ) , "G has has the wrong input/output dimensions.")

% check equivalence -- is L_inf norm of difference (error) close enough to zero?

tol = 1e-10;

assert( norm( referenceVariables.G - G , 'inf' ) < tol , "G is not correct.")

If you want to enforce the system is represented using a specific class (e.g., transfer function) then edit the class check:

% check class is tf

assert( isa(G,'tf') , "The variable G must be an transfer function.")

[BW]

Excellent! Thanks all

[Agustín]

Could we compare both numerator and denominador of the transfer function instead of evaluating the L_inf norm? This might help students to identify where the error is.