乘以 scipy.lti 传递函数

Question

我有以下 scipy.lti 对象，它基本上是表示 LTI 系统的拉普拉斯变换的对象：

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

如何将这样的传递函数与另一个传递函数相乘，例如：

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

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

我想我可以做

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

但是如果我想做怎么办：

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

有没有一种简单的方法可以用 scipy 来做到这一点？

Answer 1

有趣的是，Scipy 似乎不提供该功能。另一种方法是将 LTI 系统转换为 Sympy 有理函数。 Sympy 允许您轻松展开和取消多项式：

from IPython.display import display
from scipy import signal
import sympy as sy

sy.init_printing()  # LaTeX like pretty printing for IPython


def lti_to_sympy(lsys, symplify=True):
    """ Convert Scipy's LTI instance to Sympy expression """
    s = sy.Symbol('s')
    G = sy.Poly(lsys.num, s) / sy.Poly(lsys.den, s)
    return sy.simplify(G) if symplify else G


def sympy_to_lti(xpr, s=sy.Symbol('s')):
    """ Convert Sympy transfer function polynomial to Scipy LTI """
    num, den = sy.simplify(xpr).as_numer_denom()  # expressions
    p_num_den = sy.poly(num, s), sy.poly(den, s)  # polynomials
    c_num_den = [p.all_coeffs() for p in p_num_den]  # coefficients
    l_num, l_den = [sy.lambdify((), c)() for c in c_num_den]  # convert to floats
    return signal.lti(l_num, l_den)


pG, pH, pGH, pIGH = sy.symbols("G, H, GH, IGH")  # only needed for displaying


# Sample systems:
lti_G = signal.lti([1], [1, 2])
lti_H = signal.lti([2], [1, 0, 3])

# convert to Sympy:
Gs, Hs = lti_to_sympy(lti_G), lti_to_sympy(lti_H)


print("Converted LTI expressions:")
display(sy.Eq(pG, Gs))
display(sy.Eq(pH, Hs))

print("Multiplying Systems:")
GHs = sy.simplify(Gs*Hs).expand()  # make sure polynomials are canceled and expanded
display(sy.Eq(pGH, GHs))


print("Closing the loop:")
IGHs = sy.simplify(GHs / (1+GHs)).expand()
display(sy.Eq(pIGH, IGHs))

print("Back to LTI:")
lti_IGH = sympy_to_lti(IGHs)
print(lti_IGH)

输出为：

Answer 2

根据您对“简单”的定义，您应该考虑从

lti

派生您自己的类，在您的传递函数上实现必要的代数运算。这可能是最优雅的方法。

这是我对这个问题的看法：

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__

这定义了

ltimul

类，即

lti

加上加法、乘法、除法、减法和求反；二进制也定义为整数和浮点数作为伙伴。

我测试了它以迪特里希为例：

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

虽然

GH

很好地等于

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

GH/(1+GH) 的最终结果不太漂亮：

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

由于我对传递函数不是很熟悉，因此我不确定这是否有可能给出与基于 sympy 的解决方案相同的结果，因为该解决方案缺少一些简化。我发现可疑的是

lti

的行为已经出乎意料：

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

没有简化它的参数，尽管我怀疑这个函数是常数 1。所以我宁愿不猜测这个最终结果的正确性。

无论如何，主要信息是继承本身，因此上述实现中可能存在的错误希望只会带来轻微的不便。我对类定义也很不熟悉，所以我可能没有遵循上面的最佳实践。

---

在@ochurlaud指出后，我最终重写了上面的内容，我的原始版本仅适用于Python 2。原因是

/

操作是由Python 2中的

__div__

/

__rdiv__

实现的（并且是不明确的） “古典除法”）。然而，在 Python 3 中，

/

（真正的除法）和

//

（地板除法）之间存在区别，它们分别称为

__truediv__

和

__floordiv__

（以及它们的“右”对应项）。上面代码行中的

__future__

导入即使在 Python 2 上也会触发正确的 Python 3 行为，因此上述代码适用于两个 Python 版本。由于下限（整数）除法对我们的类没有多大意义，因此我们明确表示它不能对

//

执行任何操作（除非其他操作数实现它）。

还可以轻松地分别为

__iadd__

、

__idiv__

等定义相应的

+=

、

/=

等就地操作。