Skip to content

Fixes to margin calculation #68

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom
Merged
merged 8 commits into from
May 30, 2016
Merged

Fixes to margin calculation #68

merged 8 commits into from
May 30, 2016

Conversation

repagh
Copy link
Member

@repagh repagh commented Jul 11, 2015

It seems that back when I implemented a new margin calculation, there were some bugs left. In this pull request, I fixed the following:

  • the 'stability margin' (concept unknown in Matlab??) was fixed. We encountered a test case here in Delft that gives 3 solutions & did not match the analytical and numerical approach. Re-implemented the numerical approach, to look through the complete frequency vector and find all crossings (derivative of distance to -1 is zero), and weed out the maxima. Also corrected the analytical solution; this had a typo in there
  • checked and fixed gain margin and phase crossover frequency calculations
  • added an extra test case, and extended testing.

In an sideways related issue, I found that the "smooth" parameter was dropped from an FRD when the FRD was combined with a linear system, multiplied with a gain, etc. Fixed that to keep smooth FRD smooth, and go nonsmooth as soon as a combination with an FRD without the smooth parameter is requested.

repagh and others added 7 commits July 8, 2015 01:10
@repagh
Keep the "smooth" parameter when combining FRD functions with other s…
6523270 
…ystems
@repagh
Fix the margin calculation for FRD/data based margins.
ae09104 
TODO: verify the stability margin
@repagh
corrected the stability margin calculations
7eefe8b
@repagh
still imperfect
f309232
@repagh
further attempts at fixing
815ca2d
@repagh
corrected margin calculations; extended test to do all-margin calcula…
cd30a00 
…tions

check margin calculation on analytic tf with calculation on numeric
(interpolation with FRD) variant. Extended test cases
@repagh
Merge remote-tracking branch 'upstream/master'
04ed40d 
solved the conflict in the margin_test routine
updated the doc for margin to match order of returned values
@repagh
Copy link
Member Author

repagh commented Jul 11, 2015

I think this is OK now. I merged in roryyorke's fix for the matlab compatibility routine.
I could merge the request it in myself, but since I changed quite a bit I would appreciate another pair of eyes.

@murrayrm
Copy link
Member

I'll try to have a look through the changes and then merge. Probably won't get to it until later in the week, so if someone else has time please go ahead!

@roryyorke
Copy link
Contributor

Hi,

Not sure this is all that helpful; apologies if not.

This code [run in Ubuntu 14.04, Python 2.7.6, numpy 1.8.2, and with painyeph/master]

from __future__ import division
import numpy as np
import control

g = control.tf([-1 / 10, 1], [1, 1])
print('k=9.999 pole: {}'.format(control.feedback(9.999 * g, 1).pole()))
print('k=10.001 pole: {}'.format(control.feedback(10.001 * g, 1).pole()))
gm, _, _, _, _, _ = control.stability_margins(g)
print('gm: {}'.format(gm))

g2 = control.tf(2, [1, -1])
# g(s) + 1 = (s+1)/(s-1); |g(jw)+1| = 1 for all w
print('unstable o/l sys c/l poles: {}'.format(control.feedback(g2, 1).pole()))
w = np.logspace(-2, 2, 101)
m, p, w = g2.freqresp(w)
gjw = (m * np.exp(1j * p)).flatten()
distn1 = abs(gjw + 1)
print('max | |g(jw)+1| - 1|: {}'.format(max(abs(distn1-1))))
_, _, sm, _, _, _ = control.stability_margins(g)
print('stab margin: {}'.format(sm))

produces this result:

k=9.999 pole: [-109990.00000013]
k=10.001 pole: [ 110010.00000001]
gm: None
unstable o/l sys c/l poles: [-1.]
max | |g(jw)+1| - 1|: 4.4408920985e-16
stab margin: None

I stared at _polyimsplit for a while, but it's not obvious to me
what it does, I'm afraid. [EDIT: ah, 1j**(0+4n)=1, 1j**(1+4n)=1j,
etc. -- I like this.] I think one has to also consider the point
lim_{w->inf} g(jw) to find the above margin, but it's not obvious how
this fits in with the how the margins are currently evaluated.
Incidentally, Octave 3.8.1 does not find it either; Matlab does.

Regarding the stability margin: my possibly unrealistic test case has
constant distance from -1, and so no local minima to find.

Separately on the stability margin, I believe this can be evaluated as
1/||1/(1+g(s))||_inf [EDIT: added missing )], and painyeph submitted a
pull request to add
SLICOT's H-infinity evaulator AB13DD to jgoppert/slycot.

This idea can be tested with these two Octave scripts:

% Octave 3.8.1 doesn't have rss
% save this as rss.m
function g=rss(nx,nu,ny)

  ar= randn([nx,nx]);
  [vr,dr]= eig(ar);
  d= complex(-abs(real(dr)), imag(dr));

  a= real(vr*d/vr);
  b= randn([nx,nu]);
  c= randn([ny,nx]);
  d= randn([ny,nu]);

  g= ss(a,b,c,d);


% margstab2.m
g = tf(2,[1,-1]);

w = logspace(-3,3,501);

[gjwr,gjwi] = nyquist(g,w);
gjw = complex(gjwr,gjwi);

[mindist,minidx] = min(abs(1+gjw));
minw = w(minidx);
[rmindist_hinf,minw_hinf] = norm(1/(g+1),inf);
mindist_hinf = 1/rmindist_hinf;
fprintf('naive: mindist = %g @ %g r/s\n',mindist,minw);
fprintf('hinf:  mindist = %g @ %g r/s\n',mindist_hinf,minw_hinf);

which produce output like this:

>> margstab2
naive: mindist = 0.00425027 @ 1000 r/s
hinf:  mindist = 0.00102006 @ Inf r/s
>> margstab2
naive: mindist = 0.753816 @ 1000 r/s
hinf:  mindist = 0.753806 @ Inf r/s
>> margstab2
naive: mindist = 2.26894 @ 13.8038 r/s
hinf:  mindist = 2.27418 @ Inf r/s
>> margstab2
naive: mindist = 1.11145 @ 1000 r/s
hinf:  mindist = 1.11143 @ Inf r/s
>> margstab2
naive: mindist = 2.28208 @ 1000 r/s
hinf:  mindist = 2.28207 @ Inf r/s
>> margstab2
naive: mindist = 0.420669 @ 3.87258 r/s
hinf:  mindist = 0.420664 @ 3.94223 r/s
>> margstab2
naive: mindist = 0.52768 @ 1000 r/s
hinf:  mindist = 0.527648 @ Inf r/s

Regards,

Rory

@murrayrm murrayrm merged commit c14b880 into python-control:master May 30, 2016
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers
No reviews
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
3 participants
@repagh @murrayrm @roryyorke