L’identification de modèle par ARX (AutoRegressive with eXogenous input) est intéressante si l’on a besoin d’une estimation très rapide et discrétisée.

Elle a cependant pour grand désavantage de reposer sur l’hypothèse d’un bruit blanc gaussien comme seul perturbateur de la mesure.

Principes

On connaît l’équation aux différences d’un système linéaire invariant dans le temps :

$y(k) + a_1 y(k-1) + \>...\> + a_n y(k-n) = b_1 u(k) + b_2 u(k - 1) + \>...\> + b_m u(k - m)$

L’idée c’est de trouver tous les paramètres a et b dans l’équation pour avoir le modèle complet :

$G(z) = \frac {Y(z)}{U(z)} = \frac {B(z)} {A(z)}$

L’idée est alors de passer par les moindres carrés. Prenons un cas simple avec une équation aux différences courtes :

$y(k) + a_1 y(k-1) + a_2 y(k - 2) = b_1 u(k)$

On exprime tout en fonction de y(k) :

$y(k) = -a_1 y(k-1) - a_2 y(k - 2) + b_1 u(k)$

On pose le vecteur d’inconnus θ :

$\theta = \begin{bmatrix} a_1 \\ a_2 \\ b_1 \end{bmatrix}$

On va alors regarder notre jeu de données et constituer le vecteur Y et la matrice Φ tel que :

$Y = \Phi \times \theta$

Autrement dit :

$\begin{bmatrix} y(2) \\ y(3) \\ y(4) \\ ... \\ y(k) \end{bmatrix} = \begin{bmatrix} -y(1) & -y(0) & u(2) \\ -y(2) & -y(1) & u(3) \\ -y(3) & -y(2) & u(4) \\ & ... & \\ -y(k-1) & -y(k-2) & u(k) \end{bmatrix} \times \begin{bmatrix} a_1 \\ a_2 \\ b_1 \end{bmatrix}$

On a alors comme solution θ :

$\theta = (\Phi^T \Phi)^{-1} \Phi^T Y$

Matlab

On va dire que l’on a un fichier matlab que l’on importe dans notre environnement et qui contient toutes nos données u et y. On va alors poser :

N = length(y);
z = iddata(y, u, Ts);

On va ensuite chercher à isoler les variations. On va donc supprimer la composante continue des signaux :

z_d = detrend(z);

Puis on va séparer notre jeu de données (qui j’espère est colossale) en deux parties, une pour l’estimation du modèle, l’autre pour la validation du modèle :

N_split = floor(0.7 * N);
z_est = z_d(1:N_split);
z_val = z_d(N_split+1:end);

Ensuite, on va générer un tableau comprennant tous les systèmes cohérents possibles pour notre exemple :

NN = struc(1:4, 1:4, 1:3);

Ici, les paramètres correspondent à :

1. L’intervalle de valeurs de na : L’ordre des pôles et donc du système. Soit on regarde la réponse indicielle du système (si on la possède) et on estime à peu près l’ordre, soit on part du principe qu’un cas industriel classique ne montera pas au dessus de 4.
2. L’intervalle de valeurs de nb : Nombre de zéros du système. On prendra la même plage que na.
3. L’intervalle de valeurs de nk : C’est le retard pur. Il faut compter le nombre de pas d’échantillonnage Te entre le premier changement sur u(t) et le premier changement sur y(t). Le minimum est toujours 1, car c’est le retard causal minimum d’un échantillonneur bloqueur.

Puis on va faire tester toutes ces possibilités :

err_scores = arxstruc(z_est, z_val, NN);

Ici, err_scores contient tous les modèles qui ont été estimés à partir de z_est et qui ont ensuit été injectés dans z_val pour les scorer.

On pourra ensuite récupérer le meilleur d’entre eux avec :

best_params = selstruc(err_scores, 'mdl');

'mdl' est un paramètre qu’on passe à la fonction pour sélectionner le meilleur système. C’est l’agorithme employé pour trier les modèles. Pour faire court, il pénalise les ordres plus élevés s’ils n’apportent pas un score significativement plus élevés. On peut aussi utiliser 'aic' qui repose sur le même principe avec des pénalités un peu moins élevées. Dans les deux cas, c’est pour éviter l’overfitting où un

Ensuite, on va récupérer le modèle à partir des paramètres na, nb et nk trouvés (on peut les visualiser dans best_params) :

sys_arx = arx(z_est, best_params);

Après on peut faire de l’analyse de systèmes en obtenant son score :

figure;
compare(z_val, sys_arx);

Ou bien ses résidus, notamment pour vérifier que nous étions bien dans le cas où l’hypothèse du bruit blanc est valide :

figure;
resid(z_val, sys_arx);

Si ce test échoue, il peut être intéressant de se tourner vers d’autres modèles d’estimation comme ARMAX, ou le passage en continu avec SFRIVC. Si l’identification est simplement pour tester quelques petits comportement, OE peut être intéressant.