[SciPy-user] Fwd: signal.lti(A,B,C,D) with D=0

[josef.pktd at gmail.com]

2009/2/18 William Purcell <williamhpurcell at gmail.com>:
> I think that scipy.signal is set up to do what I need to do, but I am
> having trouble with ss2tf. Line 149 of ltisys is
>
> type_test = A[:,0] + B[:,0] + C[0,:] + D
>
> but I keep getting the error
>
> Traceback (most recent call last):
>  File "test_lti.py", line 13, in <module>
>    x = scipy.signal.ss2tf(A,B,C,D)
>  File "/usr/lib/python2.5/site-packages/scipy/signal/ltisys.py", line
> 153, in ss2tf
>    type_test = A[:,0] + B[:,0] + C[0,:] + D
> ValueError: shape mismatch: objects cannot be broadcast to a single shape
>
> This is because A is nxn, B is nxi, C is mxn, and D is mxi (I hope I
> got that right). My point is that type_test slices the 'n' dimension
> of each matrix and D doesn't have an 'n' dimension. I think that the '
> + D' needs to be removed from type_test or it needs to be padded with
> n-m elements for the test.
>
> I attached a test that reproduces the error. If I comment out '+ D' in
> ss2tf, it seems to work just fine and return what I want.
>
> One last thing, I think that signal.lti should pass an input kwarg to
> ss2zpk and ss2tf so that you don't have to always look at the first
> input (0 index).  In other words, ss2zpk and ss2tf both have a input
> kwarg to tell which input to use and I think that signal.lti should
> have the same feature.
>
> Let me know your thoughts.
>
> Thanks,
> Bill
>


I ran your test script with +D commented out as you proposed.

x = ss2tf(A,B,C,D) runs without raising an exception, but I didn't
check whether the numbers are correct. However, trying to do the
reverse operation raises different exceptions, see below.

None of the lti functions have any tests in the test suite, so it is
difficult for me to figure out what the expected behavior of these
functions is, and it makes refactoring or rewriting of the functions a
hazardous enterprise.

I'm not an expert on continuous time state space modeling, but I tried
out different versions, and my conclusion was that only single input,
single output work correctly. I attach my trying-out script.

As I see it, it is possible to use different parts of scipy.signal for
multidimensional input or output, but the conversion code that relies
on numpy polynomials cannot handle it. So part of the functionality
could be rewritten to allow for different shapes (as you did with
commenting out D in the test).
Additionally, the filters in scipy signal cannot handle multiple
signals, however the filters in nd.image can be used in an indirect
way to have multi-dimensional filters.

I looked at this more for the usage in Kalman filter, vector arma or
vector ar, and without a multidimensional polynomial class and proper
multi-dimensional filters, it is not possible to use scipy.signal for
this. So, I started to write my own VARMA filter ( for discrete time),
but without convenient conversion between the different representation
as scipy.signal has for the univariate case

For a beginning user and not knowing the jargon of the field,
signal.lti is not very approachable. A set of examples and tests would
be useful to see what the functionality and the limitations is.

If you try out different function in lti for multi-dimensional
signals, it would be useful to have a list of functions that could be
extended to the multi-dimensional case, and of functions for which
this is not possible because of the underlying limitations.

Since I'm not using scipy signal, I stopped looking into it. But,
maybe signal.ltisys needs to be adopted by someone.

Josef

-----------
>>> x = ss2tf(A,B,C,D)
>>> ss=signal.tf2ss(*x)

Warning (from warnings module):
  File "C:\Josef\_progs\building\scipy\scipy-trunk-new-r5551\dist\scipy-0.8.0.dev5551.win32\Programs\Python25\Lib\site-packages\scipy\signal\filter_design.py",
line 221
    "results may be meaningless", BadCoefficients)
BadCoefficients: Badly conditionned filter coefficients (numerator):
the results may be meaningless

Traceback (most recent call last):
  File "<pyshell#4>", line 1, in <module>
    ss=signal.tf2ss(*x)
  File "C:\Josef\_progs\building\scipy\scipy-trunk-new-r5551\dist\scipy-0.8.0.dev5551.win32\Programs\Python25\Lib\site-packages\scipy\signal\ltisys.py",
line 59, in tf2ss
    C = num[:,1:] - num[:,0] * den[1:]
ValueError: shape mismatch: objects cannot be broadcast to a single shape


>>> zpk=signal.tf2zpk(*x)

Warning (from warnings module):
  File "C:\Josef\_progs\building\scipy\scipy-trunk-new-r5551\dist\scipy-0.8.0.dev5551.win32\Programs\Python25\Lib\site-packages\scipy\signal\filter_design.py",
line 221
    "results may be meaningless", BadCoefficients)
BadCoefficients: Badly conditionned filter coefficients (numerator):
the results may be meaningless

Traceback (most recent call last):
  File "<pyshell#5>", line 1, in <module>
    zpk=signal.tf2zpk(*x)
  File "C:\Josef\_progs\building\scipy\scipy-trunk-new-r5551\dist\scipy-0.8.0.dev5551.win32\Programs\Python25\Lib\site-packages\scipy\signal\filter_design.py",
line 161, in tf2zpk
    z = roots(b)
  File "\Programs\Python25\Lib\site-packages\numpy\lib\polynomial.py",
line 133, in roots
    raise ValueError,"Input must be a rank-1 array."
ValueError: Input must be a rank-1 array.
-------------- next part --------------
import numpy as np
from scipy import signal

dvec = np.array([1,2,3,4])
A1 = np.array([-dvec,[1,0,0,0],[0,1,0,0],[0,0,1,0]])
B1 = np.array([[1,0,0,0]]).T  # wrong dimension, this requires D has one column
B1 = np.eye(4)
C1 = np.array([[1,2,3,4]])
D1 = np.zeros((1,4))
print signal.ss2tf(A1,B1,C1,D1)
#same as http://en.wikipedia.org/wiki/State_space_(controls)#Canonical_realizations

signal.ss2tf(*signal.tf2ss(*signal.ss2tf(A1,B1,C1,D1)))
np.testing.assert_almost_equal(signal.ss2tf(*signal.tf2ss(*signal.ss2tf(A1,B1,C1,D1)))[0],signal.ss2tf(A1,B1,C1,D1)[0])

'''

dx_t = A x_t + B u_t
 y_t = C x_t + D u_t


>>> dvec = np.array([1,2,3,4])
>>> A = np.array([-dvec,[1,0,0,0],[0,1,0,0],[0,0,1,0]])
>>> B = np.array([[1,0,0,0]]).T  # wrong dimension, this requires D has one column
>>> B = np.eye(4)
>>> C = np.array([[1,2,3,4]])
>>> D = np.zeros((1,4))
>>> num, den = signal.ss2tf(A,B,C,D)
>>> print num
[[ 0.  1.  2.  3.  4.]]
>>> print den
[ 1.  1.  2.  3.  4.]

>>> A1,B1,C1,D1 = signal.tf2ss(*signal.ss2tf(A,B,C,D))
>>> A1
array([[-1., -2., -3., -4.],
       [ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  0.]])
>>> B1
array([[ 1.],
       [ 0.],
       [ 0.],
       [ 0.]])
>>> C1
array([[ 1.,  2.,  3.,  4.]])
>>> D1
array([ 0.])
'''

# can only have one noise variable u_t
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
dvec = np.array([1,2,3,4])
A = np.array([-dvec,[1,0,0,0],[0,1,0,0],[0,0,1,0]])
B = np.array([[1,0,0,0]]).T  # wrong dimension, this requires D has one column
B = np.eye(4)
B[2,1] = 1
C = np.array([[1,2,3,4]])
D = np.zeros((1,4))
print signal.ss2tf(A,B,C,D)


# can only have one output variable y_t
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
dvec = np.array([1,2,3,4])
A = np.array([-dvec,[1,0,0,0],[0,1,0,0],[0,0,1,0]])
B = np.array([[1,0,0,0]]).T  # wrong dimension, this requires D has one column
B = np.eye(4)
B[2,1] = 1
C = np.array([[1,2,3,4],[1,0,0,0]])
D = np.zeros((2,4))
#print signal.ss2tf(A,B,C,D)
#this causes
##    type_test = A[:,0] + B[:,0] + C[0,:] + D
##ValueError: shape mismatch
#