Solve discrete-time algebraic ricatti equation for affine dynamical system using scipy.linalg.solve_discrete_are

control theorydynamical systemsoptimal control

I am trying to solve a discrete algebraic ricatti equation for an affine dynamical system (ADS) of the form:

$$x_{t+1} = A x_t + B u_t + b$$

where $x_t, u_t$ and $b$ are vectors and $A$ and $b$ are matrices.

As described here, I convert the ADS to a linear dynamical system (LDS) as follows:

$$x_{t+1}^{new} = A^{new} x_t^{new} + B^{new} u_t,$$

where

$$x_t^{new} = \begin{bmatrix} x_t\\ 1 \end{bmatrix}, A^{new} = \begin{bmatrix} A & b\\ 0 & 1 \end{bmatrix}, B^{new} = \begin{bmatrix} B\\ 0 \end{bmatrix},$$

However, when I try to solve the discrete-time algebraic ricatti equation for this LDS by calling scipy.linalg.solve_discrete_are in python:

import numpy as np
from scipy.linalg import solve_discrete_are, solve

x_dim = 2
u_dim = 1
A = np.ones((x_dim, x_dim))
B = np.ones((x_dim, u_dim))
b = np.ones((x_dim,))

A_new = np.zeros((x_dim + 1, x_dim + 1))
A_new[:x_dim, :x_dim] = A
A_new[:x_dim, -1] = b
A_new[-1, -1] = 1

B_new = np.zeros((x_dim + 1, u_dim))
B_new[:x_dim, :] = B
B_new[-1, :] = 0

Q = np.eye(x_dim + 1)
R = np.eye(u_dim)

S = np.zeros((x_dim + 1, u_dim))
X = solve_discrete_are(A_new, B_new, Q, R, e=None, s=S)
K = solve(B_new.T @ X @ B_new + R, B_new.T @ X @ A_new + S.T)

I get the following error message

Traceback (most recent call last):
  File ".../LQR.py", line 23, in <module>
    X = solve_discrete_are(A_new, B_new, Q, R, e=None, s=S)
  File ".../lib/python3.9/site-packages/scipy/linalg/_solvers.py", line 716, in solve_discrete_are
    raise LinAlgError('Failed to find a finite solution.')
numpy.linalg.LinAlgError: Failed to find a finite solution.

Can you someone tell me why this isn't working and how I can make it work?

Best Answer

Instead of extending the state space model with a marginally stable state (eigenvalue of one) which is uncontrollable, one can try and see if a coordinate transformation is possible, similar to the continuous time variant discussed here.

If your system is of the form

$$ x_{k+1} = A\,x_k + B\,u_k + f \\ y_k = C\,x_k + D\,u_k + g \tag{1} $$

with $x_k$, $u_k$ and $y_k$ the state, input and output of the system respectfully at time step $k$, $f$ and $g$ constant vectors. The following coordinate transformation could be used to make the affine system linear $z_k=x_k+\alpha$ and $v_k=u_k+\beta$ with $\alpha$ and $\beta$ constant vectors which satisfy

$$ \begin{bmatrix} A-I & B \\ C & D \end{bmatrix} \begin{bmatrix} \alpha \\ \beta \end{bmatrix} = \begin{bmatrix} f \\ g \end{bmatrix}. \tag{2} $$

This can always be solved if the $(A-I,B,C,D)$ matrix is full rank and if this is not the case when the vector of $(f,g)$ lies in the span of the $(A-I,B,C,D)$ matrix. After this transformation the dynamics will simply be

$$ z_{k+1} = A\,z_k + B\,v_k \\ y_k = C\,z_k + D\,v_k. \tag{3} $$

It can be noted that when applying a control technique to $(3)$, which would drive $z_k$ for example to the origin, would drive the state $x_k$ from $(1)$ to a value that is offset by $-\alpha$.

Related Question