Skip to content

[GSoC] (DO NOT MERGE) Adding Discrete-Time TransferFunction - #28115

Closed
leomanga wants to merge 33 commits into
sympy:masterfrom
leomanga:disc-time-systems
Closed

leomanga wants to merge 33 commits into
sympy:masterfrom
leomanga:disc-time-systems

Conversation

@leomanga

@leomanga leomanga commented Jun 4, 2025 •

Copy link
Copy Markdown
Contributor

References to other Issues or PRs

Brief description of what is fixed or changed

This PR lays the groundwork for extending the lti module to discrete time. It implements the TransferFunctionBase class, which is inherited by both TransferFunction and DiscreteTransferFunction. It also extends the ability to create interconnections between discrete-time systems.

Other comments

Release Notes

  • physics.control
    • Added a new class for discrete time transfer functions.
    • Added a general method for creating both discrete-time and continuous-time transfer functions

add pretty printing for discrete-time transfer functions

remove unnecessary newline in at the end of the file
@sympy-bot

sympy-bot commented Jun 4, 2025 •

Copy link
Copy Markdown

✅

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 a new class for discrete time transfer functions. (#28115 by @leomanga)

    • Added a general method for creating both discrete-time and continuous-time transfer functions (#28115 by @leomanga)

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

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

#### 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
This PR lays the groundwork for extending the lti module to discrete time. It implements the `TransferFunctionBase` class, which is inherited by both `TransferFunction` and `DiscreteTransferFunction`. It also extends the ability to create interconnections between discrete-time systems.

#### 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 a new class for discrete time transfer functions.
  * Added a general method for creating both discrete-time and continuous-time transfer functions
<!-- END RELEASE NOTES -->

Comment thread sympy/physics/control/tests/test_lti.py Outdated

def test_DTTransferFunction_construction():
z = symbols("z")
tf = DTTransferFunction(z + 1, z**2 + z + 1, z, 10)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

What do you plan for having fully symbolic transfer functions?

I've seen some recent activity about a pade approximation elsewhere in SymPy, are you going to use that?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

What do you mean?
The Padé approximation can be useful for approximating time delays. Can it be used for other purposes as well?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

I haven't taken digital control in a long long time, so I don't remember everything, but I thought that when transforming a continous transfer function into a discrete transfer function that you have a choice in the method of discretization. I thought using a pade approximation was one of them, but now I think I mis remembered.

So maybe my question should be: how are you going to manage different choices in discretization? Or is this assuming a zero order hold always?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Honestly I didn't know much about Padé approximation for discretization, but I did some research and I found helpful resources.
Anyway, I'm going to group these functions

def gbt(tf, sample_per, alpha):
r"""
Returns falling coefficients of H(z) from numerator and denominator.
Explanation
===========
Where H(z) is the corresponding discretized transfer function,
discretized with the generalised bilinear transformation method.
H(z) is obtained from the continuous transfer function H(s)
by substituting $s(z) = \frac{z-1}{T(\alpha z + (1-\alpha))}$ into H(s), where T is the
sample period.
Coefficients are falling, i.e. $H(z) = \frac{az+b}{cz+d}$ is returned
as [a, b], [c, d].
Examples
========
>>> from sympy.physics.control.lti import TransferFunction, gbt
>>> from sympy.abc import s, L, R, T
>>> tf = TransferFunction(1, s*L + R, s)
>>> numZ, denZ = gbt(tf, T, 0.5)
>>> numZ
[T/(2*(L + R*T/2)), T/(2*(L + R*T/2))]
>>> denZ
[1, (-L + R*T/2)/(L + R*T/2)]
>>> numZ, denZ = gbt(tf, T, 0)
>>> numZ
[T/L]
>>> denZ
[1, (-L + R*T)/L]
>>> numZ, denZ = gbt(tf, T, 1)
>>> numZ
[T/(L + R*T), 0]
>>> denZ
[1, -L/(L + R*T)]
>>> numZ, denZ = gbt(tf, T, 0.3)
>>> numZ
[3*T/(10*(L + 3*R*T/10)), 7*T/(10*(L + 3*R*T/10))]
>>> denZ
[1, (-L + 7*R*T/10)/(L + 3*R*T/10)]
References
==========
.. [1] https://www.polyu.edu.hk/ama/profile/gfzhang/Research/ZCC09_IJC.pdf
"""
if not tf.is_SISO:
raise NotImplementedError("Not implemented for MIMO systems.")
T = sample_per # and sample period T
s = tf.var
z = s # dummy discrete variable z
np = tf.num.as_poly(s).all_coeffs()
dp = tf.den.as_poly(s).all_coeffs()
alpha = Rational(alpha).limit_denominator(1000)
# The next line results from multiplying H(z) with z^N/z^N
N = max(len(np), len(dp)) - 1
num = Add(*[ T**(N-i) * c * (z-1)**i * (alpha * z + 1 - alpha)**(N-i) for c, i in zip(np[::-1], range(len(np))) ])
den = Add(*[ T**(N-i) * c * (z-1)**i * (alpha * z + 1 - alpha)**(N-i) for c, i in zip(dp[::-1], range(len(dp))) ])
num_coefs = num.as_poly(z).all_coeffs()
den_coefs = den.as_poly(z).all_coeffs()
para = den_coefs[0]
num_coefs = [coef/para for coef in num_coefs]
den_coefs = [coef/para for coef in den_coefs]
return num_coefs, den_coefs
def bilinear(tf, sample_per):
r"""
Returns falling coefficients of H(z) from numerator and denominator.
Explanation
===========
Where H(z) is the corresponding discretized transfer function,
discretized with the bilinear transform method.
H(z) is obtained from the continuous transfer function H(s)
by substituting $s(z) = \frac{2}{T}\frac{z-1}{z+1}$ into H(s), where T is the
sample period.
Coefficients are falling, i.e. $H(z) = \frac{az+b}{cz+d}$ is returned
as [a, b], [c, d].
Examples
========
>>> from sympy.physics.control.lti import TransferFunction, bilinear
>>> from sympy.abc import s, L, R, T
>>> tf = TransferFunction(1, s*L + R, s)
>>> numZ, denZ = bilinear(tf, T)
>>> numZ
[T/(2*(L + R*T/2)), T/(2*(L + R*T/2))]
>>> denZ
[1, (-L + R*T/2)/(L + R*T/2)]
"""
return gbt(tf, sample_per, S.Half)
def forward_diff(tf, sample_per):
r"""
Returns falling coefficients of H(z) from numerator and denominator.
Explanation
===========
Where H(z) is the corresponding discretized transfer function,
discretized with the forward difference transform method.
H(z) is obtained from the continuous transfer function H(s)
by substituting $s(z) = \frac{z-1}{T}$ into H(s), where T is the
sample period.
Coefficients are falling, i.e. $H(z) = \frac{az+b}{cz+d}$ is returned
as [a, b], [c, d].
Examples
========
>>> from sympy.physics.control.lti import TransferFunction, forward_diff
>>> from sympy.abc import s, L, R, T
>>> tf = TransferFunction(1, s*L + R, s)
>>> numZ, denZ = forward_diff(tf, T)
>>> numZ
[T/L]
>>> denZ
[1, (-L + R*T)/L]
"""
return gbt(tf, sample_per, S.Zero)
def backward_diff(tf, sample_per):
r"""
Returns falling coefficients of H(z) from numerator and denominator.
Explanation
===========
Where H(z) is the corresponding discretized transfer function,
discretized with the backward difference transform method.
H(z) is obtained from the continuous transfer function H(s)
by substituting $s(z) = \frac{z-1}{Tz}$ into H(s), where T is the
sample period.
Coefficients are falling, i.e. $H(z) = \frac{az+b}{cz+d}$ is returned
as [a, b], [c, d].
Examples
========
>>> from sympy.physics.control.lti import TransferFunction, backward_diff
>>> from sympy.abc import s, L, R, T
>>> tf = TransferFunction(1, s*L + R, s)
>>> numZ, denZ = backward_diff(tf, T)
>>> numZ
[T/(L + R*T), 0]
>>> denZ
[1, -L/(L + R*T)]
"""
return gbt(tf, sample_per, S.One)

into DTTransferFunction as classmethods and I'm going to add as many discretization methods as I can.

@leomanga
leomanga marked this pull request as ready for review June 13, 2025 17:00
leomanga added 2 commits June 15, 2025 23:20
…alization in Series divisions and returned to previous eval_frequency calculation
@moorepants

Copy link
Copy Markdown
Member

Please fill out all the top information:

References to other Issues or PRs
Brief description of what is fixed or changed
Other comments
Release Notes

NO ENTRY

Comment thread sympy/physics/control/tests/test_lti.py
raises(NotImplementedError,
lambda: TransferFunctionBase.from_coeff_lists([1], [1], s))
raises(NotImplementedError,
lambda: TransferFunctionBase.from_zpk([1], [1], 1, s))

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

If we aren't implementing them, I would simply not have them as methods.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

I've done this because these functions must be inherited by a subclass, and we need to be aware that they exist.

Comment thread sympy/physics/control/tests/test_lti.py
tf6_ = DTTransferFunction(1/2, 4, s)
assert tf6_.num == 0.5
assert tf6_.den == 4
assert tf6_.args == (0.500000000000000, 4, s, 1)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

checking floating point with == is not good practice, you should check with a tolerance, i think there is a sympy function for this

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Are you referring to this?

import math
Float(0.1, 9).is_same(Float(0.1, 10), math.isclose)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Yes, that is fine.

Comment thread sympy/physics/control/tests/test_lti.py
expect4_ = DTTransferFunction(7*p**3 + s, p - s, p)
assert SP4.subs({a0: -1, a1: -7}) == expect4_
assert SP4.subs({a0: -1, a1: -7}).evalf() == expect4
assert expect4_.evalf() == expect4

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

All of the above numerical substitution checks are ok to have, but you need to check symbolics as the primary test.

Comment thread sympy/physics/control/tests/test_lti.py Outdated
5*a2*b1*s**(S(8)/3)/3 + S(25)*b1*s/9) - \
5*3**(S(1)/3)*5**(S(2)/3)*\
b0/(9*a1*a2*s**(S(8)/3) + 15*a1*s + \
15*a2*b1*s**(S(8)/3) + 25*b1*s)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Do these expected results come from text book answers? How do you know they are correct results?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Best to avoid \ for line continuation. Just use parentheses around long lines.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Do these expected results come from text book answers? How do you know they are correct results?

These tests are the same as the ones from the continuous-time case, so I assume they hold. I just added the one at line 722, which comes from here https://it.mathworks.com/help/control/ref/dynamicsystem.evalfr.html

Comment thread sympy/physics/control/tests/test_lti.py Outdated
# Low-frequency (or DC) gain.
assert DTTransferFunction(z, z - 1, z).dc_gain() == oo
assert DTTransferFunction(z - 1, z**2 + 3*z + 1, z).dc_gain() == 0
assert DTTransferFunction(z, z + 2, z, 0.2).dc_gain() == Rational(1, 3)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Would be nice to know if these are textbook results.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

I've applied the definition of dc_gain manually. I can find some references if needed.

Comment thread sympy/physics/control/tests/test_lti.py
dtf_manual = DTTransferFunction.\
from_coeff_lists(num_manual, den_manual, z, T).expand()

assert S.Zero == (tf_test_bilinear.simplify()-tf_test_manual.simplify()).simplify().num

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

We should not change tests unless this was somehow giving an incorrect result.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

I just thought it was the wrong way to verify correctness. My changes check the same thing.

property is used for transfer functions in the ``Parallel`` class, and the
multiplicative property is used for transfer functions in the ``Series`` class.
represent Continuous-time Transfer functions in the Laplace domain, while
``DTTransferFunction`` represent Discrete-time Transfer functions in the

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

I think we should name this "DiscreteTranferFunction", i.e. use full words instead of single letters. I'm open for other options, but it is customary in Python to be descriptive with names.

Comment thread sympy/physics/control/lti.py
Comment thread sympy/physics/control/lti.py
Comment thread sympy/physics/control/lti.py
Comment thread sympy/physics/control/lti.py Outdated
@leomanga

leomanga commented Jul 4, 2025

Copy link
Copy Markdown
Contributor Author

@moorepants I added the example in the tutorial and implemented the is_stable method for the discrete transfer function.
Once I implement get_asymptotic_stability_conditions, the stability check will be updated to use the new algorithm, but it already works correctly for now.
Let me know if there’s anything else I should do.

@leomanga

Copy link
Copy Markdown
Contributor Author

@moorepants I made the new_tf() function public because I think it could be useful.
This PR should now be ready to merge.

I noticed some pre-existing style issues, but I don't think it's worth fixing them here since this PR is already quite large.
I’ve started working on the discrete counterpart of the state-space system and fixed some of those style problems there.

Comment thread doc/src/tutorials/physics/control/control_problems.rst Outdated
Comment thread doc/src/tutorials/physics/control/control_problems.rst Outdated
Comment thread doc/src/tutorials/physics/control/control_problems.rst Outdated
Comment thread doc/src/tutorials/physics/control/control_problems.rst Outdated
Comment thread doc/src/tutorials/physics/control/control_problems.rst Outdated
Comment thread doc/src/tutorials/physics/control/control_problems.rst Outdated
Comment thread sympy/physics/control/lti.py
Comment thread sympy/physics/control/lti.py Outdated

"""
m = Mul(self.num, Pow(self.den, -1, evaluate=False), evaluate=False)
return limit(m, self.var, 0)
pass

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Seems odd to have pass here.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

It's because it's an abstract method

@leomanga leomanga changed the title [GSoC] Adding Discrete-Time TransferFunction [GSoC] (DO NOT MERGE) Adding Discrete-Time TransferFunction Jul 14, 2025
- Renamed `new_tf` to `new_transfer_function` for clarity
- Adjusted LaTeX and pretty printing outputs to use "sampling time" instead of "sampling_time".
- Adjusted tutorial formatting and changed the image type to svg
@leomanga

Copy link
Copy Markdown
Contributor Author

I made the proposed changes, when we are ready I will open a new PR clone, so we won't have problems with the png images I uploaded.

Example 6
---------

.. image:: Control_Problems_Q6.svg

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

This image doesn't show correctly on a white background:

image

@moorepants

Copy link
Copy Markdown
Member

Closing in favor of #28318

@moorepants moorepants closed this Aug 20, 2025
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.

4 participants