Skip to content

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

Closed
wants to merge 27 commits into from
Closed

System norms #960

wants to merge 27 commits into from

Conversation

henriks76
Copy link

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

Henrik Sandberg added 2 commits January 8, 2024 09:43
@murrayrm
Copy link
Member

murrayrm commented Jan 8, 2024

@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 inf and nan that might be something we need to fix elsewhere (but that will be required for your function to work correctly).

* Added additional tests and warning messages for systems with poles close to stability boundary
* Added __all__
* Added more comments
@henriks76
Copy link
Author

@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 inf and nan that might be something we need to fix elsewhere (but that will be required for your function to work correctly).

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

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
Copy link
Contributor

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.

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
Copy link
Contributor

@bnavigator bnavigator Jan 9, 2024

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)

@coveralls
Copy link

coveralls commented Jan 9, 2024

Coverage Status

coverage: 94.431% (-0.4%) from 94.784%
when pulling 14fce67 on henriks76:system-norms
into a8a54d1 on python-control:main.

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.")
Copy link
Contributor

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.

@KybernetikJo
Copy link
Contributor

FWIW, there also 2 SLICOT / Slycot routines

  • ab13bd -> H2 or L2 norm of a state space system

ab13bd is part of Slycot 0.5.4, but there was a problem python-control/Slycot#199
It should be fixed in Slycot 0.6.0!

  • ab13dd -> L-infinity norm of a state space system

ab13dd is already part of python-control, called linfnorm() in
https://github.com/python-control/python-control/blob/main/control/statesp.py
but missing in the sphinx documentation
https://github.com/python-control/python-control/blob/main/doc/control.rst

henriks76 and others added 8 commits January 21, 2024 19:13
* 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
* Use of warnings package.
* Use routine statesp.linfnorm when slycot installed.
* New routine internal _h2norm_slycot when slycot is installed.
@henriks76
Copy link
Author

Thanks for the feedback @bnavigator and @KybernetikJo . Incorporated warnings and slycot routines when available now.

>>> ct.norm(Gc,'inf',tol=1e-10)
1.0000000000582077
"""
G = ct.ss(system)
Copy link
Contributor

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)


#------------------------------------------------------------------------------

def norm(system, p=2, tol=1e-10, print_warning=True, use_slycot=True):
Copy link
Contributor

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)

Copy link
Author

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.

Copy link
Contributor

@bnavigator bnavigator Jan 28, 2024

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)

bnavigator and others added 7 commits January 22, 2024 21:52
Henrik Sandberg added 5 commits January 28, 2024 16:05
* Use of warnings package.
* Use routine statesp.linfnorm when slycot installed.
* New routine internal _h2norm_slycot when slycot is installed.
* Added additional tests and warning messages for systems with poles close to stability boundary
* Added __all__
* Added more comments
* type check when calling ct.norm
* metod argument in ct.norm (slycot or scipy)
@@ -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
Copy link
Contributor

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.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I tried again.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nope.

  1. Fast forward your main branch to upstream/main
  2. 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"

Copy link
Author

@henriks76 henriks76 Jan 28, 2024

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

Copy link
Contributor

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:

image

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

Henrik Sandberg added 2 commits January 28, 2024 18:30
* Added additional tests and warning messages for systems with poles close to stability boundary
* Added __all__
* Added more comments
@henriks76 henriks76 closed this Jan 28, 2024
@bnavigator
Copy link
Contributor

@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.

@bnavigator
Copy link
Contributor

Or do git rebase -i main and use this picklist:

pick 7ace4bc New function for LTI system norm computation
pick 5103448 * Updated documentation of function norm * Added control/tests/sysnorm_test.py
pick 4fb252e Update sysnorm.py
pick 108817c Do not track changes in VS Code setup.
pick 6683eb3 Lowered tolerances in tests.
pick 32d38bf Added: * Use of warnings package. * Use routine statesp.linfnorm when slycot installed. * New routine internal _h2norm_slycot when slycot is installed.
pick 793f0d6 Added: * type check when calling ct.norm * metod argument in ct.norm (slycot or scipy)
pick fe22e33 Use function _slycot_or_scipy from ct.mateqn

# Rebase 85328d9..fe22e33 onto 85328d9 (20 commands)
#
# Commands:
# p, pick <commit> = use commit
# r, reword <commit> = use commit, but edit the commit message
# e, edit <commit> = use commit, but stop for amending
# s, squash <commit> = use commit, but meld into previous commit
# f, fixup [-C | -c] <commit> = like "squash" but keep only the previous
#                    commit's log message, unless -C is used, in which case
#                    keep only this commit's message; -c is same as -C but
#                    opens the editor
# x, exec <command> = run command (the rest of the line) using shell
# b, break = stop here (continue rebase later with 'git rebase --continue')
# d, drop <commit> = remove commit
# l, label <label> = label current HEAD with a name
# t, reset <label> = reset HEAD to a label
# m, merge [-C <commit> | -c <commit>] <label> [# <oneline>]
#         create a merge commit using the original merge commit's
#         message (or the oneline, if no original merge commit was
#         specified); use -c <commit> to reword the commit message
# u, update-ref <ref> = track a placeholder for the <ref> to be updated
#                       to this position in the new commits. The <ref> is
#                       updated at the end of the rebase
#
# These lines can be re-ordered; they are executed from top to bottom.
#
# If you remove a line here THAT COMMIT WILL BE LOST.
#
# However, if you remove everything, the rebase will be aborted.
#

@henriks76
Copy link
Author

Or do git rebase -i main and use this picklist:

pick 7ace4bc New function for LTI system norm computation
pick 5103448 * Updated documentation of function norm * Added control/tests/sysnorm_test.py
pick 4fb252e Update sysnorm.py
pick 108817c Do not track changes in VS Code setup.
pick 6683eb3 Lowered tolerances in tests.
pick 32d38bf Added: * Use of warnings package. * Use routine statesp.linfnorm when slycot installed. * New routine internal _h2norm_slycot when slycot is installed.
pick 793f0d6 Added: * type check when calling ct.norm * metod argument in ct.norm (slycot or scipy)
pick fe22e33 Use function _slycot_or_scipy from ct.mateqn

# Rebase 85328d9..fe22e33 onto 85328d9 (20 commands)
#
# Commands:
# p, pick <commit> = use commit
# r, reword <commit> = use commit, but edit the commit message
# e, edit <commit> = use commit, but stop for amending
# s, squash <commit> = use commit, but meld into previous commit
# f, fixup [-C | -c] <commit> = like "squash" but keep only the previous
#                    commit's log message, unless -C is used, in which case
#                    keep only this commit's message; -c is same as -C but
#                    opens the editor
# x, exec <command> = run command (the rest of the line) using shell
# b, break = stop here (continue rebase later with 'git rebase --continue')
# d, drop <commit> = remove commit
# l, label <label> = label current HEAD with a name
# t, reset <label> = reset HEAD to a label
# m, merge [-C <commit> | -c <commit>] <label> [# <oneline>]
#         create a merge commit using the original merge commit's
#         message (or the oneline, if no original merge commit was
#         specified); use -c <commit> to reword the commit message
# u, update-ref <ref> = track a placeholder for the <ref> to be updated
#                       to this position in the new commits. The <ref> is
#                       updated at the end of the rebase
#
# These lines can be re-ordered; they are executed from top to bottom.
#
# If you remove a line here THAT COMMIT WILL BE LOST.
#
# However, if you remove everything, the rebase will be aborted.
#

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 ?

@bnavigator
Copy link
Contributor

Both would be fine.

Note:

ben@skylab:~/src/python-control> git status
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 diff system-norms2
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

@henriks76
Copy link
Author

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.

@bnavigator
Copy link
Contributor

Remind me, what was the warning about, again?

@henriks76
Copy link
Author

warnings.warn("System has a direct feedthrough term!")
warnings.warn("System has pole(s) on the stability boundary!")
warnings.warn("System has a direct feedthrough term!")
warnings.warn("System is unstable!")
warnings.warn("Poles close to, or on, the imaginary axis. Norm value may be uncertain.")

Many of the tests gave these warnings.

@bnavigator
Copy link
Contributor

bnavigator commented Feb 4, 2024

I'd prefer to use a pytest.warns context manager or the recwarn fixture for this, if possible.

@henriks76
Copy link
Author

Thanks @bnavigator, I now use a pytest.warns context manager instead.

@bnavigator
Copy link
Contributor

Continued in #971

@bnavigator bnavigator mentioned this pull request Feb 18, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants