Re: OCST: Gain Margin

[paramaniac]

Dear Mr Favatella

I am writing to you in behalf of the Octave Control Systems Toolbox. I

wanted to ask if it is possible to add an equivalent of the Matlab function

[gamma, phi, w_gamma, w_phi] = margin(sys) in a future version of your

package.

[...]

I hope that I am writing to the right person and I am looking forward to

hearing from you.

Right person.
In the future, when contributing code to octave-forge, I think it is
good to at least cc octave-dev@....

I hope to take care of this next week.

Regards

Lukas Reichlin

[...]

Thanks for the attached function.


Regards,
Luca Favatella

Excellent :-) Below I added the discussion from help-octave@... together with a newer version of the function.

Regards,

Lukas Reichlin

---

From: lukas.reichlin@...
To: help-octave@...
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@...
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

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

------------------------------------------------------------------------------

_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Luca Favatella]

On 03/07/2009, Lukas Reichlin <lukas.reichlin@...> wrote:
[...]
> Excellent :-) Below I added the discussion from help-octave@...
> together with a newer version of the function.

Thanks.

Committed with some style fixes.

Please see
http://octave.svn.sourceforge.net/viewvc/octave?view=rev&revision=5988
and control if I forgot something.


[...]
> 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.

Mmm... It looks like you are more expert than me on this.

Can you please provide a simple test for this?
You can use this file as a reference (see the last lines)
http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/control/inst/place.m?revision=5857&view=markup


> Regards,
> Lukas

Thanks,
Luca Favatella

------------------------------------------------------------------------------
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[paramaniac]

On 06.07.2009, at 11:13, Luca Favatella wrote:

> On 03/07/2009, Lukas Reichlin <lukas.reichlin@...> wrote:
> [...]
>> Excellent :-) Below I added the discussion from help-
>> octave@...
>> together with a newer version of the function.
>
> Thanks.
>
> Committed with some style fixes.
>
> Please see
> http://octave.svn.sourceforge.net/viewvc/octave?view=rev&revision=5988
> and control if I forgot something.
>
>
> [...]
>> 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.
>
> Mmm... It looks like you are more expert than me on this.
>
> Can you please provide a simple test for this?
> You can use this file as a reference (see the last lines)
> http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/control/inst/place.m?revision=5857&view=markup
>
>
>> Regards,
>> Lukas
>
> Thanks,
> Luca Favatella

I found this hint about Matlab's implementation at the bottom of the  
following page.

http://www.mathworks.com/access/helpdesk/help/toolbox/control/index.html?/access/helpdesk/help/toolbox/control/ref/margin.html&http://www.google.com/search 
?q=%22gain%20margin%20algorithm%22&sourceid=mozilla2&ie=utf-8&oe=utf-8

On the same page, they are backing my doubts about the accuracy of the  
"bode technique" we used so far.
Is there someone skilled out there who is able to implement the real  
thing?

Regards,
Lukas


------------------------------------------------------------------------------
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev