OCST: Gain Margin

[paramaniac]

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

[Doug Stewart-4]

On Thu, Jul 2, 2009 at 5:51 PM, Lukas Reichlin <lukas.reichlin@...> wrote:

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

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

[paramaniac]

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

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

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

[dastew]

Some parts of this message have been removed. Learn more about Nabble's security policy.

---

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

## Copyright (C) 2009 Doug Stewart and Lukas

## 

## 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(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

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

[paramaniac]

---

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

[dastew]

Some parts of this message have been removed. Learn more about Nabble's security policy.

---

From: lukas.reichlin@...
To: dastew@...; help-octave@...
Subject: Re: OCST: Gain Margin
Date: Fri, 3 Jul 2009 12:19:50 +0200

---

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

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 

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

[paramaniac]

From: lukas.reichlin@...
To: dastew@...; help-octave@...
Subject: Re: OCST: Gain Margin
Date: Fri, 3 Jul 2009 12:19:50 +0200

---

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

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

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