我有一个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
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 = [sy.expand(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)
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加上加法、乘法、除法、减法和取反运算；对于整数和浮点数也定义了二元运算。我已经针对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))
当GH恰好等于
ltimul(
array([ 2.]),
array([ 1., 2., 3., 6.])
)
ltimul(
array([ 2., 4., 6., 12.]),
array([ 1., 4., 10., 26., 37., 42., 48.])
)
由于我对转移函数不是很熟悉，所以我不确定这是否会得到与基于sympy的解决方案相同的结果，因为此方法可能缺少一些简化。我觉得已经有点奇怪了，即使我想这个函数是常量1，但lti([1,2],[1,2])并没有简化它的参数。因此，我宁愿不猜测最终结果的正确性。
无论如何，主要信息是继承本身，因此上述实现中可能存在的错误希望只会造成小的不便。我对类定义也不太熟悉，因此上述实现可能并未遵循最佳实践。
在@ochurlaud指出我的原始代码仅适用于Python 2后，我最终进行了重写。原因是在Python 2中，/操作由__div__/__rdiv__实现（是模糊的"classical division"）。然而，在Python 3中，存在/（真正的除法）和//（地板除法）之间的区别，它们分别调用__truediv__和__floordiv__（以及它们各自的“右”对应项）。上述代码行中第一个__future__导入触发了适当的Python 3行为，即使在Python 2上也是如此，因此上述工作适用于两个Python版本。由于地板（整数）除法对我们的类没有多大意义，因此我们明确表示它不能处理//（除非另一个操作数实现了它）。
还可以轻松定义相应的__iadd__、__idiv__等就地操作，例如+=、/=等。
__div__运算符无效。改为__truediv__解决了问题（对于__rdiv__也是如此）。结合两个答案可以得到一个非常好的子类...我希望它已经在scipy中了... - ochurlaud