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-...@lists.sourceforge.net.
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-oct...@octave.org
together with a newer version of the function.
Regards,
Lukas Reichlin
From: lukas.reich...@swissonline.ch
To: help-oct...@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-oct...@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
_______________________________________________
Help-octave mailing list
help-oct...@octave.org
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
------------------------------------------------------------------------------
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
