Simulation Et Optimisation Des Systèmes Industriels

Universit de Lige
Facult des Sciences Appliques
SIMULATION ET OPTIMALISATION
STATIQUE DES SYSTMES INDUSTRIELS
Succession Prof. B. Kalitventzeff

G. Heyen, Matre de confrence
Laboratoire d'Analyse et de Synthse des Systmes Chimiques
Chapitre 1
CHAPITRE 1
SIMULATION STATIQUE DES SYSTEMES CHIMIQUES
Introduction
Nous allons considrer ltude du comportement dun ensemble dappareils (un systme) par
simulation. Nous considrerons que nous disposons dun modle de chacun des appareils, cest--dire
dun moyen de calculer ltat des flux sortant dun appareil lorsque ses alimentations sont connues.
La simulation statique des systmes chimiques induit systmatiquement des problmes de
convergence. Prenons comme exemple le schma trs simple ci-dessous reprsentant une petite partie dun
procd chimique. Pour pouvoir rsoudre les quations relatives une unit (un appareil), il nous faut
connatre les flux dentre de cette unit. Or, sil y a des boucles, nous ne connaissons pas ncessairement
tous les flux dentre : pour calculer lchangeur HE1, il faut connatre le flux 2. Celui-ci dpend
successivement des flux 10, 8, 5 et 6, 4. Donc il nous faut les flux 4 et 6 pour connatre le flux 2. Or le flux
4 est obtenu par la rsolution de lchangeur HE1.
Pour rsoudre ce problme, on devra raliser une coupure de flux, cest--dire que lon devra
estimer un des flux puis calculer toutes les units de la srie et enfin utiliser une mthode de promotion de
convergence afin dobtenir une meilleure estimation des grandeurs concernant le flux que lon a coup.
Figure 1.1
- 1.1 -
Chapitre 1
Par exemple, si lon coupe le flux 2, on pourra calculer lchangeur HE1, puis le racteur,
lchangeur HE2, le sparateur et le diviseur. On obtiendra ainsi de nouvelles valeurs pour le flux 2. On
peut alors effectuer la convergence sur le flux 2.
Si on considre un schma plus complexe, il nest pas vident de dterminer directement les flux
couper de manire avoir un nombre minimum de coupures (pour limiter les problmes de convergence).
Mthode de Motard1
Une mthode pour dterminer le nombre minimum de coupures a t mise au point par MOTARD.
Nous utiliserons un exemple afin dillustrer cette mthode.
y
52
(7)
y
43
(-)
y
1
12
(1)
(4)
y
y
2
23
(2)
34
(3)
y
4
45
(6)
y
5
56
(8)
y
6
(-)
(5)
y
42
Figure 1.2
La figure ci-dessus reprsente un systme compos de 6 units interconnectes entre elles

directement ou par lintermdiaire de recyclages . Nous pourrions rencontrer galement des by-pass.
Chaque flche reprsente un flux dirig dune certaine grandeur (matire, nergie,...) quil faut
dterminer. Pour chaque unit, il existe un modle pouvant tre rsolu permettant de calculer les flux de
sortie en fonction des entres. Dans le contexte du systme, elle constitue un lment connu.
Nous considrons que le systme est en rgime, cest--dire quil ny a pas daccumulation. Sil y
a une transformation chimique (racteur par ex.), il est ncessaire dintroduire un flux fictif exprimant la
quantit de produits forme et la quantit de ractifs consomme. Nous nenvisagerons pas ici ce cas de
faon explicite.
Avant de rsoudre le systme ci-dessus, nous allons tout dabord introduire la notion de graphe
dual et noncer quelques rgles systmatiques qui nous permettront de rsoudre le systme en organisant
la squence de calcul de manire optimale.
Le graphe dual
Barkley R.W., Motard R.L., Decomposition of Nets, The Chemical Engineering Journal, 1972, 3, 265-275.
- 1.2 -
Chapitre 1
Le graphe dual est construit partir du flow-sheet du systme que lon veut rsoudre. La grande
diffrence entre ces deux types de graphes se situe au niveau des significations des branches et des noeuds.
Les branches (noeuds) du flow-sheet deviennent les noeuds (branches) du graphe dual.
3
1
C
Figure 1.3 Flow-sheet
Figure 1.4 Graphe dual
Les deux figures ci-dessus montrent la relation existant entre un flow-sheet et un graphe dual.
Rgles de Motard
Une rgle gnrale consiste :
COUPER LE FLUX INTERVENANT DANS

LE PLUS GRAND NOMBRE DE CYCLES
Pour trouver ce flux de manire systmatique, nous considrons le graphe dual du flow-sheet et
nous appliquerons les rgles suivantes en revenant la premire chaque fois quune rgle aura t utilise :
1. supprimer les flux qui nont pas de flux antrieur. Ils sont ncessairement spcifis ;
2. remplacer les flux par leur antrieur, sil ny en a quun. Ds quun flux apparat deux fois, le
supprimer une fois (ou ne pas lcrire) ;
3. couper les boucles propres (ce sont les boucles sur un mme noeud) ;
Figure 1.5
Le calcul du flux 1 implique une boucle de calcul itratifs donc une coupure.
4. couper un des flux parallles de sens contraire . Couper la branche qui a le plus de flux postrieurs
(qui est le plus grand nombre de fois lantrieur dautres) ;
- 1.3 -
Chapitre 1
Figure 1.6
5. si deux boucles impliquent un mme noeud (donc un mme flux), cest celui-ci qui doit tre
coup. En effet, ce faisant on coupe les deux boucles par un seule coupure.
Figure 1.7
Remarques :
on est rarement amen appliquer la dernire rgle, car les autres suffisent gnralement,
lorsquon applique une des rgles, il faut toujours recommencer la 1re rgle pour continuer.
Ces rgles dues MOTARD ont t dictes pour permettre la dcomposition des flow-sheets
pralablement une simulation : elles permettent de dterminer et de localiser le nombre minimal de
coupures pour que le calcul des units puisse se faire de manire squentielle et itrative.
Pour connatre lendroit o nous allons couper le systme donn figure 1.2, appliquons les rgles
proposes ci-dessus.
A partir du graphe rduit, nous pouvons dresser le tableau :
FLUX
FLUX ANTERIEURS
1
2
3
4
5
6
7
8
1, 5, 7
2, 4
3
3
3
6
6
Le flux 1 nayant pas de flux antrieur est supprim (rgle 1). Do le tableau :
- 1.4 -
Chapitre 1
FLUX
FLUX ANTERIEURS
2
3
4
5
6
7
8
5, 7
2, 4
3
3
3
6
6
Remplacer les flux 4, 5, 6 par leur flux antrieur unique 3 (rgle 2). Do le tableau :
FLUX
FLUX ANTERIEURS
2
3
4
5
6
7
8
3, 7
2, 3
3
3
3
3
3
Remplacer le flux 7 par son flux antrieur 3. Ce dernier existant dj comme antrieur de 2, il est
inutile de lcrire deux fois (rgle 2). Do le tableau :
FLUX
FLUX ANTERIEURS
2
3
4
5
6
7
8
3
2, 3
3
3
3
3
3
Remplacer le flux 2 par son flux antrieur 3 et ne pas lcrire car 3 existe dj comme antrieur de
3.
Do le tableau final :
- 1.5 -
Chapitre 1
FLUX
FLUX ANTERIEURS
Il est vident quil est inutile de recopier chaque fois le tableau. Gnralement, les corrections
successives sont apportes au 1er tableau (ventuellement en couleur) jusqu obtenir ce dernier tableau.
Nous remarquons que tous les flux ont le flux 3 comme flux antrieur (sauf le flux 1 qui nen a pas).
Regardons maintenant lapplication des rgles de Motard au travers du graphe dual. Le 1er graphe
dual est obtenu partir du 1er tableau des flux (attention, ici les numros reprsentent les flux).
Figure 1.8
Grce au 2me tableau, le graphe dual devient :
Figure 1.9
- 1.6 -
Chapitre 1
Figure 1.10
Figure 1.11
Le dernier tableau donne le graphe dual :
Figure 1.12
Le dernier tableau, ainsi que ce graphe, permettent de voir que tous les flux ont le flux 3 comme
flux antrieur, mme le flux 3. Nous avons donc une boucle propre sur le flux 3.
Nous devons donc couper le flux 3.
Tous les autres flux disparaissent par application de la rgle 1 car si le flux 3 est supprim, les
autres nont plus dantrieur.
Nous pouvons alors commencer le processus itratif. Le schma suivant considre un calcul de
simulation.
- 1.7 -
Chapitre 1
Envisageons une seule variable de description des flux (dbit par ex.) et donnons au flux 3 une
premire valeur, soit a. On calcule les flux 4, 5 et 6. On calcule ensuite le flux 7 [ 7 = f6(6) ]. A ce moment,
on est en mesure de calculer le flux 2 [ 2 = f1(1, 5, 7) ] et le flux 3 [ 3 = f2(2, 4) ].
Notons a cette nouvelle valeur du flux 3.
Si a a, on recommence le calcul avec a ou plutt avec un nouveau a, fonction de a ( par ex. : a
= a, a est rebaptis a ) jusqu ce que :
| a - a | , fix.
Ayant trouv le flux 3 convergence, nous pouvons alors dterminer tous les autres flux.
On peut schmatiser la squence des calculs comme suit :
Dbut
6
a = f(a)
7
a
Satisfait
donc fin
Test : a - a<
Figure 1.13
Remarques :
Une solution moins conomique aurait, par exemple, consist couper les flux 4, 5 et 7. On aurait ainsi
obtenu un calcul squentiel, mais on aurait eu un processus itratif plus compliqu.
- 1.8 -
Chapitre 1
La thorie prsente nexclut videmment pas que lon procde dabord toutes les liminations
possibles lorsquelles sont faciles raliser.
Exemple :
1
7
5
ij
Figure 1.14
Nous pouvons maintenant traiter un cas plus complexe, soit le systme dont le flow-sheet est donn
ci-aprs :
Sortie
gaz
Cl 2
PRODUIT
11
13
12
HCl
6
10
Paraffines
9
Figure 1.15
2
3
5
Cl
11
5
18
19
14
7
11
17
12
12
16
15
10
13
13
20
10
21
Paraffines
Figure 1.16
- 1.9 -
22
HCl
Chapitre 1
Application des rgles expliques dans ce chapitre :

flux
1
1, 11
15
12
15
12
3, 9, 13
13, 15
12, 16
16
4, 14, 17
4, 14
4, 12
10
15
12
1, 11
15
12
10
15
12
10
15
12
11
15
15
12
12
3, 9, 13
13, 15
12, 16
13
16
16
16
16
14
15
15
12
15
12
12
12
16
4, 14, 17
4, 15
4, 12
17
21
18
15
12
19
20
18, 19
18, 19
18, 19
21
22
20
20
20
couper 12
couper 4
Nous aurons donc une boucle pour trouver le flux 4, accompagne dun test de convergence, ainsi
que sur le flux 12. Il sagit dun calcul o deux boucles de convergence sont imbriques, il faudra faire un
pas en avant sur 12 lorsque 4 est convergence.
Le schma des calculs est reprsent ci-aprs.
- 1.10 -
Chapitre 1
Les entres sont 1 et 21; les sorties sont 22, 6, 7.

21
17
4
b
19
16
13
14
a
12
15
1
11
20
22
b'
18
10
a'
12
9
7
Figure 1.17
Bien remarquer que tous les calculs intermdiaires qui ne sont pas ncessaires (flux 6, 19, 20, 22),
ne sont pas calculs chaque fois mais uniquement convergence.
On calcule 17 connaissant lentre 21.
On choisit un vecteur a de donnes (T, P, dbits molaires partiels) pour le flux 12, et on calcule 15
qui donne 14 et 11. Le flux 11 et lentre 1 fournissent 2 et 8, avec 2 on trouve 3 (et 18) tandis quavec 8
on trouve 10 puis 9 (et 7).
On choisit des valeurs pour les composantes du vecteur b correspondant au flux 4 (17 et 14 sont
connus), on calcule 16 (et 5) puis 13. Ensuite 9, 13 et 3 permettent de calculer 4 = b.
Test sur b = b, nouvelles valeurs pour b ... jusqu convergence de la boucle sur 4.
On calcule 12 ( partir de 13, 3, 9) = a ; test sur a - a, nouvelles valeurs pour a, etc., jusqu
convergence sur la boucle 12.
On calcule alors les flux 6, 19, 18, 20.
Choix des coupures
Lapplication des rgles de MOTARD permet de dterminer le nombre minimal de coupures
effectuer pour la rsolution du systme chimique. Cependant, parfois pour rduire le temps de calcul, il est
intressant de faire une ou deux coupures supplmentaires ou de dplacer lendroit de la coupure le long
dune boucle. Finalement, on en vient dfinir plusieurs types de coupures en fonction du nombre de
variables inconnues :
coupures totales (il faut trouver n dbits partiels, T et P),
- 1.11 -
Chapitre 1
coupures massiques (seul le dbit total est inconnu),

coupures thermiques (seul la temprature est inconnue).
Prenons le schma suivant :
Figure 1.18
Pour rsoudre ce systme, nous allons devoir effectuer une coupure. On peut reprsenter le flux en
nimporte quel point du schma, par un vecteur qui comporte les variables dtat, les spcifications du
procd, les quantits de matire : soit par exemple, un vecteur F qui contiendrait,
la temprature,
la pression,
le dbit total,
la composition,
O faut-il couper ?
Si nous coupons en 2, la temprature, les compositions et ventuellement la pression ne sont pas
connues priori.
On effectuera donc une coupure totale sur tout le flux pour rsoudre le systme.
Par contre, entre 1 et le flux dentre, le dbit total et compositions ne varient pas, la pression sera
obtenue par une perte de charge. Le seul paramtre qui va changer est la temprature et il suffira
deffectuer une coupure thermique. On ralisera ainsi une conomie de calcul car il faut itrer sur 1 seule
variable.
Ces derniers paragraphes mettent en vidence la diffrence fort importante entre les coupures
totales et partielles (compositions ou temprature, etc.).
Donc avant deffectuer une coupure, il faudra toujours examiner sil ny a pas moyen deffectuer
une coupure partielle, ce qui sera toujours intressant du point de vue temps de calcul.
- 1.12 -
Chapitre 1
Application : la synthse de lammoniac

air de
procd
CH
H 2O
Reforming
primaire
Reforming
secondaire
CC
Mthanation
Dcarbonatation
Shift
HT - BT
purge
Compression
Dssication
Synthse
NH3
Figure 1.19 Schma gnral de la synthse
La synthse du NH3 passe par diffrentes tapes : il y a dabord la prparation du gaz de synthse
comprenant notamment les reforming primaires et secondaires, la conversion du CO, la dcarbonatation et
la mthanation, et ensuite la synthse proprement dites. Dcrivons brivement lutilit des diffrentes
oprations de prparation du gaz de synthse :
Reforming primaire et secondaire : formation de lhydrogne partir de mthane et de vapeur deau (
CH4 + H2O CO + 3 H2 ).
Conversion du CO : la raction de conversion transforme le CO en CO2, avec production
complmentaire dhydrogne ( CO + H2O CO2 + H2 ).
Dcarbonatation : elle est destine liminer le CO2 par absorption et raction avec une solution de
carbonate de potassium chaude et active.
Mthanation : il sagit de purifier le gaz de synthse des traces de CO et de CO2 encore prsentes en les
transformant en CH4 ( CO + 3 H2 CH4 + H2O et CO2 + 4 H2 CH4 + 2 H2O ).
- 1.13 -
Chapitre 1
M102
3
35
E102
T101
E107
10
M101
9
8
E106
18
Figure 1.20 Schma gnral des reforming primaire et secondaire
- 1.14 -
33
35 bar
EAU
32
T102
31
30
E103
E104
20
17
7
1
19
G.N.
COMBUSTION
G.N.
PROCESS
21
E105 16
26
D101
34
14
D102
27
15
28
PURGE
23
24
22
B101
100 bar
25
M103
AIR
PROCEDE
13
E101
12
11
R101
36
AIR DE
COMBUSTION
R102
E108
29
Afin dillustrer la mthode de Motard, nous considrons uniquement la section reforming du

procd de fabrication. Cette section est schmatise ci-dessous :
Chapitre 1
Flux
1
1, 18
16
16
16
2, 11, 36
11, 16
11, 16
10, 16
30
30
3, 35
3, 35
3, 13
3, 10
30
30
30
30
5, 28
4, 28
4, 23
4, 23
30, 24
30
7, 17
16
16
16
31
31
30
30
30
30
10
8, 9
16, 31
16, 30
16, 30
30
30
11
10, 12
10, 12
10, 12
10
30
30
12
2, 11, 36
11, 16
11, 16
13
10, 12
10, 12
10, 12
10
30
30
14
13, 34
13
13
10
30
30
15
14, 25
13, 25
13, 25
10, 25
24, 30
30
16
15, 26
15, 26
15, 23
15, 23
17
16, 20
16
16
16
18
7, 17
16
16
16
19
1, 18
16
16
16
20
21
16, 20
16
16
16
22
21, 24
21, 24
16, 24
16, 24
24
23
21, 24
21, 24
16, 24
16, 24
24
24
27, 29
27, 29
27, 29
27, 29
27, 29
25
21, 24
21, 24
16, 24
16, 24
24
26
23
23
23
23
24
27
15, 26
15, 23
15, 23
15, 23
15, 24
30
28
23
23
23
23
24
29
5, 28
4, 28
4, 23
4, 23
24, 30
30
30
14, 25
13, 25
13, 25
10, 25
24, 30
30
31
30
30
30
30
30
30
32
31
31
30
30
30
30
33
32
32
30
30
30
30
34
35
13, 34
13
13
10
30
30
36
couper 12
couper 16
couper 24
couper 30
Le tableau ci-dessus montre lapplication de la mthode de Motard la section reforming du

procd de synthse de lammoniac. On voit donc quil est ncessaire deffectuer quatre coupures : flux
12, 16, 24 et 30. Les coupures 12 et 16 sont des coupures thermiques, car la composition des fumes est
connue (combustion totale).
- 1.15 -
Chapitre 2
CHAPITRE 2
MODELES MATHEMATIQUES DES SYSTEMES
CHIMIQUES ISOLES
Introduction
Le but de la simulation des systmes chimiques est de pouvoir, au moyen de modles
mathmatiques, reprsenter le comportement de ces systmes. La capacit des modles de simulation
reproduire le plus exactement possible un comportement est lie lensemble des hypothses effectues
pour tablir ces modles.
Deux types de simulation sont ralisables :
la simulation statique du systme, cest--dire lorsque ce dernier a atteint son rgime stationnaire (le
procd nvolue plus dans le temps) ;
la simulation dynamique qui permet de modliser le comportement du systme dans un tat transitoire
(dmarrage, perturbations, actions de contrle...).
Dans le cadre de ce cours, nous nous limiterons la simulation statique des systmes, la simulation
dynamique faisant lobjet dun cours option en dernire anne.
La modlisation dun systme pose la question de savoir quel modle nous allons utiliser cette fin.
En effet, un grand nombre de modles existent et donner une liste exhaustive de tous les modles
concernant lensemble des systmes chimiques nest pas lobjectif de ce cours. Nous en tablirons
toutefois quelques-uns ci-aprs, laide de modles dterministes, en explicitant la mthode gnrale des
bilans matriels, nergtiques et de quantit de mouvements.
2.1 Mthode gnrale des bilans

La modlisation dun systme fait appel des quations de bilans et des quations propres
lunit opratoire : les quations de simulation.
Pour un systme N constituants, ces bilans scrivent :
- N quations de bilan de matire (dbits molaires partiels);
- une quation de bilan nergtique (flux enthalpiques);
- une quation de bilan dimpulsion (pression).
Soit, au total, (N+2) quations de bilans.
- 2.1 -
Chapitre 2
Dautre part, un certain nombre de variables caractrisent ltat du systme. La diffrence entre le
nombre de variables du systme et le nombre total dquations est gal au nombre de degr de libert du
systme. Il est alors ncessaire dimposer la valeur dautant de variables quil y a de degrs de libert : ce
sont les spcifications du systme.
Les quations de bilans se prsentent toutes sous la forme suivante :
Accumulation
nette dans le
volume du
systme
Importation
Exportation
travers la
travers la
=
+
surface du
surface du
systme
systme
Gnration
Consommation
dans le
dans le
volume du
volume du systme
systme
Lobjectif de ce cours tant la simulation statique des systmes, le terme daccumulation sera nul
puisque le systme est suppos avoir atteint son rgime stationnaire.
La dcomposition des modles et leur rsolution sera vue plus loin (2.4). Nous ne nous tendrons
donc pas sur celle-ci par la suite.
2.2 Modles mathmatiques de diffrents appareils

Les modles que nous allons dcrire seront des modles description macroscopique, cest--dire
quils ne ncessitent pas la connaissance locale des variables. Dans les bilans, les variables dtat du
systme ont pour valeurs les moyennes correspondant lensemble du volume de ce systme.
Il existe galement des modles description microscopique dont nous verrons un exemple plus
loin dans le cours (2.3.2.1). Dans les modles examins ci-aprs, on supposera que les grandeurs
caractrisant les flux dentre sont connues : elles ne seront pas considres comme variables.
2.2.1 Le mlangeur
On ralise le mlange de plusieurs flux n constituants.
F1
T1,P1,xi1
T,P,xi
F
F2
T2,P2,xi2
Hypothse : il ny a pas de changement de phase dans le mlangeur.

Variables :
- 2.2 -
Chapitre 2
temprature T (K)
pression P (bar)
dbit molaire F (kmol/s)
fractions molaires xi (-)
du flux de sortie
idem
idem
idem
1
1
1
n
total
n+3
Equations :
bilan de matire
bilan thermique
bilan dimpulsion
F.xi = F1.xi1 + F2.xi2

F.H(T,P,xi) = F2.H2(T2,P2,xi2) + F1.H1(T1,P1,xi1)
remplac par
P = min(P1,P2)
n
normation (fractions molaires)
=1
n
1
1
1
total
n+3
Nombre de degr de libert = (n+3) variables - (n+3) quations = 0

Il ny a pas de degr de libert pour le mlangeur si nous connaissons tous les flux entrant dans le
mlangeur, cest--dire : la temprature Tk , la pression Pk , les dbits molaires totaux Fk et les fractions
molaires partielles xik de chacun des flux dentre.
Remarques :
si on fixe la temprature T, il nest plus ncessaire dcrire le bilan thermique;
de mme, si on fixe la pression P, il nest plus ncessaire de raliser un bilan dimpulsion;
on peut supprimer lquation de normation condition de remplacer les fractions molaires et les dbits
molaires totaux par les dbits molaires partiels (Xi=F.xi) et de raliser le bilan de matire sur ces
derniers. Le modle comporte alors (n+2) quations et variables et le nombre de degr de libert est
toujours nul.
2.2.2 Le diviseur
On ralise la division dun flux n constituants en plusieurs flux en sortie du diviseur. Nous
choisissons un diviseur sparant un flux en deux branches.
T1,P1,xi1
F1
T,P,xi
F
T2,P2,xi2
F2
Hypothse : il ny a pas de changement de phase dans le diviseur.
Variables :
- 2.3 -
Chapitre 2
tempratures T1,T2 (K)

pressions P1,P2 (bar)
dbits molaires F1,F2 (kmol/s)
fractions molaires xi1,xi2 (-)
2
2
2
2n
total
2n+6
Equations :
bilan de matire
bilan thermique
bilan dimpulsion
bilan dimpulsion
F.xi = F1.xi1 + F2.xi2

F.H(T,P,xi) = F2.H2(T2,P2,xi2) + F1.H1(T1,P1,xi1)
P1 = P
P2 = P
n
1
1
1
normation (fractions molaires)
i2
=1
dimensionnement
dimensionnement
n
1
xi1 = xi
T1 = T
total
2n+5
Nombre de degr de libert = (2n+6) variables - (2n+5) quations = 1

Nous disposons encore dun degr de libert, nous devons donc imposer au systme une
spcification : le taux de division par lintermdiaire dun des dbits molaires de sortie (F1 ou F2) ou du
rapport de lun de ces dbits avec celui lentre (F1/F ou F2/F). Remarques : nous pouvons faire les
mmes que dans le cas du mlangeur.
2.2.3 La vanne de dtente
On ralise la dtente dun flux n constituants dans une vanne.
T1,P1
T2,P2
F2
F1
Hypothses : il ny a pas de changement de phase dans la vanne et la dtente se fait de manire

adiabatique.
Variables :
temprature de sortie T2 (K)
pression de sortie P2 (bar)
dbits molaires totaux des flux sortants F2 (kmol/s)
fractions molaires des flux sortants xi2 (-)
total
1
1
1
n
n+3
Equations :
- 2.4 -
Chapitre 2
bilan de matire (dbits molaires partiels)

bilan de matire (dbits molaires totaux)
bilan thermique
F2.xi2 = F1.xi1
F2 = F1
F2.H2(T2,P2,xi2) = F1.H1(T1,P1,xi1)
total
n
1
1
n+2

spcification : le taux de dtente par lintermdiaire de la pression de sortie P2 ou du rapport de la pression
de sortie celle de lentre P2/P1.
2.2.4 Le sparateur liquide - vapeur
On ralise la sparation des phases vapeur et liquide dun flux n constituants.
V
zi
F
yi
T P
Q
Q
L
xi
Hypothse : lalimentation est compltement fixe.

Variables :
temprature T (K)
pression P (bar)
dbit molaire liquide L (kmol/s)
dbit molaire vapeur V (kmol/s)
fractions molaires liquide et vapeur xi, yi (-)
enthalpie transfre Q (kW)
1
1
1
1
2n
1
total
2n+5
Equations :
bilan de matire
bilan thermique
F.zi = L.xi + V.y i

F.Hf = V.HV(T,P,y i) + L.HL(T,P,xi) + Q
n
normation (vapeur)
=1
=1
y i = xi.Ki(T,P,xi,y i)
1
n
normation (liquide)
n
1
quilibre liquide-vapeur
total
2n+3
- 2.5 -
Chapitre 2
Nombre de degr de libert = (2n+5) variables - (2n+3) quations = 2

Nous disposons encore de deux degrs de libert, nous devons donc imposer au systme deux
spcifications : par exemple, la pression et la temprature dans le sparateur (quilibre P et T fixs : EPT)
ou la pression et la fraction vaporise (EPA) ou la temprature et la fraction vaporise (ETA) ou encore la
pression et la chaleur transfre.
2.2.5 Le compresseur
On ralise la compression dun flux n constituants.
W
T2,P2
F2
F1
T1,P1
Hypothses : il ny a pas de condensation et lalimentation est compltement fixe.

Variables :
dbit molaire de sortie F2 (kmol/s)
fractions molaires de sortie xi2 (-)
travail W (kW)
travail de compression isentropique Wis (kW)
temprature de compression isentropique Tis (K)
1
1
1
n
1
1
1
total
n+6
Equations :
bilan de matire
bilan thermique
F2.xi2= F1.xi1
F1.H1(T1,P1,xi1) = F2.H2(T2,P2,xi2) - W
n
normation
n
1
=1
W = W is/is
1
1
i2
dimensionnement (efficacit isentropique)

dimensionnement (travail isentropique)
dimensionnement (entropie)
W is = His(Tis,P2) - H1(T1,P1,xi1)
S(Tis,P2) = S(T1,P1)
total
n+5
- 2.6 -
Chapitre 2
spcification : le taux de compression par lintermdiaire de la pression en sortie du compresseur P2 ou du
rapport de celle-ci avec la pression lentre (P2/P1) ou bien alors le travail mcanique W fourni au
compresseur.
2.2.6 La turbine
On ralise la dtente dun flux n constituants dans une turbine afin de produire un travail
mcanique. La modlisation de la turbine se fait comme pour le compresseur, except pour la premire
quation de dimensionnement qui est modifie comme suit : W=Wis.is.
W
T1,P1
F1
T2,P2
F2
Hypothses : il ny a pas de condensation et lalimentation est compltement fixe.

Variables :
dbit molaire de sortie F2 (kmol/s)
fractions molaires de sortie xi2 (-)
travail W (kW)
travail de compression isentropique Wis (kW)
temprature de compression isentropique Tis (K)
1
1
1
n
1
1
1
total
n+6
Equations :
bilan de matire
bilan thermique
F2.xi2 = F1.xi1
F1.H1(T1,P1,xi1) = F2.H2(T2,P2,xi2) - W
n
normation
n
1
=1
W = W is.is
1
1
1
i2
dimensionnement
dimensionnement
dimensionnement (entropie)
W is = His(Tis,P2) - H1(T1,P1,xi1)
S(Tis,P2) = S(T1,P1)
total
n+5
- 2.7 -
Chapitre 2
spcification : le taux de dtente par lintermdiaire de la pression en sortie du compresseur (P2)ou du
rapport de celle-ci avec la pression lentre (P2/P1) ou bien alors le travail mcanique W fourni par la
turbine.
2.3 Modles mathmatiques des racteurs

2.3.1 Modles description macroscopique
Ce type de modle ne ncessite pas la connaissance locale des variables et donne lieu une
description mathmatique trs simplifie. Dans les bilans, les variables dtat comme la temprature et les
concentrations sont des variables dont on applique les valeurs moyennes tout le volume du systme.
Exemple : racteur continu cuve parfaitement mlange (R.C.P.M.).
2.3.1.1 R.C.P.M. isotherme
Soit la raction a A b B (ou A B avec = b/a) ralise dans un RCPM isotherme dont un
schma simplifi est donn ci-dessous :
0
T , CA , CB , F
T , V
T, C A, C B , F
Figure 2.1
avec
la temprature en K,
le volume en m3,
la concentration en A ou en B, en mole/kg,
dbit massique en kg/s plutt quen m3/s ou mole/s, car ils restent invariables mme si la
temprature et/ou le nombre de moles varie(nt),
si r est la vitesse de raction en mole/m3.s, on peut crire :

les bilans de matire :
0 = F.C0A F.C A r.V
(1)
0 = F.C 0B F.C B . r. V
(2)
- 2.8 -
Chapitre 2
le bilan nergtique :
T = T0
car le racteur est isotherme,
(3)
P = P0
si on nglige les pertes de charge,
(4)
Remarques concernant les bilans matires (quations 1 et 2) :

si la raction est du premier ordre alors, r = k.CA et on a :
0 = F.C 0A F.C A k.V.C A
do :
CA =
C 0A
k.V
1+
F
et le problme est alors trs simple ;

si la raction est du premier ordre dans un sens et du second dans lautre alors :
r = k 1 .CA k 2 . C2B
et (1) devient :
0 = F.C0A F.C A k 1 . V.C A + k 2 . V.C2B
et le problme devient dj plus complexe ;
si nous avons deux ractions parallles :
r1
A
r2
C
la premire cintique restant :

r1 = k 1 .C A k 2 .C 2B
alors que la cintique de la deuxime raction scrit :
r2 = k 3 .C A .C B
alors :
0 = F.C0A F.C A k 1 . V.C A + k 2 . V.C2B k 3 .C A .C B .V
0 = F.C0B F.C B + k 1 . V.C A k 2 . V.C 2B
0 = F.C 0C F.C C + k 3 .C A .CB . V
ce qui fait trois quations non linaires en CA, CB, et CC .
- 2.9 -
(1)
Chapitre 2
2.3.1.2 Adiabatique
Le bilan thermique scrit :
0 = F.C P .T0 F.CP .T + ( H i ) .ri . V
i
2.3.1.3 Systme non adiabatique

Un change de chaleur avec un fluide circulant dans une double enveloppe apporte encore un
terme supplmentaire au bilan thermique (3) :
0 = F.C P .T0 F.C P .T + ( H i ).ri .V - K.S.(T - TW )
i
avec
le coefficient de transfert de chaleur (kJ/m2.s.K),
la surface dchange (m2),
Tw
la temprature du fluide de la double enveloppe (K).
Il faut attirer lattention sur les H qui apparaissent dans les bilans thermiques : ce sont les
enthalpies chimiques . Elles sont calcules par rapport une temprature de rfrence T0 et si la
temprature laquelle les ractions se droulent est gale T, alors on doit appliquer la formule suivante
(cas des gaz idaux) :
H(T) = H(T ) +
0
T0
i
.Ci . dT
T0
i P
i
la temprature de rfrence (souvent celle du gaz parfait 25C et 1 atm) ;

le coefficient stoechiomtrique de la substance i.
Noublions pas que les H de raction sont ngatifs lorsque le systme chimique est exothermique,
ce qui explique les signes dans les quations des paragraphes prcdents.
2.3.2 Modles description microscopique
Ici, le systme est reprsent par un continuum. Nous distinguerons deux cas :
modles paramtres globaux (lumped parameters), cest--dire valables en tout point du systme,
modles paramtres locaux ou distribus (distributed parameters), cest--dire fonction de
grandeurs qui varient dun point lautre du systme.
2.3.2.1 Exemple : Racteur continu lit catalytique fixe
Soit le systme isol constitu par le lit catalytique dun racteur tubulaire de rayon R et de
longueur L occup par le catalyseur. Nous supposerons que le rgime dcoulement est de type piston.
- 2.10 -
Chapitre 2
Soient mc la masse (kg) et rc la masse spcifique (kg/m3) du lit catalytique qui constitue le racteur
reprsent ci-dessous :
F
dmc
dl
Figure 2.2
Dans ce cas, les bilans sont calculs sur une tranche trs petite dl. Sur celle-ci, le bilan matire du
composant A scrit :
0 = F. C A F.(C A + dC A ) r. dV. c
qui, aprs transformations, devient :

F.dCA = r.dm c
r est la somme algbrique des vitesses de raction rj , r = r j ;

dmc est la masse lmentaire de catalyseur dans dV, dmc = c.dV ;
tandis que lquation aux dimensions ne change gure par rapport celle de (1), si ce nest que la vitesse
sexprime en kmole par kg de catalyseur et par seconde.
Quant au bilan nergtique sur la tranche, il scrit :
0 = F.C p . T F.C p .( T + dT) + (- H i ). r i .dm c K. dS.( T TW )
i
On a suppos le Cp constant dans lintervalle dT, cependant, dans le cas dune raction trs
exothermique, on doit tenir compte de sa variation, on crit H(T) - H(T+dT) o H reprsente lenthalpie
"sensible" plutt que les termes avec la chaleur spcifique constante. La rfrence pour le calcul des
enthalpies sensibles ou de raction se fait toujours partir dune temprature T0 qui peut tre gale T
pour simplifier les expressions. Ajoutons que lhypothse des chaleurs spcifiques constantes entrane une
sorte dindiffrence sur le choix de la temprature de rfrence, ce qui dailleurs est source de confusion.
Si nous effectuons le changement de variable :
dl
dm c
= d ou bien
= d
L
mc
nous obtenons les quations suivantes comme bilan de matire du constituant i :
- 2.11 -
Chapitre 2
dCi
m
= c .rj
d
F j
et comme bilan nergtique :
dT
m
2 . K. L. R
= c . ( H j ). r j
.( T TW )
d
F. C p j
F. C p
Le bilan de quantit de mouvement existe si on tient compte dune perte de charge :
dC i
= J . c . F
d
J est un coefficient de perte de charge,

est un exposant agissant sur le dbit.
Remarque : Il faut vrifier sil ny a pas dquations superflues et envisager alors de rduire la dimension du
systme. Ces points font lobjet des paragraphes suivants.
2.4 Dcomposition dun modle mathmatique

Application aux quilibres physiques
2.4.1 Concepts gnraux
La simulation dune installation ou mme dun appareil nous amne couramment rsoudre des
systmes dquations de dimensions assez importantes. Plutt que de rsoudre celles-ci simultanment, on
a intrt organiser le calcul de manire rsoudre en squence une srie de systmes plus petits.
Lorganisation du calcul se dduit de lexamen dune matrice appele matrice dincidence.
Avant daborder le cas des quilibres physiques et le modle du sparateur, nous allons, partir
dun exemple trs simple, dcrire la procdure suivre :
Soit rsoudre le systme dquations suivant :
a. x 21 + b. x 2 + x 4 = k
(1)
a. x 1 + c. x 12/2 = 1
(2)
d . x 1 + e. x 3 = m
(3)
5. x 1 + 2. x 3 + x 4 = n
(4)
Construisons la matrice dincidence (ou doccurence) du systme :

Var
x1
x2
x3
x4
Eq.1
- 2.12 -
Chapitre 2
Eq.2
Eq.3
Eq.4
On peut voir que si on connat x4, on pourra trouver x1 et x2 au moyen des quations 1 et 2.
Ensuite, avec lquation 3, on obtiendra x3. Lquation 4 nous donnera x4 partir de x1 et x3. La squence
de calcul sera alors la suivante :
Estimer
x4
(1)
x1
(2)
x2
{ {
Correction
sur x4
(3) x
(4) x4
non
Test de
convergence
oui
On emploie donc un processus itratif pour rsoudre le problme. Si on ne corrige pas la valeur de
x4 donne par lquation 4 et que lon rinjecte cette valeur dans les quations 1 et 2, on emploie la
mthode de substitution simple. En gnral cependant, on corrigera cette valeur par les mthodes dcrites
au chapitre 3, ces mthodes acclrent le processus de convergence des calculs.
Si on avait eu le systme suivant :
a. x 21 + b. x 2 + x 4 = k
(1)
a. x 1 + c. x 2 = 1
(2)
d . x 1 + e. x 3 = m
(3)
5. x 1 + 2. x 2 + x 4 = n
(4)
la squence de calcul eut t fort semblable, seule lvaluation de x3 tombe hors du procd itratif et cette
inconnue est calcule lorsquon sort de la boucle.
2.4.2 Application aux quilibres physiques : quilibre L-V
2.4.2.1 Modle thermodynamique
Le modle thermodynamique est constitu :
dune part des quations :
- 2.13 -
Chapitre 2
f (variables, paramtres) = 0 ,
quations de bilans et conditions dquilibre,
dautre part, par un choix cohrent de lois thermodynamiques (quations dtat ou corrlation), fixant les
limites dapplication du modle et par lensemble des donnes ncessaires.
Exemple dquation dtat : lquation de Soave.
Lquation dtat scrit :
P=
On dfinit les variables auxiliaires Z =
R. T
a (T)
V b V.( V + b)
P. V
a. P
b. P
, A=
, les valeurs de Z tant solution de
2 et B =
R. T
R. T
( R . T)
lquation Z3 Z2 + ( A B B2 ). Z A. B = 0 . Avec ces variables, on peut calculer le coefficient de

VZ 1
A Z+B
fugacit puisque ln = Z 1 ln Z
= Z 1 ln( Z B) ln
.
V
B
Z
2.4.2.2 Examen des conditions dquilibre
Les fugacits partielles des constituants i T et P du mlange doivent tre gales dans les deux
phases :
fiL = fiV
Ces variables sexpriment gnralement en fonction du coefficient de fugacit () pour la phase
vapeur et en fonction du coefficient dactivit () et de la fugacit de rfrence ( f * ) pour la phase liquide :
i . yi . P = i . xi . f i*L
Le coefficient dquilibre du constituant est dfini par le rapport de sa fraction molaire en phase
gazeuse sa fraction molaire en phase liquide :
Ki =
yi
xi
Nous verrons plus loin limportance du rle jou par ces variables. Il est ds lors utile den
examiner les diverses expressions possibles :
Ki =
f i*L . i
P. i
Par ailleurs on sait que,

i = i (T, P, y, )
i = i (T, P, x, )
- 2.14 -
Chapitre 2
o les symboles x et y sont des vecteurs de fractions molaires tandis que indices et , sont des
vecteurs de paramtres. Lensemble de ces vecteurs sera reprsent par ci-aprs. Le nombre et les
valeurs des paramtres dpendent des choix adopts pour les lois physico-chimiques.
On peut ds lors crire en toute gnralit :
Ki = f (T, P, x, y, )
Rappelons que f i*L est la fugacit du constituant i pur liquide, la temprature et la pression du
mlange ou encore fugacit de rfrence. Ce terme ne dpend donc que de T et P. Il sagit l de la
convention symtrique. Si le constituant nexiste pas ltat liquide, cest--dire sil est incondensable, la
rfrence est ltat de dilution infinie. Cest la convention asymtrique.
A pression modre, on a :
v S .( P P S )
f i*L = f i*LS .exp i
R. T
f i*LS = PiS . *i S
o
f i*LS
*i S
PiS
fugacit du corps pur saturation la temprature du mlange,

coefficient de fugacit de la vapeur sature dans les mmes conditions,
pression de vapeur saturante du corps i (calcule laide de lquation dtat vapeur).
Si la pression du systme est choisie comme rfrence, les relations ci-dessus ne doivent pas tre
modifies.
Si lon choisit ltat de rfrence de la phase liquide comme tant celui du corps pur la pression
P et la temprature du systme, on a dans ce cas :
r
LR
i
=P
*L
i
= fi
* LS
i
LR
vSi .( P r P S )
. .exp
R. T
*S
i
v Li .( P P r )
.exp
R. T
Le volume molaire du liquide est gnralement considr comme indpendant de la pression, ce qui
justifie les approximations apparaissant dans les formules ci-dessus. De plus, cette hypothse se retrouve
implicitement quand on calcule pour le constituant i. Lavantage de choisir la pression de rfrence gale
0 est que les deux phases possdent le mme tat de rfrence.
On se rend compte que si la temprature rduite (Tr = T/Tc) du constituant i est > 1, on na plus la
possibilit de calculer sa pression de saturation, il faut ds lors prendre une valeur extrapole ou calculer la
fugacit du liquide dune autre manire. Ce cas ne sera pas trait dans le cadre de ce cours (convention
asymtrique, revoir ci-dessus).
- 2.15 -
Chapitre 2
Maintenant, nous pouvons donc supposer que nous avons accs la relation :
Ki = f(T, P, x, y, )
Avant de terminer le bref examen de cette variable, il est utile de rappeler quil existe diverses
approximations qui reprsentent un ensemble de comportements possibles pour la phase vapeur et pour la
phase liquide.
Le tableau suivant rsume ces possibilits.
Solution idale
fiL = xi .
C
o
m
p
o
r
t
e
m
e
n
t
p
h
a
s
e
v
a
p
e
u
r
l
a
n
g
e
i
d
a
l
Gaz parfait
Ki =
L
f i'
Ki =
= yi .
V
f i'
s
Pi . i
K = K = f (P ,T)
Gaz rel
V
fi
L
f Li = i . x i . f i'
P si
i* = *i = 1
Solution non idale
f *i
Ki =
*i . P
Ki =
f*i
. i
* . P
i
f 'i V
= i . P
*
Mlange
non idal
id
K = K = f(P,T)
Impossible
Ki =
f *i L . i
* . P
i
Tableau 2.1
Nous avons not f i'L ou f i'V les fugacits du constituant pur dans les conditions de T et P du
mlange, en phase liquide et en phase vapeur respectivement.
Remarquons que Kr (K Raoult) et Kid (K idal) sont indpendants des compositions des phases en
prsence et quils constituent une assez bonne approximation des valeurs de K. Nous les utiliserons par la
suite, notamment comme point de dpart des itrations. Kr peut se calculer trs simplement et Kid dune
manire un peu plus complique.
- 2.16 -
Chapitre 2
2.4.2.3 Equations du modle des quilibres liquide-vapeur - Spcifications

Les quations dquilibre scrivent :
Ki = f (T, P, x, y, )
avec i = 1, n
(1)
Les paramtres sont fixs par le choix des lois qui conditionnent la connaissance des fugacits
partielles des constituants. Les relations de dfinition des Ki sont :
y i = Ki . x i
yi
Ki
xi =
avec i = 1, n
(2)
avec i = 1, n
(2)
Les bilans de matire scrivent selon les deux quations ci-dessous :

xi =
zi
1 + .( K i 1)
(3)
yi =
Ki . z i
1 + .( K i 1)
(3)
Ces quations drivent des n quations F. z i = V. yi + L. xi avec = V / F , fraction vaporise

(ou bien 1 = L / F , fraction liquide du systme). Les variables F, V et L sont respectivement les
nombres de moles contenues dans le systme entier, dans la phase vapeur et dans la phase liquide. Les
fractions molaires correspondantes sont zi, yi et xi. Celles-ci doivent satisfaire certaines conditions de
normalit :
x
y
1 = 0
(4)
1 = 0
(4)
une de ces deux quations pouvant tre remplace par :
x y
i
=0
(4)
Les donnes du problme dquilibre comportent les variables F et z, ces dernires satisfaisant la
relation z i = 1 . Ds prsent nous considrons que F = 1 mole de mlange et que les zi auront t
norms. On aura donc V = et L = 1 - .
Un dcompte des quations et des variables conduit au tableau suivant :
Equations
Nombres
Variables
Nombre
(1)
(2)
(3)
(4)
n
n
n
2
T, P, x, y
Ki
2n+2
n
1
-
- 2.17 -
Chapitre 2
Total
3n+2
Total
3n+3
Cependant, si lon fait la somme des quations (2) et que lon tient compte de (3), on se rend
compte que les 2 quations (4) ne sont pas indpendantes. Lune delles seulement doit tre retenue ou
bien uniquement lquation (4). Nous noterons cette quation de normalit S(VAR) = 0
La notation VAR dveloppe une des variables considres prcdemment, par exemple , P,
T, etc. Lquation S(VAR) = 0 dsignera lune quelconque des quations (4) ou (4). Nous expliciterons la
fonctionnelle S(VAR) la fin de ce chapitre. Le dcompte stablit donc comme suit :
(3n+1) quations et (3n+3) variables
De plus, on peut dsirer calculer une fonction dtat quelconque - celle que lon considre le plus
souvent dans les quilibres liquide-vapeur est lenthalpie H - soit pour un constituant i dans le systme ou
dans lune des deux phases, soit pour lensemble de ceux-ci.
On aura par exemple :
H = (1 ). x i . H Li + . x i . H Vi
i
!!
(5)
Ces quations introduisent autant de variables nouvelles quil y a de relations de dfinition. On a

finalement :
(3n+2) quations et (3n+4) variables
Il sera donc ncessaire de se fixer deux spcifications supplmentaires. Selon le choix effectu, on
distingue les calculs dquilibre (ou de vaporisation partielle dsigne souvent par le terme anglais flash) :
T et P fixs,
P et fixs, cas particuliers : les calculs de T de rose ( = 1) et T de bulle ( = 0),
T et fixs, cas particuliers : les calculs de P de rose ( = 1) et P de bulle ( = 0),
P et H fixs,
P et S fixs.
2.4.2.4 Rsolution des quations du modle. Pression modre
2.4.2.4.1 Examen de la matrice doccurence du systme
Constituons cette matrice dans le cas de 2 constituants :

Inc.
Eq.
(1)
K1
x
K2
x1
x2
y1
y2
- 2.18 -
Chapitre 2
x
(2)
x
x
x
(3)
x
x
x
x
x
x
x
x
x
x
(4)
(5)
x
x
x
x
x
x
x
x
x
On remarque que seule la relation (5) permet de calculer H. Le systme rsiduel (1) (4) constitue
une matrice irrductible de (3n+1)(3n+3) (entoure en pointill). La rsolution dun tel systme
dquations non linaires posera rapidement un problme (par exemple, si n > 5).
Des mthodes existent bien sr (Newton - Raphson, Marquardt, etc.), mais elles ncessitent une
procdure itrative et sont sujettes bien des alas (dpendance par rapport lapproximation de dpart,
relaxation,...). De plus, elles demandent la connaissance des drives partielles et le remplissage des
matrices (Jacobien par exemple), en fonction du problme traiter.
2.4.2.4.2 Algorithmes de calcul
Examinons ds lors si nous ne pouvons pas rduire la dimension du problme dans les (3n+1)
quations. Pour ce faire, nous choisirons le problme T et P donns (EPT). Nous avons vu que nous
pouvons disposer dune approximation du coefficient K du constituant i par :
Kr (T, P) ou Kid (T, P)
et ds lors, si lon estime , on pourra calculer x(x) par les relations (3), ensuite y(y) par les relations (2) et
enfin les coefficients K de chaque constituant i par les relations (1). Cette procdure peut tre schmatise
par le diagramme de fluence suivant :
K
(3)
(2)
x
(1)
y
ou en simplifiant les notations :

Kr , x, y Kr
Nous avons not x, y car suivant le type de problme trait, on devra utiliser les premires ou les
deuximes expressions de (2) et (3).
Les variables K dpendant en gnral de x et y, nous devrons itrer sur K puis vrifier que la
condition de normalit est satisfaite et modifier en consquence On aura donc deux boucles de
convergence imbriques :
- 2.19 -
Chapitre 2
T et P :
x, y
K constante ?
?
S() = 0
On comprend maintenant lexpression S(), cest la valeur de :
x y
i
obtenue pour une valeur de lorsque les valeurs de K du constituant i sont convergence. On va se servir
de cette grandeur et dautres ventuellement pour dterminer la valeur corrige de avec laquelle on
espre satisfaire la condition (4). Il faut remarquer que la dimensionnalit a t rduite de manire
spectaculaire puisquon ne doit rsoudre quune quation la fois.
Diverses questions subsistent :
Comment calculer la valeur corrige de VAR ( corrig dans ce cas) en vue den promouvoir la
convergence ?
Ne peut-on promouvoir aussi la convergence de K ?
Quelles sont les mesures de "scurisation" de la mthode, cest--dire telles quelles assurent un
maximum de chances dobtenir le rsultat ?
Est-on sr par exemple que 0 < < 1 ?
La forme S(VAR) = 0 est-elle toujours la mme dans tous les cas ?
Quels sont les tests de convergence ?
Nous y rpondrons aprs avoir dress les organigrammes simplifis ditration correspondant aux
autres types de calcul dquilibre.
Examinons le cas o T et sont donns (ETA).
T et :
P K
x, y
K constante ?
On remarque sans peine la similitude des diagrammes.

On a ensuite pour P et fixs (EPA) :
- 2.20 -
?
S(P) = 0
Chapitre 2
T Kr
P et :
x, y
K constante ?
?
S(T) = 0
Le cas suivant est diffrent puisquil fait intervenir lenthalpie (EPH). Si on se reporte la matrice
doccurence, on se rend compte quil est ncessaire de passer par une variable auxiliaire. On a le choix
entre et T. Ces variables peuvent donc tre inverses ci-aprs. Dans la littrature, on prfre le choix
illustr sur le diagramme ci-dessous pour des raisons jamais voques.
H et P :
x, y
K constante ?
?
S(T) = 0
?
H() = 0
La condition H() = 0 est dcrite de manire que lenthalpie du systme soit celle qui a t
impose et appele Hd :
H() = H(T, P, x, y, ) - Hd = 0
Si lon avait rsoudre un quilibre T et H donns, la rsolution serait semblable la prcdente,
les itrations portant sur la pression ( condition bien entendu que lenthalpie dpende de la pression, ce
qui nest pas le cas des gaz idaux et parfaits).
2.4.2.4.3 Critre de convergence, rsolution et acclration
A laide du mode PT, dcouvrons tout dabord un genre de difficult, qui peut se poser. On a le
choix entre trois expressions pour S() :
S1 ( ) = x i 1 =
i
S2 ( ) =
1=
zi
1 + .(K
i
Ki . z i
1 + .(K
i
S3 ( ) = x i y i =
i
1)
1)
(1 K i ). z i
i 1)
1 + .(K
i
Pour des valeurs de Ki donnes, ces trois fonctions voluent comme lindique la figure 2.3 cidessous :
- 2.21 -
Chapitre 2
S( )
S1( )
S ( )
2
0.5
S2
0
S1
S
- 0.5
Ethane
n-butane
Z1 = Z2 = 0.5
-1
P = 5 bar
T = 280 K
Figure 2.3
Dans ce cas, seule la fonction S est unimodale et elle sera donc prfre. Si lon a choisi S1 ou S2,
il existe deux solutions et avec certaines mthodes de rsolution, on peut tre conduit vers = 0 ou vers
= 1.
Le choix du critre tant fait, comment trouver la solution de S() = 0 ? On peut songer utiliser la
mthode de Newton, la mthode regula-falsi , la mthode de Wegstein ou encore une mthode
dinterpolation parabolique... Ces mthodes seront vues chapitre 3.
Certaines mthodes de convergence peuvent conduire une valeur de en dehors du domaine 0
< < 1. Avant daccepter la valeur fournie par lalgorithme de convergence, on testera cette valeur, si
0 alors = 0 ou si 1 alors = 1. Cest une mthode de scurisation.
Les tests de convergence sont de deux types :
S( k +1 ) S( k ) < 1
k +1 k < 2
Ces deux tests devront tre satisfaits avant de terminer le calcul.
Enfin, la promotion de la convergence de la boucle interne doit porter sur chaque Ki. Lexprience
dmontre que les itrations sur cette boucle convergent rapidement par simple substitution.
- 2.22 -
Chapitre 2
Il faut tre plus prudent dans le cas des systmes fortement non idaux o des quations trs
complexes (Soave, NRTL, UNIQUAC,...) sont utilises. Il sera indispensable dutiliser des mthodes
comme celles de Marquardt ou de Broyden qui seront vues dans la suite du cours.
- 2.23 -
Chapitre 3
CHAPITRE 3
LES PRINCIPALES METHODES MATHEMATIQUES
Introduction
Si nous voulons rsoudre les modles que nous venons de dfinir au chapitre prcdent, nous
devrons disposer de mthodes de rsolution de systmes dquations simultanes, dintgration numrique
et de mthodes sappliquant au calcul itratif.
Dans ce chapitre, nous analyserons quelques-unes des plus importantes mthodes telles que celles
de Newton-Raphson, de Wegstein, de Rubin et de Runge-Kutta.
3.1 Rsolution de problmes 1 dimension

Nous allons envisager des mthodes pour rsoudre deux types dquations : les quations
explicites f ( x ) = 0 et les quations implicites x = f ( x ) . Pour le premier type, nous dcrirons les mthodes
de Newton-Raphson et de la corde tandis que pour le second type, nous dcrirons la mthode de
Wegstein.
3.1.1 Mthode de Newton-Raphson
3.1.1.1 Description
Nous devons rsoudre lquation f ( x ) = 0 , connaissant une valeur approche x0 de la solution
x* . Un dveloppement en srie de Taylor de f ( x) donne :
f ( x) = f ( x
) + ( x x ). f ( x ) +
0
(x x0 )2
. f ( x 0 )+K
2!
f , f , etc. sont les drives premires, secondes,... de la fonction f.
Lapproximation de Newton consiste ngliger les termes dordre suprieur un en supposant que
x est suffisamment prs de la solution x* et remplacer lquation initiale par lquation linaire
approche suivante :
0
f ( x 0 ) + ( x * x 0 ). f ( x 0 ) f ( x * ) = 0
On en tire une valeur approche de la racine :
x =x
1
f (x0 )
f ( x 0 )
- 3.1 -
Chapitre 3
cette valeur devant tre une meilleure estimation de la solution (ce nest pas toujours le cas !!!). Pour
trouver la solution, on utilisera donc la formule de rcurrence :
f (xn )
x n+ 1 = x n
f ( x n )
3.1.1.2 Interprtation gomtrique

y
Mn+1
x* xn+1 x n
Figure 3.1
Si on reprsente la courbe y = f ( x ) dans le systme daxes Oxy, lquation de la tangente cette

courbe au point M n dabscisse xn est :
y f ( x n ) = f ( x n ). ( x x n )
Lintersection de cette tangente avec laxe des x ( y = 0) vaut, en supposant f ( x n ) 0 :
xn
f (xn )
f ( x
= x n +1
On obtiendra de mme xn+2 , xn+3 , de proche en proche en menant les tangentes aux points Mn+1 ,
Mn+2 , ... et en cherchant leur intersection avec Ox.
Cest pourquoi la mthode de Newton est aussi appele mthode de la tangente.
3.1.1.3 Remarques
Notons que les valeurs successives des x sont trouves en tenant compte des signes de f ( x) et de
f ( x) . Dans ce cas, pour viter toute divergence, on peut utiliser une seconde formule de rcurrence
inspire de Wegstein (voir plus loin) :
~
x n+ 1 = q. ~
x n + (1 q ). x n+1
(relaxation)
Cela signifie dans lexemple donn ci-dessous, que, sur la tangente passant par M 0 , lon va sarrter en
un point dfini par la valeur de q.
- 3.2 -
Chapitre 3
La mthode de Newton-Raphson a pour inconvnient le calcul des drives.

Si la fonction f ( x) prsente des extremums, il peut arriver que la mthode soit divergente.
Exemple :
y
M1
M1
~1
x2
x*1 x 0 x*2
x
x1
x
M0
Figure 3.2
Pour acclrer la convergence, au lieu dassimiler la courbe en un point sa tangente, on peut

prendre une courbe ausculatrice du second degr : cest la mthode de Richmond.
3.1.2 Mthode de la corde ou de Regula-Falsi
Soit rsoudre lquation f ( x) = 0 ,
y
f(x k )
x k+1
0
k+1
f(x
f(x
k-1
x k-1
xk
)
)
Figure 3.3
Si on ne dispose pas de lexpression analytique de f ( x ) , on pourra choisir la mthode de la corde

qui sapparente celle de Newton. La tangente tant remplace par la corde.
Lquation de la corde est :
f ( x) f ( x k ) =
le point x
k +1
sobtient en posant f ( x
k +1
f ( x k ) f ( x k 1 )
k 1
.( x x k )
x x
) = 0 , do lquation :
k
- 3.3 -
Chapitre 3
x k +1 = x k 1. f ( x k )
avec :
=
f ( x k ) f ( x k1 )
x x
k
k 1
Cette quation diffre de celle de Newton par la valeur de qui nest quune approximation de
linverse de la drive de f au point x k .
Remarques :
Avantages : pas de calcul de drives, robustesse identique Newton-Raphson.
Inconvnients : mthode plus lente prs de la solution + inconvnients de Newton-Raphson.
3.1.3 Mthode de Wegstein
Cette mthode est trs intressante dans le cas frquent en chimie o on est amen rsoudre des
quations implicites.
Soit rsoudre lquation ( x) = x f ( x) = 0 .
Graphiquement, la solution est donne par lintersection de la bissectrice y = x avec la courbe
y = f ( x) .
La mthode est galement base sur un processus itratif. Admettons que lon ait obtenu une valeur
de x aprs i itrations. Pour obtenir la valeur suivante, on utilise la mthode dcrite ci-dessous :
y
x k+2
x k+1
x k+1
1- q
~xk+1
xk
xk
Figure 3.4
1. intersection droite x = x k avec courbe y = f ( x) ( x k , y k ) ,
2. intersection droite y = y k avec bissectrice y = x ( x k+1 , y k ) ,
3. intersection droite x = x k +1 avec courbe y = f ( x) ( x k+ 1 , y k +1 ) ,
- 3.4 -
Chapitre 3
en continuant de cette manire, on utilise la mthode dite implicite qui peut diverger. De plus, la
convergence est trs lente.
Wegstein propose une autre procdure pour viter certains dsagrments de la mthode implicite
et pour acclrer la convergence.
A cet effet, il propose de corriger la valeur de xk +1 en traant la corde entre les points k et k+1 et
en cherchant son intersection avec la droite y = x .
Pour obtenir cette nouvelle valeur de xk +1 , il faut donc rsoudre le systme suivant :
y k +1 y k
y y = k +1 ~ k . ( x ~
x k ) = .( x ~x k )
x x
y= x
k
avec , le coefficient angulaire de la corde,

y
y k+1
yk
0
x k+1
~
x k+1
xk
Figure 3.5
ce qui revient chercher q de sorte que :
~
x k +1 = q . ~x k + (1 q). x k +1
La corde passera par le point ( ~
x k +1 , ~x k+ 1 ) et comme y k = x k +1 , on aura :
~
x k +1 x k +1 = . ( ~
x k +1 ~
xk )
do,
1
. ( x k+1 . ~
xk )
1
1
par consquent, on aura q =
et 1 q =
.
1
1
~
x k +1 =
Dans le cas de la mthode implicite, suivant la valeur du coefficient angulaire de la courbe (tangente
de la courbe), on peut converger de diffrentes manires ou mme diverger.
A ce propos, jetons un rapide coup doeil sur le graphique suivant :
- 3.5 -
Chapitre 3
-1
-
x
Figure 3.6
Si , le coefficient angulaire de la corde varie de + , nous pourrons observer, pour la

mthode implicite, les phnomnes suivants (cfr. deux cas en graphique page 3.4) :
-1
divergence
oscillatoire
0
convergence
oscillatoire
1
convergence
monotone
divergence
monotone
Wegstein par contre, converge mme dans des cas o la mthode implicite diverge. Cependant,
daprs le graphique ci-aprs, on peut voir que q peut prendre des valeurs absolues trs grandes pour 0,75
< < 1,25.
On aura intrt limiter la valeur absolue de q pour restreindre lavancement et ne pas sloigner
trop de la solution.
Le graphique ci-dessous est explicite ce sujet.
-
q = --------1-
0.75
0
1 1.25
1
Figure 3.7
- 3.6 -
Chapitre 3
Procdure gnrale de calcul :

1. on choisit une valeur approche de la solution x 0 , et ensuite on calcule les deux premiers points :
x1 = f ( x 0 )
x 2 = f ( x1 )
2. on calcule q pour corriger x 2
x 2 x1
= 1
do q
x x0
3. on calcule successivement :
~
x 2 = q. x 1 + (1 q ). x 2
x 3 = f ( ~x 2 )
x3 x2
= ~2
x x1
x4 = f (~
x 3)
K
x k +1 x k
= ~ k ~ k1
x x
x k +2 = f ( ~x k+1 )
d ' o ~x 3
x k +1 x k
ou bien x k x k 1
K
et ainsi de suite jusqu la solution. Il est clair qu chaque itration un test est fait pour voir si la solution est
atteinte.
3.2 Rsolution de problmes n dimensions

Nous allons maintenant aborder la rsolution de problmes n dimensions. Pour les quations
explicites, nous envisagerons les mthodes de Newton-Raphson et de la scante gnralises, ainsi que la
mthode de Broyden. En ce qui concerne les quations implicites, nous dcrirons la mthode de Rubin.
3.2.1 Gnralisation n dimensions de Newton-Raphson
Soit rsoudre F( x ) = 0 o F est une fonction vecteur, cest--dire, selon les conventions
mathmatiques :
f1 ( x ) = f1 ( x 1 , x 2 , K, x n ) = 0
f 2 ( x) = f 2 ( x1 , x 2 ,K , x n ) = 0
f n ( x) = f n ( x 1 , x 2 ,K, x n ) = 0
- 3.7 -
Chapitre 3
On choisit cette fois, un vecteur de dpart x 0 = x 01 K x 0n

t
et par analogie avec le cas une
dimension, on limite le dveloppement de Taylor n dimensions aux termes des drives premires, ce qui
donne :
f i ( x) = f i ( x ) +
0
j =1
o x j = x j x
f i
.x j = 0 ; i = 1,K , n
x j x 0
0
j
Sous forme matricielle, ce systme dquations peut scrire :

f 1
x
1
M
f
i
x1
f n
x1
f 1
x n
0
O
M
M x 1 f1 ( x ) f 1 ( x)
L L L

f i
f i
0
L
L
. x + f ( x ) = f i ( x)
x j
x n i i

L L L
M
O
M x f ( x 0 ) f ( x)
n
n n
f n
f n
L
L
x j
x n
o les drives partielles sont calcules au point x0 ; cela scrit, avantageusement, sous la forme
condense (et conventionnalise) ci-dessous :
t
f1
x j
J .x + F( x 0 ) = F( x )
On obtient donc,
x1 = x 0 t J x10 F( x 0 )
et pour le (n+1)me terme :
x n+ 1 = x n t J x1n F( x n )
la matrice jacobienne (transpose) tant calcule en x n .
Remarques :
la mthode ncessite une inversion de matrice, il faut donc que la matrice jacobienne J soit non
singulire ;
autre inconvnient, chaque itration, toute la matrice J doit tre recalcule, ce qui fait chaque fois nn
drives, n fonctions et linverse dune matrice (nn) ;
on peut encore appliquer ici un processus de convergence sinspirant de Wegstein :
- 3.8 -
Chapitre 3
~
x n+1 = ~
x n t J x1n F( ~
x n ). (1 q )
cependant, au dbut, il y avantage choisir q grand et le diminuer par la suite. En fait, loin de la
solution, il faut avancer petits pas, tandis que prs de la solution, on peut progresser plus rapidement.
3.2.2 Mthode de la scante gnralise
La mthode de Newton dcrite au point 3.2.1 prsente le dsavantage de devoir recalculer le
jacobien chaque itration. Il est cependant possible de dvelopper un algorithme qui renouvellera le
jacobien dune manire lgrement diffrente. Nous allons dans un premier temps dcrire la mthode de la
scante gnralise avant daborder la mthode de Broyden, deux mthodes similaires.
Soit rsoudre F( x ) = 0 de dimension n. Ce systme peut contenir des quations implicites,
puisque ces dernires peuvent toujours tre ramenes sous une forme explicite. On doit trouver le vecteur
x tel que e = F( x) = 0 .
Nous pouvons tenter de rsoudre le systme en trouvant les fonctions linaires qui approximent nos
fonctions non-linaires. Nous pourrons alors rechercher les zros de ces fonctions linaires et esprer quils
sont proches des zros du systme non-linaire.
Supposons que pour les quelques essais suivant :
x0
e0 = F( x 0 )
x1
e1 = F( x1 )
x2
e2 = F( x2 )
en = F( xn )
xn
nous puissions crire un modle linaire qui permet de reproduire ces donnes :
e* = Ax + b
et donc
e0 = Ax 0 + b
e1 = Ax1 + b
M
en = Axn + b
Nous avons nn inconnues pour la matrice A et n inconnues pour le vecteur b .
- 3.9 -
Chapitre 3
Pour rsoudre ce systme, nous avons (n+1)n quations linaires, et donc la solution peut tre
trouve.
Nous pouvons en premier lieu liminer le vecteur b en crivant :
e1 = e1 e0 = A ( x1 x0 ) = A x1
e n = en e n1 = A ( x n x n1 ) = A x n
ce qui nous permet dobtenir le systme :

E = A X
Et donc
A = E X 1
nous donne A pour autant que la matrice X aie un inverse. Puisque X est entirement sous notre
contrle (nous choisissons les x ), nous pouvons garantir un inverse. Les valeurs pour b peuvent tre
trouve avec
b = en Axn
Nous pouvons maintenant estimer la valeur suivante de x comme tant celle qui rend e* = 0 .
Cette tape induit la rsolution dun systme dquations linaires :
e* = 0 = Ax n+ 1 + b
ce qui donne
x n+ 1 = A 1b = A 1 e n Ax n = x n A 1 e n
Amlioration la mthode de la scante gnralise
Lapproche reprise ci-dessus implique qua chaque tape il nous faut recalculer la matrice A
partir des donnes originales. Il est possible de dvelopper des algorithmes qui modifient A dune autre
manire.
Supposons que nous disposons dune estimation de A , appelons-la A1 . Supposons que nous
avons valu e1 pour un x1 dtermin. Nous aimerions obtenir une matrice A2 telle que
e1 = A2 x1
Pour obtenir A2 , nous pouvons modifier A1 de la faon suivante :
A2 = A1 + u1 t v1
- 3.10 -
Chapitre 3
Nous obtenons ainsi :
e1 = A1 + u1 t v1 x1
Nous pouvons choisir v 1 arbitrairement et obtenir une formule donnant u1 :
u =
1
e 1 A 1 x1
t
v 1 x 1
Tant que v 1 nest pas orthogonal x1 , nous avons un vecteur u1 bien dtermin et donc nous
pourrons obtenir une matrice A2 partir de la matrice A1 .
Tentons maintenant dobtenir un second point, e 2 pour ltape x 2 o x 2 peut tre par
( )
exemple le rsultat de : x 2 = x 1 A 2
e 1 . On a alors x 2 = x 2 x 1 . On peut alors chercher avoir :
e 2 = A 3 x 2 = A 2 + u 2 t v 2 x 2
On peut aisment trouver u 2 et v 2 comme auparavant, mais maintenant nous voudrions galement
que A 3 satisfasse notre premier point, cest--dire :
e 1 = A 3 x 1 = A 2 + u 2 t v 2 x 1
Si u 2 t v 2 x 1 tait gal au vecteur zro, alors lquation ci-dessus se rduirait lquation de

dfinition de A2 et serait donc vrifie. Nous devons donc choisir v 2 orthogonal x1 . Une solution est
de choisir un vecteur v 2 gal au vecteur x 2 orthogonalis par rapport x1 . Dans ce cas, lquation
pour obtenir v 2 est :
v 2 = x 2
v2
x1 x 2 x1
t
x x1
1
x 2
x1
Figure 3.8
Nous obtenons alors un autre point, e 3 pour ltape x 3 et imposons encore que :
- 3.11 -
Chapitre 3
e 3 = A 4 x 3
avec
e 1 = A 4 x1
e 2 = A 4 x 2
Pour A 4 = A 3 + u 3 t v 3 ,nous trouvons que v 3 devrait tre orthogonal x1 et x 2 . On peut
gnraliser ce rsultat :
i +1
=A
(e
+
A i x i
t i
v
v x
o v i est orthogonalis par rapport aux (i-1) prcdents x i .
3.2.3 Mthode de Broyden

Nous allons maintenant dvelopper lalgorithme de Broyden. Il existe une trs grande similitude
entre ce dernier et la mthode de la scante gnralise puisquune seule tape de lalgorithme les
distinguent. Cependant la mthode de Broyden est meilleure dans son utilisation du temps et de lespace
mmoire de lordinateur.
Au lieu de slectionner une nouvelle matrice de faon ce que les tapes prcdentes soient
satisfaites, lalgorithme de Broyden requiert que :
A i +1 z = A i z
pour lensemble des vecteurs z orthogonaux x i . Lavantage est que les diffrents x i prcdents ne
devront plus tre stocks puisque toute linformation sera conserve dans la matrice A . Nous avons
comme prcdemment :
e i = A i+1 x i = A i + u i t v i x i
Nous avons alors :
[A
+ ui t v i z = A i z
pour tout vecteur z 0 orthogonal x i . Donc si nous choisissons v i = x i , alors nous aurons
t i
v z = t x i z = 0 si et seulement si z 0 est orthogonal x i . Avec ce choix, nous aurons :
u =
i
e i A i x i
t
xi x i
et
(e A x )
=A +
i
i +1
x i x i
x i
- 3.12 -
Chapitre 3
Lavantage de la mthode de Broyden sur la mthode de la scante gnralise est, comme dit
prcdemment, quil nest pas ncessaire de conserver les vecteurs x i . Nanmoins le dsavantage est
que la convergence est linaire alors que la mthode de la scante converge quadratiquement.
3.2.4 Mthode de Rubin (quasi Newton)
Dans les divers modles mathmatiques examins prcdemment, nous avons vu que certains
dentre eux se prsentaient sous la forme :
x k +1 = ( x k )
o est un oprateur et pas ncessairement une fonction vecteur analytique (ici nous prendrons le cas o
cest un vecteur, donc ).
Dans les calculs de flowsheet, cette situation se prsente frquemment (voir chapitre 1).
Contentons-nous ici den donner un exemple simplifi :
7
6
x k+1
DIV
xk
MIX
3
REACTEUR
S
E
P
5
Figure 3.9
Le procd comporte un racteur R, un sparateur S, un diviseur D et une jonction J (un

mlangeur). Le flux dentre 1 est connu, les modles des appareils galement et il faut calculer les flux de
sortie 5 et 7, ce qui ncessite dailleurs de calculer tous les flux internes la boucle 2, 3, 4 et 6. On devra
donc effectuer une itration sur les valeurs des paramtres du flux 6 contenues dans le vecteur x.
Au fur et mesure quon obtient des valeurs pour x, le problme se pose de savoir comment on
peut amliorer la convergence, cest--dire (tout comme dans la procdure de Wegstein), comment arriver
plus vite (ou tout simplement arriver) la solution.
La mthode est celle de la corde gnralise, mais applique la forme particulire de la fonction
vecteur rsoudre :
x ( x) = 0
o x et sont respectivement un vecteur et une fonction vecteur. On peut sattendre devoir estimer un
Jacobien de cette fonction.
- 3.13 -
Chapitre 3
Si ( x) tait analytique, on aurait :

f ( x ) = x ( x) = 0
alors en appliquant la mthode de Newton, on a :
x k +1 = x k
{ J [f( x ) ]}
t
f (xk )
ou
]} [ x
1
x k +1 = x k E t J ( x k )
(xk )
par consquent, la recherche dune approximation du Jacobien de devient lobjectif. Ecrivons les
dveloppements de Taylor de ces fonctions de rcurrence, mais limits au premier terme :
k

1( x) = 1( x ) + 1 .( x1 x1k ) + K + 1 .( x n x kn )
x1
x n
k
L
j
j
.( x1 x1k ) + K +
.( x n x kn )
j ( x) = j ( x ) +
x1
x n
k
L
n (x) = n (x
n

.( x1 x1k ) + K + n .( x n x kn )
+
x1
x n
k dsigne le nombre ditrations,

j indique la composante,
n est la dimension du problme (nombre de variables).
Cela scrit avantageusement sous forme vectorielle,
( x) = ( x k ) +
{ J [(x )]}x
t
Au point de vue rsolution, on peut connatre les vecteurs x1 , x2 , x 3 , K , x k , aprs k itrations,
ainsi que ( x k+1 ) ( x k ) , en effet, la mthode implicite peut nous donner (on suppose tre parti dun
point x0) :
x kj = j ( x k 1 )
(composante j, itration k)
Cependant, pour un systme de n quations, on a n2 inconnues qui sont les drives de par
rapport x (pour tout x et pour tout ). Aprs une itration, on a donc n quations, aprs n itrations, on
obtient n2 quations et on peut donc obtenir lexpression du Jacobien.
Pour obtenir celui-ci, nous utiliserons deux matrices C et D pour stocker les informations donnes
par chaque itration.
- 3.14 -
Chapitre 3
Sous forme matricielle, le systme rsoudre est :

( x) ( x k ) =
{ J [(x )]}x
t
o, sous forme condense,

t
B = t J t A ou
B t A 1 = t J , d ' o J = A 1 B
avec
x11 x01
M
k
0
A = x 1 x 1
M
xn x0
1
1
L x1j x0j
O
M
L xkj x0j
M
L x x0j
n
j
1 ( x1 ) 1 ( x0 )
B = 1 ( x ) 1 ( x0 )
( x n ) ( x 0 )
1
1
L x1n x0n
M
L xkn x 0n
O
M
L xnn x0n
L j ( x 1 ) j (x 0 ) L n ( x 1 ) n ( x 0 )
O
M
M
L j ( x k ) j (x 0 ) L n ( x k ) n ( x 0 )
M
O
M
n
0
n
0
L j (x ) j (x ) L n (x ) n ( x )
Lorsque les matrices A et B seront construites, on obtiendra le Jacobien J en inversant A et en

la multipliant par B .
Dun point de vue calcul numrique, pour obtenir la matrice Jacobienne transpose, on effectue n
itrations par substitution simple et on stocke linformation dans 2 matrices C et D :
x10
1
x1
M
C= k
x 1
M
n
x1
L x 0j
L
O
x1j
M
L xkj
M
L x nj
1 ( x0 )
1
1( x )
M
D=
k
1 ( x )
M
1 (x n )
L x0n
L x1n
M
L xkn
O M
L xnn
L j ( x 0 ) L n (x 0 )
L j ( x1 ) L n ( x1 )
O
M
M
L j ( x k ) L n ( x k )
M
O
M
L j (x n ) L n ( xn )
Ces deux matrices ont chacune (n+1) lignes et n colonnes. Nous obtiendrons alors A et B en
soustrayant dans chaque matrice la premire ligne toutes les autres.
- 3.15 -
Chapitre 3
Cependant, le Jacobien ainsi calcul est celui au point x 0 , il est prfrable de calculer celui-ci aprs
litration n. Pour cette raison, les premires informations obtenues seront stockes en commenant par la
dernire ligne (cest comme si les matrices C et D taient retournes par rapport lhorizontal).
( )(
On a donc J = A1 B au point x n , la (n+1)me itration donne xn +1 , par x n+ 1 = ( x n ) , on

corrige alors cette valeur de la faon qui suit,
)} [ x
t
~
x n+ 1 = x n E A 1 B
(xn )
La mthode de Rubin est donc une mthode de Newton-Raphson qui a t modifie pour
permettre lestimation du Jacobien car dans le cas du flow-sheet, la fonction ( x) nest pas connue
analytiquement.
Cette mthode sappelle aussi mthode de la corde gnralise. En effet, le calcul du Jacobien
laide des matrices A et B revient estimer la pente de la corde entre deux points.
j( xk ) j( x n )
= pj
xkj x nj
Remarque :
Afin dacclrer encore la convergence, on appliquera la mthode de Wegstein sur la variable la
plus sensible (celle qui a t la plus modifie).
3.3 Rsolution dquations diffrentielles :

La mthode de Runge-Kutta
Le but de la mthode de Runge-Kutta est dobtenir une solution numrique approche dune
quation diffrentielle ordinaire non linaire, ou dun systme dquations diffrentielles ordinaires non
linaires.
Considrons dabord le cas dune quation, on donne :
dy
= y = f [ x , y( x) ]
dx
avec comme conditions initiales :
y( x 0 ) = y 0
x 0 , point initial.
On dispose donc de la pente de la fonction y, ce qui va nous permettre destimer la valeur de y en
un point x 0 + h , h tant lintervalle dintgration.
- 3.16 -
Chapitre 3
y
n+1
y
y
n+1
Figure 3.10
Soit y n , lestimation de la fonction y aprs n pas dintgration,
y n = y( x 0 + n. h)
et si on dveloppe en srie de Taylor, nous obtenons :
y( x n+ 1 ) = y( x n + h)
h2
h3
n
= y( x ) + h. y ( x ) +
. y ( x ) +
. y ( x n )+K
2!
3!
n
centr au point x n
(1)
Connaissant y = f ( x, y) , on va essayer dobtenir une approximation de y , y ,K , toutes ces

drives tant prises au mme point, on a :
y =
dy df f y f
=
=
. +
= fy . f + f x
dx dx y x x
y =
dy d
y
=
fy . f + f x =
f y .f + f x . +
f .f + f x
dx
dx
y
x x y
de mme, on obtiendra :
= f xx + 2. f . f xy + f 2 . f yy
(
)
+ f .( f . f + f )
y
o lon a not :
f
2f
, f xy =
,
x
x. y
dy y
et en remarquant que dans ce cas,
=
dx x
fx =
On reporte les valeurs de y (et y ventuellement) dans (1) et lon obtient :

y
n +1
=y
n +1
h2
y = h. f +
. f + fy. f
2! x
n
(1a)
On recherche alors une expression de y n+1 o ninterviennent que des estimations de la fonction f
en des points particuliers,
soit
y n +1 = y n + N 0 . k 0 + N 1 . k 1 + N 2 . k 2 +K
- 3.17 -
Chapitre 3
k 0 = h. f ( x n , y n )
o
k 1 = h. f ( x n + 1. h, y n + 10 . k 0 )
k 2 = h. f ( x n + 2 . h, y n + 20 . k 0 + 21 . k 1 )
k 3 = h. f ( x n + 3 . h, y n + 30 . k 0 + 31. k 1 + 32 . k 2 )
Si on dveloppe k 1 en laissant tomber le second ordre, on obtient :
k1 = h. f + 1. h. f x + 10 . k 0 . f y
do :
y n +1 = N0 . h. f + N1. h. f + 1. h. f x + 10 . h. f . fy
(2a)
En comparant les termes en h de (1a) et (2a), on trouve :

termes en
" "
"
"
h. f
2
h .f x
2
h .f y . f
N 0 + N1 = 1
N1 .12 =
N1 . 1 =
On a donc trois relations pour quatre paramtres. Il existe donc une infinit de solutions que lon
peut crire sous la forme :
N0 = 1
N1 =
1
2. 1
1
2. 1
10 = 1
En prenant 1 = 1 , nous aurons N1 =
, N0 =
, 10 = 1 , do :
k 0 = h. f ( x n , yn )
k1 = h. f ( xn + h , y n + k 0 )
ce qui fait que :
y n+ 1 = y n +
1
.( k 0 + k 1 )
2
Au troisime degr, nous aurons 8 paramtres pour 6 relations, au quatrime degr, 13 paramtres
pour 10 relations,... La mthode que nous utiliserons est celle qui a t dveloppe par Runge qui est une
mthode du 4me degr. On a le tableau suivant :
RUNGE
N0
N1
N2
N3
10
20
21
30
31
32
1/6
2/6
2/6
1/6
1/2
1/2
1/2
1/2
- 3.18 -
Chapitre 3
KUTTA
1/8
3/8
3/8
1/8
1/3
2/3
1/3
1/3
-1
Ce qui fait que lon a :

k 0 = h. f ( x n , y n )
k1 = h. f ( x n + h / 2, y n + k 0 / 2)
k 2 = h. f ( xn + h / 2, yn + k1 / 2 )
k 3 = h. f ( x n + h, y n + k2 )
et
y n +1 =
1
.( k 0 + 2 . k 1 + 2 . k 2 + k 3 )
6
Si lon a plusieurs quations diffrentielles, le dveloppement sera identique mais il sappliquera aux
fonctions vecteurs.
Dans ce cas cependant, on aura intrt pondrer les quations afin que les variations de chacune
delles soient du mme ordre de grandeur. Le systme rsoudre sera alors :
dy1
dx
dy 2
. f2 ( x) =
dx
dy
. f 3( x) = 3
dx
. f1( x ) =
avec , , comme coefficients de pondration.
- 3.19 -
Chapitre 3
3.4 Applications numriques

Nous allons maintenant, au moyen de quelques exemples numriques, illustrer lutilisation de certains des
modles que nous avons dcrits dans les pages prcdentes.
3.4.1 Comparaison entre la substitution simple et la mthode de Wegstein.
Premier cas : convergence monotone
x=a+bx+cx2
Substitution
iter
x
1
2
3
4
5
6
7
8
9
10
0
1,4
1,4
1,7687
1,7687
1,969
1,969
2,0959
2,0959
2,1828
2,1828
2,2454
2,2454
2,2918
2,2918
2,3271
2,3271
2,3544
2,3544
Wegstein
iter
x
1
2
3
4
5
6
0
1,4
1,4
1,9005
1,9005
2,246
2,246
2,3994
2,3994
2,4538
2,4538
2,4632
a
b
c
f(x)
1,4
0,0417
0,1583
3,5
3
1,4
1,4
1,7687
1,7687
1,969
1,969
2,0959
2,0959
2,1828
2,1828
2,2454
2,2454
2,2918
2,2918
2,3271
2,3271
2,3544
2,3544
2,3758
2,5
2
1,5
1
0,5
0
0
0,5
1,5
2,5
3,5
0,5
1,5
2,5
3,5
3,5
3
2,5
f(x)
1,4
1,4
1,7687
1,9005
2,051
2,246
2,2923
2,3994
2,4115
2,4538
2,4556
2,4632
psi
x amlior
2
0,2633
-0,357 1,9005
0,5642
-1,295
0,6982
-2,313 2,3994
0,7772
-3,488 2,4538
0,8101
-4,266 2,4632
2,246
1,5
1
0,5
0
- 3.20 -
Chapitre 3
Deuxime cas : convergence oscillatoire (uniquement pour la substitution simple)

x=a+bx+cx
a
b
c
f(x)
Substitution
iter
x
1
0
3
3
0,9
0,9
1,992
1,992
1,204
1,204
1,7247
1,7247
1,3528
1,3528
1,6074
1,6074
1,4271
1,4271
1,5521
1,5521
2
3
4
5
6
7
8
9
10
Wegstein
iter
x
1
2
3
4
5
6
3
-1,3
0,2
2,5
3
3
0,9
0,9
1,992
1,992
1,204
1,204
1,7247
1,7247
1,3528
1,3528
1,6074
1,6074
1,4271
1,4271
1,5521
1,5521
1,4641
f(x)
0
3
3
1,7647
1,7647
1,441
1,441
1,5019
1,5019
1,5
1,5
1,5
2
1,5
1
0,5
0
0,5
1,5
2,5
3,5
0,5
1,5
2,5
3,5
3
2,5
psi
3
3
0,9
1,7647
1,3287
1,441
1,542
1,5019
1,4987
1,5
1,5
1,5
-0,7
x amlior
1,5
0,4118 1,7647
1,441
-0,6588 0,3972 1,5019
0,5
-0,3471 0,2576
-0,7114 0,4157
1,5
-0,6996 0,4116
1,5
Troisime cas : divergence monotone (uniquement pour la substitution simple)

x=a+bx+cx
Substitution
iter
x
1
2
3
4
5
6
7
8
9
10
1
0,9
0,9
0,7755
0,7755
0,6219
0,6219
0,4345
0,4345
0,2091
0,2091
-0,057
-0,057
-0,366
-0,366
-0,714
-0,714
-1,096
-1,096
a
b
c
f(x)
-0,3
1,15
0,05
0,9
0,9
0,7755
0,7755
0,6219
0,6219
0,4345
0,4345
0,2091
0,2091
-0,057
-0,057
-0,366
-0,366
-0,714
-0,714
-1,096
-1,096
-1,5
0
-1,5
-1
-0,5
0,5
1,5
2,5
3,5
-1
-2
1,5
1,4
1,3
Wegstein
iter
x
1
2
3
4
5
6
1,2
f(x)
1
0,9
0,9
1,4082
1,4082
1,3691
1,3691
1,3723
1,3723
1,3723
1,3723
1,3723
0,9
0,9
0,7755
1,4082
1,4185
1,3691
1,3682
1,3723
1,3723
1,3723
1,3723
1,3723
psi
x amlior
1,1
1
1,245 5,0816 1,40816

0,9
1,2654 4,7678 1,36909
0,8
1,2889 4,4619 1,37226
0,7
0,6
1,2871 4,4835 1,37228

1,2872 4,4816 1,37228
0,5
0,75
- 3.21 -
0,9
1,05
1,2
1,35
1,5
Chapitre 3
Quatrime cas : divergence oscillatoire (uniquement pour la substitution simple)

x=a+bx+cx
Substitution
iter
x
1
1,5
2,425
2
2,425
1,2294
3 1,2294
2,7424
4 2,7424
0,7798
5 0,7798
3,2373
6 3,2373
0,0384
7 0,0384
3,9653
8 3,9653
-1,141
9 -1,141
4,8968
10 4,8968
Wegstein
iter
x
1
2
3
4
5
6
a
b
c
f(x)
2,425
2,425
1,2294
1,2294
2,7424
2,7424
0,7798
0,7798
3,2373
3,2373
0,0384
0,0384
3,9653
3,9653
-1,141
-1,141
4,8968
4,8968
-2,805
4
-0,9
-0,1
5
4
3
2
1
0
-2
-1
2,5
3,5
-1
-2
-3
4
3,5
3
1,5
2,425
2,425
1,9035
1,9035
1,9125
1,9125
1,9127
1,9127
1,9127
1,9127
1,9127
f(x)
psi
2,425
2,425
1,2294 -1,293
1,9035
1,9245 -1,333
1,9125
1,913 -1,282
1,9127
1,9127 -1,283
1,9127
1,9127 -1,283
1,9127
x amlior
0,5638 1,90349
0,5713 1,91251
2,5
2
1,5
1
0,5617 1,91271
0,5
0,5619 1,91271
0
0,5619 1,91271
0,5
1,5
3.4.2 Comparaison entre la substitution simple et la mthode de Rubin

x1=1.4 racine(1+x2)-0.4 x1
x2=1.5 racine(x1) -0.5 x2 +1.5
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
Substitution simple
It
x1
x2
f1
f2
1 3,000 2,000 1,225 3,098
2 1,225 3,098 2,344 1,611
3 2,344 1,611 1,325 2,991
4 1,325 2,991 2,267 1,731
5 2,267 1,731 1,407 2,893
6 1,407 2,893 2,200 1,833
7 2,200 1,833 1,476 2,808
8 1,476 2,808 2,142 1,918
9 2,142 1,918 1,535 2,736
10 1,535 2,736 2,092 1,990
11 2,092 1,990 1,584 2,674
12 1,584 2,674 2,050 2,051
13 2,050 2,051 1,625 2,622
14 1,625 2,622 2,014 2,101
15 2,014 2,101 1,660 2,578
16 1,660 2,578 1,984 2,143
17 1,984 2,143 1,688 2,541
18 1,688 2,541 1,959 2,178
19 1,959 2,178 1,712 2,510
20 1,712 2,510 1,938 2,208
21 1,938 2,208 1,732 2,485
22 1,732 2,485 1,921 2,232
23 1,921 2,232 1,749 2,463
24 1,749 2,463 1,906 2,252
25 1,906 2,252 1,762 2,445
- 3.22 -
1,762
1,893
1,774
1,883
1,783
1,874
1,791
1,867
1,798
1,861
1,804
1,856
1,808
1,852
1,812
1,849
1,815
1,846
1,818
1,843
1,820
1,841
1,822
1,840
1,823
1,838
1,825
1,837
1,826
1,836
2,445
2,269
2,430
2,283
2,417
2,295
2,406
2,305
2,397
2,313
2,390
2,319
2,384
2,325
2,379
2,330
2,375
2,334
2,371
2,337
2,368
2,340
2,366
2,342
2,364
2,344
2,362
2,345
2,360
2,347
1,893
1,774
1,883
1,783
1,874
1,791
1,867
1,798
1,861
1,804
1,856
1,808
1,852
1,812
1,849
1,815
1,846
1,818
1,843
1,820
1,841
1,822
1,840
1,823
1,838
1,825
1,837
1,826
1,836
1,827
2,269
2,430
2,283
2,417
2,295
2,406
2,305
2,397
2,313
2,390
2,319
2,384
2,325
2,379
2,330
2,375
2,334
2,371
2,337
2,368
2,340
2,366
2,342
2,364
2,344
2,362
2,345
2,360
2,347
2,359
Chapitre 3
Rubin
x1=1.4 racine(1+x2)-0.4 x1
x2=1.5 racine(x1) -0.5 x2 +1.5
Iter
x1
1
2
3
Iter
x2
4,0000
1,8293
1,6932
x1
4
f1
5,0000
2,0000
2,5288
x2
f2
1,8293
1,6932
1,9526
f1
2,0000
2,5288
2,1874
f2
2,3325
1,8071
1,8293
1,6932
2,3325
2,0000
2,5288
1,8329
1,6932
1,9526
2,3501
2,5288
2,1874
0,0222
-0,1139
-0,3325
0,1962
-0,1397
0,1197
0,1786
-0,1627
A-1
-5,8552
-3,3986
-9,9218
-0,6631
(A-1).B= J
-0,3697
0,3956
0,5684
-0,4992
E-J
1,3697
-0,3956
-0,5684
1,4992
(E-J)-1
0,8199
0,2163
0,3108
0,7490
Iter
x1
5
x2
1,8337
ligne modifier
Modif C
1,8071
1,8337
1,6932
1,8329
f1
2,3513
2,3501
x1-f1
x2-f2
-0,0258
-0,0176
1,8071
f2
1,8294
2,3556
x1-f1
x2-f2
0,0043
-0,0043
2
2,3325
2,3513
2,5288
modif D
- 3.23 -
1,8329
1,8294
1,9526
2,3501
2,3556
2,1874
Chapitre 3
3.4.3 Utilisation de Newton - Raphson n dimensions

f1=x1^2-x2-1
f2=x1-(x2-1)^2
Iter
Newton-Raphson
x1
x2
f1
f2
1 5,000
5,000 19,000 -11,000
J
10
1
-1
-8
J-1
0,1013
0,0127
-0,0127
-0,1266
2,937
3,367
4,257
-2,666
5,87342
1
-1
-4,7342
0,1766
0,0373
-0,0373
-0,2191
2,085
2,624
0,725
-0,552
4,17076
1
-1
-3,2481
0,2589
0,0797
-0,0797
-0,3324
1,854
2,383
0,054
-0,058
3,7075
1
-1
-2,7655
0,2989
0,1081
-0,1081
-0,4007
1,831
2,354
0,000
-0,001
3,66285
1
-1
-2,7072
0,3036
0,1122
-0,1122
-0,4108
1,831
2,353
0,000
0,000
3,66235
1
-1
-2,7064
0,3037
0,1122
-0,1122
-0,4110
1,831
2,353
0,000
0,000
3,66235
1
-1
-2,7064
0,3037
0,1122
-0,1122
-0,4110
1,831
2,353
0,000
0,000
3,66235
1
-1
-2,7064
0,3037
0,1122
-0,1122
-0,4110
1,831
2,353
0,000
0,000
3,66235
1
-1
-2,7064
0,3037
0,1122
-0,1122
-0,4110
10
1,831
2,353
0,000
0,000
3,66235
1
-1
-2,7064
0,3037
0,1122
-0,1122
-0,4110
11
1,831
2,353
0,000
0,000
3,66235
1
-1
-2,7064
0,3037
0,1122
-0,1122
-0,4110
- 3.24 -
Chapitre 3
25
20
15
10
F1
5
0
-5
4,8
4,2
3,6
2,4
1,8
1,2
0,6
00
Figure 3.11
0,6
1,2
x2
1,8
2,4
4,2
3,6
4,8
-10
x1
F1 en fonction de x1 et x2
6,00
4,00
2,00
4,8
3,2
2,4
1,6
0,8
Figure 3.12
00
0,8
1,6
2,4
3,2
x2
-8,00
-10,00
-12,00
-14,00
-16,00
4
4,8
0,00
-2,00
-4,00
-6,00
x1
F2 en fonction de x1 et x2
- 3.25 -
F2
Chapitre 4
CHAPITRE 4
OPTIMALISATION EN INGENIERIE DES PROCEDES
4.1 Introduction
Un modle de simulation statique est un ensemble dquations reliant les valeurs de toutes les
variables dtat dun procd. Le nombre dquations Ne est infrieur au nombre de variables N v . La
diffrence N f = N v N e est le nombre de degrs de libert du modle ; N f quations supplmentaires
(les spcifications) doivent tre imposes avant de tenter la rsolution du modle.
Certaines spcifications apparaissent naturellement partir de contraintes externes (par exemple, la
composition des matires premires, le minimum de puret requise par le march, la limite dmission de
polluants, etc.) Dautres spcifications peuvent apparatre arbitrairement et sont fixes soit de manire
empirique, soit par exprience.
Les rsultats de la simulation vont dpendre du choix de ces spcifications arbitraires.
Loptimalisation est un moyen de slectionner celles-ci dans le but datteindre la meilleure valeur possible
de la fonction objectif (merit function) base sur ces rsultats.
4.1.1 Glossaire
4.1.1.1 Fonction objectif
La fonction objectif doit prendre en compte de manire quantitative toutes les caractristiques du
procd qui ont un impact sur la qualit de ce dernier. Des fonctions objectifs typiques peuvent impliquer le
niveau de production, les charges en matires premires et en nergie ou les autres de cots de production
comme la main doeuvre, les catalyseurs, la dprciation de lquipement, etc.) Gnralement, les diffrents
facteurs intervenant dans la fonction objectif sont convertis en valeur montaire dans le but de comparer
plusieurs objectifs.
4.1.1.2 Variables
Les variables, dans un problme doptimalisation, correspondent tous les quantits physiques
indpendantes apparaissant dans les quations des modles. Elles peuvent tre soit continues (temprature,
pression, dbits) ou discrtes (nombre de plateaux dans une colonne distiller, vanne ouverte ou ferme,
quipement en marche ou larrt).
- 4.1 -
Chapitre 4
4.1.1.3 Contraintes
Les contraintes sont des relations qui limitent le domaine de variations des variables. Elles peuvent
tre linaires ou non linaires.
Les contraintes dgalit sont des quations reliant certaines variables entre elles : les quations du
modles sont des contraintes dgalit ; la composition des matires premires ou des variables externes
non contrlables comme la pression atmosphrique sont dautres exemples de contraintes dgalit.
Les contraintes dingalit dfinissent des limites aux valeurs de certaines variables ou de certaines
fonctions de variables. Elles apparaissent cause des caractristiques des quipements (pression ou
temprature maximum de fonctionnement dun appareil), des normes lgales (missions) ou des
spcifications commerciales (puret dun produit).
4.1.2 Types de problmes
La difficult de rsoudre les problmes doptimalisation est relie certaines caractristiques de
leur formulation. Le nec plus ultra est de progresser rgulirement avec lapparition de nouvelles mthodes
numriques et la disponibilit de calculateurs plus puissants.
4.1.2.1 Linaires - Non linaires
Quand on considre des procds rels, le comportement purement linaire est rare. Ceci ncarte
pas lutilisation dapproximations linaires mais la plus grande prudence est requise pour linterprtation des
rsultats. Les problmes linaires (fonction objectif linaire, contraintes dgalit et dingalit linaires,
limites infrieures et suprieures sur les variables) peuvent tre rsolues facilement avec les programmes
existants. Les problmes de plus grande chelle avec plus de 105 variables peuvent ltre avec les meilleurs
algorithmes. Ces algorithmes sont bass sur des oprations sur les ranges et les colonnes dune matrice de
coefficients dquations et sont similaires aux algorithmes utiliss pour rsoudre les systmes dquations
linaires. La solution est atteinte en un nombre fini doprations quand le problme est bien conditionn
(absence de singularit).
Le cas des problmes non linaires est plus compliqu. On observe frquemment la prsence de
plusieurs solutions optimales (extrema locaux). La manire dont le problme est formul (mise lchelle
des variables et des quations : scaling) et la qualit du point initial affectent lefficacit de la solution.
4.1.2.2 Fonctions convexe et concave

Une fonction est concave sur une rgion R si :
x1 , x 2 R, [ 0,1] : f x1 + (1 ) x 2 f ( x 1 ) + (1 ) f ( x 2 )
- 4.2 -
Chapitre 4
f(x)
f(x)
x2
x1
x1
x2
Fonction convexe
Fonction concave
Malheureusement, cette dfinition nest pas pratique pour tester la concavit dun systme.
Lanalyse de la drive seconde de f(x) est plus approprie pour les applications quotidiennes. La fonction
plusieurs variables f(x) est convexe quand la matrice hessienne est dfinie positive.
2 f ( x)
H ( x) = f ( x) =
x i x j
2
La matrice hessienne est dfinie positive quand pour tout x 0 : x T H x > 0 . Cette condition est
rencontre quand toutes les valeurs propres de la matrice symtrique H sont strictement positives. Quand
une fonction est convexe, le minimum non contraint est un minimum global.
4.1.2.3 Dimension du problme
Plus il y a de variables et plus le problme est difficile. Des mthodes de recherche itrative sont
disponibles pour manipuler des problmes englobant quelques centaines quelques milliers de variables.
En gnie chimique, les problmes comportent gnralement un grand nombre de contraintes et le
nombre de degrs de libert est gnralement trs infrieur au nombre de variables. Cest pourquoi les
mthodes trajectoire ralisable (feasible path methods) sont frquemment recommandes : elles
consistent identifier un ensemble minimal de variables indpendantes qui doivent tre optimises en itrant
sur une boucle externe et rsoudre les quations et contraintes dgalit du modle pour identifier les
variables dpendantes en itrant sur une boucle interne. Le principal avantage de ces mthodes est que la
srie des solutions intermdiaires correspond des problmes de simulations converges fournissant une
information intressante mme en cas dchec de la procdure globale doptimalisation : il vaut mieux une
solution sous-optimale que pas de solution du tout !
Les mthodes trajectoire irralisable (unfeasible path methods) essayent simultanment de
rechercher une meilleure valeur de la fonction objectif pendant quelles rsolvent les contraintes dgalit et
- 4.3 -
Chapitre 4
ramnent les contraintes violes lintrieur de leurs limites. Elles sont gnralement plus efficaces en temps
de calcul quand elles parviennent une solution !
4.1.2.4 Continu ou discret
Les fonctions objectifs ou les contraintes peuvent comporter des discontinuits finies. Par exemple,
le cot dun quipement dchangeur de chaleur peut ne pas voluer monotonment avec la pression, la
temprature ou la charge (quand celle-ci augmente, il peut tre ncessaire de diviser un appareil en units
parallles ou dadopter une technologie diffrente).
4.1.2.5 Objectif non contraint ; contraintes dgalit ; contraintes dingalit
Les problmes les plus simples nimpliquent pas de contraintes ou juste des limites infrieures et
suprieures sur les variables. Les contraintes dgalit et dingalit impliquent des fonctions non linaires
de variables, ce qui accrot grandement la complexit du problme. Une solution trajectoire ralisable
peut tre difficile identifier (dans certains cas, elle peut mme ne pas exister...).
4.1.3 Type de mthodes doptimalisation
De nombreux algorithmes ont t proposs pour rsoudre les problmes doptimalisation. Ils se
regroupent en trois classes.
4.1.3.1 Mthodes directes
Celles-ci sont bases sur une recherche alatoire autour du point courant dans le but destimer un
domaine ou une direction dans lesquels un progrs peut tre espr grce aux essais prcdents. Les
mthodes directes sont faciles utiliser et sont robustes, mais gnralement elles requirent plus de temps
de calcul. Elles doivent tre limites des petits problmes.
4.1.3.2 Mthodes indirectes bases sur des critres dexistence dun extremum
Elles sont bases sur des critres mathmatiques dexistence dun optimum. Elles tentent
didentifier les groupes de variables qui collent de mieux en mieux ces critres. Les mthodes indirectes
sont gnralement plus efficaces que les mthodes directes mais elles sont plus complexes et impliquent
plus de mathmatiques (par exemple, estimation ou calcul de drives).
4.1.3.3 Heuristiques
Les mthodes heuristiques sont bases sur les connaissances disponibles sur certaines classes de
problmes doptimalisation et essaient de tirer parti des motifs (patterns) dans les solutions optimales. Par
exemple, les algorithmes gntiques reprsentent les groupes de dcisions comme des squences de codes
dune manire analogue celle des gnes codant les tres vivants. Les algorithmes gntiques appliquent
- 4.4 -
Chapitre 4
des rgles de transformations pour simuler lvolution de population de gnes. La reproduction donne
plus de chance de reproduire les schmas fructueux tandis que les mutations autorisent la cration alatoire
de nouveaux motifs. Les algorithmes gntiques sont spcialement appropris pour travailler avec des
problmes impliquant des structures (patterns) ou des fonctions qui ne sont pas diffrentiables.
4.2 Optimisation pour un systme une seule variable continue

Ce problme mrite dtre tudi pour plusieurs raisons :
dans maints problmes, llimination des contraintes par insertion dans une fonction objectif rduit la
dimension du problme une variable ;
les techniques pour les optimisations dans un espace multidimensionnel gnre gnralement une
direction de recherche pour laquelle une optimisation unidimensionnelle est ralise.
Loptimisation pour un systme une seule variable implique deux tapes :
parcourir le domaine de la variable pour encadrer loptimum ;
raffiner la recherche jusqu ce que loptimum soit localis avec la prcision voulue.
Le parcours de lespace des variables est normalement fait en identifiant la direction de
dcroissance de la fonction et, ensuite, en progressant dans cette direction jusqu ce que la fonction
croisse nouveau. Les dernires valeurs recherches dfinissent ainsi une rgion o loptimum est
susceptible dtre localis (au moins si la fonction objectif est continue). Dans le but de rduire le nombre
dvaluations de la fonction objectif, la taille du pas peut tre augmenter chaque itration.
Un exemple est donn la figure suivante. Loptimum de la fonction f = 2000 + ( x 100) est
2
encadr par cet algorithme : on part de x = 1 et on se dplace vers le minimum en doublant la taille du pas
chaque essai.
- 4.5 -
Chapitre 4
40000
30000
20000
10000
0
0
50
100
150
200
250
300
4.2.1 Rduction de lintervalle

Une fois que le minimum de la fonction objectif a t encadr entre xmin et xmax, nous pouvons
calculer deux nouvelles valeurs f1 = f ( x 1 ) et f 2 = f ( x 2 ) , de telle manire que :
x min < x 1 < x 2 < x max
En comparant f1 et f2 et en supposant que la fonction est convexe, nous pouvons rejeter une partie
de lintervalle initial :
f(x)
f(x)
f1>f2
Xmin
x1
f1<f2
x2
Xmax
Xmin
x1
x2
Xmax
Si f1 > f2, on rejette [xmin,x1] et on pose xmin = x1 pour litration suivante.

Si f1 < f2, on rejette [x2,xmax] et on pose xmax = x2 pour litration suivante.
Si f1 = f2, on rejette [xmin,x1] et [x2,xmax], on pose xmin = x1 et xmax = x2 pour litration suivante.
Si xmin, x1, x2 et xmax sont rgulirement espacs, la longueur de lintervalle sera rduite dun tiers
chaque itration au prix de deux valuations de fonction. Seule lune delle est encore utilisable litration
- 4.6 -
Chapitre 4
suivante. Si L0 est la taille originale de lintervalle, elle devient Lk = (2/3)k L0 aprs k itrations (et 2k
valuations de la fonction).
4.2.2 La rgle du nombre dor
Une mthode de recherche plus efficace est obtenue en choisissant x1 et x2 dans la proportion
connue sous le nom de section du nombre dor dfinie par :
xmax - x2 = x1 - xmin
x1 x min
x x1
= max
=
x max x 1 x max x min
5 1
0,618
2
Avec cette mthode, le point lintrieur du nouvel intervalle peut tre conserv litration
suivante et un seul point doit tre rvalu chaque itration. La longueur de lintervalle devient
k
L k = ( 0,618) . L 0 aprs k itrations (k valuations de fonction).
4.2.3 Lapproximation polynomiale
Une autre de classe de mthodes de recherche unidimensionnelle procde par approximation de la
fonction objectif par une simple fonction dont loptimum est aisment identifi. Un choix possible est
linterpolation quadratique. Etant donn trois points x1, x2, x3 entourant le minimum de la fonction objectif,
nous pouvons calculer les valeurs correspondantes de la fonction objectif et ajuster les coefficients pour
une fonction quadratique interpolant les trois valeurs donnes :
f ( x) q ( x) = a + b x + c x 2
Le minimum de q(x) est directement obtenu par :
x =
b
2c
Un peu dalgbre nous amne :

2
2
2
2
2
2
1 ( x 2 x 3 ) f 1 + ( x 3 x1 ) f 2 + ( x 1 x 2 ) f 3
x =
2 ( x 2 x 3 ) f1 + ( x 3 x1 ) f 2 + ( x1 x 2 ) f 3
La fonction objectif est value en ce point et x* remplace x1 ou x3 en sassurant que les points
restants entoure toujours loptimum (la fonction objectif pour le point milieu doit tre le plus bas ce
moment).
Remarquons que la procdure ne doit pas tre applique si le minimum na pas t entour : la
fonction quadratique interpolant les trois points arbitraires peut dgnrer en une droite ou tre concave
(dans ce cas, x* localise le maximum ...).
- 4.7 -
Chapitre 4
0
0
Approximation quadratique de f(x) = 1.2 + cos(x)
4.3 Optimisation
contraintes
dun
set
de
variables
continues
non
La plupart des problmes pratiques requirent loptimisation dun modle plusieurs variables. La
plupart des mthodes examines pour les problmes unidimensionnels peuvent tre gnralises aux
problmes plus complexes.
4.3.1 Critre dexistence dun optimum
Le problme est la minimisation de f(x) pour x E n .
Pour un problme de programmation non linaire et non contraint, les conditions ncessaires pour
que x* soit un minimum local sont :
f(x) est diffrentiable en x*
f ( x *) = 0 (cest--dire, un point stationnaire existe en x*)

En plus des conditions dfinies ci-dessus, la condition suffisante pour que x* soit un minimum local
est :
2 f ( x *) > 0 (cest--dire que la matrice hessienne est dfinie positive).
Les conditions pour lexistence dun maximum sont les mmes, except le fait que la matrice
hessienne doit tre dfinie ngative.
Des points optima qui ne satisfont pas les conditions ncessaires et suffisante peuvent exister : par
exemple, pour une discontinuit ou lorsque la fonction est continue, mais que ses drives ne le sont pas.
- 4.8 -
Chapitre 4
4.3.2 Mthodes directes

Les mthodes directes ne requirent pas lutilisation de drives pour la dtermination de la
direction de recherche. Elles ont lavantage de la simplicit, mais sont gnralement moins efficaces que les
mthodes que nous verrons plus loin.
4.3.2.1 Recherche pour une seule variable
Une simple mthode doptimisation pour rsoudre un problme n dimensions est de choisir n
directions indpendantes de recherche. La fonction objectif est minimise le long de chaque direction de
recherche de manire squentielle en utilisant un algorithme de recherche unidirectionnel. Cette mthode est
inefficace moins que les directions soient bien choisies (par exemple, parallles aux crtes de la fonction
objectif).
Le but de bien des recherches qui ont t effectues tait de gnrer efficacement une direction de
recherche. Lexprience a montr que les directions conjugues sont plus efficaces que les directions
arbitraires. Un ensemble de directions linairement indpendantes s1 ,K, s n sont conjugues par rapport
une matrice carre dfinie positive Q si :
t
(s i ) Q ( s i ) = 0
pour 0 i < j n
Lorsque f(x) est une fonction quadratique, il peut tre montr que loptimum peut tre localis
aprs exactement n tapes de recherche unidirectionnelle condition que les directions de recherche soient
conjugues la matrice hessienne constante et que la minimisation soit ralise avec une prcision extrme
chaque tape. Des algorithmes ont t conus pour gnrer et mettre jour automatiquement un
ensemble de directions conjugues de recherche (Powell 1965, Brent 1973).
4.3.2.2 Simplex (Nelder, Mead)
Ces mthodes calculent la valeur de la fonction pour des points situs aux sommets dune figure
gomtrique rgulire balayant lespace des variables (un simplex). Les exemples sont le triangle pour 2
dimensions, le ttradre pour 3 dimensions... Le sommet correspondant la valeur de la fonction objectif la
plus mauvaise est supprim et remplac par sa rflexion travers le barycentre des autres sommets. Ainsi,
le simplex se dplace dans lespace des variables jusqu ce quil contienne loptimum.
Nelder et Mead ont modifi la mthode pour permettre au simplex de saccrotre ou de dcrotre
tout en se dplaant dans le but de modifier la taille du pas de recherche [voir lalgorithme et son
implmentation dans Himmelblau (1972)].
- 4.9 -
Chapitre 4
4.3.3 Mthodes indirectes du premier ordre

Les mthodes indirectes ncessitent le calcul du gradient de la fonction objectif dans le but de
choisir une direction de recherche. Une bonne direction s doit conduire dun point courant x0 un point x1
o la fonction objectif se trouve rduite :
f ( x1 ) < f ( x 0 )
Une telle direction est une direction de descente caractrise par : T f ( x) . s < 0
La plus simple est la mthode de la plus grande pente o la direction litration k est prise comme
loppos du gradient local :
s k = f ( x k )
x k +1 = x k + k s k = x k k f ( x k )
k est un scalaire gal la longueur du pas. Il peut tre soit fixer une valeur constante ou soit tre
dtermin en utilisant une recherche unidirectionnelle. La mthode de la plus grande pente est ainsi base
sur une approximation linaire de la fonction objectif.
x*
Xk
Approximation : mthode de la plus grande pente
La difficult fondamentale avec la mthode de la plus grande pente est sa grande sensibilit la
mise lchelle de f(x). La convergence peut tre trs lente et osciller autour de ce qui apparatrait comme
une direction efficace de recherche.
La mthode des gradients conjugus (Fletcher et Reeves, 1964) reprsente une nette amlioration
par rapport la mthode prcdente avec un accroissement faible des calculs ncessits. La mthode des
gradients conjugus combinent les informations sur le gradient local avec les gradients obtenus aux
itrations prcdentes. Une recherche unidirectionnelle doit tre accomplie avec prcision le long de
chaque direction de recherche.
La premire itration est similaire celle de la direction de recherche de la plus grande pente. Les n
itrations suivantes sont bases sur le schma :
x k + 2 = x k +1 + k +1s k+ 1
- 4.10 -
Chapitre 4
s k +1 = f ( x k +1 ) + s k
T f ( x k +1 ). f ( x k +1 )
T f ( x k ). f ( x k )
Aprs n itrations, la procdure est recommence. On peut montrer que les directions successives
de recherche sont conjugues quand la fonction objectif est quadratique.
4.3.4 Mthodes indirectes du second ordre
Ces mthodes tirent parti de linformation disponible grce aux drives secondes (la matrice
hessienne ou une approximation) pour amliorer la recherche en calculant la fois une direction et une
longueur de pas.
4.3.4.1 Mthode de Newton
La fonction objectif est approche au point courant xk par une srie de Taylor limite au second
ordre.
f ( x) q( x) = f ( x k ) + T f ( x k ) x k +
1
( x k ) T H ( x k ) ( x k )
2
Le minimum de cette approximation obtenu en galant zro le gradient de cette approximation

quadratique et en rsolvant le systme dquations qui en rsulte :
0 = q( x) = f ( x k ) + H( x k ) ( x k )
Ce dernier peut scrire suivant lalgorithme de Newton :
x k +1 x k = x k = [ H( x k )] 1 f ( x k )
La mthode de Newton prdit la fois une direction de recherche et une longueur pour le pas.
Quand la fonction objectif est quadratique, le minimum est localis en un seul pas. Toutefois, pour une
fonction gnrale, la mthode de Newton ne fournit quun meilleur point et une recherche par itrations est
ncessaire; la mthode peut tre raffine en appliquant une recherche le long de la direction obtenue par la
formule de Newton et en valuant numriquement la longueur du pas k dans :
x k +1 x k = x k = [ H( x k )] 1 f ( x k )
x*
Xk
Approximation quadratique : mthode de Newton
- 4.11 -
Chapitre 4
La mthode de Newton fait preuve dexcellentes proprits de convergence quand la fonction

objectif est presque quadratique. Cependant, elle requiert le calcul de la matrice hessienne qui est lourd
pour les grands systmes. La convergence peut tre amliore soit en reculant quand le pas de Newton ne
conduit pas une dcroissance de f(x), soit en contraignant le pas pour quil reste dans une rgion de
confiance pour laquelle lapproximation quadratique reste satisfaisante.
4.3.4.2 Algorithme de Levenberg-Marquardt
De plus, une prcaution doit tre prise lorsque la fonction objectif nest pas convexe (ou lorsque la
matrice hessienne nest pas dfinie positive) : le pas prdit peut ne pas correspondre la direction de
descente ou la mthode peut aboutir un point de selle. Plusieurs auteurs ont propos de remplacer la
matrice hessienne dans lalgorithme de Newton par une matrice H( x) conditionne pour rester dfinie
positive en additionnant une grande constante positive b tous les lments diagonaux de la matrice H(x) :
H = H( x ) + I
On peut montrer que pour des valeurs leves de , le pas calcul correspond un petit pas dans la
direction de la plus grande pente alors que, pour une valeur trs faible de , il y correspond le pas de
Newton. Des rgles heuristiques ont t proposes pour ajuster : il est rduit quand une itration conduit
une dcroissance de la fonction objectif et il est augment dans le cas contraire.
trust region
Xk
x*
quadratic
approximation
(forced to be
positive definite)
Illustration de lapproximation quadratique force rester dfinie positive et de lutilisation du

domaine de confiance
4.3.4.3 Approximation de la matrice hessienne : mthodes de mise jour

Nous pouvons interprter la mthode de Newton comme un changement de coordonnes qui
transforme le vecteur gradient en un pas orient vers le minimum de lapproximation quadratique locale. La
multiplication du gradient par linverse de la matrice hessienne est un exemple dune transformation dun
- 4.12 -
Chapitre 4
espace vectoriel en un autre. Nous pouvons ds lors dfinir une famille de mthodes doptimalisation
(quasi-Newton ou une mtrique variable) o le pas litration k est donn par :
x k +1 x k = k [ M k ] 1 f ( x k )
La matrice directionnelle [ M k ]
reprsente une approximation convenable de linverse de la
matrice hessienne. Mme la mthode de la plus grande pente obit ce schma en posant M = I.
La matrice M doit idalement prsenter les proprits suivantes :
mise jour chaque itration pour rendre compte des informations disponibles propos de f(x) et de
son gradient ;
ne requiert pas un calcul explicite des drives secondes de f(x) ;
rester dfinie positive ;
tre symtrique ;
converger vers la matrice hessienne relle loptimum ;
fournir de bonnes proprits de convergence lalgorithme doptimisation.
Une manire destimer la matrice hessienne est de remarquer que, pour les fonctions quadratiques,
f ( x k ) f ( x 0 ) = H( x k x 0 )
Si nous posons :
x k = x k+1 x k
et
g k = f ( x k +1 ) f ( x k )
Une proprit dsirable pour la matrice M serait de satisfaire lquation suivante chaque
itration :
M k+1 x k = g k
Cela correspond un systme n quations n2 inconnues (nombre dlments de Mk+1). Une
solution gnrale est :
M k+1 M k =
g k y T M k x k z T
y T x k
z T x k
o y et z sont des matrices telles que xk ou gk.

Mme en prenant en compte le fait que Mk+1 doive tre symtrique, une infinit de techniques
peuvent tre proposes pour calculer Mk+1 aprs litration k. Quelques algorithmes ont prouv quils
taient suprieurs en maintenant M dfinie positive, symtrique et en obtenant de bonnes proprits de
convergence. Mais des expriences numriques sont ncessaires pour identifier les meilleures techniques de
mise jour. Celle propose indpendamment par Broyden, Fletcher, Goldbarb et Shanno en 1970 est
probablement celle la plus couronne de succs :
- 4.13 -
Chapitre 4
g k ( g k ) T M k x k ( x k ) T H k
( g k ) T x k
( x k ) T H k x k
M k+1 M k =
La formule de mise jour peut tre dveloppe pour obtenir directement linverse de la matrice
M:
(M
k +1 1
k 1
(M )
[ x k - (M k ) -1 g k ] ( x k ) T + x k [ x k - (M k ) -1 g k ]T
=
( g k ) T x k
[x
- (M k ) -1 g k ] g k x k ( x k )
T
[(g ) x ][(g ) x ]
k T
k T
4.4 Optimisation dun ensemble de variables continues et

contraintes
Maints problmes industriels ne peuvent tre rsolus en appliquant aveuglement un outil
mathmatique pour optimiser le modle. En pratique, beaucoup de contraintes limitent la rgion de
recherche. Elles peuvent correspondre aux limites de validit du modle ou reprsenter les limitations
physiques de la technologie. Par consquent, une analyse numrique doit affronter ce problme et
dvelopper des outils plus labors.
4.4.1 Programmation linaire
Mme si la plupart des problmes chimiques dingnierie sont de nature non linaires, il peut tre
utile dessayer de les linariser pour pouvoir utiliser des codes doptimisation robustes et puissants
disponibles pour la rsolution de problmes de programmation linaire gnrale du type :
r
f ( x) = c x =
t
minimiser
c x
i
i =1
avec
x i 0 pour
i = 1,K , r
a x
et
ij
j = 1,K , m
= bj
i =1
a x
ij
j = m + 1,K , p
bj
i =1
Les contraintes dingalit sont converties en contraintes dgalit en y additionnant des variables
de relaxation (slack) s j 0 de manire ce que :
r
a x s
ij
= bj
i =1
- 4.14 -
Chapitre 4
G. Dantzig proposa en 1937 la mthode Simplex, un algorithme gnral pour manipuler les
problmes de programmation linaires de ce type. La procdure itrative est base sur le fait que pour tout
problme linaire impliquant r variables positives lies par un ensemble convexe de p contraintes,
loptimum est ncessairement situ lintersection des r contraintes ou limites variables. La procdure de
rsolution dmarre une telle intersection ralisable (une base) et avance chaque itration dun vertex
un autre adjacent ; elle applique des critres pour quitter les contraintes et entrer dans la base, pendant
quelle amliore la fonction objectif et reste une intersection ralisable des contraintes.
Limportance de la programmation linaire ans la rsolution des problmes pratiques est due au fait
que les logiciels disponibles sont capables de traiter les problmes de trs grandes dimensions. Mme si les
modles linaires sont seulement des approximations du monde rel et doivent tre correctement borns
pour viter les extrapolations en dehors de leurs limites de validit, ils sont encore largement utiliss pour
des problmes comme le contrle des oprations sur sites larges, les attributions de ressources et le
mlange des produits (product blending).
Les lecteurs intresss peuvent se rfrer aux ouvrages sur la programmation linaire tels que Taha
(1971). Une courte introduction peut galement tre trouve dans Edgar et Himmelblau (1988).
4.4.2 Programmation non linaire et contrainte
Un problme plus gnral, mais plus difficile peut tre dcrit de cette manire :
minimiser
f(x)
x = x 1 , x 2 , K, x n
soumis
h j ( x) = 0
j = 1,K , m
g j ( x) 0
j = m + 1,K , p
Les contraintes dingalits peuvent toujours tre transformes en contraintes dgalit par
lutilisation de variables de relaxation (slack) :
g j ( x) s 2j = 0
j = m + 1,K , p
Une procdure attrayante serait lutilisation de lensemble des contraintes dgalit pour liminer p
variables. Aprs cette substitution, on terminerait ainsi avec un problme non-contraint pour lequel des
mthodes existent.
Malheureusement, cette mthode, comme telle, ne peut tre applique en pratique : le choix des
variables qui doivent tre limines (cest--dire, la slection dune base) nest pas directe. Lensemble des
quations de contrainte peut avoir plusieurs solutions ou mme pire : pas de solution du tout pour certaines
valeurs des variables indpendantes restantes. Cest pourquoi, dautres algorithmes de rsolution ont t
proposs.
- 4.15 -
Chapitre 4
4.4.3 Formulation de Lagrange

La formulation de Lagrange permet de transformer un problme doptimisation avec contraintes en
un non-contraint de plus grande dimension. La fonction objectif non-contrainte en rsultant est appele
lagrangien et est une combinaison linaire de la fonction objectif originale et de toutes les quations de
contraintes pondres par des variables additionnelles appeles multiplicateurs de Lagrange.
Pour le problme doptimalisation donn la section prcdente, le lagrangien correspondant est :
p
h (x) + [g (x) s
m
L( x, ) = f ( x) +
j =1
2
j
j = m +1
Cette fonction doit tre minimise relativement x, et s (avec un total de n+p+(p-m) variables
indpendantes). La valeur optimum de la fonction lagrangienne est identique la valeur optimum celle de
la fonction objectif contrainte originale puisque toutes les quations de contraintes doivent tre
identiquement nulles la solution.
Les conditions ncessaire et suffisante pour obtenir un minimum en x* sont :
que f(x*) soit convexe ;
que les contraintes forment un ensemble convexe au voisinage de x* ;
et que les conditions de stationnarit soient satisfaites :
L( x )
=0
x i
L( x )
=0
j
L( x )
= 2 j sj = 0
s j
Les multiplicateurs de Lagrange expriment la pnalit apporte la valeur de la fonction objectif

suite lexistence dune contrainte. Ils sont utiles pour une analyse de sensibilit puisquils indiquent
lampleur avec laquelle la fonction objectif peut tre amliore en relaxant chaque contrainte dune unit.
Ainsi, les multiplicateurs de Lagrange des contraintes dingalit inactives sont nuls.
4.4.4 Critre dexistence dun optimum : conditions de Kuhn-Tucker
Les conditions ncessaire et suffisante dexistence dun optimum dans les problmes gnraux de
programmation non linaire sont forts complexes et leurs dmonstrations peuvent tre lues dans les
ouvrages traitant de loptimalisation (par exemple, Gill et al, 1981). Ils sont rsums dans Edgar et
Himmelblau (1972). Leur description sort du cadre de ce cours.
- 4.16 -
Chapitre 4
4.4.5 Quelques familles dalgorithmes doptimalisation avec contraintes

Etant donn la complexit et la grande varit de problmes non linaires, trouver le meilleur
algorithme doptimisation est hautement improbable. Des expriences sont ncessaires pour valuer la
valeur des codes disponibles. Plusieurs approches ont t adoptes dans les codes les plus efficaces
actuellement et sont prsentes brivement dans les sections suivantes [voir Lasdon (1983) pour une tude
plus approfondie].
4.4.5.1 Mthode du gradient rduit gnralis
Aprs lintroduction des variables de relaxation (slack variables) pour transformer toutes les
contraintes en galits h(x), les algorithmes GRG (generalized reduced gradient) utilise les p contraintes
dgalit pour rsoudre p des variables xb (appartenant lensemble de base) en fonction des n-p variables
indpendantes xnb (n'appartenant pas la base).
Les contraintes peuvent tre crites comme suit et leur jacobien peut tre partitionn :
h( x b , x nb ) = 0
h h h
=
,
= ( Bb , B nb )
x x b x nb
Lalgorithme doit partir dun point ralisable (feasible) et lensemble de base devrait tre tel que
Bb soit non singulier dans le but de permettre le calcul des variables de base en fonction de lensemble
nappartenant pas la base. La fonction objectif peut ainsi tre exprime en fonction de xnb :
F( x nb ) = f x b ( x nb ), x nb
et le problme original est transform en un problme rduit, plus simple, non contraint :
minimiser F( x nb ) soumis aux bornes sur xnb
La fonction transforme est appele fonction objectif rduite (reduced objective) et son gradient,
le gradient rduit (reduced gradient). Il est calcul chaque itration comme suit :
F
f
T
=
z Bnb
x nb x nb
o le vecteur z est obtenu en rsolvant :

Bb z =
T
f
x b
Les algorithmes GRG procdent en rsolvant une squence de problmes rduits. Les variables
n'appartenant pas la base et qui ne sont pas leur limites infrieure ou suprieure sont appeles
superbasiques. A chaque itration, le gradient rduit avec les variables n'appartenant pas la base fixes
- 4.17 -
Chapitre 4
une de leurs limites permet de vrifier si certaines de ces variables doivent tre relaxes et rejoindre
lensemble des variables superbasiques. Les directions de recherche sont gnralement obtenues en
appliquant la mthode du gradient conjugu ou une mthode mtrique variable, bases sur le gradient des
variables superbasiques Une minimisation unidimensionnelle est accomplie dans la direction de recherche,
et chaque pas, les contraintes de liaisons sont utilises pour recalculer les variables de base. Des
itrations par la mthode de Newton sont employes dans la plupart des cas pour rsoudre lensemble des
quations dans le cas o il est non linaire.
Lexprience montre que la stratgie peut tre implmente dans des algorithmes qui sont la fois
fiables et efficaces. Puisque la mthode essaye toujours de vrifier les contraintes, tous les rsultats
intermdiaires sont ralisables et peuvent fournir des informations utiles; cela mme si loptimisation choue,
puisque lon attend un progrs par rapport un point de dpart de la valeur de la fonction objectif chaque
itration importante.
4.4.5.2 Programmation linaire squentielle
Ces mthodes sont bases sur la linarisation de la fonction objectif non linaire et des fonctions de
contrainte autour du point courant de la base :
f ( x) f ( x 0 ) + f ( x 0 ) x
T
h j ( x ) h j ( x 0 ) + h j ( x 0 ) x
T
Les limites infrieure et suprieure sont imposes sur x pour tre certain que les erreurs
dapproximation sont limites. Le problme de programmation linaire rsultant est alors rsolu pour des
variables x et le nouveau point x0+x est test. Si une amlioration est observe, il est accept et devient
le point de base pour la prochaine itration. Sinon, les limites sur x sont rendues plus troites et le
problme de programmation linaire est rsolu nouveau. Les points intermdiaires ne doivent pas tre
ralisables : les mthodes SLP (successive linear programming) peuvent prendre un chemin irralisable
vers loptimum mme si le point de base initial de base est ralisable. Lavantage des mthodes SLP est la
disponibilit de codes efficaces pour les problmes de programmation linaire de grande chelle dans les
bibliothques informatiques.
4.4.5.3 Programmation quadratique squentielle
Lide sous-jacente est semblable celle des mthodes SLP, mais les mthodes SQP (successive
quadratic programming) rsolvent une srie de programmes quadratiques (cest--dire, une fonction
objectif quadratique avec des contraintes linaires). La fonction objectif quadratique est :
[ f ( x )]
0
x + ( x) Q x
T
o Q est lapproximation de la matrice hessienne de la fonction de Lagrange :
- 4.18 -
Chapitre 4
L = f T g T h
o et m sont les multiplicateurs de Lagrange
Une manire de rsoudre le problme est dappliquer la mthode de Newton lensemble des
quations exprimant le critre dexistence dun optimum (cf. les conditions de Kuhn-Tucker).
La solution de programmation quadratique nest gnralement pas adopte pour former un
nouveau point, mais bien une direction de recherche pour laquelle est accompli une minimisation de la
fonction de Lagrange dans la direction x.
Lapproximation de la matrice hessienne Q est gnralement modifie chaque itration par
lintermdiaire dune formule de mise jour qui prserve sa structure (matrice symtrique, dfinie positive).
Une autre technique pour maintenir une matrice dfinie positive Q est dutiliser un lagrangien augment
comme fonction objectif o les termes de pnalit impliquant les carrs des carts aux contraintes sont
ajouts lobjectif.
Quand la matrice Q est dfinie positive, le problme quadratique avec les contraintes linaires est
convexe et loptimum global peut tre dmontr. Toutefois, des problmes peuvent apparatre quand les
contraintes (linarises) sont telles quaucune solution ralisable nexiste pour un sous-problme
quadratique. Dautres limitations surgissent de la difficult de manipuler les matrices Q de grande taille. Le
dveloppement de codes efficaces de programmation quadratique de grande chelle et de formules
efficaces de mise jour qui conservent les proprits creuses de lestimation de la matrice hessienne font
aussi partie du dveloppement actuel.
4.4.6 Programmation en variables entires
Certains problmes doptimisation sont naturellement formuls avec certaines variables prenant
seulement des valeurs discrtes. Par exemple, les variables continues peuvent reprsenter des dbits,
tempratures, pressions, etc., alors que les variables entires des types dquipement ou de dcisions,
comme lexistence ou la disponibilit de certains quipements.
La solution de problmes de programmation non-linaire en nombres entiers et non entiers est un
domaine ardu et aucune technique gnrale doptimisation na montr tre couronne de succs dans le
traitement de tous types de problmes.
Une technique prometteuse repose sur des mthodes duales qui sont appliques de manire
cyclique. Le problme est dabord de simplifier (linarisation) pour permettre lutilisation de codes de
programmation linaire mixte (nombres entiers et non-entiers) dans le but dobtenir un ensemble dessai de
variables entires. Ensuite, lensemble de variables continues est optimis pour des valeurs fixes des
variables entires en utilisant un code NLP standard. Le processus est recommenc jusqu ce que plus
aucun progrs ne soit observ.
- 4.19 -
Chapitre 4
En conclusion, la solution de problmes d'optimalisation non-linaire de grande chelle et

numriquement mixte (entier et non-entier) est toujours un dfit et une solution pratique peut rarement tre
obtenue sans une simplification de la formulation du problme.
4.4.7 Optimalisation dans les logiciels de calcul de procds
Lexprience est ncessaire pour valuer la valeur des algorithmes disponibles et de formuler
efficacement par la bonne mthode les problmes doptimalisation. Un des problmes majeurs dans
lapplication des techniques doptimisation dans les procds chimiques dingnierie tait dadapter les
logiciels de calcul de procds existants pour permettre non seulement les simulations mais galement pour
tendre leur aptitude optimaliser.
La plupart des logiciels de simulation sont bass sur des approches squentielles modulaires. Le
procd modliser est dcompos en briques standards de construction correspondant aux principaux
types doprations unitaires. Les librairies de simulation fournissent des sous-programmes pour modliser
le comportement de ces units en se basant sur les bilans de matire, dnergie et sur quelques quations
empiriques efficaces. Les units sont rsolues de manire squentielle en suivant le trajet de la matire et de
lnergie dans le procd : les sorties des units sont calcules partir des entres et de paramtres. Les
recyclages et les boucles de contrle sont rsolus de manire itrative en commenant avec des valeurs
estimes pour les flux coups. Des bloques mathmatiques sont utiliss pour rsoudre les quations de
coupure (avec possibilit de certaines quations supplmentaires) en faisant correspondre les valeurs
estimes avec les rsultats obtenus en bouclant travers la squence dunits. Puisque les drives sont
rarement disponibles, la promotion de la convergence est gnralement base sur la substitution ou sur des
amliorations bases sur la plus grande valeur propre ou sur lextrapolation de Wegstein. Les mthodes
semblables Newton essayant de construire de manire itrative les approximations de la matrice
jacobienne en utilisant une mthode de mise jour (par exemple, Broyden) ont galement t employes
avec succs.
Cette approche souffre deux inconvnients majeurs quand on essaya de rsoudre des problmes
doptimisation : les drives ne sont pas toujours disponibles et la manipulation des contraintes nest pas
toujours directe.
La premire implmentation reposait sur des mthodes directes (Nelder-Mead par exemple). Le
module doptimisation tait trait comme une promotion de convergence ou un bloque de contrle. Il fut
complt avec des modules d'chantillonnage (pour obtenir la valeur de la fonction objectif partir des
rsultats de simulation) et des modules de consigne (pour copier la valeur des variables optimises dans les
paramtres de simulation).
Cette approche est une mthode trajectoire ralisable : chaque itration, un calcul complet du
flowsheet doit avoir converg. Puisque les mthodes directes requirent un grand nombre ditrations, cette
- 4.20 -
Chapitre 4
approche ncessite un montant excessif de ressources informatiques (plusieurs centaines dquivalent

temps de simulation) et est donc limite aux problmes bien conditionns et non contraints impliquant un
nombre faible de variables de dcision.
Un moyen plus efficace de saisir les modles de procd avec les routines doptimisation est une
approche modulaire simultane. Quand cette technique est implmente, une routine de programmation non
linaire contrainte est applique pour simultanment optimiser la fonction objectif, vrifier les contraintes
dgalits (spcifications) et dingalits, et faire converger toutes les variables de flux coupes (galits
supplmentaires). Les algorithmes SQP ont montr leur succs dans cette application et loptimalisation
peut gnralement tre effectue avec un cot numrique acceptable (3 10 quivalents temps de
simulation).
La rsolution base sur des quations est une alternative dvelopper. Le modle du procd est
reprsent par un large systme dquations impliquant des bilans de matire et dnergie, des lois de
fonctionnement des courbes opratoires des quipements, des spcifications et des proprits physiques si
possible. Lensemble des quations est rsolu comme un tout aprs un rarrangement possible. Les
drives des quations sont gnralement disponibles. Appliquer des optimisations grande chelle ces
systmes de modlisation ne requiert pas de modification importante de la structure de la routine : rsoudre
un large ensemble dquations ou les utiliser comme un ensemble de contraintes dans une optimalisation
ncessitant la mme information.
Lapproche par quation oriente conduit gnralement des calculs trs efficaces, mais le
dveloppement et lutilisation de ces routines ncessitent plus dexprience que lapproche modulaire. Pour
dire vrai, l'nonc des quations du modle et des contraintes peut tre fort dlicat : les quations doivent
tre indpendantes. Une slection impropre des variables indpendantes peut accrotre la non-linarit.
Certains phnomnes introduisent naturellement des discontinuits (changement de phase, transition dun
coulement laminaire un coulement turbulent, etc.), et peut mme ncessiter de changer le nombre
dquations. Certains modles ne se transposent pas aisment en des systmes dquations algbriques
(par exemple, les quations diffrentielles dans les racteurs o lintgration numrique pas variable peut
tre ncessaire). Une mise lchelle des variables et des quations est trs importante pour lefficacit de
la solution : linteraction de lutilisateur pour suppler lapplication de facteurs dchelle spcifiques est plus
efficace que les mthodes automatiques, mais requiert plus dexprience. Des estimations initiales doivent
tre donnes pour toutes les variables du modle, les mauvais choix pouvant entraver la convergence.
Recourir quelques sries dlimination dans les quations du modle peut tre opportun pour
linitialisation (ceci correspond en fait la premire itration qui utilise un calcul squentiel modulaire).
Lapproche modulaire simultane semble prvaloir dans les logiciels commerciaux actuels puisquils
sont plus faciles programmer et dbuguer que les ensembles dquations, et comprendre pour
lutilisateur. Mais les dveloppements dans lapproche de la rsolution dquations va srement les rendre
plus efficaces et flexibles.
- 4.21 -
Chapitre 4
4.4.8 Applications
Les exemples dapplications de loptimalisation en calcul des procds sont nombreux.
Loptimalisation est employe en conception des installations pour dimensionner correctement les
quipements et faire attention un compromis entre les investissements et les cots opratoires. Elle peut
tre aussi utilise lors du fonctionnement d'une installation pour minimiser les matires premires et les
utilitaires pour une production donne. Elle peut tre exploite en planification o les quipements et le
autres ressources doivent tre alloues plusieurs tches comptitives. On peut aussi tirer profit de
loptimisation en formulation dun produit o la spcification de qualit doit tre rencontre cot le plus
bas en mlangeant des produits intermdiaires.
4.4.8.1 Fonctionnement optimal dun racteur de synthse dammoniac
Le racteur de synthse dammoniac est exothermique, limit lquilibre est favoris par des
pressions leves. La raction a lieu sur un catalyseur solide et sa vitesse est extrmement dpendante de la
temprature. Dans une conception de racteur, le catalyseur est spar en quatre lits. Lalimentation froide
est galement divise en quatre parts. La premire est prchauffe par change thermique avec leffluent du
racteur et doit atteindre une temprature suffisamment leve pour que la raction ait lieu une vitesse
significative dans le premier lit. La vitesse de raction dcrot suite lapproche de lquilibre (la
conversion lquilibre est diminue suite llvation de la temprature). Leffluent du premier lit est
mlang avec une alimentation frache. Le mlange rsultant est plus froid et contient moins dammoniac :
ces deux facteurs permettent davantage de raction davoir lieu dans le second lit jusqu ce que la
conversion dquilibre soit presque atteinte. Un processus dinjections similaire est rpt aprs les
deuxime et troisime lits. Leffluent du dernier tage de raction quitte le racteur et prchauffe
lalimentation du premier tage.
Racteur de synthse dammoniac : 4 lits catalytiques avec injections intermdiaires
- 4.22 -
Chapitre 4
Les variables continues qui apparaissent en exploitation optimale du modle sont, pour une
composition dalimentation donne, le dbit lentre, la pression, la temprature, la temprature de
prchauffe et les fractions respectives dalimentation admises lentre de chacun des quatre lits
catalytiques. La fonction objectif est manifestement la production dammoniac et doit tre maximise.
Pour une conception optimale, le nombre de lits est une variable de dcision importante. La
quantit de catalyseur chaque tage est galement une variable cl. Elle influera sur la taille du racteur.
Le dimensionnement du racteur (rapport de la longueur au diamtre et lpaisseur des parois) sera
influence par la quantit de catalyseur, la pression et la temprature. La gomtrie du racteur influencera
galement la chute de pression et, ainsi, sur le cot de la recompression. Lconomie de la boucle de
synthse sera galement affecte par le dbit de purge et lefficacit de la section de rcupration de
lammoniac qui contrlera la quantit dinertes et dammoniac dans lalimentation du racteur. Mais ceci
implique un modlisation plus complte qui prend en compte beaucoup plus que le racteur lui-mme.
4.4.8.2 Distillation : compromis entre le reflux et la taille des colonnes
Pour laborer une colonne, la dcision fondamentale devra tre le choix du nombre dtages
dquilibre de la colonne, de la pression de fonctionnement et du taux de reflux. Le premier facteur
affectera les investissements : aprs la slection dune technologie pour la configuration interne (colonne
plateaux ou garnissages), cela fixera la hauteur de la colonne. Le taux de reflux dterminera les charges
au bouilleur et au condenseur. Cela influera non seulement sur le besoin en utilitaires (vapeur, eau de
refroidissement), mais aussi sur les investissements travers la dimension des quipements. Un plus grand
taux de reflux augmentera galement le diamtre de la colonne. Un accroissement de la pression rduira le
diamtre de la colonne mais accrotra la difficult de la sparation. Toutes les dcisions prcdentes
influenceront la puret du produit et sa rcupration. Par consquent, le problme peut tre formul comme
un problme contraint de programmation non-linaire. Le nombre de plateaux est videmment une variable
discrte, mais un modle raccourci peut tre rsolu en utilisant un nombre dtages non entier. Ltage
dalimentation doit aussi tre localis de manire optimale et dpend de la composition de lalimentation.
4.4.8.3 Fonctionnement dun rseau vapeur
Pour les grandes installations, lnergie est apporte aux diffrents procds travers un rseau
vapeur o des distributeurs offrent la vapeur diffrents niveaux de pression. La vapeur est fournie par des
gnrateurs, mais est galement fabrique avec les calories excdentaires dans des bouilleurs une vitesse
fixe par le niveau de production. La vapeur est utilise pour chauffer, mais aussi dans des turbines o elle
produit du travail mcanique en se dtendant. Loffre et la demande doivent tre gales tous les niveaux
de pression soit en utilisant des turbines contre-pression, soit en dtendant de la vapeur haute pression
vers un niveau infrieur. Les moteurs lectriques sont une alternative aux turbines un cot unitaire de
lnergie diffrent. Etant donn un niveau de production de tous les procds, le travail mcanique et
l'apport des calories excdentaires peuvent tre calculs. La fonction objectif minimiser est le cot total
- 4.23 -
Chapitre 4
de llectricit achete et du fuel brl pour la gnration de vapeur. Les contraintes sont le bilan vapeur
chaque niveau de pression et les demandes en chaleur et nergie. Des quations de modlisation sont
disponibles pour les gnrateurs de vapeur, les vannes de laminage et les turbines. Les variables binaires de
dcision sont affectes par chaque unit qui peut tre mise en marche ou arrte : ainsi, l'optimaliseur aura
dcider quels sont les gnrateurs de vapeur mettre en route et si le travail mcanique pour les diffrents
appareils sera apporte par un moteur lectrique ou une turbine vapeur.
Ce type de problme peut se rsoudre par une mthode dapproximation deux niveaux (outer
approximation method). Un ensemble initial de variables entires est choisi et le problme non-linaire
rsultant est rsolu pour des valeurs constantes des variables entires slectionnes (par exemple, en
utilisant la SQP). Les contraintes et la fonction objectif sont linarises autour de la solution sous-optimale
et la formulation rsultante est rsolue avec un algorithme MILP (par exemple, branch & bounds). Le
processus itratif est rpt jusqu ce que lensemble des variables entires ne soit plus modifi par MILP
(mixed-integer linear programming). Si la fonction de cot linarise sous-estime le cot non-linaire, la
solution NLP est utilise comme un une limite suprieure dans la stratgie branch & bounds et tous les
solutions linaires plus grandes peuvent tre limines.
4.5 Quelques rfrences utiles

D. M Himmelblau, Applied Nonlinear Programming , McGraw-Hill Book Company (1972).
T.F. Edgar, D.M. Himmelblau, Optimization of Chemical Processes , McGraw-Hill Book
Company (1988).
H.A. Taha, Operations research, an introduction , Macmillan Publishing Co (1971).
Rfrences supplmentaires
M.J.D. Powell, An efficient method for finding the minimum of a function of several variables
without calculating derivatives , Computer Journal 7 :155 (1964) ; 303 (1965).
R.P. Brent, Algorithms for Minimization Without Derivatives , Prentice-Hall, Englewood Cliffs,
New Jersey, 1973.
- 4.24 -
Chapitre 4
P.E. Gill, W. Murray, M.H. Wright, Practical Optimization , Academic Press, NY (1981).
L.S. Lasdon, A.D. Warren, Large scale nonlinear programming , Computers and Chemical
Engineering, 7(5) 595-604 (1983).
(ce sujet est entirement dvolu loptimalisation grande chelle et contient plusieurs bonnes
revues)
- 4.25 -
Chapitre 5
CHAPITRE 5
RECONCILIATION DES DONNEES : VALIDATION
5.1 Introduction
Les donnes dun procd sont le fondement sur lequel sappuient tout contrle et toute valuation
de ses performances. La fiabilit des donnes est trs importante lorsque celles-ci sont destines au
monitoring du procd (contrle, identification, ligne doptimalisation). Donc, pralablement la simulation,
loptimalisation ou au revamping dune grande usine, nous avons besoin de donnes cohrentes pour
reprsenter le procd et en identifier correctement les paramtres.
La VALIDATION ou rconciliation des donnes est une tche trs importante qui transforme
lensemble des donnes disponibles en un set cohrent dfinissant ltat du procd. Rcemment, pour
assurer le contrle des procds, on a install beaucoup dordinateurs digitaux dans les raffineries et les
complexes chimiques. Ainsi, on dispose prsent dun grand nombre de donnes, rassembles et stockes
sur ordinateur, qui peuvent tre valides systmatiquement au moyen dun programme ad hoc, qui accrot
leur prcision et assure leur cohrence.
5.1.1 Origine des erreurs de mesure
Les mesures dun procd ne sont jamais cohrentes. Les principales raisons sont reprises cidessous :
Il y a des perturbations dues linstabilit du procd mme si le systme de contrle est trs efficace.
Certaines conditions (comme les conditions atmosphriques) ne peuvent tre contrles.
Les appareils de mesure ne sont pas toujours fiables. Il se peut que les biais instrumentaux ne soient pas
compenss de manire adquate; les appareils de mesure peuvent mal fonctionner.
Les lectures des mesures et les manipulations (analyses en laboratoires) peuvent introduire des erreurs.
Le point exprimental peut tre influenc par des lments indsirables et la mesure ne correspond pas
la variable attendue (mauvaise position dun thermocouple, influence de la rpartition du dbit dans un
changeur de chaleur, effet dun condenst dans un flux de vapeur, salets sur un appareil de mesure).
Des accidents peuvent modifier les bilans attendus dun procd (pertes, clatement dun changeur de
chaleur, pertes thermiques,...).
- 5.1 -
Chapitre 5
5.1.2 Dfinition de ltat dun procd

Dfinir ltat dune installation industrielle ncessite de prendre beaucoup de mesures de variables
dcrivant le systme telles que la temprature, la pression, le dbit, la composition, etc.
Lingnieur sait trs bien que, quand le nombre de mesures est infrieur un seuil donn (que nous
appellerons nombre de spcifications), il nest pas possible de dfinir ltat du systme.
Au plus, lingnieur est-il alors capable de dfinir ltat dun sous-systme. Pour autant que ces
mesures soient choisies soigneusement et que leur nombre soit gal au nombre de spcifications, lingnieur
peut exploiter les quations de bilan matire et thermique pour calculer les autres variables du systme.
Dans de telles circonstances, une erreur systmatique de mesure, si minime soit-elle, peut biaiser,
parfois de manire considrable, les calculs des autres variables du systme (cest--dire y introduire de
grandes erreurs).
En consquence, lingnieur est forc de prendre des mesures supplmentaires pour augmenter le
niveau de confiance escompt pour le systme.
Il notera alors quil nest plus possible de satisfaire toutes les quations du systme, et se reportera
une technique personnelle danalyse des mesures.
La mthode de rconciliation des donnes dveloppe dans le logiciel BELSIM, aide rsoudre
ce problme considrable de manire univoque.
Pour valider les donnes, lingnieur doit dfinir, dans la Banque de Donnes du Procd (BDP), le
procd, les mesures et les estimations des carts standard (erreurs) affectant chaque mesure. A cette fin, il
utilise le programme BIP (Belsim Interactive Program).
5.1.3 Illustration
Dans cet exemple, nous considrerons les variables de matire et thermiques; les quations prises
en compte seront respectivement les bilans de matire et thermique. Il nest pas toujours ncessaire, ni
possible, denvisager les deux simultanment : cela dpend des donnes disponibles et de leur prcision.
Le mlange de deux flux contenant une substance pure est illustr la figure 1. La temprature et la
pression des flux sont diffrentes. Le nombre de spcifications ncessaires pour dfinir explicitement ltat
du systme est de 7.
- 5.2 -
Chapitre 5
T1 P1 D1
T3 P3 D3
T2 P2 D2
FIGURE 1
Variables :
temprature (T)
pression (P)
dbit (D)
3
3
3
total
bilan matire
bilan thermique
1
1
total
Equations :
Nombre de degr de libert = 9 variables - 7 quations = 7

Si on effectue moins de 7 mesures, il nest pas possible de dterminer compltement ltat du
systme.
1er cas : 6 variables sont mesures : P1, P2, P3, T1, T2, D3.
Dans ce cas, il nest pas possible den savoir plus sur le systme.
2me cas : 6 variables sont mesures : P1, P2, P3, T1, D2, D3.
Il est dsormais possible de calculer le dbit D1, en utilisant lquation de bilan matire. Un sous-systme
(D1) des variables non mesures est calculable. Les tempratures T2 et T3 ne sont pas calculables.
3me cas : 7 variables sont mesures : P1, P2, P3, T1, T2, T3, D1.
Elles valent : T1=300K, T2=410K, T3=310K, P1=P2=P3=1bar, D1=10kg/s.
A laide des quations de bilan matire et thermique, nous pouvons calculer les deux variables non
mesures : D2 et D3. Le systme est prsent compltement connu mais une erreur systmatique dans la
mesure de la temprature T3 peut biaiser dans une large mesure la connaissance du systme.
Supposons le Cp de la substance constant. Le tableau 1 nous fournit les valeurs calcules de D2 et
D3 pour diffrentes valeurs de T3.
- 5.3 -
Chapitre 5
Tableau 1
T3 (K)
D2 (kg/s)
D3 (kg/s)
310
305
315
1
0.48
1.58
11
10.48
11.58
Une erreur systmatique de 5 degrs sur la mesure de T3 nous donnera une trs mauvaise
valuation de D2.
5.2 Bases Theoriques

Ce chapitre a pour but de nous fournir un aperu, aussi clair que possible, des bases thoriques de
la validation. Il explique les techniques et le formalisme utiliss. Les variables utilises sont explicites, de
mme que la notion dquation de contrainte et dquation de liaison. Vous apprendrez analyser les
matrices ainsi que la mthode de rsolution.
5.2.1 Dfinitions
5.2.1.1 Variables dtat
Les variables dtat sont des variables impliques dans lcriture des quations associes aux units
physiques du procd. Les quations concernes sont principalement les quations de bilan de matire et
de chaleur.
Dans le logiciel BELSIM, les quations sont exprimes en fonction des variables dtat suivantes :
T : la temprature du mlange, si le flux matire est impliqu dans un bilan thermique.
P : la pression du mlange, si elle joue un rle dans le bilan thermique et si elle est considre comme
telle; trs souvent, il nest pas possible de valider les mesures de pression et celles-ci sont considres
constantes.
Ci (i=1,nombre de substances) : le dbit molaire partiel de la substance i dans le mlange (kmol/s).
Uj (j=1,nombre de ractions) : degr davancement de la raction j, si celle-ci intervient dans une unit
physique du systme.
FRA : la fraction darbre associe un flux, pour autant que le flux ne contienne que la fraction dun
seul arbre (dans la structure de donnes BELSIM, les dbits molaires partiels associs un flux matire
constituent un arbre) ; la variable FRA est effectivement utilise lorsquun flux est divis en plusieurs
autres de mme composition : dans ce cas, les flux de sortie du diviseur auront le mme arbre associ
que celui dentre, mais leur FRA sera infrieure 1.
- 5.4 -
Chapitre 5
NOTE : La variable temprature est avantageusement remplace par la variable enthalpie molaire
H lorsquon se trouve en prsence dun flux biphasique une seule substance. Cette dernire dfinit
univoquement ltat thermodynamique dun tel systme.
Toutes les autres grandeurs (enthalpie, fraction molaire, fraction pondrale sur extrait sec, etc.)
peuvent toujours se dduire des variables dtat (pour autant que les mthodes thermodynamiques
dcrivant le mlange soient dfinies).
Ltat du systme est connu une fois que les valeurs des variables dtat sont connues pour chaque
flux du systme.
Le flow-sheet prsent la figure 2 se compose de 4 units physiques et de 9 flux. La liste des
variables dtat associes chaque flux et chaque arbre sont donnes au tableau 2. Pour spcifier la
conversion du CH4, la variable dtat du racteur U1 est ajoute.
MELANGEUR
ECHANGEUR
DE CHALEUR
REACTEUR
F1
A1
F3
F2
A3
CH4+2O2
F4.1
CO2+2H2O
A2
F4.3
F4.2
A4
F4.4
F5.1
FLUX
ARBRES
DIVISEUR
F5.2
FIGURE 2
A5
Tableau 2
Flux
Variables
Arbres
Variables
F1
F2
F3
F4.1
F4.2
F4.3
F4.4
F5.1
F5.2
T,P
T,P
T,P
T,P
T,P
T,P,FRA
T,P,FRA
H,P
H,P
A1
A2
A3
A4
O2,N2
CH4,N2
O2,N2,CH4
O2,N2,CH4,CO2,H2O
A5
H2O
5.2.1.2 Variables mesures ou observes

Les lments en question sont ceux mesurs dans lusine. Chaque mesure est accompagne dune
estimation de lcart standard de lerreur de mesure, et est considr comme une variable du systme.
- 5.5 -
Chapitre 5
Certaines de ces variables sont dj reprises comme variables dtat du systme. Les autres,
appeles variables de liaison, sont relies aux variables dtat par des quations spcifiques : les quations
de conversion ou quations de liaison. Cest par exemple le cas du dbit molaire total DMOL qui est reli
aux dbits molaires partiels par lquation :
DMOL = C k
k
ou, de la fraction molaire FMk de la substance k :

FM k . Ci = Ck
i
Pour des applications spcifiques, lingnieur peut dcider que certains lments observs soient
considrs comme constants dans le problme de validation. Dans ce cas, les variables concernes sont
limines du systme.
5.2.1.3 Equations de contrainte
Les quations de contrainte sont : les quations de bilan matire et dnergie, les quations de
liaison, mais aussi des quations spciales (quations dquilibre liquide-vapeur, quations dgalit de
pressions ou de tempratures, etc.).
Le mlangeur illustr la figure 2 a deux flux dentre, dont les tempratures sont TF1 et TF2. Les
variables O2A1 et O2A2 sont les dbits molaires partiels du premier flux dentre (arbre A1). Les variables
N2A2 et CH4 A2 reprsentent les dbits molaires partiels du second flux dentre (arbre A2). Le premier
flux de sortie est caractris par la temprature TF3 et les dbits molaires partiels O2A3, N2A3 et CH4 A3.
Les pressions sont considres constantes.
Il faut crire trois bilans de matire :
O2A1 - O2A3 = 0
pour la substance O2
N2A1 + N2A2 - N2A3 = 0
pour la substance N2
CH4 A2 - CH4 A3 = 0 pour la substance CH4

Lquation de bilan thermique se traduit par le bilan enthalpique :
h(TF1,O2A1,N2A1) + h(TF2,N2A2,CH4 A2) - h(TF3,O2A3,N2A3,CH4 A3) = 0
avec h = fonction calculant le dbit enthalpique en fonction de la temprature et des dbits molaires
partiels.
Lquation de liaison suivante doit tre crite si la fraction molaire de la substance O2 a t mesure
dans le flux de sortie.
FMO2A3 . ( O2A3 + N2A3 + CH4 A3) - O2A3 = 0
- 5.6 -
Chapitre 5
Donc, les variables dtat du problme sont : TF1, TF2, TF3, O2A1, N2A1, N2A2, CH4 A2, O2A3,
N2A3, CH4 A3, tandis que FMO2A3 est une variable de liaison.
5.2.2 Position du problme
La technique de rconciliation se base sur lhypothse suivante. Toutes les mesures sont entaches
derreurs et les valeurs corriges ou valides diffrent des valeurs mesures. Dune part, les valeurs
valides doivent vrifier les quations de contrainte, et, dautre part, elles doivent minimiser la somme des
carrs des diffrences entre les valeurs valides et les valeurs mesures.
Cela va sans dire que de pareilles diffrences sont pondres par les carts standard
correspondants.
Dun point de vue mathmatique, ce problme est un problme de minimisation avec contraintes,
qui peut se dfinir de la manire suivante.
Soient :
y, le vecteur des valeurs mesures, de taille MES ;
Y, le vecteur des valeurs corriges ou valides ;
X, le vecteur des valeurs non mesures calculer, de taille NMES ;
F(X,Y), la fonction vectorielle dquations de bilans, de taille NEQ ;
P, la matrice de pondration qui permet de quantifier la prcision relative des mesures MES. En pratique,
cest une matrice diagonale, dont chaque lment est linverse de la variance si2 de la mesure i.
Le problme revient minimiser :
MES
i =1
(Y
yi )
sous la contrainte : F( X , Y) = 0 .
En notation matricielle (lexposant T signifie transpos ) :

Min ( Y y) P ( Y y) avec F( X, Y) = 0
X, Y
T
Ce problme de minimisation avec contraintes est rsolu par la mthode de LAGRANGE en

faisant intervenir les multiplicateurs de Lagrange pour chaque quation. Le critre scrit alors (le facteur
2 est introduit par souci de commodit, comme nous le verrons plus loin) :
Min L avec L = ( Y y ) P (Y y) + 2 F( X, Y)
X, Y,
T
- 5.7 -
Chapitre 5
La rsolution des quations dEuler fournit la solution du problme. Les quations dEuler
scrivent :
L
=0
Y i
i=1,MES
L
=0
X i
i=1,NMES
L
i
=0
i=1,NEQ
Les matrices jacobiennes des variables mesures et des variables non mesures du systme sont :
ij =
F i
Y j
ij =
F i
X j
On obtient le systme suivant dquations :
(Y y)
P + T A = 0
MES quations
T B = 0
NMES quations
F( X , Y) = 0
NEQ quations
savoir un systme de MES+NMES+NEQ quations non linaires MES+NMES+NEQ inconnues.

Une tape ncessaire toutes les mthodes de rsolution conues jusqu ce jour consiste valuer la
matrice jacobienne totale de ce systme dquations qui scrit, si on nglige la dpendance de A et de B
vis--vis de Y et X :
P
J = 0
A
0 AT
0 BT
B 0
MES lignes
NMES lignes
NEQ lignes
Considrons la validation des mesures de dbit autour du mlangeur de la figure 2. Les valeurs des
mesures et les carts standard qui leur sont associs sont donns dans le tableau 3.
Tableau 3
Variables
Y1
Y2
Y3
Y4
Y5
X1
Mesures
1
2
1
4
0.2
-
O2A1
N2A1
N2A2
CH4 A2
FMO2A3
O2A3
- 5.8 -
Ecart. Std.
0.01
0.01
0.01
0.01
0.01
-
Chapitre 5
X2
X3
N2A3
CH4 A3
( Y1 1) 2 ( Y2 2 ) 2 ( Y3 3) 2 ( Y4 4) 2 ( Y5 0.2) 2
MIN
+
+
+
+
2
2
2
2
2
0
.
01
0
.
01
0
.
01
0
.
01
0
.
01
X, Y
sous les contraintes
Y1 X1 = 0
1
0
A=
0
Y2 + Y3 X2 = 0
0
1
0
0
0
1
0
0
0
0
1
0
0
0
0
X i
Y4 X3 = 0
Y5 ( X1 + X2 + X3) X1 = 0
0
1 0
0 1 0
B=
0
0 1
Y5 Y5 Y5
5.2.3 Existence dune solution

En observant les quations de contrainte, nous constatons que :
1. Le systme accepte un nombre infini de solutions si le nombre de variables non mesures est suprieur
au nombre dquations de contrainte (NMES>NEQ). Dans ce cas particulier, la matrice B est
rectangulaire et horizontale (plus de colonnes que de lignes) et la matrice jacobienne est singulire,
puisque les lignes MES+1 MES+NMES sont systmatiquement linairement dpendantes. Le
systme nest pas soluble.
2. Si le nombre de variables non mesures est gal au nombre dquations de contrainte (NMES=NEQ),
la solution du problme est obtenue en considrant les mesures comme des constantes et en calculant
les variables non mesures laide des quations de contrainte. On dira que le systme est juste
calculable car il ny a pas assez dquations pour corriger les mesures. La matrice B est carre.
3. Si le nombre de variables non mesures est infrieur au nombre dquations de contrainte
(NMES<NEQ), le systme admet une solution unique. Les quations de contrainte sont non seulement
utilises pou calculer les variables non mesures, mais aussi pour rconcilier les mesures. La matrice B
est rectangulaire verticale (plus de lignes que de colonnes).
On peut dj affirmer quil faut, pour quil y ait une solution, que la matrice B ne soit pas
rectangulaire horizontale. En fait, il faut que la matrice jacobienne globale ne soit pas singulire, ce qui
arrive par exemple,
si les quations de contrainte sont linairement dpendantes,
si les mesures sont mal rparties et que certaines parties du procd restent indtermines,
si les variables rendues constantes sont mal choisies et crent des surspcifications.
- 5.9 -
Chapitre 5
Dans la majorit des cas, une analyse soigne de la matrice dincidence du systme permet de
dtecter ces problmes pralablement la rsolution. (La matrice dincidence dun systme dquations est
la matrice dont les lments "ij" valent 1 si la variable "j" intervient dans lquation "i", 0 sinon).
La matrice doccurrence des variables mesures et non mesures du mlangeur (figure 2, tableau
3) est :
1
0
AB =
0
0 0 0 0 1 0 0
1 1 0 0 0 1 0
0 0 1 0 0 0 1
0 0 0 1 1 1 1
ATTENTION : mme si cette analyse est satisfaisante, le systme dquations, les points de dpart et les
valeurs mesures peuvent tre tels que le programme ne converge pas vers une solution,
soit que celle-ci se trouve en dehors du domaine admissible (dfini par les intervalles de
variation physique des variables : tempratures et dbits positifs,...), soit que le systme
devienne localement singulier en cours ditrations. Ce type de problme est explicit
avec plus de dtails au chapitre traitant des problmes de convergence.
5.2.4 Analyse des matrices doccurence
5.2.4.1 Validit des quations de contrainte
Il est vident que lanalyse dun systme dquations au moyen de sa matrice dincidence est
dlicate. En effet, deux quations identiques apparatront comme deux quations faisant intervenir les mme
variables et non comme une erreur.
De mme, un sous-ensemble dquations linairement dpendantes ne saurait tre mis en vidence
dans une matrice dincidence. Les chances de rencontrer de telles singularits sont pour la plupart limines
si les quations sont gnres automatiquement par le programme. La vrification, aise, de certaines
conditions limine la gnration derreurs dans les quations. Par exemple, chaque flux ne peut tre
connect plus de deux units physiques; il est ncessaire davoir au moins un flux dentre et un flux de
sortie.
5.2.4.2 Y a-t-il assez de mesures?
Aprs gnration des mesures et des quations de liaison, les variables sont rparties en trois
catgories :
les variables spcifies constantes (appeles galement constantes),
les variables mesures,
- 5.10 -
Chapitre 5
les variables non mesures.

Pour dterminer si les mesures permettent, au moins, de dterminer ltat du systme, on analyse la
matrice dincidence quations-variables non mesures correspondant la matrice B. La technique de
dtection consiste mettre en vidence, par permutation des lignes et des colonnes, un sous-systme S2
contenant une sous-matrice B2' et une sous-matrice horizontale B2" (figure 3). La sous-matrice B2" est
horizontale. Il ny a pas assez dquations contenant les variables non mesures relatives la matrice B2".
Ds lors, les variables non mesures relatives B2" ne peuvent tre dtermines (elles interviennent dans
trop peu dquations).
Pour rendre le systme calculable, une mesure supplmentaire doit tre choisie parmi les variables
de la sous-matrice B2". Cette variable sera ajoute lensemble des variables mesures et la colonne
correspondante de la sous-matrice B2" sera limine.
Lorsquune mesure est choisie, on procde de nouveau la recherche dun sous-systme non
calculable.
S1
S2
Equations
Variables non mesures
B1
S1
sous-systme validable
(B1 vertical)
B2'
B2''
S2=B2'+B2''
sous-systme incalculable
(B2'' mat. horizontale)
variables candidates
FIGURE 3
5.2.4.3 Les mesures sont-elles validables?

Considrons maintenant que nous avons effectu un nombre de mesures supplmentaires suffisant
pour que la rsolution soit possible (redondance positive ou nulle).
Analysons de nouveau la matrice dincidence et permutons lignes et colonnes, ou vice et versa, en
essayant disoler un sous-systme dquations S2 contenant une sous-matrice B2' et une sous-matrice
carre B2" (figure 4).
Les variables non mesures associes la matrice B2" sont juste calculables puisquil y a juste
assez dquations pour les calculer.
- 5.11 -
Chapitre 5
S1
B1
S2
Equations
Variables non mesures
B2'
B2''
S1
(B1 vertical)
S2=B2'+B2''
sous-sytme juste calculable
(B2'' matrice carre)
variables juste calculables
FIGURE 4
Examinons maintenant la matrice dincidence globale correspondante quations-variables

mesures et non mesures et tassons gauche les colonnes des variables mesures intervenant dans les
quations du sous-systme S1 (figure 5).
Les variables mesures restantes ne sont pas validables car leurs valeurs peuvent tre fixes
arbitrairement leurs valeurs mesures : aucune quation ne permet de les corriger.
Equations
S1
A1
S2
Variables
mesures
A2'
Variables
non mesures
A2''
variables
invalidables
B1
B2'
S1
(B1 mat. verticale)
B2''
S2
sous-systme juste calculable
(B2'' mat. carre)
variables juste
calculables
FIGURE 5
Exemple
La figure 6 prsente un procd contenant trois units. Les trois substances sont spares dans
lunit A avant dentrer dans un racteur B; lunit C spare les produits de la raction. Dans ce cas, nous
- 5.12 -
Chapitre 5
considrons quil ny a pas assez de mesures autour des units A et C. Nous ne gnrons alors que le bilan
thermique du racteur B alors que les bilans de matire sont gnrs pour toutes les units.
F5
3 substances : a, b, c
F3
F4
REACTEUR
B
a+b -> c
F1
F6
A
F2
FIGURE 6
F3
F4
F6
F5
F4
F3
F2
a
b
c
a
b
c
a
b
c
a
b
c
U
a
b
c
a
b
c
T
T
F1
La figure 7 montre la matrice doccurrence correspondant au flow-sheet de la figure 6.
a X
X
X
A b
X
X
X
c
X
X
X
a
X
X
X
Bilan Matire
B b
X
X
X
c
X
X X
Bilan Thermique B
X X X X X X
X X
a
X
X
X
Bilan Matire
C b
X
X
X
c
X
X
X
Bilan Matire
FIGURE 7
- 5.13 -
Chapitre 5
F6
F5
F2
F4
F3
F1
VARIABLES NON MESUREES
T
T
FMa
FMa
FMb
DMOL
DMOL
FMa
FMb
a
b
c
a
b
c
a
b
c
U
a
b
c
a
b
c
a
b
c
F3
F4
F3
F4
F4
F4
F6
F5
F5
VARIABLES MESUREES
Fract.Mol
Fract.Mol
Fract.Mol
Db.Mol.Tot.
Bilan Therm.
Bilan Therm.
Bilan Therm.
Bilan Therm.
Db.Mol.Tot.
Fract.Mol
Fract.Mol
F3 a
X
F4 a
X
F4 b
X
F4
X
a
B b
c
B
X X
a
X
A b
X
c
X
a
C b
c
F6
X
F5 a
X
F5 b
X
X X X
X X X
X X X
X X X
X
X
X
X
X
X
X
X X
X X X X X X
X
X
X
X
X
X
X
X
X
X
X
X
X X X
X X X
X X X
FIGURE 8
La figure 8 montre trs clairement les mesures non validables quand les variables mesures sont :
la temprature des flux F4 et F5,
le dbit molaire total des flux F4 et F6,
la fraction molaire de la substance "a" dans les flux F3, F4 et F5,
la fraction molaire de la substance "b" dans les flux F4 et F5,
les dbits molaires partiels des substances "a", "b" et "c" dans le flux F1.
5.2.4.4 Na-t-on pas introduit des surspcifications?
Considrons maintenant que nous avons un systme dont la rsolution semble possible
(redondance positive ou nulle) et dans lequel nous avons introduit des constantes en associant certaines
mesures un cart standard trs petit.
Lors de la rsolution, les variables correspondant ces mesures seront ignores et considres
comme de relles constantes. Toutefois, lutilisateur a peut-tre mal choisi ces constantes et surspcifi le
problme.
Analysons la matrice dincidence quations-variables et permutons lignes et colonnes pour isoler
un sous-systme dquations S1 (figure 9) dont la matrice est verticale. Les variables associes cette
matrice sont surspcifies puisquil y a trop dquations pour les calculer.
- 5.14 -
Chapitre 5
S2
Equations
S1
Variables mesures et
non mesures
S1
sous-systme singulier
(AB1 mat. verticale)
surspcifi
AB1
AB2'
AB2''
FIGURE 9
Examinons la matrice dincidence globale correspondante quations - variables + constantes et

tassons gauche les colonnes des constantes intervenant dans les quations du sous-systme S1 (figure
10). Les constantes rendre variables devront tre choisies parmi ces colonnes.
S2
Equations
S1
Variables mesures
et non mesures
AB1
AB2'
Constantes
S1
sous-systme singulier
(AB1 mat. verticale)
surspcifi
AB2''
constantes candidates
FIGURE 10
5.2.4.5 Recherche des redondances triviales

Analysons la matrice dincidence quations-variables et permutons lignes et colonnes pour isoler
un sous-systme dquations S1 (figure 11) dont la matrice est carre.
- 5.15 -
Chapitre 5
S1
B2'
B2''
S1
sous-systme juste calculable
(B1 matrice carre)
Equations
B1
S2
Variables mesures et
non-mesures
variables juste calculables
FIGURE 11
Les variables associes cette matrice sont juste calculables puisquil y a juste assez dquations
pour les calculer.
Si parmi ces variables figurent des variables mesures, leurs valeurs calcules sont tout fait
indpendantes des valeurs mesures et ces mesures ne servent rien. Nous dirons quelles gnrent des
redondances triviales.
5.2.4.6 Dpendance linaire dans un sous-systme
Lorsquune partie du flow-sheet gnre un sous-systme juste calculable (sous-matrice carre, voir
figure 4), des variables non mesures, apparemment juste calculables, peuvent tre indterminables.
Cest le cas des boucles dans lesquelles aucune mesure nest donne mais pour lesquelles il est
ncessaire de mesurer un dbit de matire ou thermique pour valuer ce qui tourne dans la boucle. Une
boucle est identifie lorsquune suite de flux non orients forme une boucle. Cest le cas des recyclages
mais galement lorsquun mlange est effectu sur des flux pralablement spars.
Mathmatiquement, la matrice carre est numriquement singulire, ce qui ne peut se dtecter par
lanalyse des matrices dincidence. Cette singularit pourrait tre mise en vidence en transformant
lquation de bilan dune des units de la boucle en un bilan global autour de toutes les units de la boucle.
Exemple
- 5.16 -
Chapitre 5
x1
y1
y2
B
x2
x3
C
y3
FIGURE 12
Dbits mesurs
Dbits non mesurs
Equations de bilan de matire
- pour lunit A
- pour lunit B
- pour lunit C
y1, y2, y3
x1, x2, x3
y1 - x1 - x2 = 0
x1 + x3 - y2 = 0
x2 - x3 - y3 = 0
Matrices jacobiennes :
1 0
0
A = 0 1 0
0 0 1
1 1 0
B= 1
0
1
0
1 1
La matrice doccurrence quations-variables non mesures (matrice B) nest de toute vidence

pas singulire mais elle est numriquement singulire. En effet, une matrice doccurrence singulire est
obtenue en remplaant lquation de bilan matire de lunit C par lquation de bilan matire global.
Equation de bilan matire global :
y1-y2-y3=0
Matrices jacobiennes :
1 0
0
A = 0 1 0
1 1 1
1 1 0
B= 1
0 1
0
0 0
Une autre mthode serait de sparer la matrice carre en blocs triangulaires et de vrifier, pour
chacun deux, que toutes les variables et les quations dun bloc donn ne se rapportent pas un une
boucle : ce problme a t rsolu pour un recyclage de matire. Toutefois, il faut remarquer que dtecter
de telles singularits devient une tche trs difficile quand on est confront un recyclage de chaleur et de
matire.
- 5.17 -
Chapitre 5
Le flow-sheet de la figure 13 met en vidence un cas dindtermination matrielle et thermique. La

singularit est localise dans le bloc carr du coin droit de la figure 14. Sans nul doute, lanalyse dun tel
systme exige une certaine exprience en validation de flow-sheets.
F6'
F5'
F6
A6
F5
R2
a+b->c
A5
R1
a+b->c
F2'
F2
F1
F4
ECH1
F3
F3'
MEL
DIV
3 SUBSTANCES a b c
10 FLUX 1 2 2' 3 3' 4 5 5' 6 6'
3 ARBRES 1 5 6
ECH2
FIGURE 13
A6
R1
A5
F1
F4
F5
F5'
F6
F6'
VARIABLES NON MESUREES
a
b
c
T
T
T
T
T
T
a
b
c
U
a
b
c
U
T
T
Fr
Fr
T
T
F1
VARIABLES MESUREES
R2
F2
F3
F2
F3
F2'
F3'
A1
a X
X
R1
Bilan Matire b
X
c
X
R1
Bilan Therm.
X X X
X X
X
a
X
R2
Bilan Matire b
c
R2
Bilan Therm.
X X
X
Bilan Therm. 1
X
DIV
Bilan Therm. 2
X
Bilan Matire
ECH1 Bilan Therm.
X X X
X X
ECH2 Bilan Therm.
X X X
X X
X
MEL Bilan Therm.
X X X
X
X
X
X
X X
X X
X
X
X
X X
X
X
X X
X X X
X
X
X
X X X
X X
X X
X
X
X
X
X
X
X X X X
FIGURE 14
5.2.5 Mthode numrique utilise
Le systme est rsolu par une mthode DOGLEG. Ceci implique que la correction faite dans
lvaluation courante de la solution est une combinaison de la correction de Newton et de la direction de la
plus grande pente pour la somme des carrs des rsidus des systmes dquations.
- 5.18 -
Chapitre 5
Sans entrer dans les dtails, on peut dire que le passage vers la direction de la plus grande pente
est dautant plus important que la rduction de la somme des carrs des rsidus dans la direction de
Newton est petite, quand celle-ci nest pas totalement nulle.
Le calcul de la matrice jacobienne ncessaire pour dterminer ces directions est en partie
numrique et en partie analytique. Pour pargner du temps de calcul, cette matrice nest pas calcule
chaque itration mais selon une frquence choisie. Son calcul dpend aussi de la vitesse laquelle on
sapproche de la solution. La matrice est emmagasine dans la mmoire laide dune technique de
stockage pour les matrices parses. Seuls les lments de la matrice diffrents de zro sont stocks.
Nous vous rappelons ici que nous ngligeons les drives secondes des contraintes par rapport aux
variables mesures et aux variables non mesures. Si le nouveau point propos par la mthode se trouve
lextrieur du domaine de validit (tempratures et dbits positifs), le programme ramne le point
lintrieur du domaine de validit en relaxant le pas propos.
Selon le cas, toutes les variables sont sujettes la relaxation, ou seulement celles qui violent les
contraintes physiques.
Dans presque tous les cas, la non-convergence est attribuer la mauvaise qualit des mesures, ou
un manque de mesures qui na pas t dtect dans les analyses doccurrence.
5.3 Analyse de la sensibilit

5.3.1 Introduction
Nous avons vu que le problme de la validation peut sexprimer comme un problme de
minimisation contraint. Les dveloppements ultrieurs supposent que les contraintes sont linaires ou bien
quelles ont t linarises :
min ( Y y) P ( Y y)
X, Y
T
sous la contrainte AY + BX + C = 0
Nous avons galement vu que le problme avec contraintes peut tre transform en un problme
non contraint en utilisant la formulation de Lagrange :
min L = (Y y) P (Y y) + 2 T ( AY + BX + C)
X, Y,
T
Nous obtenons ainsi le systme dquations suivant :
- 5.19 -
Chapitre 5
+ AT
BT
PY
= Py
=0
= C
A Y + BX
On peut donc dfinir une matrice carre M et des vecteurs V et D tels que :
P
M = 0
A
0 AT
0 BT
B 0
Y

V = X

Py

D= 0
C
De cette manire, la solution du problme de validation scrit :

V = M 1 D
5.3.2 La matrice de sensibilit
La matrice M 1 est la matrice de sensibilit du systme. On remarque que les vecteurs X et Y sont
des combinaisons linaires des valeurs mesures y. La matrice de sensibilit nous permet donc dvaluer de
quelle manire la valeur valide dune variable dpend de lensemble des variables mesures et de leurs
carts standards.
En particulier :
m+ n + p
(M )
Yi =
j =1
ij
Dj
j=1
k =1
= ( M 1 ) ij Pjj y j ( M 1 ) i
m + n+ p
(M )
Xi =
j=1
n+ i j
n + m+ k
Ck
Dj
j =1
k =1
= ( M 1 ) n +i j P jj y j ( M 1 ) n+ i
n+ m + k
Ck
La variance dune combinaison linaire Z de plusieurs variables X j est calcule de la manire

suivante :
m
Z = a j Xj
j =1
( )
var ( Z) = a 2j var X j
j =1
Nous obtenons donc lestimation suivante pour la variance des variables mesures valides :
- 5.20 -
Chapitre 5
m
var (Yi ) = ( M 1 ) ij Pjj

j=1
} var(y )
2
et pour la variance des variables non mesures et estimes :

m
var ( Xi ) = ( M 1 ) n+ i j Pjj
j =1
( )
var y j
Les expressions prcdentes peuvent tre simplifies en tenant compte de :
( )
var y j = 1 P
jj
et donc :
var (Yi ) =
(M )
1 2
ij
var y
( )
j =1
var ( X i ) =
(M )
j=1
1 2
n +i j
( )
var y j
5.3.3 Illustration
Considrons un flowsheet simplifi dune boucle de synthse dammoniac (figure 1). La pression
dans le racteur passe de 280 270 bar et est compense par le compresseur de recyclage. Nous
supposerons quil ny a aucune autre perte de charge dans le systme. Les mesures sont inscrites sur le
flowsheet.
PURGE
T=20C
F=3.2t/h
0.2% CH4
4.8% Ar
7.2% NH3
64.4% H2
23.4% N2
T=20C
F=63 t/h
RECYC
FV
T=540C
RCTOUT
Q3
FLV
T=19C
T=25C
FL
F05
Q2
F06
T=40C
T=30C
W3
W=330 kW
F07
RCTIN
T=378C
F=77 t/h
0.3% CH4
3.7% Ar
5.4% NH3
67.1% H2
23.5% N2
F=12.5 t/h
0.1% CH4
0.3% Ar
97.8% NH3
1.3% H2
0.5% N2
Figure 1 : Flowsheet simplifi dune boucle de synthse dammoniac

Tableau 4 : rsultats de la rconciliation des donnes
STREAM
F05
MEASURED
40.000
STD.DEV
2.0000
- 5.21 -
ABS
RECON.VAL.
40.229
DIFF(%)
.572
PENALTY
.013
Chapitre 5
STREAM
STREAM
STREAM
STREAM
STREAM
STREAM
STREAM
STREAM
STREAM
STREAM
STREAM
STREAM
STREAM
STREAM
STREAM
STREAM
MIXTURE
MIXTURE
MIXTURE
MIXTURE
MIXTURE
MIXTURE
MIXTURE
MIXTURE
MIXTURE
MIXTURE
MIXTURE
MIXTURE
MIXTURE
MIXTURE
MIXTURE
MIXTURE
MIXTURE
F05
F06
F06
F07
F07
FLV
FLV
PURGE
PURGE
RCTIN
RCTIN
RCTOUT
RCTOUT
RECYC
RECYC
W3
FL
FL
FL
FL
FL
FL
FV
FV
FV
FV
FV
RCTIN
RCTIN
RCTIN
RCTIN
RCTIN
RCTIN
P
T
P
T
P
T
P
T
MASSF
T
P
T
P
T
MASSF
POWER
MASSF
MFCH4
MFAR
MFNH3
MFH2
MFN2
MFCH4
MFAR
MFNH3
MFH2
MFN2
MASSF
MFH2
MFCH4
MFAR
MFNH3
MFN2
bar
C
bar
C
bar
C
bar
C
t/h
C
bar
C
bar
C
t/h
kW
t/h
%
%
%
%
%
%
%
%
%
%
t/h
%
%
%
%
%
25.000
30.000
19.000
20.000
3.2000
378.00
540.00
20.000
63.000
330.00
12.500
.10000
.30000
97.800
1.3000
.50000
.20000
4.8000
7.2000
64.400
23.400
77.000
67.100
.30000
3.7000
5.4000
23.500
CONSTANT
2.0000
CONSTANT
2.0000
CONSTANT
2.0000
CONSTANT
2.0000
10.000
5.0000
CONSTANT
5.0000
CONSTANT
2.0000
10.000
30.000
5.0000
.10000
.10000
1.0000
.10000
.10000
.10000
.30000
.30000
1.0000
1.0000
10.000
1.0000
.10000
.30000
.30000
1.0000
ABS
ABS
ABS
ABS
%
ABS
ABS
ABS
%
ABS
%
ABS
ABS
ABS
ABS
ABS
ABS
ABS
ABS
ABS
ABS
%
ABS
ABS
ABS
ABS
ABS
270.00
24.631
270.00
29.309
280.00
20.435
269.00
20.544
3.2091
375.64
280.00
542.35
270.00
20.544
63.879
332.28
12.165
.18949E.17703
98.005
1.2784
.52113
.27156
4.7121
6.9320
64.871
23.214
79.253
66.849
.22207
3.8389
5.5084
23.582
-1.475
.034
-2.303
.119
7.552
.515
2.721
.286
-.624
.074
.001
.223
.436
.222
2.721
1.394
.692
-2.681
-81.051
-40.990
.209
-1.664
4.227
35.780
-1.831
-3.723
.731
-.796
2.925
-.374
-25.978
3.755
2.007
.348
.074
.019
.006
.288
.657
1.512
.042
.047
.045
.512
.086
.798
.221
.035
.086
.063
.607
.214
.130
.007
Lcart standard est fix 10% de valeurs mesures pour les dbits gazeux et 5% pour les
produits liquides. On prend = 2C pour les tempratures en-dessous de 100C et 5C au-dessus. Pour
la puissance au compresseur, on prend = 30kW. Lcart standard est pris gal 0.1 %mol pour les
fractions molaires en-dessous de 2%, 1 %mol au-dessus de 20% et 0.3 %mol autrement. Le problme de
rconciliation comprend 28 mesures, 33 variables dtat non mesures, 50 quations de contraintes et
donc 17 redondances. Dix variables sont considres constantes (toutes les pressions). Les mesures
semblent tre acceptables puisquaucune nest corrige par plus de deux fois lcart standard lui
correspondant, comme montr dans les rsultats (tableau 1).
La connaissance de la variance des variables valides nous permet de dtecter les importances
respectives de chacune des mesures pour le problme didentification de ltat dun systme. En particulier,
certaines mesures peuvent apparatre comme ayant peu dinfluence sur le rsultat, et peuvent donc tre
enleves de lanalyse. Inversement, dautres mesures peuvent avoir un impact bien plus important sur les
variables valides et sur leurs variances : ces mesures devront tre prises et manipules avec une attention
toute particulire et il peut savrer utile de multiplier les capteurs qui leur correspondent.
Le module danalyse de la sensibilit gnre deux types de rapport. Le premier type contient pour
chaque mesure la valeur mesure et la valeur valide, lcart standard suppos (Abs.Acc.) et celui estim a
- 5.22 -
Chapitre 5
posteriori (Rel.Acc.). Toutes les variables dtat dpendant dune mesure donne sont galement listes,
avec leur facteur de pondration (Contrib.) indiquant la contribution de la variance de la mesure la
variance de la valeur valide. Le second type de rapport contient la mme information mais trie par
variable dtat : pour chaque variable, on obtient la liste des mesures les plus importantes utilises pour
estimer sa valeur.
Lexemple suivant illustre le contenu du premier type de rapport. Les informations listes sont
relatives la mesure de la fraction molaire dammoniac dans le flux FL, identifi par son tag name
FL1_MFNH3. Cette valeur a t corrige de 97.800 %mol 98.005 %mol. La dviation standard de
cette valeur valide est 0.121 %mol (0.12% de la valeur valide), significativement plus faible que la
dviation standard sur la mesure qui est de 1 %mol (1.02% de la valeur exprimentale). Cette mesure
contribue faiblement (presque ngligemment) lestimation des valeurs valides de 2 autres variables dtat
mesures. Etonnamment, la variance de la fraction molaire dammoniac dans le flux FL est presque
indpendante de la mesure correspondante ( peine 1.48% de contribution).
Measurement
MFNH3
M FL
Tag Name
Reconciled
FL1_MFNH3
Variable
MFNH3
M FL
MFH2
M FL
Tag Name
FL1_MFNH3
FL1_MFH2
Value
98.005
97.800
Contrib.
1.48%
1.45%
Abs.Acc.
.12151
1.0000
Der.Val.
.14765E-01
-.97726E-02
Rel.Acc.
.12%
1.02%
P.U.
%
%
Der.Acc.
.60394E-02
-.39974E-02
P.U.
%
%
La colonne Contrib contient la contribution de la variance de la mesure k dans lestimation de la

variance de la variable dtat valide i, ce qui peut scrire :
var (Y ) var( y )
i
k
Contrib ki =
1 2
( M ) ik
Dans le tableau, on indique galement le coefficient de sensibilit reliant la valeur valide de la

variable la valeur mesure et la dviation standard.
dY
DerVal ki = k
dy i
dY
DerAcc ki = k
d i
Lexamen du rapport montre quaucune variable cl nest significativement influence par les
mesures de la composition du condensat, puisque la contribution de ces mesures la variance des valeurs
valides est plus petite que 10%. On peut en conclure que ces mesures de composition napportent que
peu dinformations. Il est vraisemblable que le cot des mesures ne se justifie pas, en comparaison de la
petite vrification permise par la rconciliation des donnes. Pour vrifier ceci, on peut se rfrer au
deuxime type de rapport de sensibilit, qui montre comment chaque variable dtat a t valide. Le
rapport contient des entres spares pour chaque variable dtat, comme illustr ci-dessous pour la
fraction molaire partielle de lazote dans le flux FL :
- 5.23 -
Chapitre 5
Variable
MFN2
M FL
Tag Name
Reconciled
FL1_MFN2
Measurement
MFH2
M FL
MFN2
M FL
T
S RECYC
T
S PURGE
MFN2
M FV
MFNH3
M FV
Tag Name
FL1_MFH2
FL1_MFN2
RECYC_T
PURGE_T
FV1_MFN2
FV1_MFNH3
Value
.52113
.50000
Contrib.
53.34%
10.35%
6.69%
6.69%
5.09%
4.19%
Abs.Acc.
Rel.Acc.
.32165E-01
6.17%
.10000
20.00%
Der.Val.
.23492
.10346
.41586E-02
.41586E-02
.72550E-02
-.21951E-01
Der.Acc.
-.10165
.43732E-01
.22633E-02
.22633E-02
-.27015E-02
.39223E-01
P.U.
%
%
P.U.
%
%
C
C
%
%
La premire ligne de ce tableau identifie la variable (fraction molaire de N2 dans le flux FL), le tag
name de la mesure correspondante, la valeur valide et sa dviation standard, enfin les units de la
variable. La valeur mesure et sa dviation standard sont indiques la ligne suivante pour comparaison :
lincertitude a t diminue dun facteur 3. Les lignes suivantes montrent quelles sont les mesures les plus
importantes utilises pour estimer la variable slectionne.
Aprs la suppression de 5 mesures de composition pour le flux FL, la redondance est rduite 12.
Lutilisation du programme de validation dmontre cependant que les variables valides ne sont pas
influence significativement par cette modification. Ceci peut galement tre illustr dune autre manire :
quand la validation est ralise aprs lajout de bruit aux mesures de composition du flux FL, ltat du
procd nest point affect et les variables perturbes sont correctement ramene leurs valeurs
normales de validation. Cet exemple nous permet de conclure que les cots dchantillonnages et
danalyse peuvent tre rduits en supprimant toutes les mesures inefficaces, tout en se concentrant sur
lamlioration des autres mesures.
Les variables non mesures sont galement dans le rapport de sensibilit, comme indiqu cidessous pour lavancement de raction :
Variable
EXTENT1 U REAC
Measurement
MASSF
M FL
T
S RCTIN
T
S RCTOUT
MFNH3
M RCTIN
MASSF
M RCTIN
S F07
Tag Name
Computed
Tag Name
FL1_MASSF
RCTIN_T
RCTOUT_T
RCTIN_MFNH3
Average
F06_MASSF
RCTIN_MASSF
F07_T
Value
.98972E-01
Contrib.
35.87%
14.91%
14.84%
7.86%
6.78%
50.00%
50.00%
5.24%
Abs.Acc.
Rel.Acc.
.41401E-02
4.18%
Der.Val.
.39674E-02
-.31975E-03
.31903E-03
-.38688E-02
.19795E-03
.98977E-04
.98977E-04
-.47368E-03
Der.Acc.
-.42549E-02
.30170E-03
.30034E-03
-.27949E-02
.16379E-03
.57909E-04
.57907E-04
.32725E-03
P.U.
kmol/s
P.U.
t/h
C
C
%
t/h
t/h
t/h
C
La premire ligne de ce tableau identifie la variable (paramtre EXTENT dans lunit REAC).
Puisque cette variable nest pas directement mesure, elle na pas de tag name. Il apparat que la plus
importante mesure est le dbit de condensat (qui est effectivement troitement li la production
dammoniac dans le racteur). Les tempratures dentre et de sortie, combine au dbit massique dans le
racteur , contribuent lestimation de ce paramtre (ils sont lis la conversion par le bilan nergtique).
Les autres variables jouent un rle moins vident. On peut remarquer que la mesure du dbit massique
- 5.24 -
Chapitre 5
dentre du racteur est en fait une moyenne pondre de deux valeurs diffrentes relies aux flux F06 et
RCTIN.
Lanalyse de sensibilit identifie correctement les variables qui sont lis par des contraintes, et sont
en fait des occurrences multiples de la mme variable dtat, comme T pour les flux FL, FV, PURGE, et
RECYC. Ceci est parfaitement illustr dans tableau de rsultats suivant, qui indique toutes les variables
influence par la mesure de la temprature de la purge :
Measurement
T
S PURGE
Tag Name
Reconciled
PURGE_T
Variable
T
T
T
T
T
MFNH3
Tag Name
RECYC_T
PURGE_T
Computed
Computed
FLV_T
FV1_MFNH3
S
S
S
S
S
M
RECYC
PURGE
FL
FV
FLV
FV
Value
293.69
20.000
Contrib.
24.38%
24.38%
24.38%
24.38%
9.31%
9.04%
Abs.Acc.
.98751
2.0000
Der.Val.
.24379
.24379
.24379
.24379
.10726
.33895E-01
Rel.Acc.
.34%
.68%
P.U.
K
C
Der.Acc.
.13268
.13268
.13268
.13268
.58378E-01
.18447E-01
P.U.
K
K
K
K
K
%
Un autre rsultat important obtenu par lanalyse de sensibilit est dindiquer exactement ce qui a t
gagn par la validation des donnes. Quand plusieurs variables sont lies par des contraintes plusieurs
mesures, la redondance peut tre exploite pour diminuer leurs incertitudes : lefficacit du processus de
validation est dmontre quand la dviation standard a posteriori est significativement infrieure la
dviation exprimentale (un facteur 2 dans lexemple prcdent).
Quelques autres variables sont limites par de fortes contraintes : leurs valeurs dpendent plus des
quations de contraintes que des valeurs directement mesures, qui donc peuvent tre limine sans perte
dinformation. Un exemple typique est une mesure de temprature pour un flux dont la pression est fixe,
ne contenant quun seul compos et tant en quilibre liquide - vapeur : la valeur valide de la temprature
sera fixe seulement par la contrainte de lquilibre liquide - vapeur, et ne sera pas influence par la valeur
mesure. Une observation similaire est faite propos de fraction molaire dammoniac dans le condensat
FL : la dviation standard a t significativement diminue par le processus de validation, mais la mesure
directe ne contribue que faiblement la valeur valide (voir le premier exemple ci-dessus).
Dun autre ct, certaines variables ne sont que peu influences par les mesures redondantes, et ne
peuvent donc tre corriges par la procdure de rconciliation des donnes. Leurs valeurs valides
dpendent principalement (ou mme seulement) de leurs propre valeur mesure : une double vrification
des mesures et une calibration minutieuse des capteurs est alors recommande. La rconciliation des
donnes peut faire des merveilles, mais pas de miracle, et quelques mesures doivent absolument tre
valable : lanalyse de sensibilit permet de les identifier. Un exemple est la puissance du compresseur, dont
lincertitude nest point diminue par la validation lorsquon utilise le groupe de mesures dont on se sert
depuis le dbut :
- 5.25 -
Chapitre 5
Measurement
POWER
S W3
Tag Name
Reconciled
W3_POWER
Variable
POWER
EFFIC
T
T
Tag Name
W3_POWER
Computed
F07_T
F06_T
S
S
S
S
W3
EFF3
F07
F06
Value
332.28
330.00
Contrib.
96.82%
87.31%
6.06%
1.44%
Abs.Acc.
29.519
30.000
Der.Val.
.96817
-.19034E-02
.96290E-02
-.40619E-02
Rel.Acc.
8.88%
9.09%
P.U.
kW
kW
Der.Acc.
.14739
-.28977E-03
.14658E-02
-.61835E-03
P.U.
kW
K
K
Comme conclusion, on peut dduire des rgles qui permettent didentifier les bonnes mesures : elles
ne devront pas tre trop modifie par la rconciliation des donnes (petite erreur de mesure) mais leur
dviation standard associ devrait diminuer grce au processus de validation (existence de redondance
amliorant la connaissance de la variable).
La connaissance de la dviation a posteriori pour toutes les variables du procd nous permet de
calculer lintervalle de confiance ou bien les marges de scurit des rsultats obtenus (lcart la limite
dexplosivit, la courbe de refoulement du compresseur). Un ensemble bien choisi de mesures
accompagn dune procdure de rconciliation des donnes permet de diminuer les marges de scurit en
ce qui concerne les conditions critiques de fonctionnement, puisque loutil de rconciliation des donnes
rduit lincertitude sur ltat actuel du procd.
5.3.4 Conclusions
Un des but de la validation est damliorer la connaissance des variables dtat du systme. Fournir
des valeurs est bien sr dune grande aide mais valuer leur fiabilit est tout aussi important, et cest dans
cette optique que des dviations standards pour les variables valides et pour celles non mesures ont t
dveloppe.
Trois types de questions peuvent tre analyses par lutilisation de ltude de la sensibilit :
Premirement, vrifier comment la prcision dun variable dtat donne est influence par le groupe de
mesures : Quelles sont les mesures qui contribuent significativement la variance du rsultat
valid pour un groupe de variables dtat ?
Le second type de problme est de dterminer les variables dtat dont la prcision est le plus influenc
par une mesure donne : Quelles sont les variables dtat dont la variance est influence
significativement par la prcision dune mesure donne ?
Le troisime type de problme est dtudier comment la valeur dune variable dtat est influence par
la valeur et la dviation standard de lensemble des mesures.
A partir de ces informations, des dcisions peuvent tre prises soit en vue de lanalyse des mesures
dun procd existant, soit pendant le design dun systme de mesure. Les analyses inutiles peuvent tre
limines, ou bien effectues moins frquemment, ce qui permet ainsi une diminution intressante des cots
opratoires. On peut galement identifier les mesures cls pour lesquels une amlioration de la prcision
permettrait un meilleur suivi du procd. On peut galement dterminer lemplacement des capteurs
- 5.26 -
Chapitre 5
permettant une bonne estimation de toutes les variables cls du procd au cot dinvestissement le plus
bas.
- 5.27 -
Chapitre 6
CHAPITRE 6
APPLICATION DES METHODES DOPTIMALISATION
AVEC CONTRAINTES AUX CALCULS DEQUILIBRES
CHIMIQUES
Introduction
Le problme choisi comme application des mthodes doptimalisation vues prcdemment, est
celui de la minimalisation dune fonction dnergie libre. Nous en profiterons pour examiner brivement de
quelle manire on peut tenir compte de contraintes.
6.1 Position du problme

Lnergie libre dun mlange gazeux est fonction de la temprature T, de la pression P et de la
composition
G = G( P, T, x)
o x est le vecteur des nombres de moles des constituants. Le nombre total de moles sera not (non
soulign) :
n
x = xi
i =1
Le potentiel chimique de chaque constituant est dfini par :

G
i =
x i P ,T
o G est lnergie libre de Gibbs (ou enthalpie libre). On sait que :

n
G = x i. i = t x
i= 1
Pour un gaz parfait, les potentiels chimiques sont calculs par les relations :
i = 0i ( T, P) + RT. ln
xi
x
0i ( T, P) = *i ( T, P 0 ) + RT.ln
- 6.1 -
P
0
P
Chapitre 6
Les potentiels chimiques standard * sont disponibles dans les tables spcifiques ou peuvent tre
obtenus en utilisant les banques de donnes : ils sont dfinis pour les gaz parfaits et sous une pression de
rfrence P 0 = 1 atm . Le problme est le suivant : calculer x T et P fixs ?
Posons :
(T, P )
Ci = i
RT
0
qui sont des constantes calculables. On a alors :

xi
G ( x) = x i . 0i ( T, P) + RT.ln
x
i =1
n
ou sous la forme adimensionnelle :

G( x )
f ( x) =
=
RT
x . (C
n
i =1
+ ln x i ln x )
La composition dquilibre est celle qui minimise f(x) en respectant la loi de conservation des
lments et le fait que les nombres de moles ne peuvent pas tre ngatifs. En termes formels :
Minimaliser f(x) par rapport x
avec x tel que g k ( x) > 0 pour k = 1, K

h j ( x ) = 0 pour j = 1, J
g sont des contraintes dingalits,

h sont des contraintes dgalit,
dans notre cas :

Min f ( x) avec x i > 0 et A x b = 0
Les bilans des lments sont crits sous forme matricielle. Un exemple concret lucidera cette
notation compacte.
Soit le cas du reforming du mthane la vapeur, le systme comporte les substances :
x = CH 4 H 2 CO H 2O CO 2
A x b = 0 reprsente les trois quations suivantes :
t
1 0 1 0 1
4 2 0 2 0
0 0 1 1 2
x1
x
2 b1
x3 = b2

x4 b 3
x 5
- 6.2 -
et
b= C H O ,
Chapitre 6
les termes de la matrice A sont les nombres datomes des lments j dans les molcules i; sil y a n
molcules et m atomes, la dimension de la matrice A est m n, les termes du vecteur b sont les nombres
datomes-grammes des lments prsents dans les systmes (au nombre de m qui est la longueur dune
colonne de la matrice A ).
Les donnes des problmes sont la temprature, la pression, les potentiels chimiques de rfrence,
les molcules prsentes et les nombres datomes-grammes des lments.
6.2 Solution
On applique la mthode de Lagrange en tenant compte des contraintes dgalit. On tient compte
des contraintes dingalit au cours de calculs itratifs en rendant positives et ngligeables les valeurs des
composantes de x qui seraient devenues ngatives; ce point sera explicit ci-aprs.
La fonction de Lagrange pour le problme pos scrit :
L( x , ) = f ( x) + t A x b
La diffrentielle totale de cette fonction est nulle loptimum :
dL = f ( x) + t A x + t A x b = 0
quelles que soient les variations x et , donc les conditions dquilibre sont :
f ( x) + t A = 0
Ax b = 0
la premire expression reprsente un systme de n quations n+m variables, les quations tant non
linaires; tandis quavec la seconde, on retrouve bien sr les contraintes, qui constituent un systme de m
quations linaires n variables. Au total, on a autant dquations que dinconnues et si elles sont
indpendantes, il existe au moins une solution.
Dun point de vue pratique, le gradient de la fonction f est calculable car on a tabli son expression
ci-dessus. La composante i scrit :
f
i
=
x i RT
puisque :
f ( x) =
G( x )
=
RT
x . (C
n
i =1
+ ln x i ln x )
- 6.3 -
Chapitre 6
Un systme dquations non linaires se rsous par itrations successives et ncessite quelques (n3)
oprations par itration. On a donc intrt diminuer, si possible, le nombre total n des quations
considrer simultanment.
6.2.1 Mthode des constantes dquilibre
Cette mthode est relativement bien connue des chimistes. Voyons rapidement comment on peut la
retrouver. Nous allons liminer les j coefficients de Lagrange et faire apparatre les avancements des r
"ractions" chimiques indpendantes (r = 1, R). Notons au passage quil peut paratre surprenant de faire
apparatre le terme ractions alors quon recherche le minimum dune fonction dtat : le chemin suivi
(impos par les ractions chimiques) ne doit pas avoir de rpercussion sur le rsultat. Les avancements de
raction (exprims dans les mmes units que les composantes de x ) sont dfinis par lquation :
x = x0 +t V 0
normalement, les degrs davancement 0 sont nuls, la matrice V contient comme lments, les
coefficients stoechiomtriques des ractions chimiques indpendantes au nombre de R. Cette relation
implique que :
x = t V
Explicitons nouveau cette notation compacte par lexemple du reforming du mthane : les
ractions dquilibre communment crites sont :
CH4 + H2O CO + 3 H2
CO + H2O CO2 + H2
qui peuvent aussi scrire :
- CH4 - H2O + CO + 3 H2 = 0
- CO - H2O + CO2 + H2 = 0
avec les bilans classiques :
[H ] = [H ]
+ 3. 1 + 2
[ CO] = [ CO] 0 + 1 2
etc.
La matrice V (nR) est dans ce cas :
1
V=
0
3
1
1 1
0
1
On pourrait crire dautres ractions et obtenir une autre matrice de coefficients. On dmontre
que (bilans sur les lments) :
- 6.4 -
Chapitre 6
V t A=t A V = 0
ici, 0 est une matrice mR et R est le nombre de ractions indpendantes dailleurs gal n (nombre de
molcules) - m (nombre datomes).
Utilisons les avancements de ractions pour liminer les coefficients de Lagrange, pour ce faire,
remplaons le vecteur x par son expression en fonction des degrs davancement dans la diffrentielle de
L:
dL= t V f ( x ) + t A + t A x b = 0
le deuxime terme entre crochets tombe vu que V A = 0 .
t
Les conditions dquilibre deviennent :

V f(x ) = 0
R = n-m quations
Ax b = 0
m quations
La dimension du problme (N) a t rduite n quations. Retrouvons maintenant les expressions

classiques des constantes dquilibres, elles ne sont rien dautre quune criture diffrente des quations :
V f ( x) = V
RT
quon explicite au maximum :

n
i =1
ri .
n
i
xi
= ri . C i + ln
RT i=1
x
avec r = 1,R
ensuite,
x
ri . C i = ln i
x
i =1
i =1
n
ri
x
= ln i
x
i =1
n
r i
ou encore
x
K = i
i =1 x
n
r
i
ri
= exp ri . C i
i =1
cette relation peut aussi scrire en termes de pressions partielles, car on sait que, si on admet la loi de
Dalton :
pi =
xi
.P
x
on obtient alors le Kp.

Il y a bien sr autant de constantes dquilibre que de ractions chimiques indpendantes.
Remarquons quon recherche souvent la solution de ces R quations explicites en fonction des
- 6.5 -
Chapitre 6
avancements de ractions. Telles quelles sont crites ci-dessus, on doit leur adjoindre les contraintes de
bilan et rsoudre ce systme directement en x.
6.2.2 Mthode de WHITE
Nous venons de voir que les n+m quations rsoudre dans le cas gnral sont :
f ( x) + t A = 0
Ax b = 0
La mthode de WHITE est base sur un examen approfondi de la structure de ces quations
linarises qui permet de rduire trs significativement la dimension du systme dquations irrductible.
Les quations (non) linaires peuvent tre rsolues par la mthode de Newton-Raphson.
Rappelons-en rapidement lessentiel.
Dans le cas dune fonction une variable f ( x) = 0 , la solution sobtient par :
k +1
=x
k
f (xk )
f ( x k )
Dans le cas dune fonction vecteur de dimension N, N variables, formellement f ( x) , la solution

sobtient par :
x k +1 = x k J x1k f ( x k )
o J est la matrice Jacobienne de la fonction vecteur, prise au point de la kme itration.
En appliquant cette dernire relation au systme des n+m quations rsoudre, on doit donc
prendre le Jacobien, de sorte quon obtient (on prend la formule sans linversion) :
H k
x
A
A x k +1 x k
=
0 k +1 k
n n|n m
m n| m m
f ( x k ) + t A k
A xk b
n 1
m1
n 1
m1
la matrice H est le Hessien de la fonction f (en effet le Hessien de L se rduit celui de f), pris au point de
la kme itration, il est issu du Jacobien du gradient de f. En ce qui concerne la fonction f que nous avons
dveloppe, le terme ij (au point de la kme itration) scrit :
H ij =
k
ij 1
k
k
xi x
avec i = 1,n et j = 1,n
o ij est la fonction de Kronecker et vaut 1 si i = j, zro sinon (si i j).
- 6.6 -
Chapitre 6
En explicitant tout fait lcriture des quations, on obtient :

1 1
x x
1
1
x
L
1 1
xi x
L
1
1
x
1
x
a ji
1 1
xn x
a ij
m
1k
+
a j1 . kj
RT j=1
x 1k+ 1 x 1k
k L
m
L
i
k
+
a
.
k+ 1
k
ji
j
x i x i
RT j=1
=
k L
L
m
n
k
x k+ 1 x k
n
RT + a jn . j
nk+ 1
k
j =1
m
j j
k
a ij . x j b j

j=1
Les n premires quations permettent dexpliciter x ki +1 en fonction de kj+1 . En effet, prenons la

ime :
m
m
ij
1 k+ 1
ki
k
k +1
k
k k .( x i x i ) + a ji . j j =
a ji . kj
x
RT j=1
i =1 x i
j =1
qui se rduit :
m
ij
1 k+1
ki
k
k+ 1
k k . ( x i x i ) + a ji . j =
x
RT
i =1 x i
j =1
n
Ceci rend facile llimination des n valeurs de x de la (k+1)me itration laide des n quations cidessus, hors des m dernires quations qui ne contiennent plus alors que des termes en de litration
(k+1). Ces m dernires sont donc rsolubles en fonction de k . En ralit, WHITE, pour des facilits
dcriture des quations, en ajoute une au systme, savoir :
x k +1 =
relation dans laquelle les x
k +1
i
x
i =1
k +1
i
sont aussi limins en fonction des kj+1 .
Il suffit ds lors de rsoudre simultanment (m+1) quations (m+1) inconnues (les m et le x cidessus) et ensuite de calculer les n valeurs de x k +1 tour de rle de manire explicite.
Lalgorithme est alors le suivant :
1. Choix de x 0 tel que A x 0 b = 0
(par exemple par tirage au hasard contraint),
2. Solution des (m+1) quations simultanes pour trouver kj+1 et x k +1
3. Tests de convergence pour les valeurs prcdentes,
4. Calcul des x ki +1 .
- 6.7 -
Chapitre 6
Le dveloppement thorique na t donn que pour des gaz idaux, mais on le gnralise
facilement au cas des gaz rels, lexpression des potentiels chimiques se compliquant videmment.
6.2.3 Contraintes dingalit et points particuliers
Dans chacune des mthodes exposes, on calcule les valeurs (k+1) des variables, partir de leurs
valeurs obtenues litration prcdente (k). On connat en fait le vecteur x k +1 x k . Avant daccepter les
valeurs x k +1 , on vrifie quelles sont toutes positives. Dans le cas contraire, on garde la direction du
vecteur, mais on en diminue la grandeur de telle manire que les x k +1 ngatifs prennent une valeur
ngligeable mais positive : il sagit dune mthode de relaxation classique. Ce test nintervient dailleurs
quaprs avoir test si le gradient de la fonction objectif le long du vecteur incrment na pas pris une valeur
positive excessive, ce qui indique que lon doit ici aussi utiliser un facteur de relaxation.
En toute gnralit, on a :
x k +1 = x k + s k
le vecteur s k est connu en direction, sens et grandeur.
est le facteur de relaxation qui est dtermin par une mthode heuristique. Notons encore qu
chaque itration, on opre une correction des erreurs darrondi; ceci se fait en ne modifiant pas les (n-m)
plus grandes valeurs des composantes de x k +1 et en corrigeant les m valeurs restantes de manire
satisfaire correctement les contraintes de bilan.
Signalons enfin que le programme relatif la mthode de White est accessible en mode
conversationnel, tant pour les gaz idaux que pour les gaz rels (programme THERMO2).
6.2.4 Exemple dapplication
Combustion de lhydrazyne : N2H4 + O2
Composition du systme de dpart :
1 atgr de N
1 atgr de H
1 atgr de O
on a le vecteur lment :
[H
N O
Composants existant 3500 K et 50 bar :

H, H2, H2O, N, N2, NH, NO, O, OH, O2
Ci-dessus, les lments sont mis dans un ordre admis pour la suite.
Donc n = 10 et m = 3.
Puisquil y a 3 lments, on a 3 quations de bilan.
- 6.8 -
Chapitre 6
On peut aussi crire 10 - 3 = 7 quations (ractions) indpendantes.

Par exemple :
H2 2 H
O2 2 O
O + H OH
N2 2 N
OH + H H2O
NO N + O
2 NH H2 + N2
On vrifiera aisment que ces ractions sont indpendantes.
Dans ces conditions, les matrices A et V scrivent :
1 2 2 0 0 1 0 0 0 1
A = 0 0 0 1 2 1 1 0 0 0
0 0 1 0 0 0 1 1 2 1
2 1
0
0
1 0
V= 0
0
0
0
1 0
1
0
2
(bilan atomique)
0
0
0
2 1 0
0
0
0 1 0
1
(coefficients
1 0
0
0
0
0 stoechiomtriques)
0
0 1 1
0
0
0
0
0
0
0 1
1
0
0
0
2 1 0
0
les matrices respectent lordre adopt ci-avant.

Recherche du point de dpart x 0 :
on fixe 7 variables x i alatoirement, de telle sorte que :
xi
0<
< 1 et
x
x
i =1
<1
Les trois quations de bilan A x b = 0 fournissent les valeurs des 3 dernires variables.
Rajustement de x (en relation avec la remarque du paragraphe prcdent).
On choisit de toujours corriger les trois composantes de x dont les valeurs sont les plus grandes.
- 6.9 -
Chapitre 6
6.3 Utilisation de diffrentes mthodes

Dans ce paragraphe, nous allons examiner assez brivement les caractristiques et efficacits des
diffrentes mthodes itratives vues au chapitre 3, en travaillant sur le problme expos prcdemment. La
formule itrative est ventuellement ramene lesprit.
Le problme gnral quil sagit de rsoudre, est :
f ( x) + t A = 0
Ax b = 0
Ce problme consiste en la recherche dun point stationnaire de la fonction Lagrangienne L.
On mentionne des temps de calcul relatifs, puisque ce sont ceux obtenus pour un mme ordinateur.
6.3.1 Steepest Descent
La recherche dun meilleur point se passe en 2 tapes :
on calcule dabord la direction du gradient quon prend comme direction dinvestigation, puis sachant que,
x k +1 = x k + L( x k )
o L est donn par :
L( x , ) = f ( x) + t A x b
et x comprend toutes les variables (les n x et les m ), il reste rechercher la distance parcourir sur
cette direction (du gradient) pour arriver un optimum.
La rsolution dun problme 13 variables a t arrte en raison de la lenteur de la convergence..
6.3.2 Newton
Cette mthode ncessite la connaissance du gradient et du Hessien de la fonction Lagrangienne. En
effet, les quations sont :
H x k h k = f ( x k )
tape (a)
x k +1 = x k + h k
tape (b)
f est la fonction optimiser (ici L), x le vecteur variable, (ici x et ).

On calcule donc :
f ( x )+ t A k
L( x , ) =
A x b
- 6.10 -
Chapitre 6
[ ]
H f ( x) t A
H L( x , ) =
0
A
et il reste alors calculer le gradient et le Hessien de la fonction f ( x)
Le problme 13 variables a t rsolu par cette mthode. Le temps de calcul est choisi comme
rfrence et nous lui affectons la valeur 100.
6.3.3 Davidon, Fletcher et Powell
Cette mthode fait appel un formalisme plus compliqu que nous ne reproduirons pas, cependant
elle offre lavantage de pas ncessiter les drives secondes, puisquelle gnre le Hessien au fur et
mesure des itrations. Le programme de loptimalisation est nettement plus lent et la rsolution a t arrte
avant lobtention dune solution converge.
6.3.4 Newton-Raphson
Rsolution complte temps quivalent = 50 (10 quations). Le gain par rapport la mthode de
Newton est d principalement au temps dinversion de la matrice. En effet, Newton exige linversion dune
matrice de 1313, tandis que Newton-Raphson exige linversion dune matrice de 1010, chaque pas (il
en faut une dizaine dans le cas trait).
6.3.5 Sequential Unconstrained Minimization Techniques (SUMT)
Cette mthode nayant pas encore t expose, nous le faisons ici.
Un problme tel quil est pos au dbut du paragraphe 6.3 peut tre rsolu en effectuant
successivement une suite doptimalisations sans contraintes. Fiacco et Mac Cornick proposent dadopter
la fonction objectif suivante minimaliser :
1 L
F( x , r ) = f ( x ) + r . g k ( x) + . h l ( x)
r l=1
k =1
K
La fonction f ( x) a t simultanment pnalise par les contraintes dgalit et dingalit.

Ayant construit la fonction F, on adopte une valeur de r, et on admet quon est litration k. On
optimalise F( x , r ) selon une des mthodes cites ci-dessus. Ensuite, on adopte une loi de dcroissance
pour les valeurs de r k , par exemple :
r k +1 =
rk
; > 1
et on rpte pour chaque valeur de r, la procdure de minimalisation choisie, accordant ainsi une
importance croissante aux contraintes dgalit.
- 6.11 -
Chapitre 6
La valeur finale de r est fonction de la prcision requise concernant ces contraintes.

La mthode permet ainsi de cheminer vers loptimum en ne respectant les contraintes quen finale,
comme on peut le voir sur la figure 6.1 :
F
Contraintes d'galits
x
0
Figure 6.1
Les points successifs sur la trajectoire doptimalisation correspondent aux diffrentes valeurs de r.
Dans le cas qui nous occupe, on a lexpression :
1
1 m t
F( x , r ) = f ( x) + r . k + k . a j x b j
r j=1
k= 1 x
k
La minimalisation de cette fonction dans un problme 13 variables par la mthode de Newton est
mene bien en un temps quivalent de 330. On a = 4 et r 0 = 0,1 ; le point de dpart ayant t pris au
hasard.
6.3.6 "Complex"
Rsolution complte en 220 (temps quivalent) (13 points). Cette mthode nutilise pas de
formalisme mathmatique, elle se contente des valeurs de la fonction optimiser. Elle est robuste mais
relativement lente. Cependant, on fait remarquer que la mthode est particulirement indique quand les
contraintes sont linaires.
6.3.7 Programmation linaire
Cette mthode est tellement gnralise que nous ne traiterons pas de la mthode comme telle,
mais bien de la transformation du problme pos, de manire le rendre apte tre rsolu par cette
mthode. Des programmes spcialiss sont disponibles et dun emploi facile.
- 6.12 -
Chapitre 6
La mthode sapplique uniquement des problmes linaires. Elle exige donc la linarisation des
contraintes et de la fonction objectif.
Linariser une fonction y = f ( x ) consiste rendre son expression linaire en toutes les n variables
du vecteur x .
Admettons des points quelconques de lhypersurface (n+1) dimensions (y et x), ces points sont
reprs par lindice suprieur k, on suppose quil y en a m. On approximera cette hypersurface par :
m
x i = n k . x ki
pour i = 1,n
k =1
y = nk . yk
(*)
k =1
o les n k sont les nouvelles variables sur lesquelles on a les contraintes :

m
n
k =1
= 1 et n k > 0
les quations plus haut constituent un systme de (n+2) quations m inconnues, (n+1) au maximum des
variables n sont nulles et elles sont ncessairement conscutives. On en dduit que m > n + 1, il faut plus
dinconnues que dquations.
Cette mthode de linarisation porte le nom de linarisation par tronons.
Exemple : Soit y = f ( x1 , x 2 ) . Cette surface 3 dimensions est reprsente sur la figure 6.2 :
y
x2
x1
Figure 6.2
La linarisation de la fonction consiste approximer la surface par des plans (ici par des triangles,
au minimum 3 variables n non nulles).
- 6.13 -
Chapitre 6
Approximons par exemple, au moyen de 4 triangles (6 points do 6 variables n) et nous obtenons

la triangulation de la figure ci-dessus (en traits pleins). Une approximation en 3 triangles (4 points do 4
variables n) donnerait la triangulation en pointill.
Application la fonction type :
n
i =1
i =1
f ( x) = x i . Ci + ln x i ln x = C i . x i + x i .ln
i =1
et on pose en premier lieu : i = x i .ln
xi
x
xi
x
Linarisons les i suivant des i (en effet nest pas du tout linaire en fonction des variables x).
Ce problme tant 2 dimensions, il exige au moins 2 variables n k ; prenons-en 3, ce qui correspond
linariser en 2 tronons comme le montre la figure 6.3 . On a ainsi :
f ( x) = Ci . x i + i = Ci . x i + ( n 1i .1i + n 2i . 2i + n 3i . 3i )
n
i =1
i =1
i =1
i =1
le dveloppement de la fonction (il y en a i), est ralis selon la formule correspondant y (se rfrer aux
3
quations * page prcdente), et on a en plus : i = n ik . ki ,

k =1
do :
3
k =1
k =1
x i = x . i = x . n ki . ik = n ki . x ik
les contraintes tant :

3
n
k =1
k
i
= 1 pour tout i = 1,n
Deux n ki au maximum (sur 3) sont nuls, n ki > 0

Choisissons les x ki tels que :
x ki
= n ki
x
Dans ces conditions, on a :
x ik
f ( x) = Ci . x i + .
i =1
i =1 k =1 x
n
x i = x ki . ki
k =1
3
x + x ki = 0
pour tout i
k =1
- 6.14 -
Chapitre 6
i1
i2
i3
Figure 6.3
et le problme revient minimaliser f ( x) suivant :

x i , x , x ik
Il comprend alors (4n+1) variables, cest--dire n pour x i , 1 pour x, 3n pour x ki . De plus, il est
soumis aux contraintes :
Ax b = 0
m quations de bilans atomiques
x + x ki = 0
n quations
k =1
3
x i = x ki . ki
n quations
k =1
x xi = 0
1 quation
i =1
toutes les variables > 0
4n + 1 quations
ce qui fait en tout (6n+m+2) quations.

Dans le cas du problme nonc prcdemment, cela fait 41 variables et 65 contraintes (41
contraintes dingalit et 24 contraintes dgalit).
La mthode est itrative, chaque pas, on diminue le domaine dinvestigation en rduisant
lintervalle de linarisation.
Le nombre lev de variables explique que la rsolution complte exige un temps quivalent de
500.
6.3.8 Quelques conclusions

Il est probablement utile de souligner ici lenseignement que lon peut tirer des rsultats exposs.
De nombreuses techniques ont t testes et on peut donc juger de leur efficacit relative. Un rsultat
intressant peut tre obtenu grce une analyse pousse de la formulation du problme, mais le travail
- 6.15 -
Chapitre 6
investi ne devient rentable que si la mthode doit tre employe un grand nombre de fois. Dans le cas
contraire, on peut se contenter dune approche mathmatique plus sommaire.
On peut galement souligner lefficacit des mthodes du second ordre, de mme que lavantage
obtenu par la rduction du rang de la matrice de travail, cest--dire le nombre de variables considres
simultanment dans les calculs.
- 6.16 -

Simulation Et Optimisation Des Systèmes Industriels

Transféré par

Droits d'auteur :

Formats disponibles

Simulation Et Optimisation Des Systèmes Industriels

Transféré par

Informations du document

Description originale:

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

Simulation Et Optimisation Des Systèmes Industriels

Transféré par

Droits d'auteur :

Formats disponibles

Universit de Lige

Facult des Sciences Appliques

Succession Prof. B. Kalitventzeff

Laboratoire d'Analyse et de Synthse des Systmes Chimiques

La figure ci-dessus reprsente un systme compos de 6 units interconnectes entre elles

Figure 1.3 Flow-sheet

Figure 1.4 Graphe dual

COUPER LE FLUX INTERVENANT DANS

Grce au 2me tableau, le graphe dual devient :

Grce au 3me tableau, le graphe dual devient :

Grce au 4me tableau, le graphe dual devient :

Le dernier tableau donne le graphe dual :

Application des rgles expliques dans ce chapitre :

Les entres sont 1 et 21; les sorties sont 22, 6, 7.

coupures massiques (seul le dbit total est inconnu),

Application : la synthse de lammoniac

Figure 1.20 Schma gnral des reforming primaire et secondaire

Afin dillustrer la mthode de Motard, nous considrons uniquement la section reforming du

Le tableau ci-dessus montre lapplication de la mthode de Motard la section reforming du

2.1 Mthode gnrale des bilans

2.2 Modles mathmatiques de diffrents appareils

Hypothse : il ny a pas de changement de phase dans le mlangeur.

F.xi = F1.xi1 + F2.xi2

normation (fractions molaires)

Nombre de degr de libert = (n+3) variables - (n+3) quations = 0

tempratures T1,T2 (K)

F.xi = F1.xi1 + F2.xi2

normation (fractions molaires)

Nombre de degr de libert = (2n+6) variables - (2n+5) quations = 1

Hypothses : il ny a pas de changement de phase dans la vanne et la dtente se fait de manire

bilan de matire (dbits molaires partiels)

Nombre de degr de libert = (n+3) variables - (n+2) quations = 1

Hypothse : lalimentation est compltement fixe.

F.zi = L.xi + V.y i

Nombre de degr de libert = (2n+5) variables - (2n+3) quations = 2

Hypothses : il ny a pas de condensation et lalimentation est compltement fixe.

dimensionnement (efficacit isentropique)

Nombre de degr de libert = (n+6) variables - (n+5) quations = 1

Hypothses : il ny a pas de condensation et lalimentation est compltement fixe.

Nombre de degr de libert = (n+6) variables - (n+5) quations = 1

2.3 Modles mathmatiques des racteurs

si r est la vitesse de raction en mole/m3.s, on peut crire :

car le racteur est isotherme,

si on nglige les pertes de charge,

Remarques concernant les bilans matires (quations 1 et 2) :

et le problme est alors trs simple ;

la premire cintique restant :

2.3.1.3 Systme non adiabatique

le coefficient de transfert de chaleur (kJ/m2.s.K),

la surface dchange (m2),

la temprature du fluide de la double enveloppe (K).

la temprature de rfrence (souvent celle du gaz parfait 25C et 1 atm) ;

qui, aprs transformations, devient :

r est la somme algbrique des vitesses de raction rj , r = r j ;

nous obtenons les quations suivantes comme bilan de matire du constituant i :

J est un coefficient de perte de charge,

2.4 Dcomposition dun modle mathmatique

Construisons la matrice dincidence (ou doccurence) du systme :

lquation Z3 Z2 + ( A B B2 ). Z A. B = 0 . Avec ces variables, on peut calculer le coefficient de

f ( x ) = 0 (cest--dire, un point stationnaire existe en x)