advanced frequency domain analysis of SISO models

Hello community,

suppose we have any given system 'G' for example:

G=tf(1, [2 0.4 1])

Suppose we search for the smallest frequency 'omega_d' at which the real part of the frequency response is zero. In other words:

real(evalfr(G,omega_d))==0

I was unable to find a function within Matlab which would help me solve this problem, except evalfr. I made an algorithm that actually does help, however it seems quite inelegant and 'brute force'. What it does is starting at omega_d=0 and increasing the frequency one digit at a time until the real part of the responce changes in sign. Then it saves the current digit and starts all over again at the next digit.

Do you know any better way to solve this problem? Just to be clear: I know that the given example here is quite easy to solve analytically. So let us also assume, that my system 'G' is a lot more complicated (including several internal delays, parallel and seriell connecting of different subsystems, etc.).

Thanks in advance...

David

PS - my code:

function omega_d=zreal(G)
% computes the smallest frequency omega_d at which the frequency response of a given system 'G' 
% is pi/4 or 3*pi/4: real(evalfr(G,omega_d))=0

% check for starting conditions
if real(evalfr(G,0))==0
  disp('Warning: real(evalfr(G,0))==0');
  omega_d = 0;
  return
elseif real(evalfr(G,0))<0
  sig=-1;
else
  sig=1;
end

% +---------------------------------+
% |search for most significant digit|
% +---------------------------------+
% increase omega_d until the sign of real(evalfr(G,omega_d)) changes,
% meaning that the previous value was the maximum value of the most significant digit
expo=11;
val=9;
err=sig;
while sig*err>0
  % compute next frequency to evaluate
  if val==9
      expo=expo-1;
      val=1;
  else
      val=val+1;
  end

    % compute frequency response
    tmp=val*10^-expo;
    err=real(evalfr(G,tmp*1i));
end

% undo last change of frequency
if val==1
  expo=expo+1;
  val=9;
else
  val=val-1;
end
omega_d=val*10^-expo; % most significant digit found

% +----------------------------+
% |search for subsequent digits|
% +----------------------------+
% increase value of current digit until
% (a) sign of real(evalfr(G,omega_d)) changes, thus last value was the maximum or
% (b) the value reaches 9, which is the maximum, too
% then change to next digit
expo=expo+1;
val=1;
maxExpo=expo+16; % maximum 17 significant digits with double precision (one digit already found)
while expo<maxExpo
  tmp=omega_d+val*10^-expo;
  err=sig*real(evalfr(G,tmp*1i));

    % check if maximum value of current digit is found
    if err<0 % sign of real(evalfr(G,omega_d)) changed, ergo...             
        omega_d=omega_d+(val-1)*10^-expo; % ...undo last change of current digit, save it and...
        expo=expo+1; % ...switch to next digit
        val=1;
    elseif val==9 % current digit is at its maximum, ergo...
        omega_d=omega_d+val*10^-expo;% ...save current digit and...
        expo=expo+1; % ...switch to next digit
        val=1;
    else % increase value of current digit
        val=val+1;
    end
end

Answer 1

[Craig]

If you have the Optimization Toolbox you might try the function FSOVLE which is used to solve equations of the form F(x) = 0.

Another option that may work is to use the bode command to do the intial gridding of the frequency range

[m,p,w] = bode(G);

Then look for a sign change in real(G) over this gridding. Then refine the grid and do a search like your performing in your code.

[David T_]

Hi,

thank you for your answer, I looked into using FSOLVE (and FZERO, too) and was quite fond of it. However I'm not sure about the following issues:

1) The systems I analyse are defined in statespace notation and usually contain internal delays. I was quite happy that most Matlab functions I use here accept the ss-class directly, but FSOLVE does not. (no wonder here, since FSOLVE is not part of the control toolbox) Redefining the systems to fit into FSOLVE seems quite tedious to me.

2) Since most of the systems I analyse have some delay in them their nyquist plot crosses the imaginary axis quite often. How do I ensure that fsolve returns the minimum(!) value at which a crossing occurs?

Using the bode command did cross my mind, too, but it feels as 'brute-force' as what I do in my code...