Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/signal/lti_conversion.py : 11%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2ltisys -- a collection of functions to convert linear time invariant systems
3from one representation to another.
4"""
5import numpy
6import numpy as np
7from numpy import (r_, eye, atleast_2d, poly, dot,
8 asarray, prod, zeros, array, outer)
9from scipy import linalg
11from .filter_design import tf2zpk, zpk2tf, normalize
14__all__ = ['tf2ss', 'abcd_normalize', 'ss2tf', 'zpk2ss', 'ss2zpk',
15 'cont2discrete']
18def tf2ss(num, den):
19 r"""Transfer function to state-space representation.
21 Parameters
22 ----------
23 num, den : array_like
24 Sequences representing the coefficients of the numerator and
25 denominator polynomials, in order of descending degree. The
26 denominator needs to be at least as long as the numerator.
28 Returns
29 -------
30 A, B, C, D : ndarray
31 State space representation of the system, in controller canonical
32 form.
34 Examples
35 --------
36 Convert the transfer function:
38 .. math:: H(s) = \frac{s^2 + 3s + 3}{s^2 + 2s + 1}
40 >>> num = [1, 3, 3]
41 >>> den = [1, 2, 1]
43 to the state-space representation:
45 .. math::
47 \dot{\textbf{x}}(t) =
48 \begin{bmatrix} -2 & -1 \\ 1 & 0 \end{bmatrix} \textbf{x}(t) +
49 \begin{bmatrix} 1 \\ 0 \end{bmatrix} \textbf{u}(t) \\
51 \textbf{y}(t) = \begin{bmatrix} 1 & 2 \end{bmatrix} \textbf{x}(t) +
52 \begin{bmatrix} 1 \end{bmatrix} \textbf{u}(t)
54 >>> from scipy.signal import tf2ss
55 >>> A, B, C, D = tf2ss(num, den)
56 >>> A
57 array([[-2., -1.],
58 [ 1., 0.]])
59 >>> B
60 array([[ 1.],
61 [ 0.]])
62 >>> C
63 array([[ 1., 2.]])
64 >>> D
65 array([[ 1.]])
66 """
67 # Controller canonical state-space representation.
68 # if M+1 = len(num) and K+1 = len(den) then we must have M <= K
69 # states are found by asserting that X(s) = U(s) / D(s)
70 # then Y(s) = N(s) * X(s)
71 #
72 # A, B, C, and D follow quite naturally.
73 #
74 num, den = normalize(num, den) # Strips zeros, checks arrays
75 nn = len(num.shape)
76 if nn == 1:
77 num = asarray([num], num.dtype)
78 M = num.shape[1]
79 K = len(den)
80 if M > K:
81 msg = "Improper transfer function. `num` is longer than `den`."
82 raise ValueError(msg)
83 if M == 0 or K == 0: # Null system
84 return (array([], float), array([], float), array([], float),
85 array([], float))
87 # pad numerator to have same number of columns has denominator
88 num = r_['-1', zeros((num.shape[0], K - M), num.dtype), num]
90 if num.shape[-1] > 0:
91 D = atleast_2d(num[:, 0])
93 else:
94 # We don't assign it an empty array because this system
95 # is not 'null'. It just doesn't have a non-zero D
96 # matrix. Thus, it should have a non-zero shape so that
97 # it can be operated on by functions like 'ss2tf'
98 D = array([[0]], float)
100 if K == 1:
101 D = D.reshape(num.shape)
103 return (zeros((1, 1)), zeros((1, D.shape[1])),
104 zeros((D.shape[0], 1)), D)
106 frow = -array([den[1:]])
107 A = r_[frow, eye(K - 2, K - 1)]
108 B = eye(K - 1, 1)
109 C = num[:, 1:] - outer(num[:, 0], den[1:])
110 D = D.reshape((C.shape[0], B.shape[1]))
112 return A, B, C, D
115def _none_to_empty_2d(arg):
116 if arg is None:
117 return zeros((0, 0))
118 else:
119 return arg
122def _atleast_2d_or_none(arg):
123 if arg is not None:
124 return atleast_2d(arg)
127def _shape_or_none(M):
128 if M is not None:
129 return M.shape
130 else:
131 return (None,) * 2
134def _choice_not_none(*args):
135 for arg in args:
136 if arg is not None:
137 return arg
140def _restore(M, shape):
141 if M.shape == (0, 0):
142 return zeros(shape)
143 else:
144 if M.shape != shape:
145 raise ValueError("The input arrays have incompatible shapes.")
146 return M
149def abcd_normalize(A=None, B=None, C=None, D=None):
150 """Check state-space matrices and ensure they are 2-D.
152 If enough information on the system is provided, that is, enough
153 properly-shaped arrays are passed to the function, the missing ones
154 are built from this information, ensuring the correct number of
155 rows and columns. Otherwise a ValueError is raised.
157 Parameters
158 ----------
159 A, B, C, D : array_like, optional
160 State-space matrices. All of them are None (missing) by default.
161 See `ss2tf` for format.
163 Returns
164 -------
165 A, B, C, D : array
166 Properly shaped state-space matrices.
168 Raises
169 ------
170 ValueError
171 If not enough information on the system was provided.
173 """
174 A, B, C, D = map(_atleast_2d_or_none, (A, B, C, D))
176 MA, NA = _shape_or_none(A)
177 MB, NB = _shape_or_none(B)
178 MC, NC = _shape_or_none(C)
179 MD, ND = _shape_or_none(D)
181 p = _choice_not_none(MA, MB, NC)
182 q = _choice_not_none(NB, ND)
183 r = _choice_not_none(MC, MD)
184 if p is None or q is None or r is None:
185 raise ValueError("Not enough information on the system.")
187 A, B, C, D = map(_none_to_empty_2d, (A, B, C, D))
188 A = _restore(A, (p, p))
189 B = _restore(B, (p, q))
190 C = _restore(C, (r, p))
191 D = _restore(D, (r, q))
193 return A, B, C, D
196def ss2tf(A, B, C, D, input=0):
197 r"""State-space to transfer function.
199 A, B, C, D defines a linear state-space system with `p` inputs,
200 `q` outputs, and `n` state variables.
202 Parameters
203 ----------
204 A : array_like
205 State (or system) matrix of shape ``(n, n)``
206 B : array_like
207 Input matrix of shape ``(n, p)``
208 C : array_like
209 Output matrix of shape ``(q, n)``
210 D : array_like
211 Feedthrough (or feedforward) matrix of shape ``(q, p)``
212 input : int, optional
213 For multiple-input systems, the index of the input to use.
215 Returns
216 -------
217 num : 2-D ndarray
218 Numerator(s) of the resulting transfer function(s). `num` has one row
219 for each of the system's outputs. Each row is a sequence representation
220 of the numerator polynomial.
221 den : 1-D ndarray
222 Denominator of the resulting transfer function(s). `den` is a sequence
223 representation of the denominator polynomial.
225 Examples
226 --------
227 Convert the state-space representation:
229 .. math::
231 \dot{\textbf{x}}(t) =
232 \begin{bmatrix} -2 & -1 \\ 1 & 0 \end{bmatrix} \textbf{x}(t) +
233 \begin{bmatrix} 1 \\ 0 \end{bmatrix} \textbf{u}(t) \\
235 \textbf{y}(t) = \begin{bmatrix} 1 & 2 \end{bmatrix} \textbf{x}(t) +
236 \begin{bmatrix} 1 \end{bmatrix} \textbf{u}(t)
238 >>> A = [[-2, -1], [1, 0]]
239 >>> B = [[1], [0]] # 2-D column vector
240 >>> C = [[1, 2]] # 2-D row vector
241 >>> D = 1
243 to the transfer function:
245 .. math:: H(s) = \frac{s^2 + 3s + 3}{s^2 + 2s + 1}
247 >>> from scipy.signal import ss2tf
248 >>> ss2tf(A, B, C, D)
249 (array([[1, 3, 3]]), array([ 1., 2., 1.]))
250 """
251 # transfer function is C (sI - A)**(-1) B + D
253 # Check consistency and make them all rank-2 arrays
254 A, B, C, D = abcd_normalize(A, B, C, D)
256 nout, nin = D.shape
257 if input >= nin:
258 raise ValueError("System does not have the input specified.")
260 # make SIMO from possibly MIMO system.
261 B = B[:, input:input + 1]
262 D = D[:, input:input + 1]
264 try:
265 den = poly(A)
266 except ValueError:
267 den = 1
269 if (prod(B.shape, axis=0) == 0) and (prod(C.shape, axis=0) == 0):
270 num = numpy.ravel(D)
271 if (prod(D.shape, axis=0) == 0) and (prod(A.shape, axis=0) == 0):
272 den = []
273 return num, den
275 num_states = A.shape[0]
276 type_test = A[:, 0] + B[:, 0] + C[0, :] + D
277 num = numpy.zeros((nout, num_states + 1), type_test.dtype)
278 for k in range(nout):
279 Ck = atleast_2d(C[k, :])
280 num[k] = poly(A - dot(B, Ck)) + (D[k] - 1) * den
282 return num, den
285def zpk2ss(z, p, k):
286 """Zero-pole-gain representation to state-space representation
288 Parameters
289 ----------
290 z, p : sequence
291 Zeros and poles.
292 k : float
293 System gain.
295 Returns
296 -------
297 A, B, C, D : ndarray
298 State space representation of the system, in controller canonical
299 form.
301 """
302 return tf2ss(*zpk2tf(z, p, k))
305def ss2zpk(A, B, C, D, input=0):
306 """State-space representation to zero-pole-gain representation.
308 A, B, C, D defines a linear state-space system with `p` inputs,
309 `q` outputs, and `n` state variables.
311 Parameters
312 ----------
313 A : array_like
314 State (or system) matrix of shape ``(n, n)``
315 B : array_like
316 Input matrix of shape ``(n, p)``
317 C : array_like
318 Output matrix of shape ``(q, n)``
319 D : array_like
320 Feedthrough (or feedforward) matrix of shape ``(q, p)``
321 input : int, optional
322 For multiple-input systems, the index of the input to use.
324 Returns
325 -------
326 z, p : sequence
327 Zeros and poles.
328 k : float
329 System gain.
331 """
332 return tf2zpk(*ss2tf(A, B, C, D, input=input))
335def cont2discrete(system, dt, method="zoh", alpha=None):
336 """
337 Transform a continuous to a discrete state-space system.
339 Parameters
340 ----------
341 system : a tuple describing the system or an instance of `lti`
342 The following gives the number of elements in the tuple and
343 the interpretation:
345 * 1: (instance of `lti`)
346 * 2: (num, den)
347 * 3: (zeros, poles, gain)
348 * 4: (A, B, C, D)
350 dt : float
351 The discretization time step.
352 method : str, optional
353 Which method to use:
355 * gbt: generalized bilinear transformation
356 * bilinear: Tustin's approximation ("gbt" with alpha=0.5)
357 * euler: Euler (or forward differencing) method ("gbt" with alpha=0)
358 * backward_diff: Backwards differencing ("gbt" with alpha=1.0)
359 * zoh: zero-order hold (default)
360 * foh: first-order hold (*versionadded: 1.3.0*)
361 * impulse: equivalent impulse response (*versionadded: 1.3.0*)
363 alpha : float within [0, 1], optional
364 The generalized bilinear transformation weighting parameter, which
365 should only be specified with method="gbt", and is ignored otherwise
367 Returns
368 -------
369 sysd : tuple containing the discrete system
370 Based on the input type, the output will be of the form
372 * (num, den, dt) for transfer function input
373 * (zeros, poles, gain, dt) for zeros-poles-gain input
374 * (A, B, C, D, dt) for state-space system input
376 Notes
377 -----
378 By default, the routine uses a Zero-Order Hold (zoh) method to perform
379 the transformation. Alternatively, a generalized bilinear transformation
380 may be used, which includes the common Tustin's bilinear approximation,
381 an Euler's method technique, or a backwards differencing technique.
383 The Zero-Order Hold (zoh) method is based on [1]_, the generalized bilinear
384 approximation is based on [2]_ and [3]_, the First-Order Hold (foh) method
385 is based on [4]_.
387 References
388 ----------
389 .. [1] https://en.wikipedia.org/wiki/Discretization#Discretization_of_linear_state_space_models
391 .. [2] http://techteach.no/publications/discretetime_signals_systems/discrete.pdf
393 .. [3] G. Zhang, X. Chen, and T. Chen, Digital redesign via the generalized
394 bilinear transformation, Int. J. Control, vol. 82, no. 4, pp. 741-754,
395 2009.
396 (https://www.mypolyuweb.hk/~magzhang/Research/ZCC09_IJC.pdf)
398 .. [4] G. F. Franklin, J. D. Powell, and M. L. Workman, Digital control
399 of dynamic systems, 3rd ed. Menlo Park, Calif: Addison-Wesley,
400 pp. 204-206, 1998.
402 """
403 if len(system) == 1:
404 return system.to_discrete()
405 if len(system) == 2:
406 sysd = cont2discrete(tf2ss(system[0], system[1]), dt, method=method,
407 alpha=alpha)
408 return ss2tf(sysd[0], sysd[1], sysd[2], sysd[3]) + (dt,)
409 elif len(system) == 3:
410 sysd = cont2discrete(zpk2ss(system[0], system[1], system[2]), dt,
411 method=method, alpha=alpha)
412 return ss2zpk(sysd[0], sysd[1], sysd[2], sysd[3]) + (dt,)
413 elif len(system) == 4:
414 a, b, c, d = system
415 else:
416 raise ValueError("First argument must either be a tuple of 2 (tf), "
417 "3 (zpk), or 4 (ss) arrays.")
419 if method == 'gbt':
420 if alpha is None:
421 raise ValueError("Alpha parameter must be specified for the "
422 "generalized bilinear transform (gbt) method")
423 elif alpha < 0 or alpha > 1:
424 raise ValueError("Alpha parameter must be within the interval "
425 "[0,1] for the gbt method")
427 if method == 'gbt':
428 # This parameter is used repeatedly - compute once here
429 ima = np.eye(a.shape[0]) - alpha*dt*a
430 ad = linalg.solve(ima, np.eye(a.shape[0]) + (1.0-alpha)*dt*a)
431 bd = linalg.solve(ima, dt*b)
433 # Similarly solve for the output equation matrices
434 cd = linalg.solve(ima.transpose(), c.transpose())
435 cd = cd.transpose()
436 dd = d + alpha*np.dot(c, bd)
438 elif method == 'bilinear' or method == 'tustin':
439 return cont2discrete(system, dt, method="gbt", alpha=0.5)
441 elif method == 'euler' or method == 'forward_diff':
442 return cont2discrete(system, dt, method="gbt", alpha=0.0)
444 elif method == 'backward_diff':
445 return cont2discrete(system, dt, method="gbt", alpha=1.0)
447 elif method == 'zoh':
448 # Build an exponential matrix
449 em_upper = np.hstack((a, b))
451 # Need to stack zeros under the a and b matrices
452 em_lower = np.hstack((np.zeros((b.shape[1], a.shape[0])),
453 np.zeros((b.shape[1], b.shape[1]))))
455 em = np.vstack((em_upper, em_lower))
456 ms = linalg.expm(dt * em)
458 # Dispose of the lower rows
459 ms = ms[:a.shape[0], :]
461 ad = ms[:, 0:a.shape[1]]
462 bd = ms[:, a.shape[1]:]
464 cd = c
465 dd = d
467 elif method == 'foh':
468 # Size parameters for convenience
469 n = a.shape[0]
470 m = b.shape[1]
472 # Build an exponential matrix similar to 'zoh' method
473 em_upper = linalg.block_diag(np.block([a, b]) * dt, np.eye(m))
474 em_lower = zeros((m, n + 2 * m))
475 em = np.block([[em_upper], [em_lower]])
477 ms = linalg.expm(em)
479 # Get the three blocks from upper rows
480 ms11 = ms[:n, 0:n]
481 ms12 = ms[:n, n:n + m]
482 ms13 = ms[:n, n + m:]
484 ad = ms11
485 bd = ms12 - ms13 + ms11 @ ms13
486 cd = c
487 dd = d + c @ ms13
489 elif method == 'impulse':
490 if not np.allclose(d, 0):
491 raise ValueError("Impulse method is only applicable"
492 "to strictly proper systems")
494 ad = linalg.expm(a * dt)
495 bd = ad @ b * dt
496 cd = c
497 dd = c @ b * dt
499 else:
500 raise ValueError("Unknown transformation method '%s'" % method)
502 return ad, bd, cd, dd, dt