Implement ERA, change api to TimeResponseData #1024

KybernetikJo · 2024-07-04T13:59:59Z

This PR implements the ERA method. Related to #1002 and #847

it uses the TimeResponseData Type.

coveralls · 2024-07-04T14:12:08Z

coverage: 94.21% (-0.3%) from 94.533%
when pulling 28744f8 on KybernetikJo:implement_era
into 4acc78b on python-control:main.

coveralls · 2024-07-04T17:30:45Z

coverage: 94.21% (-0.3%) from 94.533%
when pulling 596e5ed on KybernetikJo:implement_era
into 4acc78b on python-control:main.

coveralls · 2024-07-07T17:22:30Z

coverage: 94.374% (-0.2%) from 94.533%
when pulling 55ac0ed on KybernetikJo:implement_era
into 4acc78b on python-control:main.

…test, fix small things

coveralls · 2024-07-07T17:45:31Z

coverage: 94.533% (-0.02%) from 94.554%
when pulling 614a080 on KybernetikJo:implement_era
into 7a135d1 on python-control:main.

murrayrm

Overall this looks good. I found a few small things that you might want to update.

murrayrm · 2024-07-09T14:15:08Z

control/modelsimp.py

+        impulse-response data from which the StateSpace model is estimated,
+        1D or 3D array.


impulse → Impulse (start sentence with a capital letter) and try to fit on one line (if less than 75 characters).

murrayrm · 2024-07-09T14:15:20Z

control/modelsimp.py

+        impulse-response data from which the StateSpace model is estimated,
+        1D or 3D array.
+    data : TimeResponseData
+        impulse-response data from which the StateSpace model is estimated.


impose → Impulse

murrayrm · 2024-07-09T14:15:36Z

control/modelsimp.py

+        Number of rows in Hankel matrix.
+        Default is 2*r.


Put on a single line (here and below).

murrayrm · 2024-07-09T14:16:18Z

control/modelsimp.py

+        True indicates discrete time with unspecified sampling time,
+        positive number is discrete time with specified sampling time.
+        It can be used to scale the StateSpace model in order to match
+        the impulse response of this library.
+        Default values is True.


Can probably have fewer lines by going out to column 75 before wrapping. Default value sentence doesn't need to start on a new line.

Tried and Done!

murrayrm · 2024-07-09T14:22:40Z

control/modelsimp.py

+def era(arg, r, m=None, n=None, dt=True, transpose=False):
+    r"""era(YY, r)


Optional: to be consistent with the longer naming style we are using in python-control, I suggest we call the function eigensys_realization and then define era = eigensys_realization at the bottom of the file (see freqplot.py for some example). We should also (but not here) redefine the other functions in this file (minimal_realization, balanced_realization, etc).

Done! See #1028

coveralls · 2024-07-11T08:58:41Z

coverage: 94.608% (-0.02%) from 94.628%
when pulling 8b93264 on KybernetikJo:implement_era
into 6406868 on python-control:main.

KybernetikJo · 2024-07-11T11:01:31Z

I get a warning

../python-control/doc/generated/control.era.rst: WARNING: document isn't included in any toctree

I could not exclude the alias name era from autosummary.

__all__=[eigensys_realization, era]

How to exclude era from autosummary?

SOLUTION: The Solution is, to delete the generated folder in doc. Then the sphinx/autosummary does not use old files.

roryyorke

some minor suggestions, and a question.

roryyorke · 2024-07-14T10:31:23Z

control/modelsimp.py

+    where `response` is an `TimeResponseData` object, and `YY` is 1D or 3D
+    array and r is an integer.


Suggested change

where `response` is an `TimeResponseData` object, and `YY` is 1D or 3D

array and r is an integer.

where `data` is a `TimeResponseData` object, `YY` is a 1D or 3D

array, and r is an integer.

response is more specific than data, perhaps better?

response or data? Not sure.

For me, the function time_response_plot which uses data for time_response, is the blueprint.

The plotting functions use data for input.
bode_plot

The response data types use response for input.
NyquistResponseData

OK, in my PR #1022 I initially used response instead of data. I will go with data for now.

@roryyorke @murrayrm what should we use?

I prefer response as it's more specific, but the main thing is to be consistent within the ERA docstring(s).

roryyorke · 2024-07-14T10:33:51Z

control/modelsimp.py

+        Number of columns in Hankel matrix. Default is 2*r.
+    dt : True or float, optional
+        True indicates discrete time with unspecified sampling time,
+        positive number is discrete time with specified sampling time. It


Suggested change

positive number is discrete time with specified sampling time. It

and a positive float means discrete time with the specified sampling time. It

roryyorke · 2024-07-14T10:35:04Z

control/modelsimp.py

+        True indicates discrete time with unspecified sampling time,
+        positive number is discrete time with specified sampling time. It
+        can be used to scale the StateSpace model in order to match the
+        impulse response of this library. Default is True.


I don't understand the idea of "impulse response of this library" - what does library mean here?

This library "python-control" uses a unit area impulse "delta = 1/dt" for dt > 0 and a unit impulse "delta = 1" for dt = True.

The dt parameter is used to create a discrete time state space system (A,B,C,D,dt) to match an impulse response generated by control.impulse_response.

Done!

I change the text to " It can be used to scale the StateSpace model in order to match the unit-area impulse response of python-control. "

Is there a better way to explain that behavior. Should we add something to Library conventions in order to explain the unit-area impulse response of python-control?

Ah, thanks, now I understand.

One possibly-not-obvious side-effect of using True is that it has a value of 1, so one gets debatably natural discrete-time behaviour in that case.

The unit-area convention is documented in timeresp.impulse_resp under Notes, I suggest a reference to that.

roryyorke · 2024-07-14T10:39:06Z

control/modelsimp.py

    """
-    raise NotImplementedError('This function is not implemented yet.')
+    def block_hankel_matrix(Y, m, n):


suggest making this top-level _block_hankel.

roryyorke · 2024-07-14T10:41:01Z

control/modelsimp.py

+        A reduced order model sys=StateSpace(Ar,Br,Cr,Dr,dt).
+    S : array
+        Singular values of Hankel matrix. Can be used to choose a good r
+        value.


given this, is there an alternative interface with a minsigma argument which could be used to choose r based on the singular values instead? (or maybe tol, with minsigma = tol * maxsigma ?)

This probably a bit out-of-scope, and if the user needs this capability it is simple enough to implement in their own code.

I have thought about it, but I will not be adding this feature at this time.

I don't know the literature on this topic well enough, but I know that there are papers that study the choice of r.

Postponed. Possible improvement in the future.

roryyorke · 2024-07-14T11:05:06Z

control/modelsimp.py

+    # balanced realizations
+    Sigma_inv = np.diag(1./np.sqrt(S[0:r]))
+    Ar = Sigma_inv @ Ur.T @ Hl @ Vhr.T @ Sigma_inv
+    Br = Sigma_inv @ Ur.T @ Hf[:,0:p]*dt


It looks like the reference here is "Algorithm 1" of the arxiv paper.

Why does dt makes an appearance here? I don't see any mention of sample rate or sample period in the referenced paper. It must be needed since you have dt=0.1 in the spring-damper model unit test -- I guess I'm missing something obvious?

ERA in the paper assumes a unit impulse response "delta = 1". There is no dt in the paper.

However, python control generates a unit-area impulse response "delta = 1/dt". See also: #812, #977.

In other words, the impulse responses of StateSpace(A,B,C,D,dt=True) and StateSpace(A,B,C,D,dt=0.1) differ by a scaling.
The dt parameter in era can be used to create the right StateSpace model.

Alternative:
We could also rescale the output data by Y=Y*dt in order to get the unit impulse response instead of the unit-area impulse response.

Overall, my goal is that impulse_response, era, markov and StateSpace should work well together.

Done. I added a little comment to indicate the intentional use of Br*dt.

Is there a more elegant way to achieve this?

roryyorke

LGTM. If you think it's useful, add a reference to timeresp.impulse_response.

KybernetikJo added 5 commits July 7, 2024 19:31

Implement ERA, change api to TimeResponseData

cd87f2f

Add era example

8a9c497

Rename era example file name, small clean up

562824c

Add era example to doc

6bbad5f

Change API to work with TimeResponseData and ndarray as input, add py…

614a080

…test, fix small things

KybernetikJo force-pushed the implement_era branch from 55ac0ed to 614a080 Compare July 7, 2024 17:32

KybernetikJo marked this pull request as ready for review July 7, 2024 17:46

murrayrm requested changes Jul 9, 2024

View reviewed changes

KybernetikJo requested a review from murrayrm July 11, 2024 08:57

KybernetikJo force-pushed the implement_era branch from c5f0cbe to 4ea1d82 Compare July 11, 2024 11:24

Fix few docstring things, change name to eigensys_realization

250c448

KybernetikJo force-pushed the implement_era branch from 4ea1d82 to 250c448 Compare July 11, 2024 14:04

KybernetikJo mentioned this pull request Jul 11, 2024

Change function names in modelsimp, add aliases #1028

Merged

roryyorke reviewed Jul 14, 2024

View reviewed changes

KybernetikJo added 2 commits July 15, 2024 09:26

Change _block_hankel to top-level, update docstrings and comments

6c73062

Delete the hyphen from the Impuse response.

8b93264

KybernetikJo requested a review from roryyorke July 15, 2024 08:07

roryyorke approved these changes Jul 20, 2024

View reviewed changes

murrayrm merged commit 4c7977b into python-control:main Aug 6, 2024
23 checks passed

murrayrm added this to the 0.10.1 milestone Aug 8, 2024

murrayrm linked an issue Nov 23, 2024 that may be closed by this pull request

How to convert finite step response model to state space representation? #1002

Open

		impulse-response data from which the StateSpace model is estimated,
		1D or 3D array.

		def era(arg, r, m=None, n=None, dt=True, transpose=False):
		r"""era(YY, r)

		where `response` is an `TimeResponseData` object, and `YY` is 1D or 3D
		array and r is an integer.

	positive number is discrete time with specified sampling time. It
	and a positive float means discrete time with the specified sampling time. It

Implement ERA, change api to TimeResponseData #1024

Implement ERA, change api to TimeResponseData #1024

Conversation

KybernetikJo commented Jul 4, 2024

coveralls commented Jul 4, 2024

coveralls commented Jul 4, 2024

coveralls commented Jul 7, 2024

coveralls commented Jul 7, 2024

murrayrm left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

coveralls commented Jul 11, 2024 • edited Loading

KybernetikJo commented Jul 11, 2024 • edited Loading

roryyorke left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

KybernetikJo Jul 14, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

KybernetikJo Jul 15, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

KybernetikJo Jul 14, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

roryyorke left a comment

Choose a reason for hiding this comment

coveralls commented Jul 11, 2024 •

edited

Loading

KybernetikJo commented Jul 11, 2024 •

edited

Loading

KybernetikJo Jul 14, 2024 •

edited

Loading

KybernetikJo Jul 15, 2024 •

edited

Loading

KybernetikJo Jul 14, 2024 •

edited

Loading