OCST: Gain Margin

[Lukas Reichlin]

> 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


## 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);

% 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/7713ab33/attachment.html