NotificationsYou must be signed in to change notification settings
Fork441
Star1.8k

Implement ERA, change api to TimeResponseData#1024

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 ourterms of service andprivacy statement. We’ll occasionally send you account related emails.

Already on GitHub?Sign in to your account

Jump to bottom

Merged

murrayrm merged 8 commits intopython-control:mainfromKybernetikJo:implement_era

Aug 6, 2024

Merged

Implement ERA, change api to TimeResponseData#1024

murrayrm merged 8 commits intopython-control:mainfromKybernetikJo:implement_era

Aug 6, 2024

Conversation

Copy link

Contributor

KybernetikJo commentedJul 4, 2024

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

it uses the TimeResponseData Type.

Copy link

coveralls commentedJul 4, 2024

coverage: 94.21% (-0.3%) from 94.533%
when pulling28744f8 on KybernetikJo:implement_era
into4acc78b on python-control:main.

Copy link

coveralls commentedJul 4, 2024

coverage: 94.21% (-0.3%) from 94.533%
when pulling596e5ed on KybernetikJo:implement_era
into4acc78b on python-control:main.

Copy link

coveralls commentedJul 7, 2024

coverage: 94.374% (-0.2%) from 94.533%
when pulling55ac0ed on KybernetikJo:implement_era
into4acc78b on python-control:main.

KybernetikJo added5 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 theimplement_era branch from55ac0ed to614a080Compare

July 7, 2024 17:32

Copy link

coveralls commentedJul 7, 2024

coverage: 94.533% (-0.02%) from 94.554%
when pulling614a080 on KybernetikJo:implement_era
into7a135d1 on python-control:main.

KybernetikJo marked this pull request as ready for review

July 7, 2024 17:46

murrayrm requested changes

Jul 9, 2024

View reviewed changes

Copy link

Member

murrayrm left a comment

Choose a reason for hiding this comment

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

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

control/modelsimp.py Outdated

Comment on lines 397 to 398

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

Copy link

Member

Choose a reason for hiding this comment

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

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

Copy link

ContributorAuthor

Choose a reason for hiding this comment

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

Done.

control/modelsimp.py Outdated

		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.

Copy link

Member

Choose a reason for hiding this comment

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

impose → Impulse

Copy link

ContributorAuthor

Choose a reason for hiding this comment

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

Done.

control/modelsimp.py Outdated

Comment on lines 404 to 405

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

Copy link

Member

Choose a reason for hiding this comment

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

Put on a single line (here and below).

Copy link

ContributorAuthor

Choose a reason for hiding this comment

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

Done!

control/modelsimp.py Outdated

Comment on lines 410 to 414

		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.

Copy link

Member

Choose a reason for hiding this comment

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

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.

Copy link

ContributorAuthor

Choose a reason for hiding this comment

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

Tried and Done!

control/modelsimp.py Outdated

Comment on lines 372 to 373

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

Copy link

Member

Choose a reason for hiding this comment

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

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

Copy link

ContributorAuthor

KybernetikJoJul 15, 2024

Choose a reason for hiding this comment

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

Done! See#1028

KybernetikJo requested a review frommurrayrm

July 11, 2024 08:57

Copy link

coveralls commentedJul 11, 2024•
edited
Loading

coverage: 94.608% (-0.02%) from 94.628%
when pulling8b93264 on KybernetikJo:implement_era
into6406868 on python-control:main.

Copy link

ContributorAuthor

KybernetikJo commentedJul 11, 2024•
edited
Loading

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.

KybernetikJo force-pushed theimplement_era branch fromc5f0cbe to4ea1d82Compare

July 11, 2024 11:24

Fix few docstring things, change name to eigensys_realization

250c448

KybernetikJo force-pushed theimplement_era branch from4ea1d82 to250c448Compare

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

Copy link

Contributor

roryyorke left a comment

Choose a reason for hiding this comment

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

some minor suggestions, and a question.

control/modelsimp.py Outdated

Comment on lines 392 to 393

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

Copy link

Contributor

roryyorkeJul 14, 2024

Choose a reason for hiding this comment

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

Suggested change

	where`response`isan`TimeResponseData`object,and`YY`is1Dor3D
	arrayandrisaninteger.
	where`data`isa`TimeResponseData`object,`YY`isa1Dor3D
	array,andrisaninteger.

response is more specific than data, perhaps better?

Copy link

ContributorAuthor

KybernetikJoJul 14, 2024•
edited
Loading

Choose a reason for hiding this comment

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

response ordata? Not sure.

For me, the functiontime_response_plot which usesdata for time_response, is the blueprint.

The plotting functions usedata for input.
bode_plot

The response data types useresponse for input.
NyquistResponseData

Copy link

ContributorAuthor

KybernetikJoJul 15, 2024

Choose a reason for hiding this comment

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

Done.

Copy link

ContributorAuthor

KybernetikJoJul 15, 2024•
edited
Loading

Choose a reason for hiding this comment

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

OK, in my PR#1022 I initially usedresponse instead ofdata. I will go withdata for now.

@roryyorke @murrayrm what should we use?

Copy link

Contributor

roryyorkeJul 20, 2024

Choose a reason for hiding this comment

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

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

control/modelsimp.py Outdated

		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

Copy link

Contributor

roryyorkeJul 14, 2024

Choose a reason for hiding this comment

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

Suggested change

	positivenumberisdiscretetimewithspecifiedsamplingtime.It
	andapositivefloatmeansdiscretetimewiththespecifiedsamplingtime.It

Copy link

ContributorAuthor

KybernetikJoJul 15, 2024

Choose a reason for hiding this comment

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

Done.

control/modelsimp.py Outdated

		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.

Copy link

Contributor

roryyorkeJul 14, 2024

Choose a reason for hiding this comment

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

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

Copy link

ContributorAuthor

KybernetikJoJul 14, 2024

Choose a reason for hiding this comment

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

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

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

Copy link

ContributorAuthor

KybernetikJoJul 15, 2024

Choose a reason for hiding this comment

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

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?

Copy link

Contributor

roryyorkeJul 20, 2024

Choose a reason for hiding this comment

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

Ah, thanks, now I understand.

One possibly-not-obvious side-effect of usingTrue 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 intimeresp.impulse_resp underNotes, I suggest a reference to that.

control/modelsimp.py Outdated

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

Copy link

Contributor

roryyorkeJul 14, 2024

Choose a reason for hiding this comment

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

suggest making this top-level_block_hankel.

Copy link

ContributorAuthor

KybernetikJoJul 15, 2024

Choose a reason for hiding this comment

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

Done!

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.

Copy link

Contributor

roryyorkeJul 14, 2024

Choose a reason for hiding this comment

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

given this, is there an alternative interface with aminsigma 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.

Copy link

ContributorAuthor

KybernetikJoJul 14, 2024•
edited
Loading

Choose a reason for hiding this comment

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

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.

Copy link

ContributorAuthor

KybernetikJoJul 15, 2024

Choose a reason for hiding this comment

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

Postponed. Possible improvement in the future.

control/modelsimp.py Outdated

		# 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

Copy link

Contributor

roryyorkeJul 14, 2024

Choose a reason for hiding this comment

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

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

Why doesdt 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?

Copy link

ContributorAuthor

KybernetikJoJul 14, 2024

Choose a reason for hiding this comment

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

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.

Copy link

ContributorAuthor

KybernetikJoJul 15, 2024

Choose a reason for hiding this comment

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

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

Is there a more elegant way to achieve this?

KybernetikJo added2 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 fromroryyorke

July 15, 2024 08:07

roryyorke approved these changes

Jul 20, 2024

View reviewed changes

Copy link

Contributor

roryyorke left a comment

Choose a reason for hiding this comment

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

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

murrayrm merged commit4c7977b intopython-control:main

Aug 6, 2024

23 checks passed

murrayrm added this to the0.10.1 milestone

Aug 8, 2024

murrayrm linked an issue

Nov 23, 2024

that may beclosed by this pull request

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

Open

Labels

None yet

4 participants

Movatterモバイル変換

Implement ERA, change api to TimeResponseData#1024

Implement ERA, change api to TimeResponseData#1024

Uh oh!

Conversation

KybernetikJo commentedJul 4, 2024

Uh oh!

coveralls commentedJul 4, 2024

Uh oh!

coveralls commentedJul 4, 2024

Uh oh!

coveralls commentedJul 7, 2024

Uh oh!

coveralls commentedJul 7, 2024

Uh oh!

murrayrm left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

coveralls commentedJul 11, 2024• editedLoading Uh oh!There was an error while loading.Please reload this page.

Uh oh!

Uh oh!

KybernetikJo commentedJul 11, 2024• editedLoading Uh oh!There was an error while loading.Please reload this page.

Uh oh!

Uh oh!

roryyorke left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KybernetikJoJul 14, 2024• editedLoading Uh oh!There was an error while loading.Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KybernetikJoJul 15, 2024• editedLoading Uh oh!There was an error while loading.Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KybernetikJoJul 14, 2024• editedLoading Uh oh!There was an error while loading.Please reload this page.

coveralls commentedJul 11, 2024•
edited
Loading

KybernetikJo commentedJul 11, 2024•
edited
Loading

KybernetikJoJul 14, 2024•
edited
Loading

KybernetikJoJul 15, 2024•
edited
Loading

KybernetikJoJul 14, 2024•
edited
Loading