TD Automatique Linéaire 1A

TD3 : Identification / Système de référence / Correcteur

Etude d’un Four tubulaire

La distillation est une opération de séparation

de composés en fonction de leur différence
de volatilité. Une colonne se compose d’un
ensemble de plateau au niveau desquels se
croisent les produits « légers » montant et les
produits « lourds » descendant. Dans le cas
de la distillation du pétrole brut, appelé
distillation fractionnée, on soutire au niveau
de plusieurs plateaux où se concentrent des
produits valorisables (diesel, essence, huiles
de moteur…). L’efficacité d’une colonne à
distiller dépend de la maîtrise du gradient de
température et de pression dans le corps de la
colonne. Un des éléments critique pour cette
maîtrise est la température du brut entrant
dans la colonne. Pour optimiser le fonctionnement, le pétrole brut passe au préalable dans une
unité de préchauffage (un four tubulaire) avec lequel on va garantir une température constante
au pétrole avant son injection dans la colonne Figure 1. En fonction de l’origine du pétrole brut,
cette température d’injection va être comprise entre 350 et 450°.

Figure 1 : schéma de fonctionnement d’un four tubulaire

On se propose ici de caractériser le comportement d’un four en s’appuyant sur le schéma bloc
Figure 2.

1
 TD Automatique Linéaire 1A

Figure 2 : schéma entrées-sorties du four

𝑢(𝑡) Commande de la vanne continu d’alimentation du brûleur (en V, 𝑢 ∈ [−10V, +10V])

𝑑(𝑡) Variation de débit entrant du pétrole à chauffer (en m3/s)
𝑦(𝑡) Mesure de la température en sortie du four (en V, conversion 100°  1V)
𝐹(𝑠) Fonction de transfert du procédé par rapport à la commande T𝑢→𝑦
𝐺(𝑠) Fonction de transfert du procédé par rapport à la variation du débit T𝑑→𝑦

Dans la suite de l’étude, on travaille autour d’un point de fonctionnement. Toutes les variations
d’entrée et les observations de la réaction du système se font autour de ce point de
fonctionnement (𝑢𝑟𝑒𝑒𝑙 = 𝑢0 + 𝑢 et 𝑦𝑟𝑒𝑒𝑙 = 𝑦0 + 𝑦). Les deux transferts 𝑇𝑢→𝑦 et 𝑇𝑑→𝑦 autour
de ce point de fonctionnement sont considérés comme linéaires temps invariant (LTI).

Modélisation

On dispose d’un enregistrement de la réponse du procédé en boucle ouverte à un échelon de

commande 𝑢 unitaire avec un débit constant de pétrole (soit 𝑑 = 0)

Figure 3 : réponse à un échelon unitaire de commande

Question 1 : A partir de la Figure 3, justifier le choix d’une structure d’ordre 2, stable et sans
zéro, pour représenter l’évolution de la température en fonction de la commande de la vanne.

2
 TD Automatique Linéaire 1A

La figure présente des éléments caractéristiques dans son évolution :

• Tangente horizontale à l’origine  ordre 𝑛 = 2 ou plus
• Réponse convergente  système stable, pôles à partie réelle strictement négative
• Croissance monotone  pôles réels
• Nous n’observons pas l’influence de zéros : la tangente à l’origine est nulle, ce qui veut
dire que le degré relatif 𝑛 − 𝑚 ≥ 2 , donc nombre de zéros 𝑚 = 0. En plus, il n’y a pas
de dépassement ni de réponse indicielle « inversée » qui parte d’abord en négative.

On représente T𝑢→𝑦 par

𝑘
∀𝑠 ∈ ℂ, 𝐺(𝑠) = Eq. 1
(𝜏1 𝑠 + 1)(𝜏2 𝑠 + 1)
Avec 𝜏1 et 𝜏2 des constantes de temps

Question 2 : A partir de la Figure 3, réaliser l’identification graphique des paramètres de la

fonction de transfert 𝐺 défini Eq. 1.

Pour faire l’identification, on va commencer par passer par une représentation du type
𝜔𝑛2
∀𝑠 ∈ ℂ, 𝐺(𝑠) = 𝑘 2
𝑠 + 2𝜉𝜔𝑛 𝑠 + 𝜔𝑛2
Dont les pôles sont définis par
𝑝1 = −𝜉𝜔𝑛 + 𝑖𝜔𝑛 √1 − 𝜉 2 𝑝2 = −𝜉𝜔𝑛 − 𝑖𝜔𝑛 √1 − 𝜉 2
On peut relever le gain statique :
𝑦∞
𝑦∞ = 𝑘𝑢0 → 𝑘 = = 0.5
𝑢0

La réponse étant monotone, on a deux racines réelles, potentiellement distinctes, qui guide la
réponse indicielle du système. Cela se traduit dans l’écriture de la fonction de transfert par un
coefficient 𝜉 ≥ 1. On fixe alors arbitrairement le coefficient d’amortissement 𝜉 à 1, à l’aide de
l’abaque de temps de réponse réduit, on obtient :
5
𝑡𝑟95% 𝜔𝑛 = 5 → 𝜔𝑛 = = 0.01 𝑟𝑎𝑑. 𝑠 −1
500
On a donc deux racines réelles :
𝑝1 = 𝑝2 = −𝜉𝜔𝑛 = −0.01

Les pôles étant réels, ils sont l’inverse de constantes de temps

0.0001 1
∀𝑠 ∈ ℂ, 𝐺(𝑠) = 0.5 = 0.5
(𝑠 + 0.01)2 (100𝑠 + 1)2
De fait
𝜏1 = 𝜏2 = 100

3
 TD Automatique Linéaire 1A

Vous disposez d’un algorithme itératif d’identification permettant d’obtenir des paramètres plus
précis pour le modèle1 (TD3_identification.m et Gradient_Hessien.m).

Question 3 : Après avoir analysé le code « TD3_identification.m », initialiser l’algorithme à

l’aide de vos résultats précédents (voir la première ligne avec le commentaire « à compléter »)
et déterminer les paramètres du modèle en utilisant les 2 options proposées (variable « Choix »).
Refaites la manipulation avec une initialisation quelconque (2 constantes de temps très
distinctes, erreur de gain statique). Qu’observe-t-on ?

On constate que la méthode du gradient a « plus de mal » à converger vers une solution
satisfaisante, alors que la méthode de Newton-Raphson donne systématiquement de bon résultat
quel que soit le choix de l’initialisation. Il semble que la méthode du gradient soit plus sensible
à l’initialisation.

L’algorithme nous renvoie avec la seconde méthode le modèle suivant :

0.4996 0.5
∀𝑠 ∈ ℂ, 𝐺(𝑠) = ≈
9400 𝑠 2 + 199.7 𝑠 + 1 9400 𝑠 2 + 200 𝑠 + 1

On retiendra
0.5
∀𝑠 ∈ ℂ, 𝐺(𝑠) =
(124.5 𝑠 + 1)(75.5 𝑠 + 1)

On dispose d’un enregistrement de la réponse du procédé en boucle ouverte lors d'une variation
de débit en échelon de 0.01 m3/s (soit 𝑢 = 0)

1
Le lecteur intéressé pourra trouver en annexe une courte explication des opérations effectuées dans le programme
Gradient_Hessien.m
4
 TD Automatique Linéaire 1A

Figure 4 : réponse à un échelon de perturbation de 0.01 m3/s

Question 4 : A partir de la Figure 4, proposer une fonction de transfert 𝐺1 (𝑠) et réaliser

l’identification graphique de ses paramètres.

La figure présente des éléments caractéristiques dans son évolution :

• Tangente oblique à l’origine  ordre 𝑛 = 1
• Réponse convergente  système stable

On est sur une structure simple du type :

𝑘
∀𝑠 ∈ ℂ, 𝐺1 (𝑠) =
𝜏𝑠 + 1
On peut relever le gain statique :
𝑦∞ −0.5
𝑦∞ = 𝑘𝑢0 → 𝑘 = = = −50
𝑢0 0.01

La constante de temps est déterminée à partir du temps de réponse à 95% soit

𝑡𝑟95% 270
𝜏= ≈ = 90 𝑠𝑒𝑐
3 3
Soit
−50
∀𝑠 ∈ ℂ, 𝐺1 (𝑠) =
90 𝑠 + 1

Question 5 : En utilisant le même code « TD3_identification.m » réaliser l’optimisation du

transfert 𝑮𝟏 . Initialiser l’algorithme à l’aide de vos résultats précédents et déterminer les
paramètres du modèle.

5
 TD Automatique Linéaire 1A

On obtient
−50
∀𝑠 ∈ ℂ, 𝐺1 (𝑠) =
92.3 𝑠 + 1

Système de référence

On souhaite que la température du pétrole suive la température de consigne pour des variations
en échelon de 0.5V (soit 50°C), sans erreur, avec un dépassement inférieur ou égal à 10% et
avec un temps de réponse à 5% de 120s au maximum. De plus, elle doit rester régulée en régime
stationnaire malgré d'éventuelles variations de débit du liquide de l’ordre 0.006 m3/s. Ce débit
est mesuré.

Question 6 : Proposer un système de référence que l’on notera 𝑀𝑟𝑒𝑓 dont la réponse indicielle
respecte le cahier des charges en poursuite.

Définition de 𝑀𝑟𝑒𝑓
• Boucle fermée stable
• Suivi sans erreur de la consigne  𝑀𝑟𝑒𝑓 (0) = 1
• 𝑂𝑟𝑑𝑟𝑒(𝐺) = 2  𝑂𝑟𝑑𝑟𝑒(𝑀𝑟𝑒𝑓 ) = 2
• Dépassement max 𝐷1% = 10%  𝜉𝑟𝑒𝑓 = 0.6
1
• 𝑡𝑟95% = 120  𝑡𝑟95% 𝜔𝑟𝑒𝑓 = 𝛼(𝜉) = 5  𝜔𝑟𝑒𝑓 = 24 = 0.0417

On a alors
1 1
∀𝑠 ∈ ℂ, 𝑀𝑟𝑒𝑓 (𝑠) = =
𝑠2 𝜉 576 𝑠2 + 28.8 𝑠 + 1
2 + 2𝜔 𝑠 + 1
𝜔𝑟𝑒𝑓 𝑛

On peut en déduire
1
∀𝑠 ∈ ℂ, 𝐿𝑟𝑒𝑓 (𝑠) =
576 𝑠2 + 28.8 𝑠

Synthèse du correcteur

On décide de conserver l’instrumentation déjà présente sur le site qui permet d’implémenter un
correcteur PID mixte et un bloc d’anticipation (feed-forward) conformément à la Figure 5.

6
 TD Automatique Linéaire 1A

Figure 5 : Architecture de correction retenue

Question 7 : Le choix pour 𝐾 d’un correcteur proportionnel permettra-t-il de satisfaire aux

exigences du cahier des charges ?

Un correcteur proportionnel ne pourra pas convenir car le cahier des charges précises que l’on
ne souhaite pas d’écart statique en réponse à une variation de consigne de type échelon. Dans
ce cas, on sait que le produit de la caractéristique du correcteur 𝐾 par celle du système 𝐺 doit
porter une intégration. Comme 𝐺 n’en a pas, elle doit nécessairement être dans le correcteur 𝐾.
Ceci impose donc une préspécification au correcteur pour assurer le régime permanent désiré.

Question 8 : Un correcteur PID mixte est finalement retenu. Déterminer les paramètres du
correcteur en utilisant un calcul direct à partir du système de référence.

La synthèse par calcul direct peut être appliquée puisque les ordres de 𝐺 et de 𝑀𝑟𝑒𝑓 et de fait
𝐿𝑟𝑒𝑓 sont identiques. Le correcteur 𝐾 est alors défini par :

𝐿𝑟𝑒𝑓 (𝑠) 9400 𝑠 2 + 200 𝑠 + 1 1 9400 𝑠 2 + 200 𝑠 + 1

∀𝑠 ∈ ℂ, 𝐾(𝑠) = = =2
𝐺(𝑠) 0.5 576 𝑠 2 + 28.8 𝑠 576 𝑠 2 + 28.8 𝑠

Le correcteur obtenu est un PID série

(𝑠 + 0.01324)(𝑠 + 0.008032) 2 (124.5 𝑠 + 1)(75.5 𝑠 + 1)

∀𝑠 ∈ ℂ, 𝐾(𝑠) = 32.63 =
𝑠(𝑠 + 0.05) 28.8 𝑠 (20𝑠 + 1)

Les formes des PID série, mixte et parallèle sont formellement équivalentes. La forme mixte
est décrite par
1 𝑇𝑑 𝑠 𝑘𝑐 𝑇𝑖 (𝑇𝑑 + 𝑇𝑓 )𝑠 2 + (𝑇𝑖 + 𝑇𝑓 )𝑠 + 1
∀𝑠 ∈ ℂ, 𝐾(𝑠) = 𝑘𝑐 [1 + + ]=
𝑇𝑖 𝑠 𝑇𝑓 𝑠 + 1 𝑇𝑖 𝑠(𝑇𝑓 𝑠 + 1)

Ce qui nous donne :

𝑇𝑓 = 20 → 𝑇𝑖 + 𝑇𝑓 = 200 → 𝑇𝑖 = 180

7
 TD Automatique Linéaire 1A

2
𝑘𝑐 = × 𝑇𝑖 = 12.5
28.8
9400 9400
𝑇𝑖 (𝑇𝑑 + 𝑇𝑓 ) = 9400 → 𝑇𝑑 = − 𝑇𝑓 = − 20 = 32.2
𝑇𝑖 180

1 32.2 𝑠
∀𝑠 ∈ ℂ, 𝐾(𝑠) = 12.5 [1 + + ]
180 𝑠 20 𝑠 + 1

En fonction du besoin, on peut utiliser sous Matlab :

• la fonction « zpk(Fc) » qui donnera la forme série
• la fonction « pid(Fc) » qui donnera la forme parallèle
• la fonction « pidstd(Fc) » qui donnera la forme mixte standard ou filtrée. Avec cette
dernière fonction dans le cas filtré, la constante de temps de filtrage est fixée comme
un rapport entre 𝑇𝑑 et un nombre 𝑁. Pour respecter le correcteur 𝐾 dont vous cherchez
à la conversion, 𝑁 n’est pas nécessairement un entier.

Question 9 : Valider en simulation votre PID mixte en l’absence de perturbation, puis avec
perturbation.

La réalisation sous Simulink conduit à la simulation suivante :

La réponse en asservissement est conforme au cahier des charges. Le rejet de perturbation est
également correct au sens où elle est bien rejetée asymptotiquement. Dans ce cas particulier, le
temps de rejet est d’environ 120 secondes mais te temps n’a pas été imposé par le calcul.

8
 TD Automatique Linéaire 1A

On constate une déviation de 0.08V environ lors de l’apparition de la perturbation soit 8°C
environ.

Question 10 : Déterminer la fonction de transfert 𝐹𝐹 idéale du block d’anticipation

On a vu en cours que la fonction de transfert du feed-forward est donnée par :

𝐺1 (𝑠) 9400 𝑠 2 + 200 𝑠 + 1 −50

∀𝑠 ∈ ℂ, 𝐹𝐹(𝑠) = − =− ×
𝐺(𝑠) 0.5 90 𝑠 + 1
Soit
9400 𝑠 2 + 200 𝑠 + 1
∀𝑠 ∈ ℂ, 𝐹𝐹(𝑠) = 100
90 𝑠 + 1

Question 11 : Le feed-forward 𝐹𝐹 est-il réalisable ? Dans le cas contraire, proposer une forme
implémentable.

Le degré du numérateur de 𝐹𝐹 est supérieur au degré du dénominateur. Ce transfert est

impropre et donc non réalisable. On peut :
𝐺1 (0)
• Ce contenter d’un feed-forward statique en prenant 𝐹𝐹 = − soit 𝐹𝐹 = 100
𝐺(0)
• Rajouter un pôle négligeable, mais il faudra étudier son influence sur la dynamique de
rejet de la perturbation
9400 𝑠 2 + 200 𝑠 + 1
∀𝑠 ∈ ℂ, 𝐹𝐹(𝑠) = 100
(90 𝑠 + 1)(𝛼 𝑠 + 1)

Question 12 : Simuler la totalité de votre correcteur. Les résultats sont-ils concluants ?

• Feed-forward statique avec 𝐹𝐹 = 100

9
 TD Automatique Linéaire 1A

Les résultats sont relativement améliorés : le temps de rejet est maintenant de 80 secondes
environ et la déviation est maintenant de l’ordre de 0.06V soit 6°C environ.

• Feed-forward dynamique avec un pôle négligeable en -1 pour assurer

9400 𝑠2 +200 𝑠+1
l’implementation, ∀𝑠 ∈ ℂ, 𝐹𝐹(𝑠) = 100 (90 𝑠+1)(𝑠+1)

Cette fois, le rejet est quasi complet et le système a été insensibilisée aux variations de
perturbation.

On notera cependant que la commande présente des pics supérieurs à la valeur limite
de 10V pour la commande. Il conviendra d’en étudier les répercussions.

10
 TD Automatique Linéaire 1A

Annexe

- Le critère de performance retenue est (𝑦 la réponse enregistrée, 𝑦𝑠 la sortie simulée du

modèle)
𝑡𝑓 𝑡𝑓
1 1
𝐽= ∫ 𝑒 2 𝑑𝑡 = ∫ (𝑦 − 𝑦𝑠 )2 𝑑𝑡
2 2
𝑡0 𝑡0

- La méthode du gradient adapte le vecteur 𝑞 des paramètres à chaque itération 𝑖 avec

𝜕𝐽
𝑞𝑖+1 = 𝑞𝑖 − 𝜆−1 |
𝜕𝑞 𝑞=𝑞
𝑖

- Le gradient de 𝐽 est défini par

𝑡𝑓
𝜕𝐽 𝜕𝑦𝑠
=−∫𝑒 𝑑𝑡
𝜕𝑞 𝜕𝑞
𝑡0

𝜕𝑦𝑠
- Le calcul de , où 𝑦𝑠 = 𝑓(𝑞, •) ∗ 𝑢, s’appuie sur la propriété2 de la transformée de
𝜕𝑞
𝜕𝑓 𝜕
Laplace suivante : soit 𝑓 une fonction de paramètres 𝑞 alors ℒ [𝜕𝑞] = 𝜕𝑞 ℒ[𝑓]

𝜕𝑦𝑠 𝜕 𝜕𝑌𝑠 𝜕𝐹
∀𝑠 ∈ ℂ, ℒ[ ] (𝑠) = ℒ[𝑦𝑠 ](𝑠) = (𝑠) = (𝑠)𝑈(𝑠)
𝜕𝑞 𝜕𝑞 𝜕𝑞 𝜕𝑞

Exemple dans le cas d’un système du 1er ordre :

𝑁(𝑞, 𝑠) 𝑞2
𝐹(𝑞, 𝑠) = =
𝐷(𝑞, 𝑠) 𝑞1 𝑠 + 1
𝑞1 𝑓 ′ 𝑓 ′ 𝑔−𝑓𝑔′
Soit un vecteur de paramètre 𝑞 = [𝑞 ]. On s’appuie alors sur la propriété (𝑔) = 𝑔2
2

𝜕𝑌𝑠 𝜕𝐹 −𝑠 𝑞2 𝑠
(𝑠) = (𝑠)𝑈(𝑠) = 𝑈(𝑠) = − 𝑌 (𝑠)
𝜕𝑞1 𝜕𝑞1 (𝑞1 𝑠 + 1) 2 𝑞1 𝑠 + 1 𝑠
De la même manière, on obtient
𝜕𝑌𝑠 𝜕𝐹 1
(𝑠) = (𝑠)𝑈(𝑠) = 𝑈(𝑠)
𝜕𝑞2 𝜕𝑞2 𝑞1 𝑠 + 1

2
cf polycopié d'Analyse Appliquée de 1A ECL - proposition H5 p 45
11
 TD Automatique Linéaire 1A

𝑡𝑓

𝜕𝑌𝑠 𝜕𝑦𝑠
𝜕𝑦𝑠 ℒ −1 [ | ] 𝜕𝐽 −∫𝑒 𝑑𝑡
𝜕𝑞1 𝑞=𝑞 𝜕𝑞1
𝜕𝑦𝑠 𝜕𝑞1 𝑖
𝜕𝐽 𝜕𝑞1 𝑡0
= = → = = 𝑡𝑓
𝜕𝑞 𝜕𝑦𝑠 𝜕𝑌𝑠 𝜕𝑞 𝜕𝐽
ℒ −1 [ | ] 𝜕𝑦𝑠
[𝜕𝑞2 ] 𝜕𝑞2 𝑞=𝑞 ] [𝜕𝑞2 ] −∫𝑒 𝑑𝑡
[ 𝑖 𝜕𝑞2
[ 𝑡0 ]

𝜕𝐽
𝑞1 𝑞1 𝜕𝑞1
[𝑞 ] = [𝑞 ] − 𝜆−1
2 𝑖+1 2 𝑖 𝜕𝐽
[𝜕𝑞2 ]

12