added bilinear() to tf in control lti

[zolabar]

bilinear discretization of continuous transfer function in Control API

References to other Issues or PRs

Fixes #24544

Brief description of what is fixed or changed

Added the function bilinear() to TransferFunction objects in sympy\physics\control\lti.py

It returns falling coeffs of H(z) from numerator and denominator. H(z) is the corresponding discretized tf, discretized with bilinear transform method. Coefficients are falling, i.e. H(z) = (az+b)/(cz+d) is returned as [a, b], [c, d].

Other comments

It is a useful function for numerical implementations of discrete signal processing and control with scipy.signal, control or matlab/simulink, once the symbolic continuous transfer function design is finished.

Release Notes

• functions
  • Added bilinear, a function for the bilinear transform discretization of a continuous transfer function from Control API.

[sympy-bot]

✅

Hi, I am the SymPy bot (v169). 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:

• functions
  • Added bilinear, a function for the bilinear transform discretization of a continuous transfer function from Control API. (#24558 by @zolabar)

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

Click here to see the pull request description that was parsed.

# bilinear discretization of continuous transfer function in Control API



#### References to other Issues or PRs

Fixes #24544


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

Added the function bilinear() to TransferFunction objects in sympy\physics\control\lti.py

It returns falling coeffs of H(z) from numerator and denominator. H(z) is the corresponding discretized tf, discretized with bilinear transform method. Coefficients are falling, i.e. H(z) = (az+b)/(cz+d) is returned as [a, b], [c, d].

#### Other comments

It is a useful function for numerical implementations of discrete signal processing and control with scipy.signal, control or matlab/simulink, once the symbolic continuous transfer function design is finished.

#### Release Notes

<!-- BEGIN RELEASE NOTES -->

* functions
  * Added `bilinear`, a function for the bilinear transform discretization of a continuous transfer function from Control API.

<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[sympy-bot]

🟠

Hi, I am the SymPy bot (v169). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at sympy/sympy-bot#75.

The following commits add new files:

• fe659ef:
  • examples/notebooks/transfer_function_bilinear.ipynb
  • transfer_function_bilinear.ipynb

The following commits delete files:

• 727fd6b:
  • transfer_function_bilinear.ipynb
• 36f3651:
  • examples/notebooks/transfer_function_bilinear.ipynb

If these files were added/deleted on purpose, you can ignore this message.

[oscarbenjamin]

The following commits add new files:

• fe659ef:

  • examples/notebooks/transfer_function_bilinear.ipynb
  • transfer_function_bilinear.ipynb

The following commits delete files:

• 727fd6b:

  • transfer_function_bilinear.ipynb

These commits will need to be squashed. We can't have a commit that introduces a file followed by another that deletes it. The file must not exist in the history at all.

[faze-geek]

I see lot of blank lines in between the code which are not required.

[faze-geek]

As far as the failing test is concerned you should run python bin/strip_whitespace sympy/physics/control/lti.py. Then your code quality check will pass.

[1e9abhi1e10]

you need to add your name and email in the .mailmap file and then run python bin/mailmap_check.py for instructions read here.

[zolabar]

you need to add your name and email in the .mailmap file and then run python bin/mailmap_check.py for instructions read here.

thanks, just did that

[zolabar]

Check the .mailmap file

…

Am 21.01.2023 um 08:51 schrieb ABHISHEK PATIDAR ***@***.***>:  you need to add your name and email in the .mailmap file and then run python bin/mailmap_check.py for instructions read here. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.

[1e9abhi1e10]

Check the .mailmap file
…
Am 21.01.2023 um 08:51 schrieb ABHISHEK PATIDAR @.***>:  you need to add your name and email in the .mailmap file and then run python bin/mailmap_check.py for instructions read here. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.

After adding your name and email in the .mailmap file you need to run python bin/mailmap_check.py and then push a new commit

[github-actions]

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)

       before           after         ratio
     [41d90958]       [71817106]
     <sympy-1.11.1^0>                 
-         965±2μs          606±3μs     0.63  solve.TimeSparseSystem.time_linear_eq_to_matrix(10)
-        2.79±0ms      1.14±0.01ms     0.41  solve.TimeSparseSystem.time_linear_eq_to_matrix(20)
-     5.56±0.01ms         1.67±0ms     0.30  solve.TimeSparseSystem.time_linear_eq_to_matrix(30)

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

[zolabar]

As far as the failing test is concerned you should run python bin/strip_whitespace sympy/physics/control/lti.py. Then your code quality check will pass.

As far as the failing test is concerned you should run python bin/strip_whitespace sympy/physics/control/lti.py. Then your code quality check will pass.

Thanks, now all checks have passed :)

[faze-geek]

Thanks, now all checks have passed :)

Could you add relevant examples in form of tests in https://github.com/sympy/sympy/blob/master/sympy/physics/control/tests/test_lti.py

[zolabar]

Thanks, now all checks have passed :)

Could you add relevant examples in form of tests in https://github.com/sympy/sympy/blob/master/sympy/physics/control/tests/test_lti.py

Ok, i‘ll prepare some tests

[zolabar]

Could you add relevant examples in form of tests in https://github.com/sympy/sympy/blob/master/sympy/physics/control/tests/test_lti.py

@faze-geek, done e595caf

[hanspi42]

I think two things should be changed:

1. The method should take z and T as input arguments.
2. This can be done without creating a symbolic expression internally, as @faze-geek says: if N is the maximum of the degrees of numerator and denominator, then multiplying numerator and denominator by (z+1)**N after substitution gives polynomials in z. E.g., if the numerator polynnomial is a*s**2+b*s+c, and N=3 because we have a third-order denominator, the new numerator polynomial is a*(2/T)**2*(z-1)**2*(z+1) + b*(2/T)*(z-1)*(z+1)**2 + c*(z+1)**3. This looks simple to code.

[zolabar]

I think two things should be changed:

1. The method should take z and T as input arguments.
2. This can be done without creating a symbolic expression internally, as @faze-geek says: if N is the maximum of the degrees of numerator and denominator, then multiplying numerator and denominator by (z+1)**N after substitution gives polynomials in z. E.g., if the numerator polynnomial is a*s**2+b*s+c, and N=3 because we have a third-order denominator, the new numerator polynomial is a*(2/T)**2*(z-1)**2*(z+1) + b*(2/T)*(z-1)*(z+1)**2 + c*(z+1)**3. This looks simple to code.

@hanspi42, thank you for the feedback.

Do I understand correctly, that concerning point 1, you propose that the function should be called as (for tf as TransferFunction, discrete_var the symbolic discrete variable (e.g. z) and sample_per the symbolic sample period) :

tf.blinear(discrete_var, sample_per)
?

Concerning point 2, I don't see how to take advantage of your analysis, that is, how should I get the list of the numerator and denominator coeffcients of the discretized transfer function?

Taking into account your point 1. it could look as follows. However, I as I wrote above, I have no idea of how to avoid the defining H = self.num/self.den and then using SymPy to get the list of the numerator and denominator coeffcients of the discretized transfer function (which are needed for discrete implementation in practice then)...

`def blinear(self, discrete_var, sample_per):

    z = discrete_var
    T = sample_per
  
    H = self.num/self.den
    HZ = H.subs(self.var, (2/T)*(z-1)/(z+1))

    num, den = fraction(HZ.simplify())

    HZnum = collect(expand(num, z), z)
    HZnum = Poly(HZnum, z)

    HZden = collect(expand(den, z), z)
    HZden = Poly(HZden, z)

    num_coefs = HZnum.coeffs()
    den_coefs = HZden.coeffs()

    return num_coefs, den_coefs

`

[hanspi42]

Concerning point 2, I don't see how to take advantage of your analysis, that is, how should I get the list of the numerator and denominator coeffcients of the discretized transfer function?

For example:

from sympy import *
s, z = symbols('s z')
T, a, b, c, d = symbols('T a b c d', real=True)

n = a*s**2+b*s+c
d = b*s**3+c*s**2+d*s+a
H = n/d
display(H)

np = n.as_poly(s).all_coeffs()
dp = d.as_poly(s).all_coeffs()

# This is the suggested code starting from having numerator and denominator
# coefficients np and dp
N = max(len(np), len(dp)) - 1
n1 = Add(*[ T**(N-i)*2**i*c*(z-1)**i*(z+1)**(N-i) for c, i in zip(np[::-1], range(len(np))) ])
d1 = Add(*[ T**(N-i)*2**i*c*(z-1)**i*(z+1)**(N-i) for c, i in zip(dp[::-1], range(len(dp))) ])

# This shows the transfer function
H1 = n1/d1
display(H1)

# This is the solution you are using
H2 = H.subs(s, 2/T*(z-1)/(z+1))
display(H2)

# This shows that it is the same; must be zero.
dH = (H1-H2).simplify()
display(dH)

[hanspi42]

tf.blinear(discrete_var, sample_per)

Yes ;)

[zolabar]

Concerning point 2, I don't see how to take advantage of your analysis, that is, how should I get the list of the numerator and denominator coeffcients of the discretized transfer function?

For example:

from sympy import *
s, z = symbols('s z')
T, a, b, c, d = symbols('T a b c d', real=True)

n = a*s**2+b*s+c
d = b*s**3+c*s**2+d*s+a
H = n/d
display(H)

np = n.as_poly(s).all_coeffs()
dp = d.as_poly(s).all_coeffs()

# This is the suggested code starting from having numerator and denominator
# coefficients np and dp
N = max(len(np), len(dp)) - 1
n1 = Add(*[ T**(N-i)*2**i*c*(z-1)**i*(z+1)**(N-i) for c, i in zip(np[::-1], range(len(np))) ])
d1 = Add(*[ T**(N-i)*2**i*c*(z-1)**i*(z+1)**(N-i) for c, i in zip(dp[::-1], range(len(dp))) ])

# This shows the transfer function
H1 = n1/d1
display(H1)

# This is the solution you are using
H2 = H.subs(s, 2/T*(z-1)/(z+1))
display(H2)

# This shows that it is the same; must be zero.
dH = (H1-H2).simplify()
display(dH)

@hanspi42 , thank you for the fast reply and the elegant compact code example, i‘ll test it and think about it!

[zolabar]

tf.blinear(discrete_var, sample_per)

Yes ;)

@hanspi42 , commited the changes, thank you very much for the algebraic hints and the code prototype.
Since the discrete variable does not appear in the output coeff list I only included the sample period as input, hence the current function is used as

nz, dz = tf.bilinear(T)

[sylee957]

I think that general interface is good.
But I'm not sure it should be added under methods because transfer function basic classes only supports arithmetic now

I think that we can pull out as function outside if we are not clear.

[sylee957]

I think that adding additional features are okay if it is outside, but only want to avoid adding inside objects as much as possible because of SRP

[zolabar]

I think that general interface is good. But I'm not sure it should be added under methods because transfer function basic classes only supports arithmetic now

I think that we can pull out as function outside if we are not clear.

@sylee957 , thank you for the feedback.

For me it would be OK to inlcude bilinear as an outside function, like in

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.bilinear.html

[zolabar]

@sylee957,

in commit

36d6d4a

bilinear is a function, used as

bilinear(tf, T)

similar to its numerical counterpart from sicipy.signal

[hanspi42]

It is OK for me now, but I would suggest one final change: show the formula of the bilinear-transform substitution in the doc string.

[zolabar]

It is OK for me now, but I would suggest one final change: show the formula of the bilinear-transform substitution in the doc string.

-> done

[hanspi42]

It is OK for me now, but I would suggest one final change: show the formula of the bilinear-transform substitution in the doc string.

-> done

It looks OK, but does not appear in the documentation (preview here: https://output.circle-artifacts.com/output/job/1eb985ff-0719-4e5a-a222-ce27d1c7ee08/artifacts/0/doc/_build/html/index.html)

Unfortunately I do not know how to do this; does any of the other reviewers know?

[sylee957]

You need to add a directive .. autofunction: bilinear to doc/src/modules/physics/control/lti.rst

[zolabar]

You need to add a directive .. autofunction: bilinear to doc/src/modules/physics/control/lti.rst

-> done

[zolabar]

You need to add a directive .. autofunction: bilinear to doc/src/modules/physics/control/lti.rst

-> done

[zolabar]

It is OK for me now, but I would suggest one final change: show the formula of the bilinear-transform substitution in the doc string.

@hanspi42 , i included the docstring that @sylee957 suggested. How did you generate the test preview of the documentation? Thanks

[hanspi42]

How did you generate the test preview of the documentation?

This is done by github. Look under "All checks have passed", there is an entry "Click here to see a preview of the documentation"

[hanspi42]

@sylee957: This now looks OK to me.

[sylee957]

I think that docs should be more readable because summary is too long and there are not many examples to figure out how it works for general cases.

[sylee957]

Also, I find that the release notes are misleading to functions, I'd correct them manually