signal.lti(A,B,C,D) with D=0
I am trying to set up state space representation of a system using signal.lti. The system has no feedforward so D = 0. I've tried the three options listed in the code below to represent D. The only way I can get it to work is option 2 if C has 1 row. If C has more than 1 row it won't work. Any thoughts? -Bill Code --------------------------------------------------------------- from numpy import ones,matrix from scipy import signal r = ones(3) A = matrix([r,r,r]) B = matrix([r]).T C = matrix([r,r]) #three options of D to make it 0 #1) D=0 #2) D = matrix([0,0]).T #3) D = None D = 0 #D = None #D = matrix([0,0]).T Gss = signal.lti(A,B,C,D) ----------------------------------------------------------- Tracebacks ----------------------------------------------------------- Option 1 /usr/lib/python2.5/site-packages/scipy/signal/ltisys.py in abcd_normalize(A, B, C, D) 101 raise ValueError, "A and C must have the same number of columns." 102 if MD != MC: --> 103 raise ValueError, "C and D must have the same number of rows." 104 if ND != NB: 105 raise ValueError, "B and D must have the same number of columns." <type 'exceptions.ValueError'>: C and D must have the same number of rows. WARNING: Failure executing file: <testlti.py> Option 2 (with C as two rows...if C is a single row I do not get this traceback) /usr/lib/python2.5/site-packages/scipy/signal/ltisys.py in ss2tf(A, B, C, D, input) 147 148 num_states = A.shape[0] --> 149 type_test = A[:,0] + B[:,0] + C[0,:] + D 150 num = numpy.zeros((nout, num_states+1), type_test.dtype) 151 for k in range(nout): <type 'exceptions.ValueError'>: shape mismatch: objects cannot be broadcast to a single shape WARNING: Failure executing file: <testlti.py> Option 3 same as 1
This seems like a bug in ltisys. I think this line 149 type_test = A[:,0] + B[:,0] + C[0,:] + D should be 149 type_test = A[:,0] + B[:,0] + C[0,:] + D[0,:] for a multipe-input/multiple-output system with i inputs, n states, and m outputs, A should be n by n, B should be n by i, C should be m by n, and D should be m by i. (actually, because of this code a few lines above: # make MOSI from possibly MOMI system. if B.shape[-1] != 0: B = B[:,input] B.shape = (B.shape[0],1) if D.shape[-1] != 0: D = D[:,input] I guess line 149 should be 149 type_test = A[:,0] + B[:,0] + C[0,:] + D[0] ) But making this change exposes another problem: /usr/lib/python2.5/site-packages/scipy/signal/ltisys.py in __init__(self, *args, **kwords) 224 self.__dict__['D'] = abcd_normalize(*args) 225 self.__dict__['zeros'], self.__dict__['poles'], \ --> 226 self.__dict__['gain'] = ss2zpk(*args) 227 self.__dict__['num'], self.__dict__['den'] = ss2tf(*args) 228 self.inputs = self.B.shape[-1] /usr/lib/python2.5/site-packages/scipy/signal/ltisys.py in ss2zpk(A, B, C, D, input) 183 """ 184 Pdb().set_trace() --> 185 return tf2zpk(*ss2tf(A,B,C,D,input=input)) 186 187 class lti(object): /usr/lib/python2.5/site-packages/scipy/signal/filter_design.py in tf2zpk(b, a) 128 k = b[0] 129 b /= b[0] --> 130 z = roots(b) 131 p = roots(a) 132 return z, p, k /usr/lib/python2.5/site-packages/numpy/lib/polynomial.py in roots(p) 98 p = atleast_1d(p) 99 if len(p.shape) != 1: --> 100 raise ValueError,"Input must be a rank-1 array." 101 102 # find non-zero array entries <type 'exceptions.ValueError'>: Input must be a rank-1 array. WARNING: Failure executing file: <proj2.py> In [2]: %debug
/usr/lib/python2.5/site-packages/numpy/lib/polynomial.py(100)roots() 99 if len(p.shape) != 1: --> 100 raise ValueError,"Input must be a rank-1 array." 101
ipdb> print p [[ 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00] [ 1.81898940e-11 -3.06222400e+00 -2.79196150e+02 -5.87518033e+03 -1.59843721e+04 1.91386789e-07] [ 1.75910240e+01 1.60479206e+03 3.37698682e+04 9.12618857e+04 -1.23219975e+04 -3.29828641e+04]] ss2tf seems to correctly handle MIMO (or at least SIMO systems) correctly and returns one denominator polynomial with several (m) numerator polynomials. But tf2zpk in filter_design.py does not seem able to handle more than siso systems (which makes sense, it is expecting a transfer function which is just a SISO system). How should this be fixed? I understand why ss2tf converts a MIMO system to SIMO - trying to represent a mulitple input, multiple output system with a transfer function has some limitations. I think the real offender is in the __init__ method of signal.lti: elif N == 4: ... self.__dict__['zeros'], self.__dict__['poles'], \ self.__dict__['gain'] = ss2zpk(*args) For a MIMO system, what should the zpk representation be? For a SIMO system, I would expect a vector of poles and a matrix of zeros. This seems to be in line with what ss2tf does. It seems like the init method could find the pole by passing C[0,:] and D[0,:] to ss2zpk. The zeros and gains could then be found by calling ss2zpk in some vectorized way or simply in a for loop for each row of C and D. But there may well be a more elegant solution. Any thoughts? Ryan On Sun, Jul 13, 2008 at 2:59 PM, William Purcell <williamhpurcell@gmail.com> wrote:
I am trying to set up state space representation of a system using signal.lti. The system has no feedforward so D = 0. I've tried the three options listed in the code below to represent D.
The only way I can get it to work is option 2 if C has 1 row. If C has more than 1 row it won't work.
Any thoughts? -Bill
Code --------------------------------------------------------------- from numpy import ones,matrix from scipy import signal r = ones(3) A = matrix([r,r,r]) B = matrix([r]).T C = matrix([r,r])
#three options of D to make it 0 #1) D=0 #2) D = matrix([0,0]).T #3) D = None
D = 0 #D = None #D = matrix([0,0]).T
Gss = signal.lti(A,B,C,D) -----------------------------------------------------------
Tracebacks ----------------------------------------------------------- Option 1 /usr/lib/python2.5/site-packages/scipy/signal/ltisys.py in abcd_normalize(A, B, C, D) 101 raise ValueError, "A and C must have the same number of columns." 102 if MD != MC: --> 103 raise ValueError, "C and D must have the same number of rows." 104 if ND != NB: 105 raise ValueError, "B and D must have the same number of columns."
<type 'exceptions.ValueError'>: C and D must have the same number of rows. WARNING: Failure executing file: <testlti.py>
Option 2 (with C as two rows...if C is a single row I do not get this traceback) /usr/lib/python2.5/site-packages/scipy/signal/ltisys.py in ss2tf(A, B, C, D, input) 147 148 num_states = A.shape[0] --> 149 type_test = A[:,0] + B[:,0] + C[0,:] + D 150 num = numpy.zeros((nout, num_states+1), type_test.dtype) 151 for k in range(nout):
<type 'exceptions.ValueError'>: shape mismatch: objects cannot be broadcast to a single shape WARNING: Failure executing file: <testlti.py>
Option 3 same as 1
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
I am getting back to using signal.lti for state-space representation. This is a post and a reply that I had received a while back, but there was no follow up. Any thoughts on a solution for representing MIMO or SIMO systems with signal.lti, with the issue described below with ss2zpk? -Bill ---------- Forwarded message ---------- From: "Ryan Krauss" <ryanli...@gmail.com> Date: Jul 14 2008, 7:41 am Subject: signal.lti(A,B,C,D) with D=0 To: SciPy-user This seems like a bug in ltisys. I think this line 149 type_test = A[:,0] + B[:,0] + C[0,:] + D should be 149 type_test = A[:,0] + B[:,0] + C[0,:] + D[0,:] for a multipe-input/multiple-output system with i inputs, n states, and m outputs, A should be n by n, B should be n by i, C should be m by n, and D should be m by i. (actually, because of this code a few lines above: # make MOSI from possibly MOMI system. if B.shape[-1] != 0: B = B[:,input] B.shape = (B.shape[0],1) if D.shape[-1] != 0: D = D[:,input] I guess line 149 should be 149 type_test = A[:,0] + B[:,0] + C[0,:] + D[0] ) But making this change exposes another problem: /usr/lib/python2.5/site-packages/scipy/signal/ltisys.py in __init__(self, *args, **kwords) 224 self.__dict__['D'] = abcd_normalize(*args) 225 self.__dict__['zeros'], self.__dict__['poles'], \ --> 226 self.__dict__['gain'] = ss2zpk(*args) 227 self.__dict__['num'], self.__dict__['den'] = ss2tf (*args) 228 self.inputs = self.B.shape[-1] /usr/lib/python2.5/site-packages/scipy/signal/ltisys.py in ss2zpk(A, B, C, D, input) 183 """ 184 Pdb().set_trace() --> 185 return tf2zpk(*ss2tf(A,B,C,D,input=input)) 186 187 class lti(object): /usr/lib/python2.5/site-packages/scipy/signal/filter_design.py in tf2zpk(b, a) 128 k = b[0] 129 b /= b[0] --> 130 z = roots(b) 131 p = roots(a) 132 return z, p, k /usr/lib/python2.5/site-packages/numpy/lib/polynomial.py in roots(p) 98 p = atleast_1d(p) 99 if len(p.shape) != 1: --> 100 raise ValueError,"Input must be a rank-1 array." 101 102 # find non-zero array entries <type 'exceptions.ValueError'>: Input must be a rank-1 array. WARNING: Failure executing file: <proj2.py> In [2]: %debug> /usr/lib/python2.5/site-packages/numpy/lib/ polynomial.py(100)roots() 99 if len(p.shape) != 1: --> 100 raise ValueError,"Input must be a rank-1 array." 101 ipdb> print p [[ 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00] [ 1.81898940e-11 -3.06222400e+00 -2.79196150e+02 -5.87518033e+03 -1.59843721e+04 1.91386789e-07] [ 1.75910240e+01 1.60479206e+03 3.37698682e+04 9.12618857e+04 -1.23219975e+04 -3.29828641e+04]] ss2tf seems to correctly handle MIMO (or at least SIMO systems) correctly and returns one denominator polynomial with several (m) numerator polynomials. But tf2zpk in filter_design.py does not seem able to handle more than siso systems (which makes sense, it is expecting a transfer function which is just a SISO system). How should this be fixed? I understand why ss2tf converts a MIMO system to SIMO - trying to represent a mulitple input, multiple output system with a transfer function has some limitations. I think the real offender is in the __init__ method of signal.lti: elif N == 4: ... self.__dict__['zeros'], self.__dict__['poles'], \ self.__dict__['gain'] = ss2zpk (*args) For a MIMO system, what should the zpk representation be? For a SIMO system, I would expect a vector of poles and a matrix of zeros. This seems to be in line with what ss2tf does. It seems like the init method could find the pole by passing C[0,:] and D[0,:] to ss2zpk. The zeros and gains could then be found by calling ss2zpk in some vectorized way or simply in a for loop for each row of C and D. But there may well be a more elegant solution. Any thoughts? Ryan On Sun, Jul 13, 2008 at 2:59 PM, William Purcell <williamhpurc...@gmail.com> wrote:
I am trying to set up state space representation of a system using signal.lti. The system has no feedforward so D = 0. I've tried the three options listed in the code below to represent D.
The only way I can get it to work is option 2 if C has 1 row. If C has more than 1 row it won't work.
Any thoughts? -Bill
Code --------------------------------------------------------------- from numpy import ones,matrix from scipy import signal r = ones(3) A = matrix([r,r,r]) B = matrix([r]).T C = matrix([r,r])
#three options of D to make it 0 #1) D=0 #2) D = matrix([0,0]).T #3) D = None
D = 0 #D = None #D = matrix([0,0]).T
Gss = signal.lti(A,B,C,D) -----------------------------------------------------------
Tracebacks ----------------------------------------------------------- Option 1 /usr/lib/python2.5/site-packages/scipy/signal/ltisys.py in abcd_normalize(A, B, C, D) 101 raise ValueError, "A and C must have the same number of columns." 102 if MD != MC: --> 103 raise ValueError, "C and D must have the same number of rows." 104 if ND != NB: 105 raise ValueError, "B and D must have the same number of columns."
<type 'exceptions.ValueError'>: C and D must have the same number of rows. WARNING: Failure executing file: <testlti.py>
Option 2 (with C as two rows...if C is a single row I do not get this traceback) /usr/lib/python2.5/site-packages/scipy/signal/ltisys.py in ss2tf(A, B, C, D, input) 147 148 num_states = A.shape[0] --> 149 type_test = A[:,0] + B[:,0] + C[0,:] + D 150 num = numpy.zeros((nout, num_states+1), type_test.dtype) 151 for k in range(nout):
<type 'exceptions.ValueError'>: shape mismatch: objects cannot be broadcast to a single shape WARNING: Failure executing file: <testlti.py>
Option 3 same as 1
_______________________________________________ SciPy-user mailing list SciPy-u...@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
_______________________________________________ SciPy-user mailing list SciPy-u...@scipy.orghttp://projects.scipy.org/mailman/listinfo/scipy- user
On Tue, Feb 17, 2009 at 3:25 PM, flyeng4 <williamhpurcell@gmail.com> wrote:
I am getting back to using signal.lti for state-space representation. This is a post and a reply that I had received a while back, but there was no follow up. Any thoughts on a solution for representing MIMO or SIMO systems with signal.lti, with the issue described below with ss2zpk? -Bill
I was looking recently into scipy.signal to see whether it can be used for time series analysis. After digging around in the code and trying out some examples (I didn't find any examples in the docs), I got the impression that scipy.signal is not designed to handle more than one signal (I wanted recursive filter with interactions between several time series). I think the main problem I found was, that numpy/scipy doesn't have a multivariate/multidimensional polynomial class. Somewhere when converting from a multidimensional transfer function, scipy signal uses numpy.polynomial which only works for one-dimensional polynomials. That's when I gave up. My conclusion was, that for multidimensional signals, scipy signal would need to be rewritten quite a bit. Josef
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
2009/2/18 William Purcell <williamhpurcell@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.
another thought: Obtaining multi-dimensional output can be done by stitching together several lti systems or transfer functions by looping over the output dimension (I did something similar for multi-dimensional filters with ndimage). However, multi-dimensional inputs cannot be handled this way. I didn't see any way to merge 2 independent input signals to one output signal. Josef
I was thinking the same thing. I am working on a thin wrapper over ltisys to do tests on the dimensions of the output and then looping if there is any stitching to be done. I was also thinking that I could loop over the inputs first and then in a sub-routine, loop over the outputs to come up with a two-dimensional list of what ever I need. For example, if I am passing a MIMO system to signal.lti, each of the representation conversions would be in matrix/list form corresponding to each input/output relationship. Or lti could test if its a MIMO or even SIMO and make a matrix/list of lti instances which would each have there own alternative representations through ss2zpk/ss2tf etc (which might be the cleaner of the two alternatives.) Do you think this would be a feature for ltisys (I don't think it would take much time to implement), or do you think that it is hackish and ltisys should stick to SISO systems? On Wed, Feb 18, 2009 at 12:26 PM, <josef.pktd@gmail.com> wrote:
another thought:
Obtaining multi-dimensional output can be done by stitching together several lti systems or transfer functions by looping over the output dimension (I did something similar for multi-dimensional filters with ndimage).
However, multi-dimensional inputs cannot be handled this way. I didn't see any way to merge 2 independent input signals to one output signal.
Josef _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
On Wed, Feb 18, 2009 at 2:00 PM, William Purcell <williamhpurcell@gmail.com> wrote:
I was thinking the same thing. I am working on a thin wrapper over ltisys to do tests on the dimensions of the output and then looping if there is any stitching to be done. I was also thinking that I could loop over the inputs first and then in a sub-routine, loop over the outputs to come up with a two-dimensional list of what ever I need. For example, if I am passing a MIMO system to signal.lti, each of the representation conversions would be in matrix/list form corresponding to each input/output relationship. Or lti could test if its a MIMO or even SIMO and make a matrix/list of lti instances which would each have there own alternative representations through ss2zpk/ss2tf etc (which might be the cleaner of the two alternatives.) Do you think this would be a feature for ltisys (I don't think it would take much time to implement), or do you think that it is hackish and ltisys should stick to SISO systems?
Are you looking at parallel SISO and SIMO systems, or did you find a way to have one output depend on two inputs? Personally, I don't like to change code without having tests, to verify that what I am doing doesn't break things and delivers what I want. Since, ltisys doesn't have tests, I wouldn't want to do much surgery on it, and I would prefer a separate wrapper, or subclassing or delegation with minimal changes to the existing code. Another idea for true MIMO system would be to extract what works for that case and write a separate MIMO package. I think using only the state space representation, it shouldn't be too difficult to get simulation and similar things working correctly. The definition and representation of transfer functions and zpk would be more difficult, and I don't know much about it in the multi-dimensional input case (My intuition is more in terms of time series analysis and I haven't looked at this very closely.) Later if everything works and is tested, merging of MIMO and SISO could be considered. These are my thoughts as a (half) innocent bystander. Josef
On Wed, Feb 18, 2009 at 12:26 PM, <josef.pktd@gmail.com> wrote:
another thought:
Obtaining multi-dimensional output can be done by stitching together several lti systems or transfer functions by looping over the output dimension (I did something similar for multi-dimensional filters with ndimage).
However, multi-dimensional inputs cannot be handled this way. I didn't see any way to merge 2 independent input signals to one output signal.
Josef _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
participants (4)
-
flyeng4
-
josef.pktd@gmail.com
-
Ryan Krauss
-
William Purcell