Approximation de modèle dynamique linéaire, depuis des données fréquentielles avec MATLAB

Ici, nous générons aléatoirement un système dynamique (stable) mono-entrée mono-sortie. Ce système dynamique est décrit par une représentation d'état d'ordre 100, c'est à dire qu'il est décrit par 100 équations d'états et une équation de sortie. Remarque: du fait de l'utilisation de la commande stabsep, quelques états peuvent être éliminés.

% Create an unknown model and compute its eigenvalues
G    = stabsep(rss(200,1,1));
eigG = eig(G);

Ce système, considéré comme inconnu, mais nous considérons que nous pouvons l'évaluer. Ici, il est échantillonné en fréquence, et plus particulièremenet, 50 fréquences de 1 rad/s à 100 rad/s, avec la fonction mdspack.freqresp (extension de freqresp de MATLAB). Cette évaluation est équivalente à exciter le système par un signal temporel (riche en fréquence dans cet intervalle), puis à appliquer la transformée de Fourier.

% Compute the associated frequency domain set
w   = logspace(0,2,50);
Gjw = mdspack.freqresp(G,w);

La figure ci-dessous illustre la réponse fréquentielle (diagrame de Bode - gain et phase) du système original (inconnu), ainsi que l'évaluation (les données) dont nous disposons. De plus, la figure de droite présente les valeurs propres du système original (également non connues).

Grâce à librairie MDSPACK et la fonction mdspack.loewner, il est possible de trouver une approximation rationnelle directement depuis les données. Dans la mesure où l'approximant n'est pas nécessairement stable, il est également possible de le contraindre à l'être avec mdspack.sproj.

% Data interpolation with the Loewner method
[Gr1,info1] = mdspack.loewner(w,Gjw);
% Enforce identified model stability, and eigenvalues computation
opt.space   = 'I';
Gr1         = mdspack.sproj(Gr1,opt);
eigGr1      = eig(Gr1);

Le diagramme de Bode et les valeurs propres du modèle initial et identifié sont donnés dans la figure suivante:

A noter que l'ordre du modèle (i.e. la compléxité) a été choisie automatiquement. Dans notre cas, le nombre de données ont suffit à retrouver le modèle avec une précision machine. Afin de simplifier son analyse, il est également possible de réduire davantage la compléxité du modèle de la manière suivante (ici nous cherchons une compléxité de 20 états).

% Reduce the model complexity
[Gr2,info2] = mdspack.loewner(w,Gjw,struct('target',20));
Gr2         = mdspack.sproj(Gr2,struct('space','I'));
eigGr2      = eig(Gr2);

De la même manière, nous pouvons nous focaliser uniquement sur une zone fréquentielle, permettant de simplifier le modèle identifié, en gardant une bonne précision dans l'intervalle d'intérêt. Ici nous cherchons une bonne approximation en basses fréquences. Ceci est particulièrement intéressant pour simuler ou contrôler des sous-ensemble du système ou phénomène.

% Reduce over frequency band range
[Gr3,info3] = mdspack.loewner(w,Gjw,struct('w_band',[0 4]));
Gr3         = mdspack.sproj(Gr3,struct('space','I'));
eigGr3      = eig(Gr3);

A présent, nous pouvons comparer les réponses impulsionnelles associées. On notera que le premier modèle capture parfaitement les dynamiques du système. Le second donne une approximation qui conserve la tendance et les oscilations. Le dernier, qui se focalise seulemnt sur les basses fréquences ne restitue que la tendance et exclut les oscillation autour des 1.3 Hz (ou env. 8 rad/s).

Approximation de modèle dynamique linéaire, depuis des données fréquentielles avec Python

Ici, nous générons un système dynamique (stable) de dimension infinie, mono-entrée mono-sortie. Ce système dynamique est décrit par une fonction de transfert, c'est à dire une fonction complexe qui lie une entrée à une sortie. Remarque: la fonction de transfert est dite de dimension infinie car elle est analytique et son dénominateur à une infinité de singularités (du fait cause du terme exponentiel, caractéristique des systèmes dit à retard).

% Create an unknown model 
den = lambda s: (s**2 + s*np.exp(-4/5*s) + 1)
Hs  = lambda s: (s**1/3)/den(s)#(s**1/3)/den(s)

Ce système est considéré comme inconnu, mais nous considérons que nous pouvons l'évaluer. Ici, il est échantillonné suivant la variable "s" de Laplace, le long de l'axe imaginaire, c'est à dire en fréquence. Plus particulièremenet, 15 fréquences de 0.1 rad/s à 31 rad sont calculées. Cette évaluation est équivalente à exciter le système par un signal temporel (riche en fréquence dans cet intervalle), puis à appliquer la transformée de Fourier.

% Compute the associated frequency domain set
n       = 15
w       = list(np.logspace(-.5,1,n)) 
data    = [0 for i in range(n)]
for i in range(n):
data[i] = [[Hs(complex(0,w[i]))]]

La figure ci-dessous illustre la réponse fréquentielle (diagrame de Bode - gain et phase) du système original (inconnu), ainsi que l'évaluation (les données) dont nous disposons. Il faut noter que la réponse que nous souhaitons approcher n'est pas rationnelle, et que les points dont nous disposons dans ce sont peu nombreux.

Grâce à librairie MDSPACK et la fonction mdspacky.loewner, il est possible de trouver une approximation rationnelle directement depuis les données. Dans la mesure où l'approximant n'est pas nécessairement stable, il est également possible de le contraindre à l'être avec les options disponibles; ici stable=True. Ensuite, la réponse fréquentielle est évaluée par mdspacky.ssfr.

% Data interpolation with the Loewner method
[A1, B1, C1, D1, E1]    = mdspacky.loewner(w, data, stable=True)
Hjw1                    = mdspacky.ssfr(A1, B1, C1, D1, E1, W)
eigHr1                  = linalg.eigvals(A1, E1)

Le diagramme de Bode du modèle initial et identifié sont donnés dans la figure suivante:

A noter que l'ordre du modèle (i.e. la compléxité) a été choisie automatiquement. Dans notre cas, le nombre de données ont suffit à retrouver le modèle avec une grande précision. Dans la mesure où la fonction initiale est de dimension infinie, la précision machine serait obtenue avec plus de samples, sur une bande de fréquence plus large. Afin de simplifier son analyse, il est également possible de nous focaliser uniquement sur une zone fréquentielle, permettant de simplifier le modèle identifié, en gardant une bonne précision dans l'intervalle d'intérêt. Ici nous cherchons une bonne approximation en basses fréquences. Ceci est particulièrement intéressant pour simuler ou contrôler des sous-ensemble du système ou phénomène.

% Reduce over frequency band range
[A2, B2, C2, D2, E2]    = mdspacky.loewner(w, data, wband=[0,2], stable=True)
Hjw2                    = mdspacky.ssfr(A2, B2, C2, D2, E2, W)
eigHr2                  = linalg.eigvals(A2, E2)

Ici aussi, l'ordre est estimé automatiquement, indiquant un modèle plus simple, et pour autant, précis dans l'intervalle d'intérêt. Le diagramme de Bode du modèle initial et identifié sont donnés dans la figure suivante:

A présent, nous pouvons évaluer simplement les réponses à un échelon relative aux modèles rationnels obtenus. Le calcul de la réponse à un échelon du modèle initial est plus compliqué à obtenir car celui-ci est irrationnel et requiert un calcul dédié.