-
Notifications
You must be signed in to change notification settings - Fork 438
System norms #960
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.
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
System norms #960
Conversation
* Added control/tests/sysnorm_test.py
@henriks76 It looks like there are some problems when slycot is installed. Can you take a look at the failed unit tests above and see what is going on. From a quick look, there may be an inconsistency in |
* Added additional tests and warning messages for systems with poles close to stability boundary * Added __all__ * Added more comments
Thanks @murrayrm . Odd, I have slycot installed but did not get this fault in the tests. The problem seems to occur when the poles are very close to the stability boundary. I think the poles must first all have been evaluated as being stable, but then the Lyapunov equation is so ill-conditioned that the solution is not positive semi-definite, resulting in a negative number inside np.sqrt giving NaN. I have added more tests and also warning messages. Let's see if that solves the problem... |
control/tests/sysnorm_test.py
Outdated
s = ct.tf('s') | ||
G1 = 1/(s+1) | ||
assert np.allclose(ct.norm(G1, p='inf', tol=1e-9), 1.0, rtol=1e-09, atol=1e-08) # Comparison to norm computed in MATLAB | ||
assert np.allclose(ct.norm(G1, p=2), 0.707106781186547, rtol=1e-09, atol=1e-08) # Comparison to norm computed in MATLAB |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think these tolerances are way too tight. On a non x86_64 system, e.g. 32-bit intel, aarch64, PowerPC or even just a different LAPACK/BLAS implementation can give you differently rounded intermediate and final results.
control/sysnorm.py
Outdated
gamu = max(1.0, 2.0*gaml) # Candidate upper bound | ||
Ip = np.eye(len(D)) | ||
|
||
while any(np.isclose(la.eigvals(_Hamilton_matrix(gamu)).real, 0.0)): # Find actual upper bound |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Have you considered using SLICOT-Routines (already or to be wrapped by Slycot) optimized for certain Hamiltonian matrices instead of the generic scipy eigenvalue solver? (I didn't check if the SLICOT routines are applicable here)
control/sysnorm.py
Outdated
poles_real_part = G.poles().real | ||
if any(np.isclose(poles_real_part, 0.0)): | ||
if print_warning: | ||
print("Warning: Poles close to, or on, the imaginary axis. Norm value may be uncertain.") |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Don't use print
here and in the error cases below. For warnings use warnings
, if it is an error, please keep raising the exception. And be more specific than just Exception
.
FWIW, there also 2 SLICOT / Slycot routines
ab13bd is part of Slycot 0.5.4, but there was a problem python-control/Slycot#199
ab13dd is already part of python-control, called linfnorm() in |
* Added control/tests/sysnorm_test.py
* Added additional tests and warning messages for systems with poles close to stability boundary * Added __all__ * Added more comments
…ntrol into system-norms
Thanks for the feedback @bnavigator and @KybernetikJo . Incorporated |
>>> ct.norm(Gc,'inf',tol=1e-10) | ||
1.0000000000582077 | ||
""" | ||
G = ct.ss(system) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You might want to add a type check similar to
if not isinstance(sys, (StateSpace, TransferFunction)):
raise TypeError('Parameter ``sys``: must be a ``StateSpace`` or ``TransferFunction``)')
befor calling
G = ct.ss(sys)
control/sysnorm.py
Outdated
|
||
#------------------------------------------------------------------------------ | ||
|
||
def norm(system, p=2, tol=1e-10, print_warning=True, use_slycot=True): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The optional argument
use_slycot=True
could be replaced by an optional argument
method=None
as in https://github.com/python-control/python-control/blob/main/control/mateqn.py.
Hint: You may want to reuse the following method from mateqn.py:
def _slycot_or_scipy(method):
if method == 'slycot' or (method is None and slycot_check()):
return 'slycot'
elif method == 'scipy' or (method is None and not slycot_check()):
return 'scipy'
else:
raise ControlArgument("Unknown method %s" % method)
The argument method
has three options method=[None,"slycot","scipy"] and the following behavior:
- If method is None and slycot is installed, then use slycot.
- If method is None and slycot is not installed, then fall back on scipy. (Fallback)
- If method=="slycot" then try to use slycot or raise an Exception when slycot is not available. (No Fallback)
- If method=="scipy" then use scipy. (scipy is always available)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for the tips. It should be implemented now.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Please don't duplicate code. Reuse it, maybe move it to a common place, but don't repeat yourself (or previous contributors)
* Added control/tests/sysnorm_test.py
* Added additional tests and warning messages for systems with poles close to stability boundary * Added __all__ * Added more comments
* Added additional tests and warning messages for systems with poles close to stability boundary * Added __all__ * Added more comments
…ntrol into system-norms
@@ -5,7 +5,7 @@ | |||
import os | |||
import matplotlib.pyplot as plt # Grab MATLAB plotting functions | |||
from control.matlab import * # MATLAB-like functions | |||
from scipy import pi | |||
from numpy import pi |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This has already been resolved in #965. Please rebase.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ok, I tried again.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Nope.
- Fast forward your main branch to upstream/main
- In system-norms use
git rebase -i main
https://github.com/python-control/python-control/wiki/How-to-contribute-with-a-pull-request
This gets rid of the extra commits which are already in the main branch. While you are at it, you can also squash some of the commits and use better messages than "Update sysnorm.py"
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I did steps 7-9:
... `git checkout main; git fetch --all; git merge origin/main
git push
and then bring those changes into your branch:
git checkout ; git rebase main`
Is it step 12 I should do?
git checkout main; git fetch --all; git merge upstream/main
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I am afraid your multiple merges and manual fixes like 7a5af50 have created a situation which git cannot resolve automatically:
[ben@skylab:~/src/python-control]% git status [0]
On branch system-norms
Your branch is up to date with 'henriks76/system-norms'.
nothing to commit, working tree clean
[ben@skylab:~/src/python-control]% git rebase main [0]
Auto-merging control/sysnorm.py
CONFLICT (add/add): Merge conflict in control/sysnorm.py
error: could not apply fd076c5... New function for LTI system norm computation
hint: Resolve all conflicts manually, mark them as resolved with
hint: "git add/rm <conflicted_files>", then run "git rebase --continue".
hint: You can instead skip this commit: run "git rebase --skip".
hint: To abort and get back to the state before "git rebase", run "git rebase --abort".
Could not apply fd076c5... New function for LTI system norm computation
[ben@skylab:~/src/python-control]% git status [1]
interactive rebase in progress; onto 85328d9
Last commands done (10 commands done):
pick 49f7e5f Update sysnorm.py
pick fd076c5 New function for LTI system norm computation
(see more in file .git/rebase-merge/done)
Next commands to do (9 remaining commands):
pick 0fe2c57 * Updated documentation of function norm * Added control/tests/sysnorm_test.py
pick cdf4bab Update sysnorm.py
(use "git rebase --edit-todo" to view and edit)
You are currently rebasing branch 'system-norms' on '85328d9'.
(fix conflicts and then run "git rebase --continue")
(use "git rebase --skip" to skip this patch)
(use "git rebase --abort" to check out the original branch)
Unmerged paths:
(use "git restore --staged <file>..." to unstage)
(use "git add <file>..." to mark resolution)
both added: control/sysnorm.py
no changes added to commit (use "git add" and/or "git commit -a")
You have to go through it step by step, because only you know what is relevant.
* Added additional tests and warning messages for systems with poles close to stability boundary * Added __all__ * Added more comments
…ntrol into system-norms
@henriks76, you can use https://github.com/bnavigator/python-control/tree/system-norms2 and re-open if you want. I basically picked the relevant commits in order, without the merge and duplication stuff. |
Or do
|
Thanks @bnavigator, this was very helpful! I squashed some commits in https://github.com/henriks76/python-control as well. Or is it better to continue on https://github.com/bnavigator/python-control/tree/system-norms2 ? |
Both would be fine. Note:
diff --git a/control/tests/sysnorm_test.py b/control/tests/sysnorm_test.py
index 917e98d..74e8de0 100644
--- a/control/tests/sysnorm_test.py
+++ b/control/tests/sysnorm_test.py
@@ -13,35 +13,35 @@ def test_norm_1st_order_stable_system():
"""First-order stable continuous-time system"""
s = ct.tf('s')
G1 = 1/(s+1)
- assert np.allclose(ct.norm(G1, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB
- assert np.allclose(ct.norm(G1, p=2), 0.707106781186547) # Comparison to norm computed in MATLAB
+ assert np.allclose(ct.norm(G1, p='inf', tol=1e-9, print_warning=False), 1.0) # Comparison to norm computed in MATLAB
+ assert np.allclose(ct.norm(G1, p=2, print_warning=False), 0.707106781186547) # Comparison to norm computed in MATLAB
Gd1 = ct.sample_system(G1, 0.1)
- assert np.allclose(ct.norm(Gd1, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB
- assert np.allclose(ct.norm(Gd1, p=2), 0.223513699524858) # Comparison to norm computed in MATLAB
+ assert np.allclose(ct.norm(Gd1, p='inf', tol=1e-9, print_warning=False), 1.0) # Comparison to norm computed in MATLAB
+ assert np.allclose(ct.norm(Gd1, p=2, print_warning=False), 0.223513699524858) # Comparison to norm computed in MATLAB
def test_norm_1st_order_unstable_system():
"""First-order unstable continuous-time system"""
s = ct.tf('s')
G2 = 1/(1-s)
- assert np.allclose(ct.norm(G2, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB
- assert ct.norm(G2, p=2) == float('inf') # Comparison to norm computed in MATLAB
+ assert np.allclose(ct.norm(G2, p='inf', tol=1e-9, print_warning=False), 1.0) # Comparison to norm computed in MATLAB
+ assert ct.norm(G2, p=2, print_warning=False) == float('inf') # Comparison to norm computed in MATLAB
Gd2 = ct.sample_system(G2, 0.1)
- assert np.allclose(ct.norm(Gd2, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB
- assert ct.norm(Gd2, p=2) == float('inf') # Comparison to norm computed in MATLAB
+ assert np.allclose(ct.norm(Gd2, p='inf', tol=1e-9, print_warning=False), 1.0) # Comparison to norm computed in MATLAB
+ assert ct.norm(Gd2, p=2, print_warning=False) == float('inf') # Comparison to norm computed in MATLAB
def test_norm_2nd_order_system_imag_poles():
"""Second-order continuous-time system with poles on imaginary axis"""
s = ct.tf('s')
G3 = 1/(s**2+1)
- assert ct.norm(G3, p='inf') == float('inf') # Comparison to norm computed in MATLAB
- assert ct.norm(G3, p=2) == float('inf') # Comparison to norm computed in MATLAB
+ assert ct.norm(G3, p='inf', print_warning=False) == float('inf') # Comparison to norm computed in MATLAB
+ assert ct.norm(G3, p=2, print_warning=False) == float('inf') # Comparison to norm computed in MATLAB
Gd3 = ct.sample_system(G3, 0.1)
- assert ct.norm(Gd3, p='inf') == float('inf') # Comparison to norm computed in MATLAB
- assert ct.norm(Gd3, p=2) == float('inf') # Comparison to norm computed in MATLAB
+ assert ct.norm(Gd3, p='inf', print_warning=False) == float('inf') # Comparison to norm computed in MATLAB
+ assert ct.norm(Gd3, p=2, print_warning=False) == float('inf') # Comparison to norm computed in MATLAB
def test_norm_3rd_order_mimo_system():
"""Third-order stable MIMO continuous-time system"""
@@ -55,9 +55,9 @@ def test_norm_3rd_order_mimo_system():
[-0.863652821988714, -1.214117043615409, -0.006849328103348]])
D = np.zeros((2,2))
G4 = ct.ss(A,B,C,D) # Random system generated in MATLAB
- assert np.allclose(ct.norm(G4, p='inf', tol=1e-9), 4.276759162964244) # Comparison to norm computed in MATLAB
- assert np.allclose(ct.norm(G4, p=2), 2.237461821810309) # Comparison to norm computed in MATLAB
+ assert np.allclose(ct.norm(G4, p='inf', tol=1e-9, print_warning=False), 4.276759162964244) # Comparison to norm computed in MATLAB
+ assert np.allclose(ct.norm(G4, p=2, print_warning=False), 2.237461821810309) # Comparison to norm computed in MATLAB
Gd4 = ct.sample_system(G4, 0.1)
- assert np.allclose(ct.norm(Gd4, p='inf', tol=1e-9), 4.276759162964228) # Comparison to norm computed in MATLAB
- assert np.allclose(ct.norm(Gd4, p=2), 0.707434962289554) # Comparison to norm computed in MATLAB
+ assert np.allclose(ct.norm(Gd4, p='inf', tol=1e-9, print_warning=False), 4.276759162964228) # Comparison to norm computed in MATLAB
+ assert np.allclose(ct.norm(Gd4, p=2, print_warning=False), 0.707434962289554) # Comparison to norm computed in MATLAB
|
Yes, I stopped the printout of warnings to not clutter the output of pytest. If it is better to put it back in I can do it, of course. |
Remind me, what was the warning about, again? |
Many of the tests gave these warnings. |
I'd prefer to use a |
Thanks @bnavigator, I now use a |
Continued in #971 |
Hello,
I have written a new function norm to compute H_2 and L_infinity system norms of LTI objects.
Please do let me know if you have any suggested improvements.
Best regards,
Henrik