GeezerVison is showing github.com — Observability Gramian for discrete-time systems · Issue #967 · python-control/python-control
Observability Gramian for discrete-time systems #967DescriptionI noticed in the source code for the gram function that it says only continuous-time systems are supported for now:
|
# TODO: Check for continuous or discrete, only continuous supported for now |
|
# if isCont(): |
|
# dico = 'C' |
|
# elif isDisc(): |
|
# dico = 'D' |
|
# else: |
I checked and it doesn't raise an exception when you pass a discrete time system but I don't think this means it is calculating a gramian for a discrete-time system:
A = [[-0.31, 0.21], [-0.68, -0.57]]
B = [[1.23], [1.42]]
C = [[ 1.32, -0.55]]
D = [[0.]]
sysd = ct.ss(A, B, C, D, dt=1)
assert sysd.isdtime()
Wo = ct.gram(sysd, 'o')
print(Wo)
array([[ 3.24632551, -0.19876604],
[-0.19876604, 0.19212128]])
What would it take to get it working for discrete time systems? Is it simply a different call to the same slycot module, sb03md, or is it more involved?
Since I need it I'd be happy to try implementing it. Activitybilltubbs commented on Feb 8, 2024 ContributorAuthor More actions Expected output should be this I think:
array([[ 2.06193690, -0.57482614],
[-0.57482614, 0.78661650]])
billtubbs commented on Feb 10, 2024 ContributorAuthor More actions Where is the documentation on slycot? I tried googling and found the API but couldn't find anything specific on sb03md.
From wikipedia, the two equations that need to be solved are quite similar:
Observability Gramian for a continuous-time system:
$$A^T W_o + W_o A = -C^T C$$
Observability Gramian for a discrete-time system:
$$A^T W_{o} A - W_{o} = -C^T C$$
So I guess I just need to figure out how to do the second one.
Any advice/pointers would be appreciated. bnavigator commented on Feb 11, 2024 There is no hosted doc of Slycot (yet) (python-control/Slycot#202).
- You can look at the documentation of SLICOT's SB03MD
- There are the docstrings, which are used by IDEs like spyder to render documentation on-the-fly.
You are probably looking for the dico='D' parameter of slycot.sb03md57.
billtubbs commented on Feb 11, 2024 ContributorAuthor More actions Okay well that was easy. If I simply allow dico='D' option it calculates the correct Gramians.
The code is all there, just commented-out:
# TODO: Check for continuous or discrete, only continuous supported for now
# if isCont():
# dico = 'C'
# elif isDisc():
# dico = 'D'
# else:I simply replaced above with:
if sys.isctime():
dico = 'C'
elif sys.isdtime():
dico = 'D'
else:
raise ValueError("sys")Is there a reason this wasn't implemented? Maybe at the time of writing isCont and isDisc didn't exist and the author was planning to come back and finish it but didn't perhaps? billtubbs commented on Feb 11, 2024 ContributorAuthor More actions Unless anyone knows of a problem with this, I'll just go ahead and implement it, including updated unit tests. billtubbs commented on Feb 11, 2024 ContributorAuthor More actions I noticed this check for stability in the continuous-time case:
# Check for continuous or discrete
if sys.isctime():
dico = 'C'
# TODO: Check system is stable, perhaps a utility in ctrlutil.py
# or a method of the StateSpace class?
if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
raise ValueError("Oops, the system is unstable!")Is this the right way to do this for discrete-time systems?
elif sys.isdtime():
dico = 'D'
if np.any(np.abs(sys.poles()) >= 1.):
raise ValueError("Oops, the system is unstable!")Correct me if I'm wrong but we don't have an sys.isstable method yet. billtubbs commented on Feb 11, 2024 ContributorAuthor More actions I've submitted the pull request.
One other thing I noticed, the gram function has an argument type which is a reserved Python keyword.
Should we change this to something else?
Or we could replace gram with obsv_gram and ctrl_gram. murrayrm commented on Mar 30, 2024
MetadataMetadataDevelopmentNo branches or pull requestsIssue actions
|
