Skip to content
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

[GSoC] Added methods to StateSpace class to calculate state and output vectors. #26685

Merged
merged 1 commit into from
Jun 19, 2024

Conversation

abhiphile
Copy link
Contributor

Added methods to StateSpace class to calculate state and output vectors by using the ODE module.

  1. State vector can be found by solving the equation $\dot x= Ax + Bu$ by using the ODE module. linodesolve() function type2 has been used to find the solution.
  2. Output vector is calculated by using the state_vector function and by using the equation $y = Cx + Du$.

References to other Issues or PRs

Brief description of what is fixed or changed

Other comments

Release Notes

  • physics.control
    • Added methods to find state and output vectors for StateSpace.

@sympy-bot
Copy link

sympy-bot commented Jun 7, 2024

Hi, I am the SymPy bot. I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

  • physics.control
    • Added methods to find state and output vectors for StateSpace. (#26685 by @abhiphile)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.14.

Click here to see the pull request description that was parsed.
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->
Added methods to `StateSpace` class to calculate state and output vectors by using the ODE module.
1. State vector can be found by solving the equation $\dot x= Ax + Bu$ by using the ODE module. `linodesolve()` function `type2` has been used to find the solution.
2. Output vector is calculated by using the `state_vector` function and by using the equation $y = Cx + Du$.

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->


#### Brief description of what is fixed or changed


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* physics.control
  * Added methods to find state and output vectors for StateSpace.
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

@abhiphile abhiphile changed the title Added methods to StateSpace class to calculate state and output vectors. [GSoC] Added methods to StateSpace class to calculate state and output vectors. Jun 7, 2024
@abhiphile
Copy link
Contributor Author

Doctest

For doctest I've used this resource https://web.mit.edu/2.14/www/Handouts/StateSpaceResponse.pdf.
Example 1 for state_vector
Example 2 for output_vector

Unit Test

For unit test the following resources are used

  1. Test 1 : https://lpsa.swarthmore.edu/Transient/TransMethSS.html
  2. Test 2 : https://web.mit.edu/2.14/www/Handouts/StateSpaceResponse.pdf (Example 11)

Copy link

github-actions bot commented Jun 7, 2024

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)

Significantly changed benchmark results (master vs previous release)

Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

@abhiphile
Copy link
Contributor Author

@Sc0rpi0n101 @faze-geek Please take a look at this PR.

@Sc0rpi0n101 Sc0rpi0n101 merged commit 84b9fc7 into sympy:master Jun 19, 2024
48 checks passed
t = symbols('t')
ss1 = StateSpace(A1, B1, C1, D1)
assert ss1.state_vector(initial_conditions=i1, var=t) == Matrix([[ 4*exp(-t) - 3*exp(-2*t)],
[-4*exp(-t) + 6*exp(-2*t)]])
Copy link
Member

Choose a reason for hiding this comment

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

There are no tests with symbolic matrices. Will this work with such matrices?

Copy link
Contributor

Choose a reason for hiding this comment

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

I don't think the StateSpace model supports symbolic matrices as of now. Kindly raise an issue and we can discuss it there.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It may be possible that some functions in StateSpace model doesn't support symbolic matrices but I've tested that most of the methods in StateSpace support symbolic.

>>> a11, a12, a21, a22 = symbols('a11 a12 a21 a22')
>>> b1, b2 = symbols('b1 b2')
>>> c1, c2 = symbols('c1 c2')
>>> d1 = symbols('d1')
>>> A = Matrix([[a11, a12], [a21, a22]])
>>> B = Matrix([[b1], [b2]])
>>> C = Matrix([[c1, c2]])
>>> D = Matrix([d1])
>>> ss = StateSpace(A, B, C, D)
>>> ss
StateSpace(
Matrix([
[a11, a12],
[a21, a22]]),

Matrix([
[b1],
[b2]]),

Matrix([[c1, c2]]),

Matrix([[d1]]))
>>> ss.rewrite(TransferFunction)
[[TransferFunction(-c1*(a12*b2 - a22*b1 + b1*s) - c2*(-a11*b2 + a21*b1 + b2*s) + d1*(-a11*a22 + a11*s + a12*a21 + a22*s - s**2), -a11*a22 + a11*s + a12*a21 + a22*s - s**2, s)]]

Copy link
Contributor

Choose a reason for hiding this comment

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

Ohh yes, I recall this. Thanks for letting me know.

input_vector = zeros(self._B.shape[1], 1)
sv = self.state_vector(initial_conditions, input_vector, var)
res = self._C*sv + self._D*input_vector
return res.simplify()
Copy link
Member

Choose a reason for hiding this comment

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

In general, simplify should not be called in library code due to large potential performance costs.

@@ -3838,6 +3840,108 @@ def num_outputs(self):
"""
return self._D.rows

def state_vector(self, initial_conditions=None, input_vector=None, var=None):
Copy link
Member

Choose a reason for hiding this comment

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

Given these are all kwargs, what happens if you call .state_vector(). I don't see a test case for that.


initial_conditions : Matrix
The initial conditions of `x` state vector.
input_vector : Matrix
Copy link
Member

Choose a reason for hiding this comment

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

I don't see any test cases for when you pass something like u = Matrix([Heaviside(t)]). Does something like this work? It is super common to use a step function. What about u = Matrix([sin(t)])? Or any other common input function to test?

Copy link
Contributor

Choose a reason for hiding this comment

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

@abhiphile It is your responsibility to make such tests. If you stick to MATLAB examples then only numerical tests will be supported and not symbolic ones. Let's make a separate PR only for such tests (it would cover all the StateSpace methods not just this one).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes @faze-geek I'm working on a new PR which resolves the issues of this PR.

"""

if not var:
var = Symbol('t')
Copy link
Member

Choose a reason for hiding this comment

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

Would this not be better in the function signature? def state_vector(var=Symbol('t')) then the user will know what symbol is used if they don't pass one.

mat1 = mat1.replace(dummy_symbols[ind], v)
return mat1

def output_vector(self, initial_conditions=None, input_vector=None, var=None):
Copy link
Member

@moorepants moorepants Jun 19, 2024

Choose a reason for hiding this comment

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

I'm not a huge fan of the names of these new methods. You are solving the linear ordinary differential equations initial value problem. initial_value_problem_solution or ivp_sol or solve_ivp (which is same as scipy's numerical method) seem more aligned with what the methods actually do.

Copy link
Contributor

Choose a reason for hiding this comment

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

+1

Copy link
Member

Choose a reason for hiding this comment

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

We also already have dsolve() terminology in SymPy, so copying what we already have for these names is also probably a good idea.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Do you mean that the output_vector() should be renamed as dsolve() or ivp_sol() ?

Copy link
Member

Choose a reason for hiding this comment

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

Yes, both state_vector and output_vector coudl have better names. You also do not need a state_vector method because you can always pass C = eye(len(A)) to make y = x. I would only expose the single method to the users and name it something that aligns with the naming already in sympy or the name in scipy.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

6 participants