OCST: Gain Margin

[Lukas Reichlin]

> 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

## Copyright (C) 2004 Doug Stewart
##
## 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     bode1w (n,d)
%n= numerator
%d=denominator
%
% This was written to help my students with Gain Margin and Phase Margin
%GM and PM lines are added to the Bode plot and
%GM and PM values are printed on the screen.


## bode1w

## Author: doug
##
## 2004-09-27 doug
## * Initial revision
## Rev #1  Nov 29 2004 Doug Stewart
## Rev #2 Sept 2005 Doug Stewart
##              - Added the If for negative phase steps.
##              - and the DO loop
## Rev #3 Sept 2005 Doug Stewart
##              - used unwrap.m


function [gamma, phi, w_gamma, w_phi] = 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(ph);
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/af76ffe1/attachment.html