Main Content

La traduction de cette page n'est pas à jour. Cliquez ici pour voir la dernière version en anglais.

Identification de modèles linéaires avec l'interface en lignes de commande

Introduction

Objectifs

Estimer et valider des modèles linéaires à partir de données à entrées multiples/sortie unique (MISO) pour trouver celui qui décrit le mieux la dynamique du système.

À la fin de ce tutoriel, vous serez capable d’accomplir les tâches suivantes à l’aide de la ligne de commande :

  • Créer des objets de type data pour représenter les données.

  • Tracer les données.

  • Traiter les données en supprimant les décalages des signaux d’entrée et de sortie.

  • Estimer et valider des modèles linéaires à partir de données.

  • Simuler et prédire la sortie des modèles.

Remarque

Ce tutoriel utilise des données dans le domaine temporel pour montrer comment estimer des modèles linéaires. Le même workflow s’applique à l’ajustement de données dans le domaine fréquentiel.

Description des données

Ce tutoriel utilise le fichier de données co2data.mat, qui contient les données dans le domaine temporel de deux expérimentations sur un système à deux entrées et une sortie unique (MISO) perturbé par l’opérateur depuis un état stationnaire autour de ses valeurs d’équilibre.

Dans la première expérimentation, l’opérateur a introduit une onde pulsée sur les deux entrées. Dans la deuxième expérimentation, l’opérateur a introduit une onde pulsée sur la première entrée et un signal indiciel sur la deuxième entrée.

Préparation des données

Chargement des données dans l'espace de travail MATLAB

Chargez les données.

load co2data

Cette commande charge les cinq variables suivantes dans l’espace de travail MATLAB :

  • Input_exp1 et Output_exp1 sont respectivement les données d’entrée et de sortie de la première expérimentation.

  • Input_exp2 et Output_exp2 sont respectivement les données d’entrée et de sortie de la deuxième expérimentation.

  • Time est le vecteur temps de 0 à 1 000 minutes, augmentant par incréments égaux de 0.5 min.

Pour les deux expérimentations, les données d’entrée consistent en deux colonnes de valeurs. La première colonne de valeurs est le taux de consommation de produits chimiques (en kilogrammes par minute) et la deuxième colonne de valeurs est l'intensité du courant (en ampères). Les données de sortie sont représentées dans une colonne unique indiquant le taux de production de dioxyde de carbone (en milligrammes par minute).

Tracer les données d’entrée/sortie

Tracez les données d’entrée et de sortie des deux expérimentations.

plot(Time,Input_exp1,Time,Output_exp1)
legend('Input 1','Input 2','Output 1')
figure
plot(Time,Input_exp2,Time,Output_exp2)
legend('Input 1','Input 2','Output 1')

Le premier graphique montre la première expérimentation, où l’opérateur applique une onde pulsée à chaque entrée pour perturber le système de son équilibre en état stationnaire.

Le deuxième graphique montre la deuxième expérimentation, où l’opérateur applique une onde pulsée à la première entrée et un signal indiciel à la deuxième entrée.

Supprimer les valeurs d’équilibre des données

Tracer les données, comme décrit dans Tracer les données d’entrée/sortie, montre que les entrées et les sorties ont des valeurs d’équilibre non nulles. Dans cette partie du tutoriel, vous allez soustraire les valeurs d’équilibre aux données.

La raison pour laquelle on soustrait les valeurs moyennes de chaque signal est que, généralement, on construit des modèles linéaires qui décrivent les réponses aux déviations d'un équilibre physique. Avec des données à l’état stationnaire, il est raisonnable de supposer que les niveaux moyens des signaux correspondent à cet équilibre. Ainsi, vous pouvez rechercher des modèles autour de zéro sans modéliser les niveaux absolus d’équilibre en unités physiques.

Faites un zoom sur les graphiques pour voir que le premier moment où l’opérateur applique une perturbation aux entrées se produit après 25 minutes de conditions d’état stationnaire (ou après les 50 premiers échantillons). Ainsi, la valeur moyenne des 50 premiers échantillons représente les conditions d’équilibre.

Utilisez les commandes suivantes pour supprimer les valeurs d’équilibre des entrées et des sorties des deux expérimentations :

Input_exp1 = Input_exp1-...
   ones(size(Input_exp1,1),1)*mean(Input_exp1(1:50,:));
Output_exp1 = Output_exp1-...
   mean(Output_exp1(1:50,:));
Input_exp2 = Input_exp2-...
   ones(size(Input_exp2,1),1)*mean(Input_exp2(1:50,:));
Output_exp2 = Output_exp2-...
   mean(Output_exp2(1:50,:));

Utiliser des objets pour représenter les données pour l'identification de système

Les objets pour la représentation de données de System Identification Toolbox™, iddata et idfrd encapsulent les valeurs des données et leurs propriétés dans une seule entité. Vous pouvez utiliser les commandes de System Identification Toolbox pour manipuler aisément ces objets de données comme des entités uniques.

Dans cette partie du tutoriel, vous allez créer deux objets iddata, un pour chacune des deux expérimentations. Vous allez utiliser les données de l’expérimentation 1 pour l’estimation du modèle et les données de l’expérimentation 2 pour la validation du modèle. Vous allez travailler avec deux jeux de données indépendants, car vous utilisez un jeu de données pour l’estimation du modèle et l’autre pour la validation du modèle.

Remarque

S’il n’y a pas deux jeux de données indépendants disponibles, vous pouvez diviser un jeu de données en deux parties, à condition que chaque partie contienne suffisamment d’informations pour représenter de manière adéquate la dynamique du système.

Création d’objets iddata

Vous devez avoir déjà chargé les données d’échantillon dans l’espace de travail MATLAB®, comme décrit dans Chargement des données dans l'espace de travail MATLAB.

Utilisez ces commandes pour créer deux objets de données, ze et zv :

Ts = 0.5; % Sample time is 0.5 min
ze = iddata(Output_exp1,Input_exp1,Ts);
zv = iddata(Output_exp2,Input_exp2,Ts);

ze contient les données de l’expérimentation 1 et zv contient les données de l’expérimentation 2. Ts est le pas d'échantillonnage.

Le constructeur de iddata nécessite trois arguments d'entrée pour des données dans le domaine temporel et a la syntaxe suivante :

data_obj = iddata(output,input,sampling_interval);

Pour visualiser les propriétés d’un objet iddata, utilisez la commande get. Par exemple, saisissez cette commande pour obtenir les propriétés des données d’estimation :

get(ze)
ans = 

  struct with fields:

              Domain: 'Time'
                Name: ''
          OutputData: [2001x1 double]
                   y: 'Same as OutputData'
          OutputName: {'y1'}
          OutputUnit: {''}
           InputData: [2001x2 double]
                   u: 'Same as InputData'
           InputName: {2x1 cell}
           InputUnit: {2x1 cell}
              Period: [2x1 double]
         InterSample: {2x1 cell}
                  Ts: 0.5000
              Tstart: 0.5000
    SamplingInstants: [2001x1 double]
            TimeUnit: 'seconds'
      ExperimentName: 'Exp1'
               Notes: {}
            UserData: []

Pour en savoir plus sur les propriétés de cet objet de données, consultez la page de référence iddata.

Pour modifier les propriétés de données, vous pouvez utiliser la notation avec le point « . ») ou la commande set. Par exemple, pour attribuer aux canaux des noms et des unités qui étiquettent les axes du graphique, saisissez la syntaxe suivante dans la fenêtre de commande MATLAB :

% Set time units to minutes
ze.TimeUnit = 'min';
% Set names of input channels
ze.InputName = {'ConsumptionRate','Current'};
% Set units for input variables
ze.InputUnit = {'kg/min','A'};
% Set name of output channel
ze.OutputName = 'Production';
% Set unit of output channel
ze.OutputUnit = 'mg/min';

% Set validation data properties
zv.TimeUnit = 'min';
zv.InputName = {'ConsumptionRate','Current'};
zv.InputUnit = {'kg/min','A'};
zv.OutputName = 'Production';
zv.OutputUnit = 'mg/min';

Vous pouvez vérifier que la propriété InputName de ze est modifiée ou « indexer » dans cette propriété :

ze.inputname;

Conseil

Les noms de propriétés, comme InputUnit, ne sont pas sensibles à la casse. Vous pouvez également abréger les noms de propriétés qui commencent par Input ou Output en remplaçant u par Input et y par Output dans le nom de la propriété. Par exemple, OutputUnit est équivalent à yunit.

Tracer des données représentées dans un objet de données

Vous pouvez tracer des objets iddata avec la commande plot.

Tracez l'erreur d’estimation.

plot(ze)

Les axes des graphiques du bas montrent les entrées ConsumptionRate et Current et l' axe du graphique du haut montrent la sortie ProductionRate.

Tracez les données de validation dans une nouvelle fenêtre Figure.

figure   % Open a new MATLAB Figure window
plot(zv) % Plot the validation data

Faites un zoom sur les graphiques pour voir que le processus expérimental amplifie la première entrée (ConsumptionRate) par un facteur 2 et amplifie la deuxième entrée (Current) par un facteur 10.

Sélectionner un sous-ensemble de données

Avant de commencer, créez un nouveau jeu de données contenant uniquement les 1 000 premiers échantillons des jeux de données d’estimation et de validation originaux pour accélérer les calculs.

Ze1 = ze(1:1000);
Zv1 = zv(1:1000);

Pour plus d'informations sur l’indexation dans des objets iddata, consultez la page de référence correspondante.

Estimer des modèles de réponse impulsionnelle

Pourquoi estimer des modèles de réponse indicielle et de réponse fréquentielle ?

Les réponses fréquentielles et les réponses indicielles sont des modèles non paramétriques qui peuvent vous aider à comprendre les caractéristiques de la dynamique de votre système. Ces modèles ne sont pas représentés par une formule mathématique compacte avec des paramètres ajustables. À la place, ils se composent de tables de données.

Dans cette partie du tutoriel, vous estimez ces modèles à l’aide du jeu de données ze. Vous devriez avoir déjà créé ze, comme décrit dans Création d’objets iddata.

Les graphiques de réponse de ces modèles montrent les caractéristiques suivantes du système :

  • La réponse de la première entrée à la sortie pourrait être une fonction de second ordre.

  • La réponse de la deuxième entrée à la sortie pourrait être une fonction de premier ordre ou une fonction suramortie.

Estimer la réponse fréquentielle

System Identification Toolbox propose trois fonctions pour estimer la réponse en fréquence :

  • etfe calcule la fonction de transfert empirique avec une analyse de Fourier.

  • spa estime la fonction de transfert à l’aide d’une analyse spectrale pour une résolution de fréquence fixe.

  • spafdr vous permet de spécifier une résolution de fréquence variable pour estimer la réponse en fréquence.

Utilisez spa pour estimer la réponse en fréquence.

Ge = spa(ze);

Tracez la réponse en fréquence sous la forme d’un diagramme de Bode.

bode(Ge)

L’amplitude est maximale à la fréquence de 0,54 rad/min, ce qui suggère un possible comportement de résonance (pôles complexes) pour la première combinaison entrée-sortie - ConsumptionRate par rapport Production.

Sur les deux graphiques, la phase décroît rapidement, ce qui suggère un retard pour les deux combinaisons entrée/sortie.

Estimer la réponse indicielle empirique

Pour estimer la réponse indicielle à partir des données, estimez d’abord un modèle de réponse impulsionnelle non paramétrique (filtre FIR) à partir des données, puis tracez sa réponse indicielle.

% model estimation
Mimp = impulseest(Ze1,60);

% step response
step(Mimp)

La réponse indicielle pour la première combinaison entrée/sortie suggère un dépassement, ce qui indique la présence d’un mode sous-amorti (pôles complexes) dans le système physique.

La réponse indicielle entre la deuxième entrée et la sortie ne montre pas de dépassement, ce qui indique une réponse de premier ordre ou une réponse d’ordre supérieur avec des pôles réels (réponse sur-amortie).

Le tracé de la réponse indicielle indique un retard non nul dans le système, ce qui est cohérent avec l'amortissement rapide de la phase constatée dans le diagramme de Bode que vous avez créé dans Estimer la réponse indicielle empirique.

Estimer des retards dans un système à entrées multiples

Pourquoi estimer des retards ?

Pour identifier des modèles boîte noire paramétriques, vous devez spécifier le retard en entrée/sortie pour la définition de l’ordre du modèle.

Si l’expérimentation ne vous permet pas de connaître les retards en entrée/sortie de votre système, vous pouvez utiliser System Identification Toolbox pour estimer le retard.

Estimer des retards à l’aide de la structure de modèle ARX

Dans le cas de systèmes à entrée unique, vous pouvez lire le retard sur le graphique de la réponse impulsionnelle. Cependant, dans le cas de systèmes à entrées multiples comme celui de ce tutoriel, il se peut que vous ne puissiez pas savoir quelle entrée a provoqué le changement initial de sortie et que vous deviez utiliser la commande delayest à la place.

La commande estime le retard dans un système dynamique en estimant un modèle ARX à temps discret d’ordre inférieur avec un ensemble de retards, puis en choisissant le retard qui correspond au meilleur ajustement.

La structure de modèle ARX est l’une des structures paramétriques boîte noire la plus simple. En temps discret, la structure ARX est une équation aux différences de la forme suivante :

y(t)+a1y(t1)++anay(tna)=         b1u(tnk)++bnbu(tnknb+1)+e(t)

y(t) représente la sortie à l’instant t, u(t) représente l’entrée à l’instant t, na est le nombre de pôles, nb est le nombre de paramètres b (égal au nombre de zéros plus 1), nk est le nombre d’échantillons avant que l’entrée affecte la sortie du système (appelé retard ou temps mort du modèle) et e(t) est la perturbation du bruit blanc.

delayest suppose que na=nb=2 et que le bruit e est un bruit blanc ou négligeable, puis estime nk.

Pour estimer le retard dans ce système, saisissez :

delayest(ze)
ans =

     5    10

Ce résultat inclut deux nombres, car il y a deux entrées : le retard estimé pour la première entrée est 5 échantillons de données et le retard estimé pour la deuxième entrée est 10 échantillons de données. Comme le pas d’échantillonnage des expérimentations est de 0.5 min, cela correspond à un retard de 2.5 min avant que la première entrée affecte la sortie et à un retard de 5.0 min avant que la deuxième entrée affecte la sortie.

Estimer des retards à l’aide d’autres méthodes

Il existe deux autres méthodes d’estimation du retard dans le système :

  • Tracez le graphique en fonction du temps des données d’entrée et de sortie et lisez la différence de temps entre le premier changement en entrée et le premier changement en sortie. Cette méthode est pratique uniquement pour un système à entrée unique/sortie unique ; dans le cas de systèmes à entrées multiples, il se peut que vous ne puissiez pas savoir quelle entrée a provoqué le changement initial de sortie.

  • Tracez la réponse impulsionnelle des données avec une région de confiance de 1 écart-type. Vous pouvez estimer le retard en utilisant le moment où la réponse impulsionnelle sort pour la première fois de la région de confiance.

Estimer des ordres des modèles à l’aide d’une structure de modèle ARX

Pourquoi estimer l’ordre d'un modèle ?

L’ordre d'un modèle correspond à un ou plusieurs nombres entiers qui définissent la complexité du modèle. En général, l’ordre d'un modèle est lié au nombre de pôles, au nombre de zéros et au retard (temps correspondant au nombre d’échantillons avant que la sortie réponde à l’entrée). La signification spécifique de l’ordre d'un modèle dépend de la structure du modèle.

Pour identifier des modèles boîte noire paramétriques, vous devez spécifier l’ordre du modèle en entrée. Si vous ne connaissez pas l’ordre de votre système, vous pouvez l’estimer.

Après avoir effectué les étapes de cette section, vous obtenez les résultats suivants :

  • Pour la première combinaison entrée/sortie : na=2, nb=2 et le retard nk=5.

  • Pour la deuxième combinaison entrée/sortie : na=1, nb=1 et le retard nk=10.

Ensuite, vous examinez différentes structures de modèle en spécifiant des valeurs d’ordre de modèle qui sont de légères variations autour de ces estimations initiales.

Commandes pour l’estimation de l’ordre du modèle

Dans cette partie du tutoriel, vous allez utiliser struc, arxstruc et selstruc pour estimer et comparer des modèles polynomiaux simples (ARX) pour une plage d’ordres de modèle et de retards, puis sélectionnez les meilleurs ordres en fonction de la qualité du modèle.

La liste suivante décrit les résultats de l’utilisation de chaque commande :

  • struc crée une matrice de combinaisons modèle-ordre pour une plage spécifiée des valeurs de na, nb et nk.

  • arxstruc prend la sortie de struc, estime systématiquement un modèle ARX pour chaque ordre du modèle et compare la sortie du modèle à la sortie mesurée ; arxstruc renvoie la fonction de perte pour chaque modèle, qui est la somme normalisée des erreurs de prédiction au carré.

  • selstruc prend la sortie de arxstruc et ouvre la fenêtre ARX Model Structure Selection qui ressemble à la figure suivante pour vous aider à choisir l’ordre du modèle.

    Vous pouvez utiliser ce graphique pour sélectionner le modèle avec le meilleur ajustement.

    • L’axe horizontal est le nombre total de paramètres — na + nb.

    • L’axe vertical, appelé Unexplained output variance (in %), est la partie de la sortie non expliquée par le modèle, à savoir l’erreur de prédiction du modèle ARX pour le nombre de paramètres indiqué sur l’axe horizontal.

      L’erreur de prédiction est la somme des carrés des différences entre la sortie des données de validation et la sortie prédite du modèle un pas avant.

    • nk est le retard.

    Trois rectangles sont surlignés sur le graphique en vert, bleu et rouge. Chaque couleur indique un type de critère de meilleur ajustement, comme suit :

    • Rouge — Le meilleur ajustement minimise la somme des carrés des différences entre la sortie des données de validation et la sortie du modèle. Ce rectangle indique le meilleur ajustement global.

    • Vert — Le meilleur ajustement qui minimise le critère MDL de Rissanen.

    • Bleu — Le meilleur ajustement qui minimise le critère AIC d’Akaike.

    Dans ce tutoriel, la valeur Unexplained output variance (in %) reste approximativement constante pour le nombre combiné de paramètres de 4 à 20. Une telle constance indique que les performances du modèle ne s’améliorent pas aux ordres supérieurs. Ainsi, les modèles d’ordre inférieur pourraient correspondre tout aussi bien aux données.

    Remarque

    Lorsque vous utilisez les mêmes données pour l’estimation et la validation, utilisez les critères MDL et AIC pour sélectionner les ordres des modèles. Ces critères compensent le surajustement (overfitting) qui résulte de l’utilisation de trop nombreux paramètres. Pour plus d'informations sur ces critères, consultez la page de référence selstruc.

Ordre du modèle pour la première combinaison entrée-sortie

Dans ce tutoriel, le système a deux entrées et une sortie et vous estimez les ordres du modèle pour chaque combinaison entrée/sortie indépendamment. Vous pouvez estimer les retards des deux entrées simultanément ou d’une entrée à la fois.

Il est logique d’essayer les combinaisons d’ordres suivantes pour la première combinaison entrée/sortie :

  • na=2:5

  • nb=1:5

  • nk=5

En effet, les modèles non paramétriques que vous avez estimés dans Estimer des modèles de réponse impulsionnelle montrent que la réponse pour la première combinaison entrée/sortie peut avoir une réponse de deuxième ordre. De même, dans Estimer des retards dans un système à entrées multiples, le retard estimé pour cette combinaison entrée/sortie était égal à 5.

Pour estimer l’ordre du modèle pour la première combinaison entrée/sortie :

  1. Utilisez struc pour créer une matrice des ordres possibles du modèle.

    NN1 = struc(2:5,1:5,5);
  2. Utilisez selstruc pour calculer les fonctions de perte pour les modèles ARX avec les ordres de NN1.

    selstruc(arxstruc(ze(:,:,1),zv(:,:,1),NN1))

    Remarque

    ze(:,:,1) sélectionne la première entrée dans les données.

    Cette commande ouvre la fenêtre ARX Model Structure Selection.

    Remarque

    Le critère MDL de Rissanen et le critère AIC d’Akaike produisent des résultats équivalents et sont tous les deux indiqués par un rectangle bleu sur le graphique.

    Le rectangle rouge représente le meilleur ajustement global, qui se produit pour na=2, nb=3 et nk=5. La différence de hauteur entre les rectangles rouge et bleu est négligeable. Par conséquent, vous pouvez choisir la combinaison de paramètres qui correspond à l’ordre le plus faible et au modèle le plus simple.

  3. Cliquez sur rectangle bleu, puis cliquez sur Select pour choisir la combinaison d’ordres :

    na=2

    nb=2

    nk=5

  4. Pour continuer, appuyez sur n’importe quelle touche lorsque vous êtes dans la fenêtre de commande MATLAB.

Ordre du modèle pour la deuxième combinaison entrée-sortie

Il est logique d’essayer les combinaisons d’ordres suivantes pour la deuxième combinaison entrée/sortie :

  • na=1:3

  • nb=1:3

  • nk=10

En effet, les modèles non paramétriques que vous avez estimés dans Estimer des modèles de réponse impulsionnelle montrent que la réponse pour la deuxième combinaison entrée/sortie peut avoir une réponse de premier ordre. De même, dans Estimer des retards dans un système à entrées multiples, le retard estimé pour cette combinaison entrée/sortie était égal à 10.

Pour estimer l’ordre du modèle pour la deuxième combinaison entrée/sortie :

  1. Utilisez struc pour créer une matrice des ordres possibles du modèle.

    NN2 = struc(1:3,1:3,10);
  2. Utilisez selstruc pour calculer les fonctions de perte pour les modèles ARX avec les ordres de NN2.

    selstruc(arxstruc(ze(:,:,2),zv(:,:,2),NN2))

    Cette commande ouvre la fenêtre ARX Model Structure Selection.

    Remarque

    Le critère AIC d’Akaike et le critère de meilleur ajustement global produisent des résultats équivalents. Les deux sont indiqués par le même rectangle rouge sur le graphique.

    La différence de hauteur entre tous les rectangles est négligeable et tous ces ordres de modèle donnent lieu à des performances de modèle similaires. Par conséquent, vous pouvez choisir la combinaison de paramètres qui correspond à l’ordre le plus faible et au modèle le plus simple.

  3. Cliquez sur le rectangle jaune à l’extrême gauche, puis cliquez sur Select pour choisir l’ordre le plus faible : na=1, nb=1 et nk=10.

  4. Pour continuer, appuyez sur n’importe quelle touche lorsque vous êtes dans la fenêtre de commande MATLAB.

Estimation des fonctions de transfert

Spécification de la structure de la fonction de transfert

Dans cette partie du tutoriel, vous estimez une fonction de transfert à temps continu. Vous devez avoir déjà préparé vos données, comme décrit dans Préparation des données. Vous pouvez utiliser les résultats suivants d’estimation des ordres du modèle pour les spécifier :

  • Pour la première combinaison entrée/sortie, utilisez :

    • Deux pôles, correspondant à na=2 dans le modèle ARX.

    • Un retard de 5, correspondant à nk=5 échantillons (ou 2,5 minutes) dans le modèle ARX.

  • Pour la deuxième combinaison entrée/sortie, utilisez :

    • Un pôle, correspondant à na=1 dans le modèle ARX

    • Un retard de 10, correspondant à nk=10 échantillons (ou 5 minutes) dans le modèle ARX.

Vous pouvez estimer une fonction de transfert de ces ordres à l’aide de la commande tfest. Pour l’estimation, vous pouvez également choisir de visualiser un rapport d’avancement en réglant l’option Display sur on dans l’ensemble d’options créé par la commande tfestOptions.

Opt = tfestOptions('Display','on');

Collectez les ordres du modèle et les retards dans des variables à transmettre à tfest.

np = [2 1];
ioDelay = [2.5 5];

Estimez la fonction de transfert.

mtf = tfest(Ze1,np,[],ioDelay,Opt);

Visualisez les coefficients du modèle.

mtf
mtf =
  From input "ConsumptionRate" to output "Production":
                  -1.382 s + 0.0008305
  exp(-2.5*s) * -------------------------
                s^2 + 1.014 s + 5.444e-12
 
  From input "Current" to output "Production":
                 5.93
  exp(-5*s) * ----------
              s + 0.5858
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: [2 1]   Number of zeros: [1 0]
   Number of free coefficients: 6
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                         
Estimated using TFEST on time domain data "Ze1".
Fit to estimation data: 78.92%                  
FPE: 14.22, MSE: 13.97                          
 

L’affichage du modèle montre un ajustement aux données d’estimation supérieur à 85 %.

Valider le modèle

Dans cette partie du tutoriel, vous créez un graphique qui compare la sortie réelle et la sortie du modèle à l’aide de la commande compare.

compare(Zv1,mtf)

La comparaison montre un ajustement d’environ 60 %.

Analyse des résidus

Utilisez la commande resid pour évaluer les caractéristiques des résidus.

resid(Zv1,mtf)

Les résidus montrent un degré élevé d’autocorrélation. Cela n’est pas surprenant, car le modèle mtf n’a pas de composante décrivant séparément la nature du bruit. Pour modéliser la dynamique d’entrée-sortie mesurée et la dynamique des perturbations, vous devez utiliser une structure de modèle qui contient des éléments décrivant la composante de bruit. Vous pouvez utiliser les commandes bj, ssest et procest qui créent respectivement des modèles avec des structures polynomiales, de représentation d’état et de processus. Ces structures contiennent notamment des éléments rendant compte du comportement de bruit.

Estimer des modèles de processus

Spécifier la structure du modèle de processus

Dans cette partie du tutoriel, vous estimerez une fonction de transfert à temps continu d’ordre inférieur (modèle de processus). System Identification Toolbox supporte des modèles à temps continu avec au maximum trois pôles (qui peuvent contenir des pôles sous-amortis), un zéro, un élément de retard et un intégrateur.

Vous devez avoir déjà préparé vos données, comme décrit dans Préparation des données.

Vous pouvez utiliser les résultats suivants d’estimation des ordres du modèle pour les spécifier :

  • Pour la première combinaison entrée/sortie, utilisez :

    • Deux pôles, correspondant à na=2 dans le modèle ARX.

    • Un retard de 5, correspondant à nk=5 échantillons (ou 2,5 minutes) dans le modèle ARX.

  • Pour la deuxième combinaison entrée/sortie, utilisez :

    • Un pôle, correspondant à na=1 dans le modèle ARX.

    • Un retard de 10, correspondant à nk=10 échantillons (ou 5 minutes) dans le modèle ARX.

Remarque

Comme il n’y a pas de relation entre le nombre de zéros estimé par le modèle ARX à temps discret et son homologue à temps continu, vous n’avez pas d’estimation du nombre de zéros. Dans ce tutoriel, vous pouvez spécifier un zéro pour la première combinaison entrée/sortie et aucun zéro pour la deuxième combinaison.

Utilisez la commande idproc pour créer deux structures de modèle, une pour chaque combinaison entrée/sortie :

midproc0 = idproc({'P2ZUD','P1D'}, 'TimeUnit', 'minutes');

Le cell array {'P2ZUD','P1D'} spécifie la structure du modèle pour chaque combinaison entrée/sortie :

  • 'P2ZUD' représente une fonction de transfert avec deux pôles (P2), un zéro (Z), des pôles sous-amortis (conjugué complexe) (U) et un retard (D).

  • 'P1D' représente une fonction de transfert avec un pôle (P1) et un retard (D).

L’exemple utilise également le paramètre TimeUnit pour spécifier l’unité de temps employée.

Visualisation de la structure du modèle et des valeurs des paramètres

Visualisez les deux modèles résultants.

midproc0
midproc0 =

Process model with 2 inputs: y = G11(s)u1 + G12(s)u2
  From input 1 to output 1:                         
                       1+Tz*s                       
  G11(s) = Kp * ---------------------- * exp(-Td*s) 
                1+2*Zeta*Tw*s+(Tw*s)^2              
                                                    
         Kp = NaN                                   
         Tw = NaN                                   
       Zeta = NaN                                   
         Td = NaN                                   
         Tz = NaN                                   
                                                    
  From input 2 to output 1:                         
             Kp                                     
  G12(s) = ---------- * exp(-Td*s)                  
            1+Tp1*s                                 
                                                    
        Kp = NaN                                    
       Tp1 = NaN                                    
        Td = NaN                                    
                                                    
Parameterization:
    {'P2DUZ'}    {'P1D'}
   Number of free coefficients: 8
   Use "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                         
Created by direct construction or transformation. Not estimated.
 

Les valeurs des paramètres sont définies par NaN car elles ne sont pas encore estimées.

Spécification des estimations initiales pour les retards

Réglez la propriété de retard de l’objet de modèle sur 2.5 min et 5 min pour chaque paire d’entrée/sortie comme estimations initiales. Définissez également une limite supérieure pour le retard, car les estimations disponibles sont bonnes.

midproc0.Structure(1,1).Td.Value = 2.5;
midproc0.Structure(1,2).Td.Value = 5;
midproc0.Structure(1,1).Td.Maximum = 3;
midproc0.Structure(1,2).Td.Maximum = 7;

Remarque

Lors du réglage des contraintes de retard, vous devez spécifier les retards en unités de temps réelles (minutes, dans ce cas) et non en nombre d’échantillons.

Estimation des paramètres du modèle à l’aide de procest

procest est un estimateur itératif de modèles de processus, ce qui signifie qu’il utilise un algorithme des moindres carrés non linéaire itératif pour minimiser une fonction de coût. La fonction de coût est la somme pondérée des carrés des erreurs.

Selon ses arguments, procest estime différents modèles polynomiaux de type boîte noire. Vous pouvez utiliser procest, par exemple, pour estimer les paramètres des structures de modèle linéaire à fonction de transfert en temps continu, de modèle de représentation d’état ou polynomial. Pour utiliser procest, vous devez fournir une structure de modèle avec des paramètres inconnus et les données d’estimation comme arguments d’entrée.

Pour cette partie du tutoriel, vous devez avoir déjà défini la structure du modèle, comme décrit dans Spécifier la structure du modèle de processus. Utilisez midproc0 comme structure de modèle et Ze1 comme données d’estimation :

midproc = procest(Ze1,midproc0);
present(midproc)
                                                                  
midproc =                                                         
                                                                  
Process model with 2 inputs: y = G11(s)u1 + G12(s)u2              
  From input "ConsumptionRate" to output "Production":            
                       1+Tz*s                                     
  G11(s) = Kp * ---------------------- * exp(-Td*s)               
                1+2*Zeta*Tw*s+(Tw*s)^2                            
                                                                  
         Kp = -1.1807 +/- 0.29986                                 
         Tw = 1.6437 +/- 714.6                                    
       Zeta = 16.036 +/- 6958.9                                   
         Td = 2.426 +/- 64.276                                    
         Tz = -109.19 +/- 63.733                                  
                                                                  
  From input "Current" to output "Production":                    
             Kp                                                   
  G12(s) = ---------- * exp(-Td*s)                                
            1+Tp1*s                                               
                                                                  
        Kp = 10.264 +/- 0.048404                                  
       Tp1 = 2.049 +/- 0.054901                                   
        Td = 4.9175 +/- 0.034374                                  
                                                                  
Parameterization:                                                 
    {'P2DUZ'}    {'P1D'}                                          
   Number of free coefficients: 8                                 
   Use "getpvec", "getcov" for parameters and their uncertainties.
                                                                  
Status:                                                           
Termination condition: Maximum number of iterations reached..     
Number of iterations: 20, Number of function evaluations: 279     
                                                                  
Estimated using PROCEST on time domain data "Ze1".                
Fit to estimation data: 86.2%                                     
FPE: 6.081, MSE: 5.984                                            
More information in model's "Report" property.                    
                                                                  

Contrairement aux modèles polynomiaux à temps discret, les modèles à temps continu vous permettent d’estimer les retards. Dans ce cas, les valeurs de retard estimées sont différentes des estimations initiales spécifiées respectivement à 2.5 et 5. Les grandes incertitudes des valeurs estimées des paramètres de G_1(s) indiquent que la dynamique entre l’entrée 1 et la sortie n’est pas correctement représentée.

Pour en savoir plus sur l’estimation des modèles, consultez Process Models.

Valider le modèle

Dans cette partie du tutoriel, vous créez un graphique qui compare la sortie réelle et la sortie du modèle.

compare(Zv1,midproc)

Le graphique montre que la sortie du modèle correspond raisonnablement à la sortie mesurée : il y a une correspondance de 60 % entre le modèle et les données de validation.

Utilisez resid pour effectuer une analyse des résidus.

resid(Zv1,midproc)

La corrélation croisée entre la deuxième entrée et les erreurs résiduelles est significative. Le graphique d’autocorrélation montre des valeurs à l’extérieur de la région de confiance et indique que les résidus sont corrélés.

Remplacez l’algorithme d’estimation des paramètres itérative par un algorithme de Levenberg-Marquardt.

Opt = procestOptions;
Opt.SearchMethod = 'lm';
midproc1 = procest(Ze1,midproc,Opt);

Pour essayer d’améliorer les résultats de la simulation, vous pouvez modifier les propriétés de l’algorithme ou spécifier des estimations de paramètres initiales, puis exécuter une nouvelle fois l'estimation. L’ajout d’un modèle de bruit peut améliorer les résultats de la prédiction, mais pas forcément les résultats de la simulation.

Estimation d’un modèle de processus avec modèle de bruit

Cette partie du tutoriel montre comment estimer un modèle de processus et inclure un modèle de bruit dans l’estimation. L’ajout d’un modèle de bruit améliore généralement les résultats de la prédiction, mais pas forcément les résultats de la simulation.

Utilisez les commandes suivantes pour spécifier un bruit ARMA de premier ordre :

Opt = procestOptions;
Opt.DisturbanceModel = 'ARMA1';
midproc2 = procest(Ze1,midproc0,Opt);

Vous pouvez saisir 'dist' au lieu de 'DisturbanceModel'. Les noms des propriétés ne sont pas sensibles à la casse et vous devez seulement inclure la partie du nom qui identifie de manière unique la propriété.

Comparez le modèle résultant aux données mesurées.

compare(Zv1,midproc2)

Le graphique montre que la sortie du modèle correspond raisonnablement à la sortie des données de validation.

Effectuer une analyse résiduelle.

resid(Zv1,midproc2)

Le graphique des résidus montre que les valeurs d’autocorrélation sont à l’intérieur des limites de confiance. Ainsi, l’ajout d’un modèle de bruit produit des résidus non corrélés. Cela indique que le modèle est plus précis.

Estimation de modèles polynomiaux de type boîte noire

Ordres des modèles pour l’estimation de modèles polynomiaux

Dans cette partie du tutoriel, vous allez estimer différents types de modèles polynomiaux d’entrée-sortie de type boîte noire.

Vous devez avoir déjà préparé vos données, comme décrit dans Préparation des données.

Vous pouvez utiliser les résultats suivants d'estimation des ordres de modèle pour spécifier les ordres du modèle polynomial :

  • Pour la première combinaison entrée/sortie, utilisez :

    • Deux pôles, correspondant à na=2 dans le modèle ARX.

    • Un zéro, correspondant à nb=2 dans le modèle ARX.

    • Un retard de 5, correspondant à nk=5 échantillons (ou 2,5 minutes) dans le modèle ARX.

  • Pour la deuxième combinaison entrée/sortie, utilisez :

    • Un pôle, correspondant à na=1 dans le modèle ARX.

    • Aucun zéro, correspondant à nb=1 dans le modèle ARX.

    • Un retard de 10, correspondant à nk=10 échantillons (ou 5 minutes) dans le modèle ARX.

Estimation d’un modèle ARX linéaire

À propos des modèles ARX.  Pour un système à entrée unique/sortie unique (SISO), la structure de modèle ARX est :

y(t)+a1y(t1)++anay(tna)=         b1u(tnk)++bnbu(tnknb+1)+e(t)

y(t) représente la sortie à l’instant t, u(t) représente l’entrée à l’instant t, na est le nombre de pôles, nb est le nombre de zéros plus 1, nk est le nombre d’échantillons avant que l’entrée affecte la sortie du système (appelé retard ou temps mort du modèle) et e(t) est la perturbation du bruit blanc.

La structure de modèle ARX ne fait pas de distinction entre les pôles des chemins entrée/sortie individuels : la division de l’équation ARX par A, qui contient les pôles, montre que A apparaît dans le dénominateur des deux entrées. Par conséquent, vous pouvez définir na comme la somme des pôles de chaque paire entrée/sortie, ce qui est égal à 3 dans ce cas.

System Identification Toolbox estime les paramètres a1an et b1bn à l’aide des données et des ordres du modèle que vous avez spécifiés.

Estimation de modèles ARX avec arx.  Utilisez arx pour calculer les coefficients polynomiaux avec la méthode non itérative rapide arx :

marx = arx(Ze1,'na',3,'nb',[2 1],'nk',[5 10]);
present(marx) % Displays model parameters
              % with uncertainty information
                                                                              
marx =                                                                        
Discrete-time ARX model: A(z)y(t) = B(z)u(t) + e(t)                           
                                                                              
  A(z) = 1 - 1.027 (+/- 0.02907) z^-1 + 0.1678 (+/- 0.042) z^-2 + 0.01289 (   
                                                         +/- 0.02583) z^-3    
                                                                              
  B1(z) = 1.86 (+/- 0.189) z^-5 - 1.608 (+/- 0.1888) z^-6                     
                                                                              
  B2(z) = 1.612 (+/- 0.07392) z^-10                                           
                                                                              
Sample time: 0.5 minutes                                                      
                                                                              
Parameterization:                                                             
   Polynomial orders:   na=3   nb=[2 1]   nk=[5 10]                           
   Number of free coefficients: 6                                             
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
                                                                              
Status:                                                                       
Estimated using ARX on time domain data "Ze1".                                
Fit to estimation data: 90.7% (prediction focus)                              
FPE: 2.768, MSE: 2.719                                                        
More information in model's "Report" property.                                
                                                                              

MATLAB estime les polynômes A, B1 et B2. L’incertitude de chaque paramètre du modèle est calculée à 1 écart-type et apparaît entre parenthèses à côté de chaque valeur de paramètre.

Sinon, vous pouvez utiliser la syntaxe abrégée suivante et spécifier les ordres du modèle sous forme d’un vecteur unique :

marx = arx(Ze1,[3 2 1 5 10]);

Accès aux données du modèle.  Le modèle que vous avez estimé, marx, est un objet idpoly à temps discret. Pour obtenir les propriétés de cet objet de modèle, vous pouvez utiliser la fonction get :

get(marx)
                 A: [1 -1.0267 0.1678 0.0129]
                 B: {[0 0 0 0 0 1.8599 -1.6084]  [0 0 0 0 0 0 0 0 0 0 1.6118]}
                 C: 1
                 D: 1
                 F: {[1]  [1]}
    IntegrateNoise: 0
          Variable: 'z^-1'
           IODelay: [0 0]
         Structure: [1x1 pmodel.polynomial]
     NoiseVariance: 2.7436
        InputDelay: [2x1 double]
       OutputDelay: 0
         InputName: {2x1 cell}
         InputUnit: {2x1 cell}
        InputGroup: [1x1 struct]
        OutputName: {'Production'}
        OutputUnit: {'mg/min'}
       OutputGroup: [1x1 struct]
             Notes: [0x1 string]
          UserData: []
              Name: ''
                Ts: 0.5000
          TimeUnit: 'minutes'
      SamplingGrid: [1x1 struct]
            Report: [1x1 idresults.arx]

Vous pouvez accéder aux informations enregistrées par ces propriétés à l’aide de la notation par points. Par exemple, vous pouvez calculer les pôles discrets du modèle en trouvant les racines du polynôme A.

marx_poles = roots(marx.a)
marx_poles =

    0.7953
    0.2877
   -0.0564

Dans ce cas, vous accédez au polynôme A à l’aide de marx.a.

Le modèle marx décrit la dynamique du système avec trois pôles discrets.

Conseil

Vous pouvez également utiliser pole pour calculer directement les pôles d’un modèle.

En savoir plus.  Pour en savoir plus sur l’estimation des modèles polynomiaux, consultez Input-Output Polynomial Models.

Pour plus d'informations sur l’accès aux données du modèle, consultez Data Extraction.

Estimation des modèles de représentation d'état

À propos des modèles de représentation d'état.  La structure générale des modèles de représentation d’état est :

dxdt=Ax(t)+Bu(t)+Ke(t)y(t)=Cx(t)+Du(t)+e(t)

y(t) représente la sortie à l’instant t, u(t) représente l’entrée à l’instant t, x(t) est le vecteur d’état à l’instant t et e(t) est la perturbation du bruit blanc.

Vous devez spécifier un nombre entier unique comme ordre du modèle (dimension du vecteur d’état) pour estimer un modèle de représentation d’état. Par défaut, le retard est égal à 1.

System Identification Toolbox estime les matrices de représentation d'état A, B, C, D et K à l’aide de l’ordre du modèle et des données que vous spécifiez.

La structure de modèle de représentation d'état est un bon choix pour une estimation rapide, car elle contient seulement deux paramètres : n est le nombre de pôles (la taille de la matrice A) et nk est le retard.

Estimation des modèles de représentation d’état avec n4sid.  Utilisez la commande n4sid pour spécifier une plage d’ordres des modèles et évaluer les performances de plusieurs modèles de représentation d’état (ordres 2 à 8) :

mn4sid = n4sid(Ze1,2:8,'InputDelay',[4 9]);

Cette commande utilise la méthode non itérative rapide (méthode de sous-espace) et ouvre le graphique suivant. Vous utilisez ce graphique pour décider quels états contribuent de manière significative au comportement entrée/sortie et quels états y contribuent le moins.

L’axe vertical est une mesure relative de la manière dont chaque état contribue au comportement d'entrée/sortie du modèle (log de valeurs singulières de la matrice de covariance). L’axe horizontal correspond à l’ordre du modèle n. Ce graphique recommande n=3, indiqué par un rectangle rouge.

Le champ Chosen Order affiche par défaut l’ordre recommandé du modèle, 3 dans ce cas. Vous pouvez changer la sélection de l’ordre avec le menu déroulant Chosen Order. Appliquez la valeur dans le champ Chosen Order et fermez la fenêtre de sélection de l’ordre en cliquant sur Apply.

Par défaut, n4sid utilise un paramétrage libre de la forme de la représentation d'état. Pour estimer une forme canonique à la place, réglez la valeur de la propriété SSParameterization sur 'Canonical'. Vous pouvez également spécifier le retard entre l’entrée et la sortie (en échantillons) avec la propriété 'InputDelay'.

mCanonical = n4sid(Ze1,3,'SSParameterization','canonical','InputDelay',[4 9]);
present(mCanonical);  %  Display model properties
                                                                              
mCanonical =                                                                  
  Discrete-time identified state-space model:                                 
    x(t+Ts) = A x(t) + B u(t) + K e(t)                                        
       y(t) = C x(t) + D u(t) + e(t)                                          
                                                                              
  A =                                                                         
                       x1                  x2                  x3             
   x1                   0                   1                   0             
   x2                   0                   0                   1             
   x3  0.0737 +/- 0.05919  -0.6093 +/- 0.1626    1.446 +/- 0.1287             
                                                                              
  B =                                                                         
             ConsumptionR             Current                                 
   x1   1.844 +/-   0.175   0.5633 +/-  0.122                                 
   x2   1.063 +/-  0.1673    2.308 +/- 0.1222                                 
   x3  0.2779 +/- 0.09615    1.878 +/- 0.1058                                 
                                                                              
  C =                                                                         
               x1  x2  x3                                                     
   Production   1   0   0                                                     
                                                                              
  D =                                                                         
               ConsumptionR       Current                                     
   Production             0             0                                     
                                                                              
  K =                                                                         
               Production                                                     
   x1  0.8674 +/- 0.03169                                                     
   x2  0.6849 +/- 0.04145                                                     
   x3  0.5105 +/- 0.04352                                                     
                                                                              
  Input delays (sampling periods): 4  9                                       
                                                                              
Sample time: 0.5 minutes                                                      
                                                                              
Parameterization:                                                             
   CANONICAL form with indices: 3.                                            
   Feedthrough: none                                                          
   Disturbance component: estimate                                            
   Number of free coefficients: 12                                            
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
                                                                              
Status:                                                                       
Estimated using N4SID on time domain data "Ze1".                              
Fit to estimation data: 91.39% (prediction focus)                             
FPE: 2.402, MSE: 2.331                                                        
More information in model's "Report" property.                                
                                                                              

Remarque

mn4sid et mCanonical sont des modèles à temps discret. Pour estimer un modèle à temps continu, réglez la propriété 'Ts' à 0 dans la commande d’estimation ou utilisez la commande ssest :

mCT1 = n4sid(Ze1, 3, 'Ts', 0, 'InputDelay', [2.5 5])
mCT2 = ssest(Ze1, 3,'InputDelay', [2.5 5])

En savoir plus.  Pour en savoir plus sur l’estimation de modèles de représentation d’état, consultez State-Space Models.

Estimation d’un modèle Box-Jenkins

À propos des modèles Box-Jenkins.  La structure générale de Box-Jenkins (BJ) est :

y(t)=i=1nuBi(q)Fi(q)ui(tnki)+C(q)D(q)e(t)

Pour estimer un modèle BJ, vous devez spécifier les paramètres nb, nf, nc, nd et nk.

Tandis que la structure de modèle ARX ne fait pas de distinction entre les pôles des chemins d'entrée/sortie individuels, le modèle BJ offre davantage de flexibilité en modélisant les pôles et les zéros des perturbations séparément des pôles et des zéros de la dynamique du système.

Estimation d’un modèle BJ avec polyest.  Vous pouvez utiliser polyest pour estimer le modèle BJ. polyest est une méthode itérative et a la syntaxe générale suivante :

polyest(data,[na nb nc nd nf nk]);

Pour estimer le modèle BJ, saisissez :

na = 0;
nb = [ 2 1 ];
nc = 1;
nd = 1;
nf = [ 2 1 ];
nk = [ 5 10];
mbj = polyest(Ze1,[na nb nc nd nf nk]);

Cette commande spécifie nf=2, nb=2, nk=5 pour la première entrée et nf=nb=1 et nk=10 pour la deuxième entrée.

Affichez les informations du modèle.

present(mbj)
                                                                              
mbj =                                                                         
Discrete-time BJ model: y(t) = [B(z)/F(z)]u(t) + [C(z)/D(z)]e(t)              
  B1(z) = 1.823 (+/- 0.1792) z^-5 - 1.315 (+/- 0.2367) z^-6                   
                                                                              
  B2(z) = 1.791 (+/- 0.06431) z^-10                                           
                                                                              
  C(z) = 1 + 0.1068 (+/- 0.04009) z^-1                                        
                                                                              
  D(z) = 1 - 0.7452 (+/- 0.02694) z^-1                                        
                                                                              
  F1(z) = 1 - 1.321 (+/- 0.06936) z^-1 + 0.5911 (+/- 0.05514) z^-2            
                                                                              
  F2(z) = 1 - 0.8314 (+/- 0.006441) z^-1                                      
                                                                              
Sample time: 0.5 minutes                                                      
                                                                              
Parameterization:                                                             
   Polynomial orders:   nb=[2 1]   nc=1   nd=1   nf=[2 1]                     
   nk=[5 10]                                                                  
   Number of free coefficients: 8                                             
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
                                                                              
Status:                                                                       
Termination condition: Near (local) minimum, (norm(g) < tol)..                
Number of iterations: 6, Number of function evaluations: 13                   
                                                                              
Estimated using POLYEST on time domain data "Ze1".                            
Fit to estimation data: 90.75% (prediction focus)                             
FPE: 2.733, MSE: 2.689                                                        
More information in model's "Report" property.                                
                                                                              

L’incertitude de chaque paramètre du modèle est calculée à 1 écart-type et apparaît entre parenthèses à côté de chaque valeur de paramètre.

Les polynômes C et D donnent respectivement le numérateur et le dénominateur du modèle de bruit.

Conseil

Sinon, vous pouvez utiliser la syntaxe abrégée suivante qui spécifie les ordres sous forme d’un vecteur unique :

mbj = bj(Ze1,[2 1 1 1 2 1 5 10]);

bj est une version de polyest qui estime spécifiquement la structure de modèle BJ.

En savoir plus.  Pour en savoir plus sur l’identification de modèles polynomiaux d’entrée-sortie comme BJ, consultez Input-Output Polynomial Models.

Comparaison de la sortie du modèle et de la sortie mesurée

Comparez la sortie des modèles ARX, Box-Jenkins et de représentation d’état à la sortie mesurée.

compare(Zv1,marx,mn4sid,mbj)

compare trace la sortie mesurée du jeu de données de validation par rapport à la sortie simulée à partir des modèles. Les données d’entrée du jeu de données de validation servent d’entrée aux modèles.

Effectuez une analyse des résidus sur les modèles ARX, Box-Jenkins et de représentation d’état.

resid(Zv1,marx,mn4sid,mbj)

Les trois modèles simulent aussi bien la sortie et ont des résidus non corrélés. Par conséquent, choisissez le modèle ARX, car il est le plus simple des trois modèles polynomiaux d’entrée-sortie et rend compte de la dynamique du processus de manière adéquate.

Simulation et prédiction de la sortie du modèle

Simulation de la sortie du modèle

Dans cette partie du tutoriel, vous allez estimer la sortie du modèle. Vous devez avoir déjà créé le modèle à temps continu midproc2, comme décrit dans Estimer des modèles de processus.

La simulation de la sortie du modèle exige les informations suivantes :

  • Valeurs d’entrée du modèle

  • Conditions initiales de la simulation (également appelées états initiaux)

Par exemple, les commandes suivantes utilisent les commandes iddata et idinput pour construire un jeu de données d’entrée et utilisent sim pour simuler la sortie du modèle :

% Create input for simulation
U = iddata([],idinput([200 2]),'Ts',0.5,'TimeUnit','min');
% Simulate the response setting initial conditions equal to zero
ysim_1 = sim(midproc2,U);

Pour maximiser l'ajustement ou correspondance entre la réponse simulée d’un modèle à la sortie mesurée pour la même entrée, vous pouvez calculer les conditions initiales correspondant aux données mesurées. Les conditions initiales du meilleur ajustement peuvent être obtenues avec findstates sur la version de représentation d’état du modèle estimé. Les commandes suivantes estiment les états initiaux X0est à partir du jeu de données Zv1 :

% State-space version of the model needed for computing initial states
midproc2_ss = idss(midproc2);
X0est = findstates(midproc2_ss,Zv1);

Ensuite, simulez le modèle à l’aide des états initiaux estimés à partir des données.

% Simulation input
Usim = Zv1(:,[],:);
Opt = simOptions('InitialCondition',X0est);
ysim_2 = sim(midproc2_ss,Usim,Opt);

Comparez la sortie simulée et la sortie mesurée sur un graphique.

figure
plot([ysim_2.y, Zv1.y])
legend({'model output','measured'})
xlabel('time'), ylabel('Output')

Prédiction de la future sortie

De nombreuses applications de design de système de contrôle exigent que vous prédisiez les futures sorties d’un système dynamique en utilisant des données d’entrée/sortie passées.

Par exemple, utilisez predict pour prédire la réponse du modèle cinq pas à l’avance :

predict(midproc2,Ze1,5)

Comparez les valeurs de sortie prédites aux valeurs de sortie mesurées. Le troisième argument de compare spécifie une prédiction cinq pas à l’avance. Lorsque vous ne spécifiez pas de troisième argument, compare suppose un horizon de prédiction infini et simule la sortie du modèle à la place.

compare(Ze1,midproc2,5)

Utilisez pe pour calculer l’erreur de prédiction Err entre la sortie prédite de midproc2 et la sortie mesurée. Ensuite, tracez le spectre d’erreur avec la commande spectrum.

[Err] = pe(midproc2,Zv1);
spectrum(spa(Err,[],logspace(-2,2,200)))

Les erreurs de prédiction sont tracées avec un intervalle de confiance de 1 écart-type. Les erreurs sont supérieures aux hautes fréquences en raison de la nature haute fréquence des perturbations.