Asservissement par LQR

On a un système sous forme de représentation d’états. On veut asservir ce système par placement de pôles, mais cet asservissement doit répondre à un cahier des charges (contraintes physiques, optimisations). Comment optimiser l’asservissement en utilisant un Linear Quadratic Regulator (LQR) ?

Principes

L’asservissement LQR peut être démontré par l’équation HJB (Hamilton-Jacobi-Bellman) ou par le principe du maximum de Pontryagin. Le but est de déterminer les composantes du gain K pour une loi de commande par retour d’état simple :

$u(t) = -Kx(t)$

sur un horizon de temps infini (cas le plus courant et le plus pratique), tel que l’on minimise un critère J défini comme suit :

$J = \int_{0}^{+\infty} x(t)^T Q_c x(t) + u(t)^T R_c u(t) dt$

où Q_c et R_c sont des matrices choisies préalablement tel que : $Q_c \geq 0 \quad et \quad R_c \gt 0$

Dans le cas où x est estimé par un observateur, et que cet observateur est un filtre de Kalman, alors le processus est appelé LQG pour Linear Quadratic Gaussian.

Observations

On remarque alors que :

Le système est linéaire : $\begin{cases} \dot{x}(t) &= A x(t) + B u(t) \\ y(t) &= C x(t) + D u(t) \end{cases}$
L’optimisation est quadratique dans le critère J.
On effectue une régulation car l’on est en train de déterminer les paramètres de la loi de commande du système ; loi qui elle même est linéaire.

Résolution

Sans rentrer dans les détails de la résolution de l’HJB pour le cas LQR (ce qui est long, chiant et plus mathématique que de l’ingénierie à proprement parler), trouver K qui minimise J passe par la résolution de l’équation algébrique de Riccati (ARE) :

$A^T P + P A - P B R_c^{-1} B^T P + Q_c = 0$

Lorsque l’on a trouvé P, on peut alors déterminer K :

$K = R_c^{-1} B^T P$

Optimisation

Si on est passé par le LQR, c’est pour éviter de choisir de mauvaises composantes à K lors du placement sur un retour d’état classique. Il ne faut donc pas perdre le sens physique des choses. Ainsi :

Plus l’on augmente les composantes de Q_c par rapport à celles de R_c, plus le système suivra fidélement la consigne, mais plus la consommation d’énergie sur le moteur sera élevée.
Plus l’on augmente les composantes de R_c par rapport à celles de Q_c, plus le système sera économe en énergie, mais il sera beaucoup plus lent pour atteindre une consigne.

Aussi, on peut se fixer comme heuristique la règle de Bryson :

$Q_{c,ii} = \frac {1} {max \> acceptable \> pour \> x_i^2(t)}, \quad R_{c,jj} = \frac {1}{max \> acceptable \> pour \> u_j^2(t)}$

Implémentations

Implémenter un LQR, ça devient vite compliqué. Il faut savoir dans quel cas l’on se trouve. J’ai énuméré la résolution en milieu continu car elle nous permet de garder le sens physique du gain K et nous assurer qu’il n’y a pas d’aberration. Mais lorsqu’il faut implémenter sur une MPU en temps réel ou une SoC, il faudra forcément passer par de la discrétisation.

À cette première difficulté s’en ajoute une autre : quid de la variation de Q_c et R_c ? À première vue, la question passerait pour une excentricité académique, mais elle a un véritable cas d’usage lorsque, par exemple, pour le décollage d’une fusée, on souhaite adopter des stratégies différentes en termes de consommation d’ergol entre le début de l’atterrissage et la fin de la descente.

Nous avons donc deux cas à traiter.

N’oublions cependant pas que, quelque soit le cas, il est nécessaire que :

$Q_c(t) = Q_c^T(t) \geq 0 \quad \forall t$ $R_c(t) = R_c^T(t) \gt 0 \quad \forall t$

Q et R ne varient pas dans le temps

On garde l’ARE. Puisque Q et R sont invariants dans le temps, K le sera aussi. Nous pouvons donc suivre la méthodologie suivante :

Identifier le système linéaire en milieu continu (A, B, C & D), voir qu’il est instable en boucle ouverte et proposer une loi de commande linéaire suivant un critère d’optimalité quadratique.
Résoudre l’ARE en milieu continu et déterminer K. Analyser physiquement K en cas d’éventuelles aberrations.
Discrétiser le système (F, G, H & J) et résoudre l’ARE en milieu discret. On obtient alors K_d, le gain pour le système discret.

Donc, une fois le système discrétisé, il faut résoudre l’ARE discrète :

$P = F^T P F - (F^T P G) \left( R_c + G^T P G \right)^{-1} (G^T P F) + Q_c$

Pour trouver K_d :

$K_d = \left( R_c + G^T P G \right)^{-1} G^T P F$

Répondant à notre loi de commande en discret :

$u_k = -K_d x_k$

Matlab

Pour calculer K_d en vue de l’implémenter, on peut utiliser Matlab.

On commence par définir le système :

A = [0 1; -2 -3];
B = [0; 1];
C = [1 0];
D = 0;

sys_c = ss(A, B, C, D);

Et les paramètres de notre critère :

Q_c = [diag([1 2])];
R_c = 2;

Et on peut ensuite obtenir différentes matrices très intéressantes :

[K, P, A_bf_eig] = lqr(sys, Q_c, R_c)

Où :

K est le gain pour notre retour d’état
P est la matrice solution de l’ARE
A_bf_eig est un vecteur contenant les valeurs propres de la matrice d’état en boucle fermée.

On peut aussi passer par le système discrétiser :

Te = 0.01;
[Kd, P, F_bf_eig] = lqrd(sys_c, Q_c, R_c, Te);

Où :

K_d est le gain du système discret dans le régime stationnaire
P est la matrice solution de l’ARE en discret
F_bf_eig est un vecteur contenant les valeurs propres de la matrice d’état en boucle fermée du système discret.

Implémentation

L’implémentation d’une telle loi de commande est alors très simple :

/// Calcule la loi de commande u = -K * x
/// - `M` : dimension de la commande (taille de u, nombre de lignes de K)
/// - `N` : dimension de l'état (taille de x, nombre de colonnes de K)
pub fn cmd_law<const M: usize, const N: usize>(
    k: &[[f32; N]; M],
    x: &[f32; N],
) -> [f32; M] {
    let mut u = [0.0f32; M];

    for i in 0..M {
        let mut sum = 0.0f32;
        for j in 0..N {
            sum += k[i][j] * x[j];
        }
        u[i] = -sum; 
    }

    u
}

On peut alors hardcode la matrice K comme une constante :

let k: [[f32; 2]; 1] = [[1.5, 0.8]];

Aussi, il est possible d’ajouter une précompensation statique N à la loi de commande pour permettre le suivi d’une consigne r(t).

Cette implémentation nous fait gagné numériquement en calcul, mais nous perdons en efficacité lors du régime transitoire.

Cas où Q(t) et R(t) dépendent du temps

On abandonne l’ARE pour lui préférer l’équation différentielle de Riccati (DRE) :

$-\dot{P}(t) = A^T P(t) + P(t) A - P(t) B R_c(t)^{-1} B^T P(t) + Q_c(t)$

Le truc, c’est que l’on a peu d’intérêt de la résoudre. Q_c(t) et R_c(t) varient dans le temps. Il va donc falloir constamment recalculer K sur le micro-contrôleur. Donc pas de résolution temporelle et passage direct par la discrétisation.

Il ne faudra pas oublier de discrétiser également Q_c et R_c.

On cherche alors le gain K_k pour notre loi de commande :

$u_k = -K_k x_k$

Il est donné par la formule en temps discret :

$K_k = \left( R_{c,k} + G^T P_{k+1} G \right)^{-1} G^T P_{k+1} F$

Tel que P est régi par la DRE en milieu discret :

$P_k = F^T P_{k+1} F - (F^T P_{k+1} G) \left( R_{c,k} + G^T P_{k+1} G \right)^{-1} (G^T P_{k+1} F) + Q_{c,k}$

En regardant de plus près cette équation, on se rend compte d’un problème fondamental : l’équation aux différences de Riccati s’intégre à rebours dans le temps. Notre malheur vient du fait que pour obtenir P(k) je dois connaître P(k+1).

En soi, la situation se résout d’elle même si nous connaissons à l’avance toute la trajectoire. Ainsi, l’implémentation d’une trajectoire pour un bras robotique sur un laps de temps de 10 secondes ne poserait pas de problème : on calcule toutes les occurrences de P (et donc par extension, de K) sur Matlab à l’avance, et nous serons rendus au cas où Q_c et R_c ne varient pas, avec la subtilité que K serait plutôt K_tab, une liste de vecteurs K répondant au problème d’optimalité LQR.

En revanche, si la trajectoire est indéterminée, il faudra passer par des modèles d’implémentation plus sophistiqués comme le MPC.