Problem 668561 · medium · Level 06 Heuristics & Optimization

The Crane Model Back as a Difference Equation

state space · transfer function · characteristic polynomial · impulse response · Markov parameters · numpy

An overhead crane's sway model was built in state space, x[k+1] = A x[k] + B u[k], y[k] = C x[k] (one input, one output). The cheap microcontroller on the hook only runs difference equations, so you need the same system as

y[k] + a[1]*y[k-1] + ... + a[n]*y[k-n] = b[1]*u[k-1] + ... + b[n]*u[k-n]

with a = [1, a[1], ..., a[n]] and b = [0, b[1], ..., b[n]] (coefficient lists in z^-1, constant term first; b[0] = 0 because u[k] first reaches y one step later).

  • The denominator is the characteristic polynomial of A: numpy.poly(A) returns [1, a[1], ..., a[n]], the coefficients of det(z I - A), highest power first. Its roots are the eigenvalues of A, which are the poles. (Take numpy.real(...): the coefficients of a real matrix are real.)
  • The numerator comes from the impulse response. Start from rest and give one input kick u[0] = 1: then x[1] = B, x[2] = A B, ..., so h[0] = 0 and h[m] = C A^(m-1) B for m >= 1 (these numbers are also called the Markov parameters). The difference equation says that a convolved with h is b:
b[j] = a[0]*h[j] + a[1]*h[j-1] + ... + a[j]*h[0]        for j = 1, ..., n

Write to_difference(A, B, C) that returns (a, b), both lists of n + 1 floats.

Examples

Input:  A = [[1.5, -0.7], [1, 0]], B = [1, 0], C = [0.4, 0.2]
Output: ([1.0, -1.5, 0.7000000000000001], [0.0, 0.4, 0.19999999999999996])
Explanation: a model in companion form gives back its own equation:
y[k] = 1.5 y[k-1] - 0.7 y[k-2] + 0.4 u[k-1] + 0.2 u[k-2].

Input:  A = [[1, 0.5], [0, 1]], B = [0.125, 0.5], C = [1, 0]
Output: ([1.0, -2.0, 1.0], [0.0, 0.125, 0.125])
Explanation: a lift's height for an acceleration held for 0.5 s. h[1] = 0.125, h[2] = 0.375,
so b[2] = 0.375 - 2 * 0.125 = 0.125.

Input:  A = [[0.8, 0], [0, 0.5]], B = [1, 1], C = [1, 0]
Output: ([1.0, -1.3, 0.4], [0.0, 1.0, -0.5])
Explanation: the sensor never sees the second state. b = 1 - 0.5 z^-1 cancels the pole at 0.5:
the output behaves like y[k] = 0.8 y[k-1] + u[k-1].

Constraints

  • 1 <= n <= 6
  • answers are compared with a tolerance of 1e-6; numpy arrays are accepted in place of the lists

Goals

  • Get a state-space model's denominator from the characteristic polynomial of A
  • Compute the impulse response C A^(m-1) B of a state-space model
  • Recover the numerator by convolving the denominator with the impulse response
Starting Python…