将 scipy.lti 传递函数相乘 [英] Multiply scipy.lti transfer functions
问题描述
我有以下 scipy.lti 对象，它基本上是一个表示 LTI 系统拉普拉斯变换的对象:
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_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
有没有一种简单的方法可以用 scipy 做到这一点?
Is there an easy way to do this with scipy?
推荐答案
根据你对简单"的定义，你应该考虑从 lti 派生你自己的类，在你的代码上实现必要的代数运算传递函数.这可能是最优雅的方法.
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__
这里定义了ltimul类，就是lti加上加法、乘法、除法、减法、求反；二进制的也定义为整数和浮点数作为伙伴.
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.
我测试了它以 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.])
)
GH/(1+GH) 的最终结果不太漂亮:
the final result for GH/(1+GH) is less pretty:
ltimul(
array([ 2., 4., 6., 12.]),
array([ 1., 4., 10., 26., 37., 42., 48.])
)
由于我对传递函数不是很熟悉，因此我不确定这与基于 sympy 的解决方案得到相同结果的可能性有多大，因为该解决方案缺少一些简化.我发现 lti 的行为已经出乎意料了: lti([1,2],[1,2]) 并没有简化它的参数，即使我会怀疑这个函数是常数 1.所以我宁愿不去猜测这个最终结果的正确性.
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.
在@ochurlaud指出后，我最终重写了上述内容，我原来只适用于 Python 2.原因是 / 操作是由 Python 2 中的 __div__/__rdiv__ 实现的(并且是模棱两可的经典划分").然而，在 Python 3 中，/(真除法)和 //(底除法)是有区别的，它们调用 __truediv__ 和__floordiv__(以及它们的正确"对应物).上面代码行中的 __future__ import first 会触发正确的 Python 3 行为，即使在 Python 2 上也是如此，因此上述内容适用于两个 Python 版本.由于地板(整数)除法对我们的类没有多大意义，我们明确表示它不能用 // 做任何事情(除非另一个操作数实现它).
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).
还可以轻松地为 +=、/=<定义相应的 __iadd__、__idiv__ 等就地操作/code> 等，分别.
One could also easily define the respective __iadd__, __idiv__ etc. in-place operations for +=, /= etc., respectively.
这篇关于将 scipy.lti 传递函数相乘的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
