What is Riccati Algebra Equation (ARE)?

The following formula. But $ A, Q, P \ in \ mathbb R ^ {n \ times n}, B \ in \ mathbb R ^ {n \ times m}, R \ in \ mathbb R ^ {m \ times m} $ , Q, R are definite matrix symmetric matrices, and $ A, B, Q, R $ are given to find $ P $.

A P + P A^T + Q - P B R^{-1} B^T P = 0

Scene field where ARE appears, example

It appears when finding the gain with the LQR of modern control theory. LQR is an optimal control theory, and the system $ \frac{dx}{dt} = Ax+Bu $ When can be expressed as, the control law that satisfies $ \ min J = \ int (x ^ TQx + u ^ TRu) dt $ uses $ P $ that satisfies the Riccati algebra equation above. $ u=-R^{-1}B^TPx $ The story that you can write.

solve

policy

We will implement and solve the Arimoto-Potter method (explained in the next section). At that time, the eigenvalue problem of the matrix will be solved, but this uses numpy or the like.

Arimoto / Potter method

The procedure will be explained. See books for proof.

1. Place a matrix

The matrix $ H \ in \ mathbb R ^ {2n \ times 2n} $ called the Hamilton matrix

H=\left[ \begin{array}{cc} A^T & -B R^{-1} B^T \newline -Q & -A \end{array} \right]

far.

2. Eigenvalue decomposition of Hamilton matrix

There are n eigenvalues of H whose real part is negative. Place this as $ \ lambda_1, \ lambda_2, ..., \ lambda_n $ and the corresponding eigenvectors are $ \ vec w_1, \ vec w_2, ..., \ vec w_n $.

3. Place a matrix that divides the eigenvectors side by side

Use this to get $ Y, Z \ in \ mathbb R ^ {n \ times n} $ $ \left[ \vec w_1 \vec w_2 ... \vec w_n \right] = \left[ \begin{array}{c} Y \newline Z \end{array} \right] $ far.

4. P can be found

P=ZY^{-1}

Implementation

import numpy as np
import numpy.linalg as LA

def solve_are(A, B, Q, R):
    # 1.Put the Hamilton Matrix
    H = np.block([[A.T, -B @ LA.inv(R) @ B.T],
                  [-Q , -A]])
    # 2.Eigenvalue decomposition
    eigenvalue, w = LA.eig(H)
    # 3.Put an auxiliary matrix
    Y_, Z_ = [], []
    n = len(w[0])//2
    for i in range(2*n):
        if eigenvalue[i].real < 0.0:
            Y_.append(w.T[i][:n])
            Z_.append(w.T[i][n:])
    Y = np.array(Y_).T
    Z = np.array(Z_).T
    # 4.P is sought
    return Z @ LA.inv(Y)

test

Test code

A = np.array([[3., 1.],[0., 1.]])
B = np.array([[1.2], [1.]])
Q = np.array([[1., 0.2], [0.2, 1.0]])
R = np.array([[1.]])
P = solve_are(A, B, Q, R)
print("P")
print(P)
print("Left side of Riccati algebraic equation")
print(A@P + [email protected] + Q - P@[email protected](R)@B.T@P)

result

P
[[ 69.20010326 -66.19334596]
 [-66.19334596  67.7487967 ]]
Left side of Riccati algebraic equation
[[-1.13686838e-13  5.11590770e-13]
 [ 2.55795385e-13 -5.68434189e-13]]

I solved it.

Implement the solution of Riccati algebraic equations in Python