Behavior of isPassive and hinfnorm.

For the transfer function HzPw attached, I am not getting the behavior I expect from the function isPassive and hinfnorm. Any suggestions on what I am interpreting incorrectly would be much appreciated.

 load HzPw.mat;
 
 nopt = nyquistoptions;
 nopt.XLim = [-5e-3 2e-2];
 nyquistplot(HzPw, nopt);

The Nyquist plot of this transfer function is entirely in the right half plane. So, I expect this transfer function to be passive. However,

 [pf,R] = isPassive(HzPw)

pf = logical
   0

R = Inf

Then, I look at a scattering function corresponding to this transfer function, and the Bode plot of this scattering function.

 Hscat = (1-HzPw)/(1+HzPw);
 bode(Hscat);

I expect the norm of this to be 1. Indeed, by my thinking, the transfer function is passive, and so the scattering function norm should be ≤ 1. However,

 hinfnorm(Hscat)

ans = Inf

I am wondering what I am not thinking about correctly.

Thank you.

Answer 1

[Jon]

I'm thinking that you may have some numerical issues with poles very close to the jw axis. Your transfer function uses quite high order polynomials. Your numerator is a 9th order polynomial and your denominator is a 10th order polynomial. For such high order polynomials very small changes in the coefficient values may result in large changes in the pole and zero location. They are therefore generally to be avoided in favor of other representations.

In particular I see using MATLAB's pole command

p = pole(HzPw)

p =

   1.0e+02 *

  -9.0685 + 0.0000i
  -0.9015 + 3.3574i
  -0.9015 - 3.3574i
  -0.3718 + 0.2699i
  -0.3718 - 0.2699i
  -0.0566 + 0.1218i
  -0.0566 - 0.1218i
  -0.0031 + 0.0036i
  -0.0031 - 0.0036i
   0.0036 + 0.0000i

So you in fact have an unstable pole in the right half plane at 0.0036.

I'm not an expert in the relation between stability and passivity, but my intuition is that an unstable system is not passive.

[Jon]

Following up on @Paul answer, note that you have a close but not exact pole zero cancellation of that unstable pole

>> zero(HzPw)

ans =

   1.0e+02 *

  -9.3875 + 0.0000i
  -0.8611 + 3.3513i
  -0.8611 - 3.3513i
  -0.3093 + 0.2775i
  -0.3093 - 0.2775i
  -0.0031 + 0.0036i
  -0.0031 - 0.0036i
   0.0035 + 0.0000i
  -0.0000 + 0.0000i

>> pole(HzPw)

ans =

   1.0e+02 *

  -9.0685 + 0.0000i
  -0.9015 + 3.3574i
  -0.9015 - 3.3574i
  -0.3718 + 0.2699i
  -0.3718 - 0.2699i
  -0.0566 + 0.1218i
  -0.0566 - 0.1218i
  -0.0031 + 0.0036i
  -0.0031 - 0.0036i
   0.0036 + 0.0000i

[Paul]

Actually, that unstable pole is at s = 0.36, not 0.0036.

Looks like there are some other near cancellations as well, like at -0.31 + 1j*.36.

I wonder how HzPw was constructed, myabe by "transfer function math" that can introduce some additional poles that aren't really in the system .. more info needed from @Siva

[Siva]

Thank you @Jon and @Paul. I had missed HzPw having an unstable pole, and the near pole-zero cancellation.

Following up, I did

load HzPw.mat;
HzPw_reduced = minreal(HzPw,3e-2);
bode(HzPw, HzPw_reduced)

The frequency responses are nearly identical, but the reduced order transfer function is stable (the near pole-zero cancellation dynamics has been discarded?).

zero(HzPw_reduced)

ans = 
   1.0e+02 *

  -9.3875 + 0.0000i
  -0.8611 + 3.3513i
  -0.8611 - 3.3513i
  -0.3093 + 0.2775i
  -0.3093 - 0.2775i
  -0.0000 + 0.0000i

pole(HzPw_reduced)

ans = 
   1.0e+02 *

  -9.0685 + 0.0000i
  -0.9015 + 3.3574i
  -0.9015 - 3.3574i
  -0.3718 + 0.2699i
  -0.3718 - 0.2699i
  -0.0566 + 0.1218i
  -0.0566 - 0.1218i

The reduced transfer function is indeed passive.

[pf,R] = isPassive(HzPw_reduced)

pf = logical
   1

R = 1.0000

And the corresponding scattering function has norm = 1.

Hscat_reduced = (1-HzPw_reduced)/(1+HzPw_reduced);
hinfnorm(Hscat_reduced)

ans = 1

HzPw is a product of using systune. I am doing an optimal control calculation with a tuanbleTF in the model. I was trying to apply a TuningGoal.Passivity constraint on HzPw, but was unable to, even though it ended up looking passive. That is what led me to this line of exploration. I can post a minimal version of the systune calculation, or perhaps that would better belong in a different thread?

I am guessing I have to scale my problem better before passing it off to systune to avoid issues like near cancellations of unstable poles and zeros?

In fact, I can reduce the order even further without changing the dynamics much.

HzPw_reduced2 = minreal(HzPw,5e-2);
bode(HzPw, HzPw_reduced2)

Maybe this does point to poor conditioning?

[Paul]

@Siva

Does "product of using systune" mean HzPw is the output, CL, of systune() evaluated at the tuned parameters?

How was CL0, the input to systune() constructed? Using interconnection functions, like feedback(G,1), or using transfer function algebra, like G/(1 + G)?

[Siva]

@Paul Yes, HzPw is part of the output, CL, of systune evaluated at the tuned parameters. I get HzPw as

HzPw = tf(getIOtransfer(CL,'w','zP'));

CL0 was constructed using transfer function blocks, sumblks and the connect() function.

I am not sure if this uses transfer function algebra.

One thing I realized in the meanwhile is that my tunableTF had degree 5 numerator and denominator polynomials. If I instead used degree 2 polynomials, I get the same same HzPw as if I used the (5,5) tuableTF and then constructed a minimum realization. The additional degrees of freedom I had for tuning probably led to the near pole-zero cancellations.

[Paul]

Using connect() is the right approach. Of course, there could still be pole/zero cancellations depending on the blocks in the model, but that's different than what can happen using transfer function algebra.

[Paul]

Another thing to keep in mind is that it looks like isPassive only checks the poles, and does not account for pole/zero cancellations

isPassive(zpk([],-1,1))

ans = logical
   1

isPassive(zpk(1,[-1 1],1))

ans = logical
   0

Of course, checking for pole/zero cancelations can be a tricky business, so it's understandable tht isPassive() takes no action on its own in this regard, but it's something to bear in mind as the user.

[Siva]

@Paul Thank you.

I am curious how transfer function algebra and functions like connect() and feedback() are implemented differently in MATLAB.

It makes sense that isPassive() returns false for the transfer function with the unstable pole-zero cancellation, since it likely uses a state space based algorithm and the state space representations of the two transfer functions are different. However, this raises the following question for me. In the documentation for isPassive, it says, "Equivalently, a system is passive if its frequency response is positive real, ...". Should that be, "... if it is internally stable and its frequency response is positive real ..."?

[Paul]

I don't know anymore how connect() is implemented and so would only be able to speculate, but I will speculate that at least it's using state space math, since it returns an ss object for numeric LTI inputs.

C = pid(2,1); 
G = zpk([],[-1,-1],1);
blksys = append(C,G);
connections = [2 1; 1 -2];
T = connect(blksys,connections,1,2)

T =
 
  A = 
           x1      x2      x3
   x1       0  -1.414       0
   x2       0      -1       1
   x3  0.7071      -2      -1
 
  B = 
          u1
   x1  1.414
   x2      0
   x3      2
 
  C = 
       x1  x2  x3
   y1   0   1   0
 
  D = 
       u1
   y1   0
 
Continuous-time state-space model.

zpk(T)

ans =
 
             2 (s+0.5)
  --------------------------------
  (s+0.4302) (s^2 + 1.57s + 2.325)
 
Continuous-time zero/pole/gain model.

As expected, T is third order.

But doing transfer function / zpk algebra we get

zpk(C*G/(1+C*G))

ans =
 
             2 s (s+0.5) (s+1)^2
  ------------------------------------------
  s (s+1)^2 (s+0.4302) (s^2 + 1.57s + 2.325)
 
Continuous-time zero/pole/gain model.

which would be the same as zpk(T) after cancelling the common poles/zerost that arise from essentially doing the polynomial math.

Those common poles/zeros can be eliminated using minreal(), but minreal() would also eliminate common poles/zeros that are actually inherent to the system, that may actually be of interest. For example

C = zpk(-1,3,1); 
G = zpk([],[-1,-1],1);
blksys = append(C,G);
connections = [2 1; 1 -2];
T = connect(blksys,connections,1,2);
zpk(T)

ans =
 
            (s+1)
  --------------------------
  (s-2.732) (s+1) (s+0.7321)
 
Continuous-time zero/pole/gain model.

Here, T maintains the zero from C that cancels the pole in G. But if we do

C*G/(1 + C*G)

ans =
 
            (s+1)^3 (s-3)
  ----------------------------------
  (s-3) (s-2.732) (s+1)^3 (s+0.7321)
 
Continuous-time zero/pole/gain model.

then it would be difficult to know, in general, that only two of the poles/zeros at s = -1 are phantom that should be eliminated.

Moreal of the story is that it's generally best to use interconnection functions like series, ,parallel, connect, etc. The doc pages talk about this as well.

Should that be, "... if it is internally stable and its frequency response is positive real ..."?

That's probably correct wrt to how isPassive is actually implemented, but I thought passivity is really an input/output characteristic of the transfer function as discussed on isPassive , so then going on to say anything about internal stability without furhter explanation would be odd. But I agree with your sentiment that the doc page should do a better job of explaining how isPassive() handles systems with RHP poles..

[Siva]

@Paul That is a great example. That clarifies all of my current questions. Thank you @Paul and @Jon.

[Siva]

Just ran into this, related to the discussion in this thread.

[Paul]

This link might also be of interest, though I disagree with the "guarantee" in the very last sentence on that page, which is clearly incorrect as shown above.

Answer 2

[Paul]

Hi Siva,

I'm not able to load the .mat file, so can't really look at the problem.

load('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1020880/HzPw.mat');

Error using load
Unable to read file 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1020880/HzPw.mat'. If it is a Version 7 or earlier MAT-file, consider saving your data afresh in Version 7.3
MAT-files to access it from a remote location.

Can you save in v7.3 format and upload, or, if simple enough, just show what it is otherwise, like displaying the A,B,C,D or zpkdata, or the tfdata?

If not, can you post the output of the following commands:

pole(HzPw)
zero(HzPw)
pole(Hscat)

I'm asking because this doc page states "As a result, passive systems are intrinsically stable."

So maybe HzPw is not BIBO stable, but has a Nyquist plot all in the RHP, like this

h = zpk(2,1,1)
nyquist(h)

Similary, the H-infinity norm, IIRC, is not defined for unstable systems, so if Hscat is not BIBO stable, then hinfnorm will return inf as stated on hinfnorm

[Jon]

The OP's transfer function has a pole in the right half plane, please see my answer above

[Paul]

Yes, I saw your Answer after I clicked Submit.

I guess it's safe to assume that HzPw doesn't have a zero that exactly cancels that pole, though isPassive doesn't seem to care about that.

I assume Hscat also has pole in the RHP?

[Jon]

Yes, in fact a double one

>>  Hscat = (1-HzPw)/(1+HzPw);
>> pole(Hscat)

ans =

   1.0e+02 *

  -9.0685 + 0.0000i
  -9.0684 + 0.0000i
  -0.9015 + 3.3574i
  -0.9015 - 3.3574i
  -0.9015 + 3.3573i
  -0.9015 - 3.3573i
  -0.3718 + 0.2699i
  -0.3718 - 0.2699i
  -0.3720 + 0.2695i
  -0.3720 - 0.2695i
  -0.0566 + 0.1218i
  -0.0566 - 0.1218i
  -0.0576 + 0.1213i
  -0.0576 - 0.1213i
  -0.0031 + 0.0036i
  -0.0031 - 0.0036i
  -0.0031 + 0.0036i
  -0.0031 - 0.0036i
   0.0036 + 0.0000i
   0.0036 + 0.0000i

[Jon]

And here is the transfer function given for HzPw

HzPw =
 
  From input "w" to output "zP":
                                                                                                    
  0.2415 s^9 + 283.4 s^8 + 8.506e04 s^7 + 3.184e07 s^6 + 1.805e09 s^5 + 4.735e10 s^4 + 1.238e10 s^3 
                                                                                                    
                                                          + 9.535e07 s^2 - 3.645e09 s - 1.714e04    
                                                                                                    
  --------------------------------------------------------------------------------------------------
                                                                                                  
    s^10 + 1173 s^9 + 3.809e05 s^8 + 1.375e08 s^7 + 1.036e10 s^6 + 3.57e11 s^5 + 4.286e12 s^4     
                                                                                                  
                                             + 4.277e13 s^3 + 1.068e13 s^2 - 3.252e11 s - 3.334e12
                                                                                                  
 
Continuous-time transfer function.

[Siva]

I think the mat file was v7.3 already, but just in case I saved explicitly with the -v7.3 option.

[Paul]

load('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1021010/HzPw.mat');

Error using load
Unable to read file 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1021010/HzPw.mat'. If it is a Version 7 or earlier MAT-file, consider saving your data afresh in Version 7.3
MAT-files to access it from a remote location.

Hmm, same problem. I'm probably not using load() correctly for a remote file.

[Siva]

@Paul From this thread, this seems to work:

load(websave('HzPw', 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1020880/HzPw.mat'));
pole(HzPw)

ans = 
   1.0e+02 *

  -9.0685 + 0.0000i
  -0.9015 + 3.3574i
  -0.9015 - 3.3574i
  -0.3718 + 0.2699i
  -0.3718 - 0.2699i
  -0.0566 + 0.1218i
  -0.0566 - 0.1218i
  -0.0031 + 0.0036i
  -0.0031 - 0.0036i
   0.0036 + 0.0000i