OCST: Gain Margin

[Lukas Reichlin]

> From: lukas.reichlin at swissonline.ch
> To: dastew at sympatico.ca; help-octave at octave.org
> Subject: Re: OCST: Gain Margin
> Date: Fri, 3 Jul 2009 12:19:50 +0200
>
>
> From: lukas.reichlin at swissonline.ch
> To: help-octave at octave.org
> Subject: Re: OCST: Gain Margin
> Date: Fri, 3 Jul 2009 09:48:16 +0200
>
> Dear Octave Community,
>
> I wrote an objective function to optimize an Aström/Hägglund PID
> Controller numerically by fminsearch. Versions for Octave and Matlab/
> SImulink can be found here:
>
> http://n.ethz.ch/student/lukasre/download/optiPID/
>
> The Octave control package is quite limited. However, I managed to get
> along quite easily except for one thing: the gain margin of a system.
> In Matlab, there's [gamma, phi, w_gamma, w_phi] = margin(sys). I
> couldn't think of a way to calculate the gain margin numerically. Is
> there any control systems engineer out there who knows how to
> implement an algorithm for the problem? Help would be very  
> appreciated.
>
> Regards,
> Lukas
>
> BTW: Despite I implemented the routine quite differently, the
> conformity of the results between Matlab and Octave is simply
> fantastic :-)
> _______________________________________________
> Help-octave mailing list
> Help-octave at octave.org
> https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
>
>
> If you look at what I did you might get some ideas.
>
> http://dougs.homeip.net/qtoctave/bode1w.m
>
> Doug Stewart
>
> Thanks for your help! If I understood your code correctly, the  
> following code should do the job:
>
> Regards,
> Lukas
>
>
> I cleaned this up and added some comments. If it works the way you  
> want, then
> we should put it in the controls tool box.
>
>
> Doug Stewart
>
> Thank you. I corrected some minor typos and added my family name :-)  
> I think this hack needs some serious testing before it's added to  
> anything ;-) Especially because the accuracy relies on how tightly  
> spaced the frequency vector of bode(sys) is.
>
> Regards,
> Lukas
>
>
> Thanks Lucas
>
> If you think it is important to use a finer resolution for the  
> frequencies then we could:
> 1) find the frequencies the way we are now.
> 2) call bode again with a frequency vector that is just in a narrow  
> band around these frequencies;
> 3) recalculate the GM and PM
>
> Do you think this is a good plan?
> Doug


I've done it for the gain margin. I don't know what happens if w(ix-1)  
or w(ix+1) doesn't exist ...

Another idea would be to calculate the roots of a polynomial. I have  
to think about that.

Regards,
Lukas



## Copyright (C) 2009 Doug Stewart and Lukas Reichlin
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  
USA

%
%function     function [gamma, phi, w_gamma, w_phi] = margin(sys);
%   gamma	=  Gain Margin  ( as gain not dbs)
%   w_gamma	=  radian frequency for the Gain Margin
%   phi 	=  Phase Margin ( in degrees)
%   w_phi	=  radian frequency for the Phase margin
%
% This function calls bode - see  help bode for details about bode.
% This Function was written for continuous time systems

function [gamma, phi, w_gamma, w_phi] = margin(sys);

%if(! isstruct(sys) )
%	error("The input must be a system type variable.	margin(sys)")

% Get Frequency Response
[mag, pha, w] = bode(sys);

% fix the phase wrap around
phw = unwrap(pha * pi/180);
pha = phw * 180/pi;

% find the Gain Margin and its location
[x, ix] = min(abs(pha + 180));

if (x > 1)
	gamma = "INF";
	ix = length(pha);
else
	gamma = 1/mag(ix);
endif

w_gamma = w(ix);


% Fine Tuning of Gain Margin Result
w_low = w(ix-1);
w_high = w(ix+1);
n_steps = 10000;
dw = (w_high - w_low) / n_steps;

w_fine = w_low : dw : w_high;

[mag_fine, pha_fine, w_fine] = bode(sys, w_fine);

[x_fine, ix_fine] = min(abs(pha_fine + 180));

gamma = 1/mag_fine(ix_fine);
w_gamma = w_fine(ix_fine)


% find the Phase Margin and its location
[pmx, ipmx] = min(abs(mag - 1));

phi = 180 + pha(ipmx);
w_phi = w(ipmx);

endfunction













-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://www-old.cae.wisc.edu/pipermail/help-octave/attachments/20090703/38059872/attachment.html