Re: OCST: Gain Margin

[Lukas Reichlin]

From: address@hidden
To: address@hidden; address@hidden
Subject: Re: OCST: Gain Margin
Date: Fri, 3 Jul 2009 12:19:50 +0200

---

From: address@hidden
To: address@hidden
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
address@hidden
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