Multiply scipy.lti transfer functions

Question

I have the following scipy.lti object that is basically an object representing a Laplace transform of an LTI system:

G_s = lti([1], [1, 2])

How to multiply such a transfer function with another one, such as i.e.:

H_s = lti([2], [1, 2])

#I_s = G_s * H_s <---- How to multiply this properly?

I guess I could do

I_s = lti(np.polymul([1], [2]), np.polymul([1, 2], [1, 2]))

But what if I want to do:

#I_s = H_s / (1 + H_s) <---- Does not work since H_s is an lti object

Is there an easy way to do this with scipy?

Answer 1

Depending on your definition of "easy", you should consider deriving your own class from lti, implementing the necessary algebraic operations on your transfer functions. This is probably the most elegant approach.

Here's my take on the subject:

from __future__ import division

from scipy.signal.ltisys import TransferFunction as TransFun
from numpy import polymul,polyadd

class ltimul(TransFun):
    def __neg__(self):
        return ltimul(-self.num,self.den)

    def __floordiv__(self,other):
        # can't make sense of integer division right now
        return NotImplemented

    def __mul__(self,other):
        if type(other) in [int, float]:
            return ltimul(self.num*other,self.den)
        elif type(other) in [TransFun, ltimul]:
            numer = polymul(self.num,other.num)
            denom = polymul(self.den,other.den)
            return ltimul(numer,denom)

    def __truediv__(self,other):
        if type(other) in [int, float]:
            return ltimul(self.num,self.den*other)
        if type(other) in [TransFun, ltimul]:
            numer = polymul(self.num,other.den)
            denom = polymul(self.den,other.num)
            return ltimul(numer,denom)

    def __rtruediv__(self,other):
        if type(other) in [int, float]:
            return ltimul(other*self.den,self.num)
        if type(other) in [TransFun, ltimul]:
            numer = polymul(self.den,other.num)
            denom = polymul(self.num,other.den)
            return ltimul(numer,denom)

    def __add__(self,other):
        if type(other) in [int, float]:
            return ltimul(polyadd(self.num,self.den*other),self.den)
        if type(other) in [TransFun, type(self)]:
            numer = polyadd(polymul(self.num,other.den),polymul(self.den,other.num))
            denom = polymul(self.den,other.den)
            return ltimul(numer,denom)

    def __sub__(self,other):
        if type(other) in [int, float]:
            return ltimul(polyadd(self.num,-self.den*other),self.den)
        if type(other) in [TransFun, type(self)]:
            numer = polyadd(polymul(self.num,other.den),-polymul(self.den,other.num))
            denom = polymul(self.den,other.den)
            return ltimul(numer,denom)

    def __rsub__(self,other):
        if type(other) in [int, float]:
            return ltimul(polyadd(-self.num,self.den*other),self.den)
        if type(other) in [TransFun, type(self)]:
            numer = polyadd(polymul(other.num,self.den),-polymul(other.den,self.num))
            denom = polymul(self.den,other.den)
            return ltimul(numer,denom)

    # sheer laziness: symmetric behaviour for commutative operators
    __rmul__ = __mul__
    __radd__ = __add__

This defines the ltimul class, which is lti plus addition, multiplication, division, subtraction, and negation; binary ones also defined for integers and floats as partners.

I tested it for the example of Dietrich:

G_s = ltimul([1], [1, 2])
H_s = ltimul([2], [1, 0, 3])
print(G_s*H_s)
print(G_s*H_s/(1+G_s*H_s))

While GH is nicely equal to

ltimul(
array([ 2.]),
array([ 1.,  2.,  3.,  6.])
)

the final result for GH/(1+GH) is less pretty:

ltimul(
array([  2.,   4.,   6.,  12.]),
array([  1.,   4.,  10.,  26.,  37.,  42.,  48.])
)

Since I'm not very familiar with transfer functions, I'm not sure how likely it is that this gives the same result as the sympy-based solution due to some simplifications missing from this one. I find it suspicious that already lti behaves unexpectedly: lti([1,2],[1,2]) doesn't simplify its arguments, even though I'd suspect this function to be constant 1. So I'd rather not guess the correctness of this final result.

Anyway, the main message is inheritance itself, so possible bugs in the above implementation hopefully pose only a minor inconvenience. I'm also quite unfamiliar with class definitions, so it's possible that I didn't follow best practices in the above.

---

I eventually rewrote the above after @ochurlaud pointed out, that my original only worked for Python 2. The reason is that the / operation is implemented by __div__/__rdiv__ in Python 2 (and is the ambiguous "classical division"). In Python 3, however, there is a distinction between / (true division) and // (floor division), and they call __truediv__ and __floordiv__ (and their "right" counterparts), respectively. The __future__ import first in the line of the above code triggers the proper Python 3 behaviour even on Python 2, so the above works on both Python versions. Since floor (integer) division doesn't make much sense for our class, we explicitly signal that it can't do anything with // (unless the other operand implements it).

One could also easily define the respective __iadd__, __idiv__ etc. in-place operations for +=, /= etc., respectively.