advanced frequency domain analysis of SISO models

2 次查看（过去 30 天）
David T_
David T_ 2012-5-30
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

请先登录，再回答此问题。

采纳的回答

Craig
Craig 2012-7-18
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.
  1 个评论
David T_
David T_ 2012-8-30
编辑：David T_ 2012-8-30
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...

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Frequency-Domain Analysis 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by