Implémenter la solution de l'algèbre de Riccati en Python

Qu'est-ce que l'équation algébrique de Riccati (ARE)?

La formule suivante. Cependant, $ 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 sont des matrices symétriques à valeur constante positive, et $ A, B, Q, R $ sont donnés pour trouver $ P $.

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

Champ de scène où ARE apparaît, exemple

Il apparaît lors de la recherche du gain avec LQR de la théorie du contrôle moderne. LQR est la théorie du contrôle optimal, et le système $ \frac{dx}{dt} = Ax+Bu $ Quand peut être exprimé comme, la loi de contrôle qui satisfait $ \ min J = \ int (x ^ TQx + u ^ TRu) dt $ utilise $ P $ qui satisfait l'équation d'algèbre de Ricati ci-dessus. $ u=-R^{-1}B^TPx $ L'histoire que vous pouvez écrire.

résoudre

politique

La méthode d'Arimoto-Potter (expliquée dans la section suivante) sera implémentée et résolue. A ce moment, le problème des valeurs propres de la matrice sera résolu, mais cela utilise numpy ou similaire.

Méthode Arimoto / Potter

La procédure sera expliquée. Voir les livres, etc. pour la preuve.

1. Placez une file d'attente

Matrice $ H \ in \ mathbb R ^ {2n \ fois 2n} $ appelée Hamilton Matrix

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

loin.

2. Décomposition de la matrice de Hamilton par valeur propre

Il existe n valeurs propres de H dont la partie réelle est négative. Placez ceci comme $ \ lambda_1, \ lambda_2, ..., \ lambda_n $ et les vecteurs propres correspondants comme $ \ vec w_1, \ vec w_2, ..., \ vec w_n $.

3. Placer une matrice divisée en disposant les vecteurs propres

Utilisez ceci pour obtenir $ 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] $ loin.

4. P peut être trouvé

P=ZY^{-1}

la mise en oeuvre

import numpy as np
import numpy.linalg as LA

def solve_are(A, B, Q, R):
    # 1.Mettez la matrice de Hamilton
    H = np.block([[A.T, -B @ LA.inv(R) @ B.T],
                  [-Q , -A]])
    # 2.Décomposer les valeurs propres
    eigenvalue, w = LA.eig(H)
    # 3.Mettre une matrice auxiliaire
    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 est recherché
    return Z @ LA.inv(Y)

tester

Code de test

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("Côté gauche de l'équation algébrique de Riccati")
print(A@P + [email protected] + Q - P@[email protected](R)@B.T@P)

résultat

P
[[ 69.20010326 -66.19334596]
 [-66.19334596  67.7487967 ]]
Côté gauche de l'équation algébrique de Riccati
[[-1.13686838e-13  5.11590770e-13]
 [ 2.55795385e-13 -5.68434189e-13]]

Je l'ai résolu.