Acknowledgement

I would like first to thank my supervisors Mr. Didier Dumur, Mr. Pedro Rodriguez-Ayerbe, Mr. Fayçal Lawayeb, and Mr. Alain Mouchette who have guided and encouraged me since the beginning of my doctoral program. I am grateful for their help, their patience, their availability to correct my last minute papers and presentations.

I am very thankful to Mr. Edouard Laroche and Mr. Pierre Riedinger for reviewing my thesis and to Ms. Dorothée Normand-Cyrot, Ms. Estelle Courtial and Mr. Mohammed M'saad for being part of the examination committee of my thesis.

I also would like to thank all professors and personnel of the Automatic Control Department of CentraleSupélec, especially to Ms. Josiane Dartron, Ms. Pascale Lepeltier and Mr. Léon Marquet for their help during my three years at CentraleSupélec. I will not forget my colleagues who supported me with their pleasant presence during coffee breaks and lunch time, specially: Dung, Tri, Ngoc Anh, Binh, Djawad, Imad, Séif and Quang.

It is also a pleasure for me to thank the personnel of the Measurement and Control Department of ArcelorMittal R&D Center, especially to Ms. Laurette Duarte Domingues for her administrative support during my stays at the company.

I am particularly thankful to the colleagues of automatic control team from whom I have learned a lot about technical issues, about organizational and professional skills.

Finally, I am grateful to my girlfriend, my parents, my sister and my friends for their support and their encouragements.

Introduction

La sidérurgie a en particulier pour but de transformer le minerai et la ferraille de fer en acier. Pour cela, le fer doit subir plusieurs processus thermiques, mécaniques et thermomécaniques avant de devenir adaptés à d'autres applications industrielles telles que la construction, l'automobile ou l'industrie spatiale. Ces processus peuvent être regroupés en quatre étapes : la fabrication du fer, la fabrication de l'acier, la coulée continue et le laminage à chaud. Au cours de la fabrication du fer, le liquide de fer est produit à partir du minerai et de la ferraille. A ce liquide est ensuite ajouté différents éléments comme le nickel, ou le chrome pour fabriquer l'acier liquide. La coulée continue permet ensuite d'obtenir à partir de l'acier liquide des brames d'acier qui seront transformées selon les besoins du client. Enfin, la fonction principale du processus de laminage à chaud est de réchauffer les brames d'acier à une température proche de leur point de fusion, puis de les laminer par passage dans plusieurs laminoirs successifs et enfin d'enrouler les tôles d'acier sous forme de bobines, voir Fig. 1.

Fig. 1 : Laminage à chaud Les fours de réchauffage sont positionnés juste avant les laminoirs à chaud, et utilisés pour le réchauffage, ainsi que pour le traitement métallurgique des brames. Les brames d'acier sont habituellement chauffées de la température ambiante jusqu'à une température visée de 1200°C. La performance d'un four de réchauffage est évaluée par la précision de la température de défournement, par l'homogénéité de la température à l'intérieur des brames et par le rendement énergétique du four. Certains phénomènes tels que la formation de calamine, la décarburisation, la fonte de surface ou la distorsion géométrique des brames doivent être réduits autant que possible.

Dans une usine sidérurgique, les fours de réchauffages sont classés deuxième consommateurs d'énergie avec une consommation représentant 17% de l'énergie utilisée dans la production de l'acier. La consommation de gaz du four considéré dans l'étude est équivalente à la consommation d'électricité de 100,000 foyers. De fait, l'efficacité énergétique du four de réchauffage a une influence directe sur le coût de production. Par ailleurs, une réduction de la consommation de gaz entraîne une diminution des émissions de gaz à effet de serre ainsi que du CO 2 , N 2 O, etc, émissions objet de régulations de plus en plus restrictives. Plusieurs facteurs de conception affectent l'efficacité énergétique du four, tels que la géométrie du four, le système d'isolation thermique, et en particulier, la régulation de la température dans le four.

De nos jours, la diversification des produits fait intervenir une grande variété de brames en termes de matériaux, de dimensions, et de températures visées. Par conséquent, une mauvaise régulation de la température du four oblige les opérateurs à imposer de grandes marges de température afin d'assurer un fonctionnement satisfaisant du four. Cela induit une augmentation importante de la consommation en termes d'énergie.

Pour toutes ces raisons, l'objectif des travaux de cette thèse est d'améliorer la commande de la température du four afin d'en augmenter son efficacité énergétique et la qualité de réchauffage des produits.

Pour réaliser la tâche de réchauffage, les fours sont habituellement construits en plusieurs zones équipées de bruleurs à gaz comme indiqué Fig. 2 (typiquement récupération, préchauffage, chauffage, égalisation). Les longerons permettent de déplacer et donc de transporter les brames à travers le four. Le mur réfractaire est l'élément d'isolation thermique de four. La fumée issue de la combustion des gaz est évacuée à l'entrée du four, cela permet de récupérer une partie de la chaleur afin de chauffer les brames situées dans la zone de récupération. Le chauffage est réalisé dans les zones de préchauffage et de chauffage principal. La température à l'intérieur des brames est finalement homogénéisée dans la zone dite d'égalisation. Fig. 2 : Vue latérale du four de réchauffage Le but de la régulation de la température des fours de réchauffage peut être décrit de façon simplifiée comme étant la stratégie à développer pour assurer la meilleure distribution de l'énergie thermique le long du four afin de chauffer les brames à la bonne température de défournement. Les variables disponibles pour effectuer cette régulation de la température sont les débits de gaz et d'air des bruleurs et, plus indirectement, les consignes de température de zones. Dans certains cas, il est également possible d'ajuster la cadence des brames dans le four à condition que cela ne dépasse pas par la suite la capacité des laminoirs. Pour gérer la complexité du système, la structure actuelle de commande du four exploite l'architecture cascade hiérarchisée représentée Fig. 3. Au niveau 2, la température visée de brames est contrôlée à l'aide d'un modèle de transfert thermique du four, à ce niveau, les consignes de température de zones sont calculées et envoyées ensuite à niveau 1. A partir des consignes calculées par le niveau 2, le niveau 1 détermine les puissances nécessaires de zones permettant d'atteindre ces consignes.

Dans le four considéré, le niveau 2 du système de commande actuel est fondé sur un modèle de transfert thermique stationnaire. Ce modèle ne permet pas d'estimer la consommation du four. A cet effet, un choix lié à la précision de la température de défournement est réalisé, et aucun critère d'efficacité énergétique n'est considéré. Au niveau 1, des régulateurs PID traditionnels sont mis en oeuvre pour le suivi des températures des zones. Cependant, ces régulateurs PID ne peuvent pas prendre en compte l'interaction entre les zones du four, qui représente pourtant un phénomène non négligeable à l'intérieur du four. Par conséquent, le chauffage des brames n'est pas optimisé, ce qui induit un faible rendement énergétique du four. Afin d'améliorer la qualité de chauffage et de réduire la consommation d'énergie, ce travail propose une nouvelle structure de commande illustrée Fig. 3. Les régulateurs PID de niveau 1 y sont remplacés par une stratégie de commande prédictive distribuée permettant la prise en compte de l'interaction entre les zones du four. En outre, au niveau 2, le modèle d'état stationnaire est remplacé par un modèle dynamique non linéaire ; ce modèle peut prédire le comportement transitoire du four, permettant une meilleure prévision et l'estimation de la température des brames. Le travail vise principalement la qualité de chauffage et la réduction de la consommation énergétique du four. Toutefois, une réduction de la calamine et de la décarburisation peut être également espérée. De plus, la flexibilité de production est étudiée afin d'augmenter la productivité du four. Fig. 3 : Approche proposée pour la commande du four

Modélisation du four

Au sein du four, les variables à contrôler sont principalement les profils de températures de brames. Ces variables ne peuvent être mesurées directement. La modélisation du four est une étape fondamentale permettant par la suite une estimation la meilleure possible des températures de brames. Un modèle de transfert thermique est donc nécessaire, se fondant sur la connaissance des phénomènes de conduction, de rayonnement et de convection dans le four.

) (t q i  et ) (t q i
 désignent les débits thermiques absorbés par les surfaces supérieure et inférieure de la brame, respectivement. Le débit thermique ) , ( t y q i à l'intérieur de la brame dépend de la capacité thermique du matériau, de la conductivité thermique, du gradient de température à la position y, et des conditions aux limites aux deux surfaces de la brame. Le problème de conduction suit la loi de Fourier et la loi de diffusion thermique :

y t y T t y q i i i     ) , ( ) , (  (1) 0 ) , ( ) , (  
   y t y T c t y q i i i i  (2)
A partir des équations (1) et (2), nous obtenons l'équation décrivant la conduction de la chaleur dans l'épaisseur de la brame : . La méthode aux différences finies peut être ensuite utilisée pour déterminer H en fonction du temps et la variable de Kirchhoff.

              
Pour résoudre l'équation de conduction dans l'épaisseur de la brame, on peut également utiliser la méthode des volumes finis en intégrant l'équation (3) dans l'espace et dans le temps (Fig. 5), à savoir :

                     

Le rayonnement thermique

Le rayonnement est le phénomène de transfert thermique dominant dans le four de réchauffage, comme indiqué Fig. 6. Environ 90% du débit thermique absorbé par les surfaces des brames est issu du rayonnement. Actuellement, le modèle du four considéré utilise une méthode de calcul stationnaire pour le rayonnement thermique. Les hypothèses de cette méthode sont :  seules les surfaces supérieure et inférieure sont prises en compte dans l'échange thermique  l'épaisseur des brames n'est pas prise en compte  une surface d'échange thermique est homogène en température, et réflecteur idéal  la température de fumée est homogène le long du four  tout échange radiatif se produit entre la fumée, les brames et la zone correspondante  débit issu du mur reçu par la brame i w q ,  débit issu de la fumée, réfléchi sur le mur et reçu par la brame i gw q ,  débit issu de la brame, réfléchi sur le mur et reçu par la brame elle-même i sw q ,  débit issu de la fumée, reçu directement par la brame i g q , , , , , i i i g i w i gw i sw i q q q q q q q       

Le rayonnement suit la loi de Stefan-Boltzmann appliquée aux corps gris. Les éléments de (11) sont calculés par : 4 , , , , T , sont les température de la fumée, du mur de la zone j et de la brame i.

n s g g j g i g i k s g k w s g ε ε σT qq εA a ε ε α           2 1 1 1 1 (1 )(1 )(1 )
s g w w j w i s w w s g ε α ε σT q εA a ε ε α         (12) 4 , , 2 (1 ) 1 1 1 (1 )(1 )(1 )
) 1 1 1 (1 )(1 ) (13) 4 , , , 2 (1 )(1
s g w g g j gw i s gw i w s g ε α ε ε σT q εA a ε ε α        (14) 24 , , (1 )
La démarche tenant compte du comportement transitoire considère le rayonnement comme un phénomène global. L'intérieur du four est composé de plusieurs surfaces élémentaires grises et diffusantes S j de température T j et de coefficient d'émissivité j  (Fig. 8). Cette méthode nécessite des calculs numériques, à l'opposé des formules analytiques données par une approche stationnaire. Fig. 8 : Rayonnement thermique transitoire dans un four de réchauffage

La convection thermique

La convection représente environ 5% des échanges thermiques entre les brames et le four. Elle se produit quand il y a mouvement de la fumée à l'intérieur du four. La loi de refroidissement de Newton décrit la convection comme dépendant de la vitesse de la fumée, et de la différence entre la température de la fumée et la température de la brame, à savoir :

  , , , , ( ) ( ) ( 2,
)

c i c j g j s i j q t h T t T D t       (16)
où les coefficients de convection  j c h , sont proportionnels à la vitesse de la fumée.

Conception de la commande prédictive

La commande prédictive, souvent appelée MPC, est une méthode de commande avancée qui suscite un engouement fort à la fois du point de vue théorique et par les nombreuses applications industrielles. L'idée de base du MPC consiste à utiliser un modèle du système pour prédire le comportement des sorties à contrôler. En particulier, le MPC est capable de gérer des contraintes du système et de traiter les problèmes de commande multivariable. Le problème de commande est très souvent (en particulier en présence de contraintes) transformé en un problème d'optimisation, faisant intervenir pour notre application des critères économiques comme la réduction de la consommation énergétique ou l'amélioration de la qualité des produits.

Dans le cas de la commande des fours de réchauffage et avec l'approche proposée dans les travaux de thèse, le MPC est utilisé pour optimiser la commande des deux niveaux principaux. Au niveau 1, une structure MPC distribuée est construite comme le montre la Fig. 9. Dans cette structure, chaque zone du four est pilotée par un régulateur MPC local. Malgré tout, les zones interagissent du fait du mouvement de la fumée et de l'effet du rayonnement thermique. En conséquence, à l'opposé des régulateurs PID traditionnels, la structure MPC distribuée est développée pour prendre en compte les interactions entre les zones du four. Deux approches sont envisagées : une structure MPC distribuée coopérative et une stratégie MPC non-coopérative. Dans l'approche non-coopérative, chaque régulateur MPC local élabore la puissance de zone nécessaire afin d'optimiser la température de la zone en prenant en compte les influences des autres zones. Dans l'approche coopérative, chaque régulateur MPC local vise à optimiser non seulement la performance de la zone correspondante mais aussi celles des autres zones, particulièrement les zones voisines ayant des couplages directs avec la zone considérée. Au niveau 2, le contrôle de température des brames est réalisé par l'optimisation des consignes de température des zones. Différents critères sont pris en compte, comme la précision de la température de sortie des brames, l'homogénéité de la température à l'intérieur des brames, la consommation énergétique. Pour cela, le modèle précédent de rayonnement transitoire est utilisé afin de mieux capturer la dynamique du four et d'être capable également d'estimer la consommation du four, ce qui n'est pas permis par le modèle stationnaire de rayonnement thermique. Les contraintes du four sont prises en compte dans le problème d'optimisation sous forme de pénalités.

Les critères correspondant à chaque brame y sont pondérés en fonction de l'ordre de sortie de la brame. Les brames situées dans les zones de fin du four sont considérées plus importantes lors de l'optimisation parce qu'il reste moins de temps pour les chauffer à la bonne température. Par ailleurs, les brames situées dans la zone d'égalisation sont prioritaires pour atteindre l'homogénéité en température. Les puissances consommées des zones sont également pondérées. Dans une simulation du four, plus on avance dans l'horizon de prédiction, moins la puissance estimée par le modèle sera pondérée dans les critères d'optimisation.

Les paragraphes suivants détaillent la construction du MPC pour les deux niveaux.

Régulation de la température de zone

Des essais expérimentaux réalisés sur le four montrent que l'évolution de la température des zones du four se comporte comme un système du premier ordre retardé comme le montre la Fig. 10 sur un cas particulier. Les puissances u 1 , u 3 y sont maintenues constantes pendant que u 2 est alimentée par un échelon. Le système global peut dès lors être identifié par des fonctions de transfert du premier ordre avec retard comme donné ci-dessous :

M i i i i i ij j j j i i i i x k A x k B u k B y k y k C x k            (18) ( 1) ( ) ( ) ( ) ( ) ( )
où M est le nombre de zones du four. La relation (18) permet alors de prédire les sorties du système sur l'horizon de prédiction N P à partir d'une séquence de commande sur l'horizon de commande N u , comme suit :

y k s k C A x k C A B u k p k C A B y k p k s N

               (19)
Dans une approche non-coopérative, le régulateur local de la zone i résout le problème de minimisation suivant :

1 2 2 , , , 1 ( )
i i i u i i i u i i i p u u k p u p N u u k p u p N y y k s k y s N               (20)
Pour cela, un algorithme itératif est mis en oeuvre, décrit par l'organigramme de la Fig. 11. La communication entre régulateurs locaux est décrite Fig. 12. A chaque instant k, chaque régulateur communique les prédictions de la température de la zone correspondante et reçoit celles envoyées par les autres régulateurs. A partir des informations reçues, chaque régulateur MPC résout son problème d'optimisation locale, et ensuite met à jour les prédictions de température. La première valeur de la séquence de commande après convergence est ensuite appliquée à la zone correspondante. Fig. 12: Communication au sein de la structure MPC distribuée non-coopérative Dans une approche coopérative, le régulateur de la zone i résout le problème de minimisation suivant :

1 2 2 , () 1 ( ) ( ) ( ) ( ) , min i i j j ref i R Q i U k k j i ji J k Y k k Y k U k k       
                (21) Avec ' '' ˆˆ( ) ( 1 )
,..., ( ) ,

i i i p Y k k y k k y k N k      ' ''
( ) ( ),..., ( 1 ) . Le régulateur de la zone i calcule ainsi la commande U i (k|k) en optimisant non seulement la performance de la zone i mais aussi celle des zones voisines i+1 et i-1. Pour cela, les régulateurs locaux ont besoin de plus de communication entre eux. Un régulateur communique à la fois ses prédictions de température, sa séquence de commande et l'état de sa zone, et reçoit les informations similaires émises par les autres régulateurs (Fig. 13).

i i i u U k k u k k u k N k      (22)
Les deux approches précédentes ont été testées sur un système couplé défini par les relations :  

yk x k x k u k yk y k x k                                   
yk x k x k u k yk y k x k                             1
yk x k x k u k yk y k x k                                            
Fig. 14 : Résultats issus de la commande MPC distribuée Les résultats obtenus par les deux approches MPC distribuées (coopérative et non-coopérative) sont donnés Fig. 14. Nous constatons que les deux approches procurent de bons résultats, en particulier l'algorithme non-coopératif même en présence de couplage. Du fait du couplage peu important entre les zones de réchauffage, la stratégie MPC distribuée non-coopérative a été privilégiée pour la comparaison avec la structure PID et l'implantation sur site. Le train à chaud du site de Florange d'ArcelorMittal possède deux fours de réchauffage ayant la même technologie et chauffant souvent les mêmes types de produits. A titre de comparaison, le MPC distribué non-coopératif a été implanté au sein du niveau 1 du four A tandis que les régulateurs de type PID étaient conservés dans le four B. La similarité de fonctionnement des deux fours permet de comparer efficacement les comportements obtenus avec le MPC distribué non-coopératif et les régulateurs PID, comme le montre la Fig. 15. On constate alors clairement que la structure MPC procure moins de fluctuations de température et de dépassements de consignes.

Régulation de la température de brame

Cette partie est dédiée à la commande de niveau 2 dans le cas de cadence prévisionnelle fixe. Les données de cadencement des brames sont envoyées par le niveau 3, cela signifie que le temps de séjour de chaque brame et le contenu du four à un instant donné sont connus à l'avance.

La boucle fermée cascade de niveau 1 et niveau 2 est illustrée Fig. 19. La température visée de la brame j, notée , j end T , est reçue depuis le niveau 3. La température réalisée de la brame j, end j T , est estimée via le modèle de transfert thermique transitoire. Ce modèle utilise les mesures de puissance P i et celles de température des zones T z,i afin de calculer la température de brames. Le niveau 2 optimise les consignes de température des zones i z T , ~ via un modèle simplifié du four et la stratégie MPC. Ce modèle simplifié est construit sur le même principe que le modèle transitoire, la différence concerne le pas de discrétisation spatiale et échantillonnage temporel. Le MPC de niveau 1 vu précédemment assure le suivi de consigne de température de zone en optimisant la puissance de zones. Fig. 19 : Boucle fermée cascade du système de commande niveaux 1 et 2 avec stratégie MPC Les objectifs de la commande de niveau 2 sont : erreur de température de sortie de brames minimale, inhomogénéité de la température à l'intérieur de la brame minimale, perte de matière issue de la formation de calamine minimale, et décarburation ayant une influence directe sur la qualité du produit minimale. Les contraintes de la commande sont liées aux limites des températures de zones, aux puissances maximales de zones et à la température des brames. Les objectifs et contraintes sur la température des brames sont illustrés Fig. 20, Fig. 21 et Fig. 22. La Fig. 20 montre que la température de brames est contrainte par T j,max pour éviter les pertes d'énergie et de matière. La température de défournement des brames est contrainte dans un intervalle [T j,end,min , T j,end,max ], voir Fig. 21. En effet, cette marge est imposée pour prendre en compte l'erreur issue de la régulation de température. Plus la commande est précise, plus cette marge peut être réduite. Le gradient de température est aussi objet de contraintes comme l'illustre la Fig. 22. La limite supérieure du gradient ΔT j,max est imposée pour éviter la déformation de la brame. A la sortie du four, les brames doivent avoir une bonne homogénéité définie par ΔT j,end,max pour être acceptées dans les laminoirs. ,max ,0

z i z i z i z i z i z i T t T min , , ,max , ,min , , ,max ( ) ( ) ( ) ( ) ( ) ( )
, .

j j j j j j j j T y t T y D D t t t               (25) , ,min ,1 , ,max ( , ) , 2, 2
j end j j j j end j j j T T y t T y D D        (26)
,max ,0

              (28)
Ces contraintes de température de brames sont récapitulées dans le Tab. 2 ci-dessous.

Variables Description

Valeurs nominales

,max j T Borne supérieure de température de défournement de brame 1350 °C Le problème de la régulation de la température de brames peut être formalisé comme le problème d'asservissement de la température de défournement moyenne (et de son gradient) définie par :

  ,1 ,1 ,1 ,1 ( , ) , Si , () 0,
                  (29)
où S est l'ensemble des brames à chauffer par le four. Les consignes de température à suivre et contraintes associées sont alors définies par :

, , 1 , ,min ,1 ,min , ,max ,1 ,max , Si , () 0, autrement , Si , () 0, autrement
, Si , () 0, autrement j end j s j end j s j end j s

T j S t t Tt T j S t t Tt T j S t t Tt

                           (30) , ,max ,1 max

Δ

, Si , () 0, autrement

j end j T j S t t gt          (31)
Pour les contraintes de température de brames, celles-ci s'expriment en fonction des variables d'intérêt x 1,

x 2 :      
              0 1 1 2 2 3 3 4 4 5 5 () ( ) ( ) ( ) ( ) ( ) ( ) min U τ JU τ α J U τ α J U τ α J U τ α J U τ α J U τ      (33)
sous les contraintes des températures de zones, des température de sortie, des gradients de température de sortie définies comme suit : , ,min

                          
    1 1 ,1 ,1 ,1 () ( ) ( ) ( ) ( ) , , Δ j s j s j j p jSτ JUτ w z τ T t T t t τ τ t          (35)     2 2 ,1 ,1 () ( ) ( ) ( ), , Δ j j j p jSτ JU τ w z τ g t t τ τ t         (36)
où w 1 (z j ) et w 2 (z j ) sont des coefficients de pondération de la brame j, dépendant de la position z j de la brame. De même les autres termes du critère sont donnés par les relations :

  Δ 33 1 ( ) ( ) ( ) τt p M i i τ JU τ w t P t dt      (37)   Δ 44 1 ( ) ( ) ( ) τt p M i i τ JU τ w t u t dt      (38)   Δ 5 1 5 1 ( ) ( ) ( ) ( ) τt p T τ JU τ x t w t x t dt        (39)
Les contraintes (34) peuvent être exprimées sous forme de pénalités :

    Δ , 1 , ,max , ,min 1 (
) max 0, ( ) ( ) min 0, ( ) ( ) τt c M c i z i i z i i τ J τ u t T t u t T t dt           (40)     Δ , 2 , ,max , ,min 1 (
) max 0, ( ) ( ) min 0, ( ) ( ) τt c M c i z i i z i i τ J τ u t T t u t T t dt            (41)     Δ ,3 max ,1 , ,max () ( ) max 0, ( ) ( ) max 0, ( ) Δ τt p c j j end jSτ τ J τ g t g t dt g t T         (42)         Δ , 4 ,max ,min ,1 ,max ,1 ,1 ,min ,1 () (
) max 0, ( ) ( ) min 0, ( ) ( ) max 0, ( ) ( ) min 0, ( ) ( )

τt p c s s s s τ s j s j s j s j jSτ J τ T t T t T t T t dt T t T t T t T t                 (43)             Δ ,5 1 max 1 max Δ ,6 2 max 2 max ( ) max 0, ( ) ( ) ( ) ( ) ( ) max 0, ( ) Δ ( ) ( ) Δ ( ) τt p T c τ τt p T c τ J τ x t T t x t T t dt J τ x t T t x t T t dt           (44)
Par conséquent, la minimisation de (33) sous les contraintes (34) devient une minimisation sans contrainte :

      56 0, () 11 ( ) ( ) ( ) min i i c i U τ ii JU τ α J U τ μ J U τ         (45)
La discrétisation temporelle est alors réalisée par substitution aux variables continues des variables discrètes récapitulées dans le Tab.

    1 1 ,1 ,1 ,1 () ( ) ( ) ( ) ( ) , k , j s j s j j p j S k J U k w z k T k T k k k N          (46)     2 2 ,1 ,1 () ( ) ( ) ( ), , j j j p j S k J U k w z k g k k k k N         (47)   33 11 ( ) ( ) ( ) N p M i ji J U k w k j P k j      (48)   44 11 ( ) ( ) ( ) N M c i ji J U k w k j u k j      (49)   5 1 5 1 1 ( ) ( ) ( ) ( ) N p T i J U k x k i w k i x k i         (50)     ,1 , ,max , ,min 11 ( ) max 0, ( ) ( ) min 0, ( ) ( ) N M c c j z j j z j ij J k u k i T k i u k i T k i             (51)     ,2 , ,max , ,min 11 ( ) max 0, ( ) ( ) min 0, ( ) ( ) N M c c j z j j z j ij J k u k i T k i u k i T k i              (52)   ,3 ,1 , ,max ()
( ) max 0, ( ) Δ c j j end j S k J k g k T    (53)     ,4 ,1 ,max ,1 ,1 ,min ,1 ()
( ) max 0, ( ) ( ) min 0, ( ) ( ) c s j s j s j s j j S k J k T k T k T k T k         (54)             ,5 1 max 1 max 1 ,6 2 max 2 max 1 ( ) max 0, ( ) ( ) ( ) ( )
( ) max 0, ( ) Δ ( ) ( ) Δ ( ) N p T c i N p T c i J k x k i T k i x k i T k i J τ x k i T k i x k i T k i                   (55)
De sorte que la solution recherchée sera issue de la minimisation du critère suivant sous forme discrète :

      56 0 , () 11 ( ) ( ) ( ) min i i c i Uk ii J U k α J U k μ J U k         (56)
Le coefficient associé aux termes de pénalités,  , est significativement plus grand que les autres coefficients. Les coefficients 𝛼 1 , 𝛼 2 associés aux objectifs sur la température de défournement des brames sont plus forts (car objectifs plus prioritaires) que ceux liés aux objectifs de consommation. Nous avons donc la relation entre les coefficients :

𝜇 ≫ 1 2 3 4 5 , , , α α α α α  (57)
L'horizon de prédiction doit être choisi de sorte que l'optimisation prenne en compte toutes les brames situées dans les zones actives du four. Le temps de séjour d'une brame dans les zones actives représente environ 2/3 du temps de séjour total de la brame dans le four, par conséquent nous avons la contrainte sur l'horizon de prédiction sur la sortie :

  ,1 ,0 2 Δ 3 p j j t t t  (58)
Le temps d'échantillonnage du niveau 2 correspond à l'intervalle temporel de défournement des brames qui est de l'ordre de 5 min.

La résolution de problème de minimisation sans contrainte (54) est fondée sur la simulation du four comme le montre la Fig. 16. La simulation de four est réalisée via un programme numérique complexe, par conséquent ne permettant pas d'obtenir une forme analytique de la fonction de cout (54) ni de ses dérivées. Il est alors intéressant d'utiliser l'algorithme d'optimisation du simplex Nelder-Mead, qui a uniquement besoins de la valeur de la fonction de coût. Cet algorithme est connu pour sa faible charge de calcul tout en fournissant une solution satisfaisante dans les applications industrielles. Fig. 23 : Optimisation de la stratégie MPC fondée sur la simulation du four La simulation de tout cet ensemble pour la validation de la commande (Fig. 24) est effectuée en Fortran, utilisant le simulateur du four fondé sur le modèle transitoire de four. Dans le scénario de test, des données de trois jours du four sont utilisées ; environ 900 brames sont chauffées durant cette période de trois jours. Les consignes de température réelles appliquées dans les fours sont réutilisées en parallèle avec celles données par le régulateur MPC. Les résultats sont ensuite comparés en termes de température finale de brames, de gradient de température, et de consommation du four. Les consignes de température de la zone d'égalisation sont données sur la Fig. 30. On constate que les consignes calculées par la structure MPC sont inférieures à celles fournies par le régulateur existant. En revanche, concernant les consignes des zones de chauffe 2 et de chauffe 1 des Fig. 31 et Fig. 32, on remarque que les consignes données par la loi MPC sont en général supérieures à celles données par la commande existante. En zone de préchauffage, les consignes de température données par la loi MPC Fig. 33 La Fig. 35 fournit l'histogramme de l'erreur de la température de défournement des brames par rapport à la température visée. Le nombre de brames défournées ayant une erreur de moins de 5°C a augmenté de 10%. La surchauffe est réduite, ce qui peut avoir un effet positif sur la consommation du four.

Augmentation de la productivité

La régulation de la température de la brame peut être effectuée par la régulation de la température proprement dite ou par la régulation du cadencement des brames, ou par un mixte des deux solutions. Dans la section précédente, nous avons mis l'accent sur le problème de l'optimisation de la consigne de température de zone du niveau 2 et de la régulation de la température de ces zones au niveau 1, avec un cadencement fixe des brames. Or le cadencement des brames est connu à l'avance, via une table fournissant l'intervalle de temps entre les défournements consécutifs des brames comme le montre la Fig. 45. Dans la pratique, lorsque les laminoirs à chaud ne constituent pas un bouchon de la chaîne de production, la priorité du four de réchauffage est de maximiser sa productivité, et de toute évidence, dans le même temps, d'assurer l'exigence de qualité de chauffage. Par conséquent, nous allons reprendre la solution précédemment construite pour l'optimisation des consignes de température et l'adapter au problème de maximisation de la productivité.

Identifiant de brame

Intervalle

Les variables de commande sont maintenant l'intervalle moyen de défournement et les consignes de températures de zones :

,1 ,2 , ( ) ( )
, ( ), ( ),..., ( )

sch z z z M uk δt k T k T k T k    (57)
Le problème de minimisation est réécrit comme suit :

      67 0 , () 11 ( ) ( ) ( ) min i i c i Uk ii J U k α J U k μ J U k         (58)
où les termes J 1 ,…,J 5 et J c,1 ,…, J c,6 sont ceux définis dans la section précédente. J 6 représente l'objectif visant à réduire l'intervalle de temps entre défournements et est déterminé par :

  6 6 1 1 ( ) ( ) ( ) N c i J U k w k i u k i       (59)
Le terme de pénalité J c,7 représente la contrainte de limite inférieure de l'intervalle de temps 𝛿𝑡 𝑠𝑐ℎ,𝑚𝑖𝑛 qui définit la capacité physique du laminoir à chaud :

    ,7
,min 1 ( ) min 0, ( ) Le but de ce travail était de développer une stratégie de commande prédictive MPC d'un four de réchauffage de brames au sein d'une architecture distribuée hiérarchisée. La fonction première d'un four de réchauffage de brames est de ramener la température des brames d'acier, partant d'une température initiale donnée, jusqu'à environ 1200°C avec une bonne homogénéité avant que les brames n'entrent dans les laminoirs à chaud. La modélisation du four s'avère indispensable pour l'observation des variables telles que la température des brames, la température des parois, etc., et pour la régulation de la température du four. Un modèle permettant de reproduire le comportement transitoire du four a dès lors été élaboré, venant remplacer le modèle stationnaire actuel du four. Ce nouveau modèle permettra en outre de prendre en compte la diversification croissante des brames en termes de température visée, de dimensions, de propriétés des matériaux, etc. Le modèle transitoire a été validé sur des données industrielles et est ensuite utilisé dans la conception de la commande de four.

N c c sch sch i J U k δt k i δt      (
La tâche de commande des fours de réchauffage concerne principalement deux niveaux. Le niveau 1 a pour objectif la commande des températures de zone du four en agissant sur les débits de gaz des zones (boucle interne). Le niveau 2 quant à lui doit assurer la maitrise de la température des brames en faisant varier les consignes de température de zones (boucle externe).

Le niveau 1 doit permettre de suivre avec précision les consignes de température de zones du four, tout en tenant compte de l'objectif de minimisation de la consommation d'énergie. Pour résoudre ce problème de commande, une stratégie MPC distribuée a été élaborée en se fondant sur un modèle de zones du four identifié et tenant compte de l'interaction entre les zones, déduit expérimentalement sous la forme de fonctions de transfert du premier ordre avec retard pour une entrée en puissances et des sorties en température. Deux approches ont été mises en oeuvre pour capturer cet effet de couplage entre les zones du four. La première approche considère que la température d'une zone est une fonction non seulement de l'entrée de cette zone, mais aussi de celles de ses zones voisines. Cela conduit à ce que l'on appelle un modèle de couplage par les entrées. La seconde approche est plus directe; considérant que la température d'une zone est fonction de la puissance de cette zone et de la température de ses zones voisines. Cela conduit à ce que l'on appelle un modèle de couplage par les sorties. Les résultats des validations ayant montré que le modèle de couplage par les sorties s'avérait meilleur que le couplage par les entrées, ce modèle de couplage par les sorties a été retenu pour la conception de la structure MPC distribuée.

En simulation dans l'environnement Matlab dans un premier temps, les algorithmes MPC distribués proposés ont été testés et validés à partir d'un scénario type issu du fonctionnement réel du four. Les résultats avec un modèle validé de zones de four illustrent un meilleur suivi et une consommation d'énergie plus faible avec la loi MPC distribuée par rapport aux régulateurs PID actuels. La méthode proposée a enfin été mise en oeuvre expérimentalement au sein d'un niveau 1 d'un four réel. Corroborant les observations en simulation, les résultats obtenus montrent que la commande MPC distribuée permet d'économiser jusqu'à 3% de la consommation d'énergie tout en assurant un bon suivi des consignes de température de zone. En outre, le taux d'oxygène dans le four est réduit, ce qui conduit à moins de formation de calamine et de décarburation, et donc moins de pertes de matière. La fluctuation de demande globale de gaz est réduite, ce qui induit une pression plus stable de la station d'alimentation en gaz. Par ailleurs, l'utilisation du mode automatique est augmentée de 10% dans la zone d'égalisation, ce qui signifie que l'opérateur du four a moins d'intervention manuelle à réaliser pendant le fonctionnement du four.

En termes d'algorithme d'optimisation, des algorithmes de recherche directs ont été étudiés, et afin d'éviter d'avoir recours à une forme analytique du gradient de la fonction de coût ou de son estimation, l'algorithme de Nelder-Mead a été privilégié et mis en oeuvre en raison de ses avantages, notamment, en plus du précédent, une mise en oeuvre simple et une décroissance rapide de la fonction de coût.

En simulation dans un environnement Fortran, la stratégie et l'optimisation MPC proposées pour le niveau 2 ont été testées à partir d'un scénario réel de four. Dans ce scénario, 900 brames sont traitées sur trois jours avec un temps de séjour moyen de 180 min. La production est de 281 tonnes/h. La simulation réalisée avec le simulateur validé du four permet une comparaison de la stratégie MPC et de la stratégie existante de niveau 2. Les résultats obtenus montrent que l'utilisation de la loi MPC peut économiser 5% de la consommation d'énergie du four tout en assurant la température visée et l'homogénéité de la température de brames. De surcroit, malgré un modèle transitoire coûteux utilisé lors de l'optimisation, l'algorithme du simplexe Nelder-Mead est en mesure de fournir les consignes de température dans la période d'échantillonnage imposée, ce qui signifie que le régulateur MPC proposé pour le niveau 2 peut être mis en oeuvre en temps réel.

Une extension de la démarche précédente visant la maximisation de la productivité du four a également été effectuée. Dans cette approche, le cadencement des brames et l'optimisation de la température de brame sont considérés simultanément. En plus de la fonction de coût utilisée précédemment lors de l'optimisation des consignes de température de zones à cadence fixe, l'intervalle de temps moyen entre les défournements est ajouté cette fois comme variable de contrôle. En tant que nouvelle entrée de commande à côté des consignes de température de zones, l'intervalle de temps moyen entre les défournements est soumis à une limite inférieure en raison des contraintes des laminoirs à chaud.

Une originalité dans la conception du MPC visant l'augmentation de la capacité du four est liée au fait que le nombre de brames prises en compte à chaque exécution de l'optimisation est constant alors que l'horizon temporel de prédiction est variable d'une simulation du four à une autre.

La simulation dans l'environnement Fortran montre alors que la production du four peut être augmentée avec un surcoût très modéré de la consommation d'énergie. En outre, la précision de la température de sortie de brames et son homogénéité peut être significativement améliorée. Les résultats obtenus montrent que la productivité de four peut être augmentée jusqu'à 15 tonnes/h. Dans la pratique, ces résultats prouvent que le MPC peut être utilisé pour informer l'opérateur du four, lui permettant d'ajuster ainsi la productivité instantanée du four. Cette information aidera l'opérateur du four pour ajuster le cadencement des brames afin que le four puisse fonctionner à sa productivité maximale.

Perspectives

La stratégie MPC distribuée proposée a été mise en place et fonctionne actuellement à temps plein dans les fours de l'usine de laminage à chaud de l'usine ArcelorMittal de Florange. Les résultats du MPC niveau 2 obtenus en simulation à partir des données du système réel sont convaincants, de sorte que cette stratégie MPC sera implantée à court terme sur le four pour sa validation expérimentale. Sur le même laminoir à chaud se trouve un autre four de réchauffage de conception similaire mais avec un dimensionnement des zones différent. Une adaptation à court terme des stratégies développées est donc à prévoir pour ce nouveau four, qui ne concernera en fait que la modification des paramètres du modèle du four.

Dans le développement du MPC pour le niveau 2, nous avons utilisé des bornes constantes de températures de zone lors de la formulation des contraintes du problème d'optimisation. Dans la mise en oeuvre expérimentale sur four réel, ces limites seront amenées à varier dans le temps. En effet, à chaque instant, les limites de températures de zone dépendent non seulement de la puissance disponible dans les zones du four, mais aussi de la charge de brames située dans chaque zone. Pour résoudre ce problème, le modèle de four pourra effectivement donner des estimations de ces limites qui pourraient ensuite être utilisées pour l'optimisation.

Enfin, la grande perspective de ce travail est à plus long terme la mise en oeuvre du MPC niveau 2 pour tous les fours de réchauffage de brames des usines sidérurgiques d'ArcelorMittal, qui devrait avoir des impacts économiques et environnementaux significatifs. Plus largement donc, il conviendra de généraliser la démarche à d'autres types de fours de réchauffage dédiés à différents types de produits tels que les blooms, les billettes, les tubes, afin de leur appliquer cette stratégie MPC distribuée et hiérarchisée.

List of Tables

Part I

Process description and modeling

Chapter 1

Introduction

In this chapter, a brief introduction to steel industry is firstly given. Next, the context of slab reheating furnaces is introduced. The challenge of control task and its influence on energy performance of the furnaces are then presented. Finally, an overview of the dissertation is given summarizing the flow of the work through the chapters.

Context and motivation of the work

In steel industry, iron ore (F e 2 O 3 ) and ferrous scrap (F e) are transformed into steel products. Iron must undergo several thermal, mechanical, or thermo-mechanical processes before it is suitable for other industrial applications such as construction, automotive, or space industry, etc. Those processes can be grouped into four stages: iron making, steel making, continuous casting, and hot rolling. In iron making, molten metal is produced by melting ferrous scrap or burning the mixture of coke and iron ore at 1450 • C in blast furnaces [START_REF] Trinks | Industrial furnaces[END_REF]. In steel making, impurities such as nitrogen, silicon and excess carbon are removed from raw iron liquid. Then, alloying elements such as manganese, nickel, chromium and vanadium are added to produce different grades of steel. In the step of continuous casting, the obtained steel liquid in steel making is cooled down and formed into solid blocks.

These steel blocks are semi-finished products that have form of slabs with dimensions of 0.25m x1.5m x10m for example. These slabs are then stocked in a slab yard or transported directly to hot rolling process when they are needed to fill client orders. The primary function of hot rolling process is to heat steel slabs nearly to their melting point, then roll them thinner and longer through successive rolling mill stands driven by motors and finally coil up the lengthened steel sheet for transport to the next process, see Figure 1.1.

Reheating furnaces are positioned just before the hot rolling mills, utilized for reheating, and also for chemical and metallurgical treatment of steel slabs. The role of these furnaces is to provide sufficient hot slabs for rolling mills. The furnace design is based on capacity of hot rolling mill. Steel slabs are usually heated from ambient temperature to 1250 • C; this case is called cold charging. In an integrated plant which has both continuous casting and hot rolling processes, slabs might be loaded into reheating furnaces at 800 • C [START_REF] Mårtensson | Energy efficiency improvement by measurement and control: A case study of reheating furnaces in the steel industry[END_REF]; this case is called To realize the task, reheating furnaces are constructed with several zones which can be thermally isolated or coupled one from others. Transportation mechanism is a system of mobile and fixed skids. In order to move the slabs forward, the mobile skids are firstly lifted up, next translated, then brought down to the same level of the fixed skids, and finally brought back to the initial position. The presented mechanism is used in walking-beam reheating furnace type which differs from pusher-type reheating furnace where slab rows are pushed by machine at the entry of furnace and slabs slide on the support skids. For instance, the considered furnace in our study is a walking-beam type furnace, it has dimensions of 6m x12m x60m and a production rate that can come up to 400 tons of steel per hour.

In terms of energy consumption, about 20 percent of steel production cost is from energy use [START_REF] Lankford | United[END_REF][START_REF] Sakamoto | Estimation of energy consumption for each process in the japanese steel industry: a process analysis[END_REF].

In a steel plant, the largest energy consumers are blast furnaces in ironmaking process with 44 percent of total energy use and reheating furnaces of hot rolling process are second energy consumers with 17 percent of total energy use [START_REF] Sakamoto | Estimation of energy consumption for each process in the japanese steel industry: a process analysis[END_REF]. Particularly, annual energy consumption of the considered slab reheating furnace is equivalent to the consumption of 100,000 houses [START_REF] Poquet | Pour les ménages, la recherche du confort prime encore, sur les économies dénergie[END_REF]. Therefore, the energy efficiency of reheating furnaces influences significantly the energy performance of the whole production chain. The aftereffects of energy savings are the lowering of greenhouse gas (CO 2 , N 2 O, etc.) emissions and the reduction of environmental footprint which are subjects to ever more restrictive regulations.

Many design factors affect the energy efficiency of reheating furnaces like furnace geometry, number of zones, wall and roof insulation, and design of transport mechanism [START_REF] Mårtensson | Energy efficiency improvement by measurement and control: A case study of reheating furnaces in the steel industry[END_REF]. Charging temperature of slabs may vary from ambient temperature to 800 • C as mentioned before, which is a significant condition affecting energy efficiency of the furnace. These are prefixed conditions of a furnace and can not be intervened in temperature control system of furnace.

Furthermore, operational practices are also crucial factors for the energy efficiency. Operational conditions may change from mill to mill but also for the normal operation of each furnace in the same mill. Different disturbances and deviations occur in a furnace operation like delays in the rolling mill, diversity of slab dimensions, fuel composition, charging temperature. For these reasons, optimal manual control is in general an overwhelming task. Even skilled operators generally tend to run the furnace with large safety margins to assure the required discharge temperature of slabs. This practice often leads to excess use of energy and to quality problem due to overheating of slabs.

To increase energy efficiency, more efficient control of the furnace is required. More accurate control will improve accuracy of slab temperature and combustion process. Therefore, safety margins on slab temperatures will be reduced and the overall temperature of the furnace will be lowered.

The control of a slab reheating furnace is an intricate task. From theoretical point of view, reheating furnaces are non-linear systems with dynamic interaction of multiple subsystems which refer to furnace zones. The thermal behavior of the system is rather slow due to high thermal inertia of slabs and furnace enclosure. Efficient operation of the furnace essentially needs knowledge of certain physical quantities such as slab temperature profile, furnace gas and wall temperatures, fuel and combustion air flows, oxygen content of exhaust gas, position of each slab inside the furnace [START_REF] Mårtensson | Energy efficiency improvement by measurement and control: A case study of reheating furnaces in the steel industry[END_REF]. Unfortunately, these measurements may be difficult, unreliable, or even in-feasible as a consequence of harsh environmental conditions of the furnace like high temperature, dirt, dust, electrical noise, wetness [START_REF] Steinböck | Model-based control and optimization of a continuous slab reheating furnace[END_REF]. Since slab temperature profile can not be directly measured; it is usually estimated by a heat transfer model. In addition, homogeneity of slab temperature profile is a critical quantity that needs to meet requirement to prevent rejection of slabs in rolling mills. Process variables such as refractory walls and roof temperature, fuel and combustion air flow rates are also subject to constraints. Furthermore, the diversity of product is continuously increasing in terms of dimensions, steel grade, desired discharge temperature. Changes of product types in reheating furnace are more and more frequent. Therefore, heating process is usually a non-steady state operation.

The purpose of temperature control of reheating furnaces can be simplified as regulation of thermal energy injection and distribution along the furnace length to heat slabs to correct temperature. Control variables available for temperature control are burner fuel and air flows or more indirectly zone temperature set-points for burner regulators. In some cases, it is also possible to adjust the throughput rate of furnace with condition that it does not exceed the rolling mill capacity. Dealing with complexity of the system, the current control structure of the furnace exploits the hierarchical cascade control concept, as shown in Figure 1.2. At level 2, desired discharge temperatures of slabs are tracked via a heat transfer model of furnace. The control variables are temperature of furnace zones. Level 2 then communicates with level 1 the required zone set-point temperature to obtain the target temperatures of slabs. At level 1, necessary heat power to achieve set-point temperatures of each zone is calculated and then the corresponding fuel and combustion air flow rates are derived. In the considered furnace, current level 2 of control system is designed based on a steady-state model. This model can estimate slab temperature profiles using radiations between slabs, furnace enclosure and gas. The estimations of temperatures are used to decide the temperature set points of each zone. For this purpose, a performance-based choice is used, and no energy efficiency criterion is considered. At level 1, traditional PID controllers are used to control zone temperatures. In fact, PID controllers can not consider the interaction of furnace zones which is an important phenomenon inside furnace. Therefore, fuel combustion process is not optimized, which leads to low energy efficiency. The work mainly aims for heating quality and energy efficiency improvements. However, the reduction of scale formation, decarburization of slabs are also expected to come with those improvements. In addition, flexibility of production rate is studied as well to boost furnace productivity.

Structure of this dissertation

After the Abstract and the Introduction, Chapter 2 presents in details the important operational aspects of walkingbeam slab reheating furnace such as combustion system, mechanical system, energy flows inside the furnace, etc.

In Chapter 3, a brief literature review on reheating furnace modeling is given. The steady-state mathematical model that gives the basis for current control system of furnace is described. Afterwards, a dynamical non-steadystate numerical model is presented. The advantages, disadvantages and utilization of two models are discussed.

At the end of the chapter, measurement data from the considered walking-beam slab reheating furnace is used to validate the numerical model. Chapter 4 will give the overview of significant works on control system of reheating furnaces. The existing control strategy of considered furnace is described. And the proposed new solution is then presented.

In Chapter 5, an algorithm of distributed MPC is developed for level 1. First, modeling and identification of furnace zones are addressed. Second, formulation of optimization and iterative algorithm of distributed MPC is presented. Finally, simulation results and industrial validation are discussed at the end of the chapter.

Chapter 6 is dedicated to the control at level 2 where set points of zone temperature are optimized. Firstly, the optimization is formulated. Then, potential dynamic optimization methods are discussed. Finally, simulation results based on industrial scenario are presented and discussed.

Chapter 7 studies the improvement of production flexibility and boosting furnace production rate. Literature review on important studies is given. And the optimization problem is formulated. Simulation results are discussed at the end of the chapter showing possibility for an increase of furnace throughput.

In Chapter 8, an effort of generalization of control solution is given for different types of furnaces with considerations of their distinctive characteristics.

Conclusions and perspectives of future work are given in Chapter 9.

Chapter 2

Slab walking-beam reheating furnace

The context and motivation of this work have been clarified in the Introduction. This chapter describes more in detail the considered walking-beam slab reheating furnace by firstly distinguishing it with other types of reheating furnace. Second, the physical structure and transportation mechanism of the furnace is detailed. The combustion, insulation and cooling systems are then presented. Flows of energy inside the furnace are also discussed.

Transportation of slabs

Steel slabs are conveyed from slab yard to the entry of the furnace. An automated mechanism is programmed to take the slab standing at the furnace door when there is enough space and predetermined loading event of the slab happens. Between slabs, some distances are imposed to facilitate gas circulation and therefore heat exchange inside the furnace, see Figure 2.1.

Inside the furnace, slabs are placed on skid beam support system. These beams are subject to heat corrosion because the operating temperature are usually around 1250 • C. Nowadays, skid beams may consist of watercooled steel members and topped with refractories which is the only parts exposed directly to the heat of the furnace.

The movement of slabs can be executed in the two directions of furnace length, which can help modifying the slab sequencing. In the furnace, two type of skid beams are utilized, one is mobile also called walking beam and the other is stationary. To move slabs forward, the mobile skid beams are lifted up above the stationary ones, which also brings up the slabs. The whole set mobile skids and slabs is then translated forward through a given step length. Next, mobile skids are pushed down so that slabs return to the level of fixed skids. Finally, mobile skids come back to the initial position ready for a new cycle of movement. Slabs are now in new positions inside the furnace.

Slab movement is a discontinuous process. At one position, a slab will stay for an amount of time so that it can absorb enough heat before being discharged from the furnace, see Figure 2.3. The scheduling of slab movements is programmed in advance. The loading and discharging events of a slab are referenced relatively to those of the previous slab. For instance, a slab could be loaded 5 minutes and discharged 4 minutes after the previous slab.

Therefore, the scheduling task can be simplified into imposing interval time of loading events and discharging events of consecutive slabs.

In general, a walking beam movement happens following four configurations. The first one is when a new slab enters, and a slab is discharged from the furnace. Second, no slabs are removed from furnace but a new slab is loaded into it. In the third case, a slab is discharged but no slabs enter the furnace. And in final case, all slabs are moved forward but no slabs are removed from or loaded into the furnace.

The contact between slabs and fixed skids causes marques on the lower surface of slabs. These marques indicate lower temperature of steel compared to other regions of slab surface. Therefore, as shown in Figure 2.2, fixed skids of soaking zone are interspersed between those of other zones to reduce the apparition of marques.

Combustion and heating of the furnace

Thermal energy is supplied into furnace by gas burners which are installed in the refractory wall of the furnace.

Each zone of the furnace is equipped with a certain amount of burners, and its input power is limited. In addition, between upper and lower sub-zones of the same zone, burners are coupled with a certain ratio so that heat flux absorbed into upper and lower surface are as similar as possible.

In terms of combustion process, gas and combustion air must be appropriately mixed. With a required heat power for a zone, the gas flow rate of each burner is determined based on the gas heating capacity which defines the thermal energy released as the gas undergoes complete combustion. The flow rate of combustion air is then computed accordingly to gas flow rate. In practice, air/gas ratio is controlled to avoid excessive rich burning which results in complete combustion with partially burned or unburned gas escaping from the furnace without releasing heat where it can be used effectively [START_REF] Trinks | Industrial furnaces[END_REF]. On the other hand, air/fuel ratio control also prevents excessive lean burning, in which extra unused air passes through the furnace, absorbs heat, and carries that heat out of the furnace, therefore unabsorbed by the loads [START_REF] Trinks | Industrial furnaces[END_REF].

Heating process of a slab starts by the exchange of energy with the exhaust gas in the recuperation zone. This zone is not provided with burners. At preheating and heating zones, main heat inputs are released to bring the average temperature of slabs to the desired ones. However, internal temperature gradient of slabs are still high at the end of heating zone. In fact, slab surface temperature can rise rapidly but due to thermal inertia of slab the core slab temperature takes more time to increase. At the end of heating zone, slabs are generally not ready for rolling. Therefore, soaking zone is dedicated for equalization of internal temperature of slabs. In general, soaking zone is isolated from other zones to preserve heat of slabs. Surface temperatures of slabs are maintained at desired discharge temperature so that homogenization will take place during residence period inside the zone.

Temperature evolution of a slab inside the furnace is shown in Figure 2.4. A phenomenon that may significantly affect heating process is the scale formation. The slabs inside furnace are exposed to hot-oxidizing combustion gas that contains H 2 0, O 2 , C0 2 and N 2 . These hot gases transfer also heat into slabs. However, reactions between slab surfaces and furnace gas happen as slab temperatures increase, which causes the formation of an iron oxide film usually called scale. The thickness of the scale depends on residence time in the furnace of slabs, their surface temperature, temperature and composition of the furnace gas.

The scale formation causes material losses of the slabs. Moreover, because thermal conductivity of the scale is much smaller than that of steel, the existence of the scale can have high effect on heat transfer behavior of slabs.

On the other hand, decarburization is a chemical process that affects quality of slabs. Decarburization is the change of structure and content of steel slabs in which some carbon in the slab surfaces is lost due to the reaction with O 2 . Since carbon content determines the physical performance of steel, the decarburization phenomenon of the furnace has to be reduced.

Temperature measurement

Temperature control can not have better precision than the sensors upon which it relies [START_REF] Trinks | Industrial furnaces[END_REF]. For temperature control of walking-beam slab reheating furnace, it is necessary to know the temperature distribution in the slabs, in the refractory furnace walls and in the atmosphere of the furnace. In practice, the use of thermal sensors is a challenging and costly task due to high temperature and harsh environment of the furnace. Thermocouples are frequently utilized for measuring wall temperatures and gas temperature in reheating furnaces.

For instance, several thermocouples measure gas temperatures at different positions of a zone, an average temperature is then used for zone temperature regulation. The method is relatively robust, and sufficiently accurate for most control purposes. In contrast to pyrometers, themocouples require thermal contact, may give significant time delays, and can only make local measurements.

On the other hand, thermocouples can be implemented in holes drilled at several locations in in a product. During the heating time, the thermocouples signals were recorded on an acquisition system travelling with the product inside the furnace. Because of the extreme thermal conditions in the furnace, the acquisition system is placed inside a water box to keep its temperature lower than 100 • C. The thermal box is also covered with additional insulation as it shown in Figure 2.5 to protect it from radiation. In the example, the dimensions of the instrumented bloom were 160 mm in width, 128 mm in thickness and 10 m in length. The cooling chamber was place at as shown in Figure 2.5. This acquisition system is utilized only for the validation of furnace model which needed exact measurement of slab temperature through the length of the furnace.

Energy flows of reheating furnace

As previously mentioned, in a hot rolling mill, approximately 80 percent of energy consumption is due to reheating furnaces [START_REF] Steinböck | Model-based control and optimization of a continuous slab reheating furnace[END_REF]. Therefore, energy efficiency is a key performance value that directly affects the operational costs. The most significant energy flows of a walking-beam reheating furnace are shown in the Figure 2.6. The heat released by gas combustion is distributed to all furnace interior. A considerable portion of the total energy is stored in hot exhaust gas that leaves the furnace. Energy losses are minimized by a recuperation unit that reuses hot exhaust gas for preheating the combustion air and gas. Effectively, the recuperation process cannot totally recirculate the heat back into furnace, which means that some thermal energy will always be lost. In some cases, the heat losses to the flue is observed as 30% of the total supplied energy [START_REF] Chen | Analysis on energy consumption and performance of reheating furnaces in a hot strip mill[END_REF][START_REF] Zhang | Optimization of combustion control based on fuzzy logic[END_REF]. Without the recuperation unit, heat losses due to flue gas may increase to 50% [START_REF] Steinböck | Model-based control and optimization of a continuous slab reheating furnace[END_REF].

The other losses of energy are caused by door openings or expansion of walking beam joints. In addition, refractory walls and roof conduct heat from furnace interior to the outer surface. Heat is subsequently emitted to the external environment by radiation and convection. Energy that heats walking beams or skids may account for more than 10% of the total energy losses [START_REF] Chen | Analysis on energy consumption and performance of reheating furnaces in a hot strip mill[END_REF]. Cooling system equipped for skids or walking beams is connected to the heat recovery system improving the overall energy efficiency of furnace.

The heat absorbed by the slabs is considered as helpful energy. Additionally, its absolute value may be predetermined by the material properties, the initial temperatures, and the desired discharge temperatures of slabs.

According to [START_REF] Chen | Analysis on energy consumption and performance of reheating furnaces in a hot strip mill[END_REF], from 40% to 60% energy of the heat supplied to a reheating furnace is effectively absorbed by the slabs. Energy losses can be reduced by insulation of the furnace, recuperation of flue gas, heat recovery from cooling water or steam. On the other hand, energy savings can be obtained by reducing consumption of raw material while volume of heated products remains the same. This method is associated with improved heating quality, e.g., reduced material loss through scale formation, decarburization, reduced steel melting, and reduced need for reworking.

Chapter 2. Furnace presentation 25 A better precision of final slab temperature may lead to energy savings. If the slabs temperature can be controlled more accurately, the desired final temperatures can be reduced because the safety margins between the controlled final slab temperature and the temperature at least required for rolling mill can be smaller.

Hot charging may constitute also a method for energy savings. Instead of storing slabs after continuous casting process, which cools slabs to ambient temperature, they may be immediately transported to the reheating furnace. This method can save energy and operational costs. It reduces residence time of slabs, gas combustion and therefore flue gas. Higher furnace throughput, and increase furnace capacity is also obtained with hot charging compared to cold charging. However, integrated conveyance system is required to handle hot slabs from continuous casting stage. In addition, it happens practically that not all the slabs are loaded into reheating furnace at high temperature. Instead, a mixture of hot and cold slabs is fed to the furnace, which may make temperature control more complicated and more considerations will be needed.

After all, it is important to improve the strategy of temperature control, for which a model with high precision is necessary. As a result, the modeling of the furnace will be discussed in the next chapter.

Chapter 3

Modeling of slab walking-beam reheating furnaces

The presentation of the furnace has given an overview of its functional aspects and the main challenges for furnace control task. This chapter deals with furnace modeling. Firstly, a non exhaustive literature review of furnace models that are usually used is considered. Secondly, a steady-state model of the considered furnace describing the heat transfer phenomena inside furnace is presented. Finally, a dynamical model is introduced. This model is developed and tailored for use in transient phase of the furnace. The advantages and disadvantages of these two models are also discussed.

Literature review

A model of reheating furnace can be used for energy consumption analysis, estimation of physical quantities, simulation of control design, or analysis of furnace design. The specifications are detailed as follows [START_REF] Steinböck | Model-based control and optimization of a continuous slab reheating furnace[END_REF]:

• Energy consumption analysis of influence of insulation, oxygen enrichment, excess air coefficient, temperature of supplied combustion air and fuel, slab hot charging, distribution of heat input through burners along the furnace • Estimation of quantities that cannot be measured such as: slab profile temperatures, walls and roof temperatures, scale formation on slab surfaces, composition, temperature, and flow rates of flue gases.

• Simulation for control design which needs analysis on control strategies, controller gains, fuel consumption, influence of slab characteristics (dimensions, initial temperature, desired final temperature profile, material, etc.), and influence of residence time and time intervals between discharging events.

• Analysis of furnace design like: furnace geometry, shape and cross sectional area of flue gas ducts, shading effects, influence of refractory wall materials, location, size and type of burners, location and type of temperature sensors, modifications of existing furnaces.

Mathematical models are crucial for estimating or predicting slab temperature profiles by an observer algorithm or by a simulator. In order to compute slab temperatures, a model has to consider the heat transfer into surface of the slabs, and then the conduction through the slab thickness. Sometimes, temperature of slab surfaces can be measured by radiation pyrometry but there is still a need for a model that calculates the temperature distribution inside the slabs. In addition, mathematical models of furnaces can provide statistics about the gas flow, pressure inside the furnace, and energy consumption of the furnace [FRLD06a, GSR + 98a, HCK10a, HHLW08, KHK00, KKY + 00, ZIS97], excess oxygen ratio [KKY + 00, ZIS97], scale formation [JLKK10, Led86b, RSCG86, Sch96], and decarburization of the material [START_REF] Scholey | A transient, three-dimensional, thermal model of a billet reheating furnace[END_REF].

In real-time applications, computational constraints demand a reduction of mathematical furnace models to less expensive models while sufficient accuracy is still maintained [GSR + 98b, vDSK02, HVS04, KvD94, ESYS97, KKY + 00, Led86a]. Mathematical furnace models can be simplified by means of reduction techniques. However, for a highly nonlinear model, the application of reduction techniques is a rather challenging task.

In general, models of reheating furnaces are classified based on the integrated physical theory. White box models are derived and created from fundamental equations of chemistry and physics. Therefore, structure and parameters of white-box models have an explicit connections to the phenomena and quantities observed in the real system. On the contrary, black-box models or grey-box models utilize generic structures of system identification theory which do not allow direct physical interpretations [START_REF] Kusters | Mimo system identification of a slab reheating furnace[END_REF][START_REF] Van Ditzhuijzen | Identification and model predictive control of a slab reheating furnace[END_REF]. However, black-box models are able to simulate unknown disturbances which affect physical effects.

Essentiallly, there are three methods for modeling heat transfer and the temperature distribution inside the furnace enclosure [JLKK10] [START_REF] Kim | Three-dimensional analysis of the walking-beam-type slab reheating furnace in hot strip mills[END_REF]. First, the furnace enclosure is divided into distinct zones which are assumed to be a stirred vessel [ESYS97, GSR + 98a, HBK09, KMP94, KKY + 00, KvD94, PW98, YNY91, ZCX + 02]. Second, the furnace volume may be considered as a tube with a 1-dimensional plug flow in the longitudinal direction of the furnaces [START_REF] Steinböck | Model-based control and optimization of a continuous slab reheating furnace[END_REF]. Then the furnace temperature profile is determined along this direction [BBG08, GSR + 98b, ]. An alternative approach is adopted in [START_REF] Yoshitani | Optimal slab heating control with temperature trajectory optimization[END_REF], such as the method of weighted residuals with up to 5 polynomial trial functions. Some works account for temperature dependence of material parameters [WK99, SWKK11, SWKK10, Led86a, JLKK10, JK07].

Most of the models reported in the literature are based on the Fourier's law for heat conduction problem inside the slabs. Nevertheless, there are different approaches of reflecting the heat exchange outside the slabs. Three physical principles contribute to the heat exchange between the slabs and the environment:

• Heat conduction from skids and hearths • Convective heat transfer from flue gas • Radiative heat transfer from flue gas and furnace walls Because a furnace operates in relatively high temperature, radiative heat transfer is usually the dominant mode of heat exchange [HCK10a, HHLW08, JK07, Kim07, KKY + 00, Mår92, ZIS97, Li86]. The radiative effect of reheating furnace is characterized by Stefan-Boltzmann radiation law, which contains 4 th power of gas and surface temperatures. In the first category, furnace models do not take into account the radiative interaction between the flue gas and the surrounding surfaces. The low accuracy of these models can be reduced if local temperature of flue gases and furnace walls are averaged. The resulting mean temperatures are sometimes referred to as zone temperatures. Methods that compute the important effect of gas radiation like the net radiation method [START_REF] Frank | Fundamentals of heat and mass transfer[END_REF] and the zone method [START_REF] Noble | The zone method: explicit matrix relations for total exchange areas[END_REF], constitute the second category. Some models neglect the radiative heat exchange in the longitudinal direction of the furnace [vDSK02, YUU94, WCGS99, ESYS97, HVS04, JK07, KKY + 00]. Then, if the gas flow inside the furnace is not taken into account, there is no themodynamic interaction between slabs. This approach leads to particularly simple mathematical formulations as presented in Section 3.2.2. Linearized dynamic model is also a possibility of furnace modeling, it could be sufficiently accurate for the respective purpose [vDSK02, ESYS97, KKY + 00, KvD94, WK99]. However, it shows some limits since the 4 th power non-linearity associated with radiative heat transfer is significant. In addition, the temperature of slabs during heating process varies in a wide range, which means that there is no steady operating point for linearization.

A more elaborate furnace model is constructed if a full energy balance supplements the computation of radiative heat fluxes. Energy balance method helps evaluate the system in terms of energy consumption and efficiency As discussed in [START_REF] Laurinen | An adaptive neural network model for predicting the post roughing mill temperature of steel slabs in the reheating furnace[END_REF], the parameterization of a model can be a laborious task with numerous measurement series. However, it is a necessity to have accurate estimation of slab temperature, and model parameterizing is generally rewarded once achieved with the capacity of predicting online tracking quantities that cannot be measured due to consequent costs. This argument especially applies to the slab temperature profiles. The effort put into parameterizing physical models may be reduced if a stochastic model is used as in [START_REF] Honner | Exodus stochastic method application in the continuous reheating furnace control system[END_REF].

Black box methods can avoid physical modeling as used in [KKY + 00, KvD94].

Step response experiments are carried out on the real furnace system to obtain data for non parametric system identification with a discrete-time autoregressive model with exogenous input enhanced by time-lag behavior. Furnace models are utilized as comprehensive model and reduced model. The comprehensive model accounts for detailed physical effects and therefore captures the system dynamic as accurately as possible. This model is intended for furnace simulation and estimation of physical quantities that cannot be measured in the real system.

Inputs to this model can be flow rates of fuel and combustion air. Possible outputs are energy flows within the furnace, temperature profiles of the slabs, of furnace walls, as well as composition, temperature, velocity of flue gases. Comprehensive model might constitute computational burden. Therefore, reduced model may be constructed by neglecting convective heat transfer and simplifying furnace operation conditions. Optimization and control applications usually need reduced model to fit real-time operation requirements.

Nowadays, fast development of computer hardware offers more and more computation power as well as data storage for industrial applications. Therefore, in our case, detailed furnace models can be implemented providing high precision estimation or prediction of slab temperatures and other furnace variables. Furthermore, based on these detailed models, advanced control technologies can be designed for furnace control and optimization to improve furnace performance. In the following section, a steady-state model which is currently implemented in the considered furnace will be described. After that, a dynamical model is presented, this model is utilized in the design of a new furnace controller to increase the furnace performance.

Steady-state model

For the development of a model, the geometry of the furnace is simplified as described in Figure 3.1. The upper and lower half of the furnace are separated by the slab row. The furnace enclosure is partitioned into N z zones, this partition serves as a discretization in furnace length. The slab i last is the last loaded slab into furnace, and the slab i f irst will be the next slab discharged from the furnace.

In the frame reference of the furnace, the center of a slab i has the current z-position z i as indicated in Figure 3.2. In general, slabs are moved in positive z-direction which means that z i is a monotonously non-decreasing T i ([x, y, z], t) > 0 is the absolute temperature profile of slab i. The terms q + i (t) and q - i (t) are heat fluxes received by the upper and lower surface of the furnace, respectively. Heat flux q i ([x, y, z], t)(W.m -2 ) inside the slab depends on specific heat capacity c i (J.kg -1 .K -1 ), thermal conductivity λ i (W.m -1 K -1 ), temperature gradient, and boundary conditions at the surface of the slab. The parameters of a steel grade c i and λ i are usually temperature-dependent variables.

Unsteady one-dimensional heat conduction inside slabs

Heat conduction problem inside slab follows Fourier's law. Let ▽ be the gradient operator with respect to the spatial coordinates [x, y, z]. Knowing that heat flux is a vector quantity, the conduction rate equation from Fourier's law is expressed as follows:

q i ([x, y, z], t) = -λ i ▽ T i (3.1a) = -λ i (u x ∂T i ∂x + u y ∂T i ∂y + u z ∂T i ∂z ), (3.1b)
with u x , u y , u z are unit vectors of three axes OX, OY, OZ, respectively. The thermal conduction inside slab is therefore described by heat diffusion law as presented in [START_REF] Frank | Fundamentals of heat and mass transfer[END_REF]:

ρ i c i ∂T i ∂t = -▽ q i (3.2a) = ▽[λ i ▽ T i ], (3.2b)
where ρ i represents the mass density of the slab material.

Since the slab dimension ratios (cf.Figure 3.2) are frequently characterized by

L i ≫ D i and W i ≫ D i , therefore T i ([x, y, z], t
) is reasonably independent from both x and z, and therefore can be approximated by

T i (y, t) with y ∈ [-D i /2, D i /2].
Consequently, the heat flux q i is considered as function of y and time t, and is determined by one-dimensional version of Fourier's law through the y-direction, which is expressed as follows:

q i (y, t) = -λ i ∂T i (y, t) ∂y . ( 3.3)
In this cas, Equation 3.2 is simplified to:

ρ i c i ∂T i ∂t = ∂ ∂y (λ i ∂T i ∂y ). (3.4)
The initial conditions concerns temperature profile of the slab when it is loaded into the furnace, which is described as:

T i (y, t i,0 ) = T i,0 (y), (3.5
)

y ∈ [-D i /2, D i /2] ,
where t i,0 is time when slab i is loaded into the furnace and T i,0 (y) is temperature at position y inside the slab when it is charged into the furnace.

The boundary conditions of Equation 3.3 is given as:

q + i (t) = -λ i ∂T i ∂y | y=D i /2 , (3.6a) q - i (t) = -λ i ∂T i ∂y | y=-D i /2 . (3.6b)
In summary, one-dimensional heat conduction problem inside slab i is now defined by one-dimensional heat diffusion law given in Equation 3.4, initial condition in Equation 3.6, and boundary condition described in Equation 3.6.

Kirchhoff transformation

In practice, specific heat capacity c i is assumed to be independent of temperature. However, the dependence of thermal conductivity λ i on temperature has to be considered. This leads to non-linearity of Equation 3.4 which can be solved by Kirchhoff transformation as described in the following:

We first introduce the Kirchhoff variable

θ(T i ) = ∫ T i 0 λ i (T )dT. (3.7)
It should be noted that if λ i (T ) is specified, Equation 3.7 gives the relationship between θ and T i . To transform Equation 3.4 in terms of the variable θ, dθ/dT i is obtained by:

dθ dT i = λ i (3.8)
And ∂T i /∂t and ∂T i /∂y are then constructed as:

∂T i ∂t = dT i dθ ∂θ ∂t = 1 λ i ∂θ ∂t , (3.9)
and

∂T i ∂y = dT i dθ ∂θ ∂y = 1 λ i ∂θ ∂y . ( 3.10)
Substituting Equations 3.9 and 3.10 into Equation 3.4, we have

∂ 2 θ ∂y 2 = 1 α i ∂θ ∂t , (3.11)
where α i is the thermal diffusivity of slab i, defined as

α i = λ i ρ i c i . (3.12)
Since ρ i and c i are temperature-dependent, the diffusivity α i is also a function of temperature. It shall be noted that in steady heat conduction problem, ∂θ/∂t = 0. Then, Equation 3.11 becomes linear. To continue the development, we introduce the temperature-dependent enthalpy H(T i ) which is defined by ∂H/∂t = c i ∂T i /∂t into Equation 3.4, then we have

ρ i ∂H ∂t = ∂(λ i (∂T i /∂y) ∂y . (3.13)
In order complete the Kirrchhoff transformation, the boundary conditions in Equation 3.6 and enthalpy in Equation 3.13 are also expressed in term of θ i . By substituting Equation 3.9 into Equation 3.6 and Equation 3.10 into 3.13, we have:

q + i (t) = ∂θ ∂y | y=D i /2 , (3.14) q - i (t) = ∂θ ∂y | y=-D i /2 ,
and

ρ i ∂H ∂t = ▽ 2 θ.
(3.15)

Finite difference method

The Equations 3.14 and 3.15 constitute the problem of one-dimensional heat conduction inside a slab of continuoustime and continuous space domains. However, a discrete representation is required because the model has to be implemented on a computer. The discretization is needed for both time and space domain. For instance, temperatures at 2n + 1 points through slab thickness will be computed. First, the sampling distance is △ y = DX = D i /2n, and as shown in Figure 3.3, the temperature at position y j = j× △ y inside the slab i will be denoted T i,j . Sampling instants t k , ∀k ∈ N are discrete time, and △ t k = t k+1 -t k defines the corresponding sampling period. The main assumptions are:

• At a fixed position z i of slab, the heat flow input q ± i (t) are constant within the period t

∈ [t k , t k+1 [
• Conductivity of steel is piesewise linear between two temperatures T i,j (t k ) and T i,j (t k+1 )

• where θ j represents the Kirchhoff variable corresponding to temperature T i,j . Similarly, second order difference is expressed as:

∂ 2 θ ∂y 2 | y j = θ j+1 + θ j-1 -2θ j (△y) 2 + O((△y) 2 ), (3.17)
when number of points is sufficiently large, i.e, △y is significantly small, O((△y) 2 ) can be neglected. Therefore, from Equation 3.17, the problem defined by Equations 3.14 and 3.15 is rewritten as

ρ i H j (t k+1 ) -H j (t k ) △ t k = (θ j+1 -θ j ) -(θ j -θ j-1 ) (△ y) 2 , ( 3.18)
where j ∈ [-n + 1, n -1] indicates the position inside slab and H j represents the enthalpy corresponding to the temperature T i,j . Thus, the evolution of enthalpy H j is defined by as follows:

H j (t k+1 ) = H j (t k ) + △ t k (θ j+1 + θ j-1 -2θ j ) ρ i (△ y) 2 . (3.19)
To compute boundary values of H j with j = ±n, from Equation 3.16, Equation 3.19 is rewritten as

H j (t k+1 ) = H j (t k ) + △t k (2(θ j+1 -θ j ) -(θ j+1 -θ j-1 )) ρ i (△y) 2 (3.20a) = H j (t k ) + 2△t k ( (θ j+1 -θ j ) -△y ∂θ i ∂y | y j ) ρ i (△y) 2 . (3.20b)
Substituting Equation 3.14 into Equation 3.20 The heat flows received by slab surfaces q ± i (t k ) are subsequently determined by steady-state radiative heat transfer model based on position of slabs inside furnace.

H -n (t k+1 ) = H -n (t k ) + 2△t k ( (θ -n+1 -θ -n ) -△yq - i (t k ) ) ρ i (△y) 2 , (3.21a) H n (t k+1 ) = H n (t k ) + 2△t k ( (θ n+1 -θ n ) -△yq + i (t k ) ) ρ i (△y) 2 . ( 3

Steady-state radiative heat transfer in the furnace

The heat fluxes q ± i (t k ) received by surfaces of slab i introduced in the previous section as boundary conditions are necessary to compute temperature profile of the slab, and they represent radiative thermal fluxes inside the furnace. This section presents a simple method for calculation of q ± i (t k ). In general, radiative heat transfer couples the dynamic of all furnace zones, and all slabs inside the furnace because radiation is, in contrast to heat conduction, not just a local phenomenon. Inside the furnace, the heat exchange through radiation is determined by:

• The areas, shape, and position of the participating surfaces such as: walls, roof, skids, slabs • Temperature of surfaces that participate to the radiation • Emissive, absorptive, and reflective properties of the surfaces • The transmissive properties of the medium residing between the surfaces, which is, for instance, the exhaust gas circulating inside the furnace.

The general problem of radiative heat transfer is described with fundamental concepts. After that, the assumptions for steady-state operation of the furnace will be given. The radiative heat transfer is then computed based on the simplifying assumptions.

Assumptions

• Thermal radiation is the only considered mode of heat exchange between the slabs and their surrounding as shown in Figure 3.4. Moreover, in high temperature environment of the furnace, the dominant mode of heat exchange is radiation because radiative heat transfer is proportional to the difference of the 4 th powers of the temperatures whereas most other types of heat exchange depend only on the difference or the gradient of the temperatures themselves. Therefore, heat conduction between the slabs and the skids, or convection between the slabs and the flue gas, are considered negligible. However, heat convection will be examined in the next section.

• Only the bottom and the top surface of a slab serve as interfaces for heat exchange between the slab and its environment. Lateral, front, and back faces of the slab do not participate in the radiative heat exchange process. This assumption is inline with one-dimensional heat conduction inside the slabs as considered in Section 3.2.1

• For the computation of radiative heat exchange, the thickness of the slabs is not taken into account, i.e., all top surfaces of the slabs are assumed to be at the same level. The heat fluxes received by two surfaces are assumed to be equal:

q + i (t k ) = q - i (t k ) = q i (t k )
• At any time, the temperature of a participating surface is homogeneously distributed. This assumption seems adequate for slab surfaces. For furnace wall surface, however, its validity depends largely on the surface size, the distribution of burners, and the gaseous flow inside the furnace.

• All participating surfaces are ideal diffuse reflectors which means that an incoming ray of light is equally reflected in any direction over the hemisphere covering the surface.

• The furnace is divided into zones. It is assumed that the temperature and properties of the gas are homogeneous within each volume zone.

• When considering slab i residing in zone j, furnace wall and slab are combined into a closed space filled with flue gas as shown in Figure 3.5. That is all radiation exchanges are between furnace wall, gas of zone j, and slab i.

Computation of heat radiation

Let T g,j with ∀j ∈ 1, 2, ..., N z be the local gas temperatures which are usually measured by thermocouples, T s,i and T w,j are respectively the surface temperature of slab i and wall temperature of zone j. These temperatures are homogeneous as indicated with the assumptions. Radiative fluxes absorbed by slabs come mainly from furnace walls, and gas. Between slab surface and walls of the corresponding zone, the radiative heat flows are absorbed by the flue gas, continuously reflected and absorbed by solid surfaces. Radiative properties of walls and slabs are respectively characterized by their emissivities ε w , ε s , reflectivities

τ w = 1 -ε w , τ s = 1 -ε s , absorbtivities α w = ε w , α s = ε s .
Flue gas radiative characteristics are defined by its emissivity ε g and its absorbtivity α g . For instance, four radiative fluxes have to be considered in order to compute radiative heat flux absorbed by the slab as shown in Figure 3.4

• Heat flows received by slab i from wall: q w,i

• Heat radiation emitted by the gas to wall, reflected by the wall and finally absorbed by the slab i: q gw,i

• Heat flows emitted from the slab to wall and reflected back to the slab i: q sw,i

• Heat radiation emitted to slab i from flue gas: q g,i

The total heat flux absorbed by surfaces of a slabs is computed as

q + i = q - i = q i = q g,i + q w,i + q gw,i + q sw,i . (3.22)
Notice that the subscripts w, g, s stand for wall, gas and slab, respectively.

Radiation from furnace walls absorbed by a slab.

The radiation of furnace walls may be described by the Stefan-Boltzmann law applied for a grey body

q w = ε w σT 4 w,j . (3.23)
where ε w ∈ [0, 1] is the emissivity of furnace walls as previously mentioned, σ is the Stefan-Boltzmann constant.

In order to reach slab surface, radiative fluxes from furnace walls have to travel through the gas. Flue gas will absorb a portion of radiative energy, which is characterized by absorbtivity coefficient α g . Once the radiation has reached slab surface, it is generally reflected back through the gas and then comes back to furnace walls. The radiation reflection continues until its energy is attenuated. The reflectivity τ s on slab surface is characterized by its emissivity ε s by the relation τ s + ε s = 1. The heat flux absorbed by slab from the initial radiation (indicated by the subscript 1) is computed as

q w,i,1 = ε s (1 -α g )q w,0 = ε s A w , ( 3.24)
with A w = (1 -α g )q w,0 . Then, the reflection of the radiation to furnace walls through the gas is determined as

q * w,i,1 = (1 -ε s )A w . (3.25)
The reflective radiation is then absorbed partially by furnace walls and then reflected back to slab surface as

q * * w,i,1 = (1 -ε w )(1 -α g ) 2 q * w,i,1 = (1 -ε w )(1 -ε s )(1 -α g ) 2 A w (3.26)
The radiative heat absorbed by slab surface from the second radiation is calculated as

q w,i,2 = ε s q * * w,i,1 = ε s (1 -ε w )(1 -ε s )(1 -α g ) 2 A w . (3.27)
From Equations 3.24 and 3.27 , we can generalize the heat radiation of wall absorbed by the slab surfaces after n th time reflected back to slab surface as

q w,i,n = ε s [ (1 -ε w )(1 -ε s )(1 -α g ) 2 ] n-1 A w . (3.28)
The total radiative heat flux received by slab surface originated from furnace walls is

q w,i = n ∑ k=1 q w,i,k = ε s A w ( 1 + a + a 2 + ... + a n-1 ) = ε s A w 1 -a n 1 -a (3.29) where a = (1 -ε w )(1 -ε s )(1 -α g ) 2 ;
and it is obvious that a < 1. Therefore, when radiation reflects many time, i.e.,n → ∞, a n → 0. Then, Equation 3.29 can be simplified as

q w,i = ε s A w 1 1 -a = ε s (1 -α g )ε w σT 4 w,j 1 -(1 -ε w )(1 -ε s )(1 -α g ) 2 .
(3.30)

Direct radiation from gas to slab surface.

The gas is assumed to radiate directly on the slab surface like a normal solid object. Consequently, the initial heat flux coming from gas absorbed by the slab is computed as

q g,i,1 = ε s ε g σT 4 g,j . (3.31)
Following the Equation 3.28, the energy from the consequent reflective radiation on furnace walls after k time is computed as

q g,i,k = ε s [ (1 -ε w )(1 -ε s )(1 -α g ) 2 ] k-1 A g (3.32)
where

A g = ε g σT 4 g,j .
Similar to the Equations 3.30, the sum of direct radiation flux from gas to slab surface and its reflections can be calculated as

q g,i = n ∑ k=1 q g,i,k = ε s A g 1 1 -a = ε s ε g σT 4 g,j 1 -(1 -ε w )(1 -ε s )(1 -α g ) 2 .
(3.33)

Radiation from gas to furnace walls reflected to slab surface.

The direct radiation reflected by furnace walls from gas, and comes to slab surface is computed as

A gw,i = (1 -α g )(1 -ε w )ε g σT 4 g,j . (3.34)
Heat flux is then absorbed by the slab as the given formula

q gw,i,1 = ε s A g,w,i . (3.35)
Radiative heat flux after k reflections

q gw,i,k = ε s [ (1 -ε w )(1 -ε s )(1 -α g ) 2 ] k-1 A gw,i . (3.36)
The total radiative thermal flux absorbed by slab surfaces in this case is constructed as

q gw,i = ε s A gw,i 1 1 -a = ε s (1 -α g )(1 -ε w )ε g σT 4 g,j 1 -(1 -ε w )(1 -ε s )(1 -α g ) 2 .
(3.37)

Radiation issued from the slab, reflected on furnace walls, then back to slab surface.

The radiation flux arriving at furnace walls from slab surface across gas is constructed as

q i,w = (1 -α g )ε s σT 4 s,i . (3.38)
The radiation is then reflected by the walls, and partially absorbed by gas as:

q sw,i = (1 -ε w )(1 -α g ) 2 σT 4 s,i = A sws . (3.39)
Absorption by the slab of this first reflected radiation is calculated by

q s,i,1 = ε s A sws . (3.40)
The heat flux absorbed by slab surface after k time reflected on furnace walls is given as

q s,i,k = ε s [ (1 -ε w )(1 -ε s )(1 -α g ) 2 ] k-1 A s,w,s . (3.41)
Total heat flux absorbed by slab surface issued from its own radiation

q sw,i = ε s A sw,i 1 1 -a = ε s (1 -ε w )(1 -α g ) 2 σT 4 s,i 1 -(1 -ε w )(1 -ε s )(1 -α g ) 2 . (3.42)
In summary, solution of radiative heat transfer is determined by Equations 3.30, 3.33, 3.37, 3.42.

Dynamical model

Nowadays, furnace production has to respond rapidly to the client demand. significantly enhances the competitivity of the final product. New product with special characteristic requires more and more strict control on heating quality of the furnace. Therefore, the furnace operators would want to have a flexibility on production rate, being able to integrate an optimization of energy consumption and at the same time improve the heating process. For these reasons, furnace model needs to be able to monitor and predict gas, wall, slab temperature, energy consumption, etc., by taking into account available process variables of furnace conditions such as excess air, input power to each zone, or heat loss to skids.

In the following a dynamical model of the furnace is described

Heat conduction

In Section 3.2.1, one-dimensional heat conduction problem has been defined by Equations 3.4 and 3.6. The problem is then solved by introducing enthalpy variable specified by Equation 3.13, utilizing finite difference method and Kirchhoff transformation. In dynamical model of the furnace, finite volume method is used to construct onedimensional heat conduction model inside slabs. Consider again the unsteady one-dimensional heat conduction equation for a slab i assuming that specific heat capacity c i is temperature independent, we have

ρ i c i ∂T i ∂t = ∂ ∂y ( λ i (T i ) ∂T i ∂y ) . (3.43)
To find discretized equation, we use the grid-point cluster shown in Figure 3.7. We focus on the grid point P , which has grid points E and W as its neighbors. The dashed lines show the faces of a so-called control volume; Temperature profile of a slab at each time instant is defined by Equation 3.43 and its initial distribution of temperature. Thus, in a time step ∆t, the problem to be solved is: Given the grid-point values of T at time t, find the values of T at time t + ∆t. The given values of T at time t are denoted by T 0 P , T 0 E , T 0 W , and the new values at time t + ∆t by T 1 P , T 1 E , T 1 W . Now, discretization equation is described by integrating Equation 3.43 over the control volume shown in Figure 3.7, and over the interval from t to t + ∆t. Thereby,

ρ i c i ∫ e w ∫ t+∆t t ∂T i ∂t dtdy = ∫ t+∆t t ∫ e w ∂ ∂y ( λ i (T i ) ∂T i ∂y ) dydt, (3.44)
where the order of the integration is chosen according to the nature of the term. Assuming that the whole control volume has the same temperature with the grid point P , then ∂T /∂t = dT /dt. Then,

ρ i c i ∫ e w ∫ t+∆t t ∂T i ∂t dtdy = ρ i c i ∆y(T 1 P -T 0 P ). (3.45)
At each time instant, heat conduction inside slab is considered as in steady-state condition, that is ∂T /∂y = dT /dy; therefore,

ρ i c i ∆y(T 1 P -T 0 P ) = ∫ t+∆t t [ λ e (T E -T P ) (δy) e - λ w (T P -T W ) (δy) w ] dt. (3.46)
At this point, we need to make assumption about how T P , T E , T W evolve with time from t to t + ∆t. Many choices are possible, some of them can be given by the following formula

∫ t+∆t t T P dt = [ f T 1 P + (1 -f )T 0 P ] ∆t, (3.47)
where f is a weighting factor in [0, 1]. By assuming the same formulas for the integrals of T E and T W , Equation 3.46 becomes

ρ i c i ∆y ∆t (T 1 P -T 0 P ) = f [ λ e (T 1 E -T 1 P ) (δy) e - λ w (T 1 P -T 1 W ) (δy) w ] +(1-f ) [ λ e (T 0 E -T 0 P ) (δy) e - λ w (T 0 P -T 0 W ) (δy) w ] (3.48)
To obtain a more general form, we shall drop the superscript 1, and remember that T P , T E , T W from now on stand for the new values of T at time t + ∆t. The resulting equation is

a P T P = a E [ f T E + (1 -f )T 0 E ] + a W [ f T W S + (1 -f )T 0 W ] + [ a 0 P -(1 -f )a E -(1 -f )a W ] T 0 P , (3.49)
where

a E = λ e (δy) e , ( 3.50a
)

a W = λ w (δy) w , ( 3.50b
)

a 0 P = ρ i c i ∆y ∆t , ( 3.50c
)

a P = f a E + f a W + a 0 P . ( 3.50d)
The Equations 3.49 and 3.50 constitutes the discritization equation of one-dimensional conduction inside a slab i. At each time instant, temperature profile of a slab is determined if the temperature of lower and upper surfaces is known.

Now, an important question remains, that is: how the weighting factor f is chosen so that the discretization equation is physically realistic. The answer is f = 1 as well investigated in [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF]. As written in the Equation 3.47, f indicates how temperature evolve from T 0 P at time t to T 1 P at time t + ∆t. In Figure 3.8, explicit scheme corresponds to f = 0, that is the old value T 0 P endures throughout the entire time, except at time t + ∆t. Crank-Nicolson scheme defined by f = 0.5 assumes the linear variation of T P . Fully implicit scheme determined by f = 1 supposes that, at time t, T P suddenly jumps from T 0 P to T 1 P then stays at T 1 P over the whole time step. As reasoned in [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF], most realistic situations will be such that temperature at one point in the grid is influenced by the temperature of neighboring points only through process of conduction and diffusion. Then, decrease of temperature at one point, without changing of other conditions, leads to decrease (and not increase) of temperature at neighboring points. Applying this reasoning to the discretization equation 3.49, if T P decreases then T W and T E also decrease. It implies that the coefficients a P ,a W , and a E must have the same sign, and for that only f = 1 is satisfactory.

Radiative heat transfer

The furnace enclosure is considered as consisting of n grey and diffusing surface S i with emissivity ε i and homogeneous temperature T i as shown in Figure 3.9. The gas has a uniform temperature T g . These surfaces can be slab surfaces, wall surfaces, roof surfaces, etc,. The radiosity of the surface S i is the total heat flux coming out from S i . A part of the radiosity of surface S i is its own radiative emission S i ε i σT 4 i . A second part is the radiative fluxes transmitted through furnace gas, reflected by S i coming from all other surfaces of furnace enclosure including S i if it is concave. The final part is the heat flux of furnace gas emitted to S i and reflected on S i with a reflectivity τ i . The total radiative flux leaving S i is therefore given as follows:

S i J i = S i ε i σT 4 i + τ i n ∑ j=1 S j F j,i β j,i J j + τ i S i ε g,i σT 4 g . (3.51)
In this expression: β j,i is the total transmissivity of the gas separating surfaces S j and S i , F i,j and F j,i are the view factors between two surfaces, J i is the radiative flux of surface S i , ε g , i is the total emissivity of the gas to

S i .
Knowing that S j F j,i = S i F i,j , also τ i = 1 -ε i , and dividing both side of Equation 3.51 by S i we obtain the following equation:

J i = ε i σT 4 i + (1 -ε i ) n ∑ j=1 F i,j β j,i J j + (1 -ε i )ε g,i σT 4 g . (3.52)
Using Kronecker operator: δ i,j = 1 if i = j; δ i,j = 0 if i ̸ = j, we have the following form of Equation 3.52:

n ∑ j=1 [δ i,j -(1 -ε i )F i,j β j,i ]J j = ε i σT 4 i + (1 -ε i )ε g,i σT 4 g . (3.53)
This equation is only for the surfaces with imposed temperatures. For those with imposed flux, we have to establish another equation by using the net flux leaving S i . This latter net flux is equal to total flux leaving S i reduced with flux arriving on this surface, as given in the following equation:

q i = J i - n ∑ j=1 F i,j β j,i J j + (1 -ε i )ε g,i σT 4 g .
Substituting with Equation 3.53, we have the following form:

n ∑ j=1 [δ i,j -F i,j β j,i ]J j = q net,i + ε g,i σT 4 g . (3.54)
The radiative heat transfer inside furnace is then defined by the following system of equations: The resolution of system of equations in Table 3.1 provides the radiosity J i of a surface S i . The calculation of the unknown temperature T i or unknown net flux q i is defined by the following relation:

∑ n j=1 A i,j J j = B i , for i = 1 : n S i imposed with temperature T i imposed with net flux q net,i A i,j = δ i,j -(1 -ε i )F i,j β j,i δ i,j -F i,j β j,i B i = ε i σT 4 i + (1 -ε i )ε g,i σT 4 g q i + ε g,i σT
q i = ε i 1 -ε i (σT 4 i -J i ) (3.55)
For the furnace, different surface is defined based on division into zones. All slab surfaces participate to the overall heat radiation inside the furnace. In this dynamical model, local assumptions are eliminated.

Convective heat transfer

In continuous reheating furnace, typically less than 5% of the total heat input into the slabs comes from heat convection. This may explain why convection can be neglected in most of the models. However, if a model needs to be refined, convection is surely to be considered.

Fundamental differences between convection and radiation transfer which is dominant mode of heat exchange in a slab reheating furnace are:

• Convection can only occur between the exhaust gas and an adjacent surface of slabs, furnace walls/roof.

• Heat convection is due to the gas motion which may be generated by external source like pumps, fans, suction devices, or only by the density differences in the gas occurring due to temperature gradients.

• Convection is a local phenomenon, as opposed to radiation which transfers heat over large distances in the furnace Typically, the convective heat flux is modeled with linear dependence on the temperature difference. Thus, the convective heat flux into the surface of slab i which resides currently in zone j ∈ 1, ..., N z is computed as

q ± c,i (t) = h ± c,j ( T ± g,j (t) -T s,i (±D j /2, t) ) (3.56)
where T s,i (±D j /2, t), T ± g,j (t) are slab surface temperatures, and gas temperature, respectively. This expression is often known as Newton's law of cooling. The convection heat transfer coefficient h ± c,j is generally constant for a zone and is proportional to velocity of gas. The convective heat fluxes q ± c,i (t) are only added to the radiative heat fluxes.

Discussion

Now, we discuss the differences between steady-state model and dynamical model. About heat conduction problem, steady-state model utilizes the initial derivative form of conduction Equation 3.4. Finite difference method assumes the polynomial form of Kirchhoff variable and uses Taylor series expansion to approximate the temperature at a number of nodes inside slab. As investigated in [START_REF] Botte | Comparison of finite difference and control volume methods for solving differential equations[END_REF], this method is not mass conservative because the boundary conditions are not included in the governing equation. However, using a large number of nodes for the temperature profile can solve this problem.

The dynamical model, on the other hand, utilizes a integral form of heat conduction problem, which is the basic idea of finite volume method defined by Equation 3.44. The integral form assures the conservation of mass for any number of nodes in temperature profile. Therefore, with a profile consisting of small number of nodes, the finite volume method is superior to finite difference method. In fact, furnace control only needs a representative temperature profile of slabs. Thus, temperature estimation at several points through slab thickness shall be satisfactory. This validates the finite volume method and therefore validates the new approach to conduction problem adopted in the dynamical model.

In steady-state model, it is assumed that radiation in a zone is a local phenomenon. Only the radiative heat exchanges between slabs, exhaust fume, and wall of the zone are considered. The dynamical model considers heat radiation in the furnace as a global phenomenon, therefore heat radiation from zone to zone and slabs to slabs are taken into account. In practice, the developed dynamical model has been validated experimentally. In Figure 3.10, the measurement of slab average temperature trajectory is compared with simulated values. We see that the model is highly precise, which means it can be utilized to replace the measurement.

Chapter 4

Control system of walking beam reheating furnace General structure of reheating furnace was introduced in Chapter 1 and Chapter 2. Chapter 3 has described the currently used heat transfer model of the considered furnace. A dynamical model which is able to deal with unsteady-state furnace operation mode was also presented. This chapter will firstly give an overview of past works on the control structure of reheating furnaces. Then, the existing control strategy of the considered furnace will be presented. Finally, basic ideas of proposed control method and structure will be clarified.

Literature review

Main performance indicators of reheating furnaces are energy consumption, production rate, accuracy of the controlled temperatures, and operational costs. These indicators are influenced by construction factors, characteristic of charge, operational factors and control system of the furnace. Construction factors can be listed as: furnace geometry, refractory of furnace enclosure, measurement system, e.g., thermocouples, construction of heat sinks and heat sources, e.g., cooling system of skids and fuel burners, composition of fuel and limitation on fuel flow rates. Charge factors concern initial temperatures of slabs, desired final temperatures of slabs, material properties, sequence, geometry, geometric arrangement of slab inside the furnace, and metallurgical constraints on the slab temperature trajectories. Operational and system control factors involve the following aspects: In reheating furnace control, slab temperatures are the most relevant process variables. These variables are controlled by varying furnace zone temperatures. Subsequently, zone temperatures are controlled by varying the fuel flow to the burners of the zones. In a particular case, production rate can be regulated to control slab temperature while furnace zone temperatures are maintained constant [START_REF] Knoop | Nonlinear pi-controller design for a continuous-flow furnace via continuous gain scheduling[END_REF]. Therefore, it may be distinguished between controllers that regulate only temperatures or just the movement speed of slabs through the furnace [START_REF] Knoop | Nonlinear pi-controller design for a continuous-flow furnace via continuous gain scheduling[END_REF] or both [START_REF] Balbis | Nonlinear predictive control of steel slab reheating furnace[END_REF]. In this part of the work, like in most control strategies, only the temperature is controlled, the scheduling and so velocity of the slabs are forecast parameters. It shall be noted that temperature control can refer to slab temperature or zone temperature control or both.

Most furnace control strategies reported in the literature belong to one of the following strategies:

• Open-loop control • State-feedback control • Two-degree-of-freedom control

• Control strategies based on mathematical optimization

Open-loop control schemes are only suitable when furnace operations are stable or do not have considerable uncertainties. The references [START_REF] Ezure | Development of a simulator to calculate an optimal slab heating pattern for reheat furnaces[END_REF][START_REF] Chen | Modification of reheat furnace practices through comprehensive process modeling[END_REF] report on open-loop control of slab temperature by adjusting zone temperatures, i. e., there is no feedback from slab temperatures whether they are measured or estimated.

The set-points of zone temperatures are sent to subordinate control device, which regulates the supply of fuel and combustion air. The major drawback of this type of control is the lack of feedback [START_REF] Schurko | Computer control of reheat furnace: A comparison of strategies and application[END_REF]. In [START_REF] Hollander | Design, development and performance of online computer control in 3-zone reheating furnace[END_REF], the In two-degree-of-freedom control structure, the output of feedback controllers are added to optimal preplanned zone temperature trajectories [DBS88, Led86a, Pri80, RB92, Sta04, SWK11, VJKM08, WCGS99, WWC04] or optimal preplanned fuel flow rate [DCM + 92, Ves86]. Here, the motivation is to compensate for remaining deviations between the real furnace and the mathematical model used in feedforward control. Usually, the feedback law is based on the control error of the slab temperature trajectory and therefore requires the slab temperature to be estimated by some observer. For control of zone temperatures, a nonlinear feedback law can be developed and guarantees asymptotic stability [START_REF] Pedersen | On the reheat furnace control problem[END_REF]. Before looking into model predictive control strategy, we will now discuss on some particular aspects of furnace control.

In case where steady-state furnace models are applied for non-steady-state operating conditions like diversity of slabs, averaging methods are utilized. To obtain the control inputs, these methods use the averaged quantities of all considered slabs that are situating currently in the same zone. A feedback control law may be applied either before or after averaging. Thus, quantities that can be averaged are desired final values of slab temperatures, control errors of slab temperatures, heat inputs demanded by individual slabs, furnace temperatures demanded by individual slabs, fuel flow rates demanded by individual slabs, etc. In some cases, only one representative slab may be used, the approach is considered as weighted averaging technique. Prior knowledge of the future configuration and performance objectives of the furnace about slab sequencing, desired final temperatures, etc., can be utilized by some control strategies to improve furnace operating efficiency.

Among these control strategies are many optimization methods. Optimization techniques are good choices for furnace control because they can cope well with situations where the number of control objectives or controlled variables (slab temperature profiles) significantly exceeds the number of controllable inputs (furnace zone temperatures or fuel supplies). Furthermore, these techniques allow an explicit consideration of integral control objectives like minimum total energy consumption or maximum throughput. The following classes of optimization strategies are relevant to slab reheating furnaces:

• Static optimization under steady-state operating conditions • Heuristic improvements based on results of simulation models

• Dynamic optimization and model predictive control or also called receding horizon control In furnace control case where steady-state conditions are assumed, static optimization is usually efficient and especially less computational demanding than dynamic optimization. Therefore, highly detailed and accurate Chapter 4. Control system of the furnace 55 nonlinear models can be utilized in this case. The results for different operating scenarios may be stored in databases and commanded to the furnace upon demand. In published studies, optimal slab temperature trajectories can be derived as presented in [ESYS97, Ste11, FBF + 83], whereas optimum zone temperature set-points can be calculated as provided in [CA08, RB92, WCGS99, WWC04, Wö, YL88]. Fuel flow rates can also be used as input variables for optimization of steady-state furnace operation [START_REF] Icev | Supervisory-plus-regulatory control design for efficient operation of industrial furnaces[END_REF]. A general method of optimizing the distribution of local fuel inputs and the shape of counterflow reheating furnaces is presented in [START_REF] Kang | Constructal architecture for heating a stream by convection[END_REF].

An example of heuristic algorithms for furnace control can be found in [START_REF] Yang | Development of a computer control model for slab reheating furnaces[END_REF]. In general, heuristic methods are undemanding in terms of development and computational costs. These criteria are also important for the fuzzy logic control approach discussed in [START_REF] Zhang | Optimization of combustion control based on fuzzy logic[END_REF]. Optimal slab reheating trajectories are empirically designed based on a set of rules obtained from the operators' expertise. Fuzzy logic control is used for defining setpoint values of furnace zone temperatures depending on control errors of slab temperatures [START_REF] Keplinger | Automation system for reheating furnaces at hutta sendzimira[END_REF]. Another heuristic furnace control system based on genetic algorithm is proposed in [LGC].

On the contrary, dynamic optimization methods usually demand more computational power and can be executed offline or online. An offline preplanning of furnace fuel flow rates are published in [START_REF] Pike | Optimization studies of a slab reheating furnace[END_REF]. An optimal timevariant feedback law suitable for online implementation and direct computation of fuel flow rates is presented in [START_REF] Banerjee | A methodology to control direct-fired furnaces[END_REF]. These optimization methods need estimations of slab temperatures based on mathematical models which take fuel flow rates as inputs. The furnace throughput rate and fuel flow rates can be simultaneously optimized by means of mixed integer linear programming as provided in [FUK + 10].

If optimization algorithms are implemented online, the furnace models must be kept at an average level in term of computation demand. Therefore, the complexity of underlying furnace models is a constraint for control problem.

In [MPS04, SII + 99, Ste11], a quadratic objective function is minimized to find optimal furnace zone temperature set-points. This objective function contains weighted average of deviations between the realized and the desired heat inputs into the slabs. The heat inputs are calculated based on the desired slab temperatures at the end of each furnace zone.

A notable control strategy that uses dynamic optimization in closed loop is model predictive control. This method has been applied to slab reheating furnaces as reported in [BBG08, CXZ + 03, IZS04, KKY + 00, NvWS, ALRS05, YNY91, YUU94, Ste11, SWK13].

A hierarchical cascaded control structure where model predictive control is used in each loop is proposed in [START_REF] Balbis | Nonlinear predictive control of steel slab reheating furnace[END_REF]. The primary control loop calculates optimal furnace zone temperatures assuming a nonlinear furnace model subject to steady-state conditions. A linear model predictive control routine which is implemented in the secondary loop, provides dynamic optimization of both slab advancing rate and fuel supplies. Similar model predictive control system is suggested in [KKY + 00]. Fuel flow rates are controlled based on linearized MIMO models of the furnace zone temperatures. Energy consumption and furnace throughput can be both optimized at the same time with model predictive control strategies presented in [START_REF] Auzinger | Process for the optimized operation of a reheating furnace[END_REF]. There, the desired final slab temperature is only accounted as constraints of the optimization problem. A highly elaborated dynamic optimization scheme for a slab reheating furnace is reported in [START_REF] Yoshitani | Optimal slab heating control for reheating furnaces[END_REF][START_REF] Yoshitani | Optimal slab heating control with temperature trajectory optimization[END_REF].

Initially, a nonlinear optimization problem is solved to obtain reference trajectories of slab temperatures. They are designed as piece-wise polynomial functions, which allow the formulation of a low-dimensional and unbounded optimization problem. Moreover, the method of weighted residuals for computing the one-dimensional temperature profile inside the slabs is used to keep an acceptable complexity level of the optimization problem.

A second optimization problem based on a linearized model is finally solved for selecting fuel flow rates which ensure that the slabs follow their respective optimized temperature trajectories. Here, a generalized predictive control method suggested in [START_REF] Clarke | Generalized predictive control-part i. the basic algorithm[END_REF] is applied. Both optimization procedures are executed in real time.

Recent development on model predictive control for slab reheating furnace is reported in [START_REF] Steinböck | Model-based control and optimization of a continuous slab reheating furnace[END_REF][START_REF] Steinboeck | Nonlinear model predictive control of a continuous slab reheating furnace[END_REF]. In these publications, a nonlinear model predictive controller is designed based on a first-principle mathematical model, the controller defines local furnace temperatures so that the slabs reach their desired final temperatures.

The controller is able to cope with non-steady-state operating conditions of the furnace and reaches the desired slab temperature profiles. The nonlinear unconstrained dynamic optimization problem is solved by the quasi-Newton method.

In [START_REF] Zhang | Networked model predictive control based on neighbourhood optimization for serially connected large-scale processes[END_REF], a distributed model predictive control based on a neighborhood optimization is demonstrated for slab reheating furnace case. The furnace zones are considered to be serially connected (cascaded processes), that is, interactions between furnace zones are taken into account for furnace control.

In terms of algorithms, there are several standard methods for solving mathematically formulated optimization problems. The best choice among these algorithms for the general reheating furnace problem remains absent.

Actually, it depends on the available computing power and the mathematical structure of both the model and the

Existing and proposed control design of the considered furnace

In the considered reheating furnace, temperature control is designed based on hierarchical cascaded structure as shown in Figure 4.2. The scheduling and calculation of target temperature for each slab are at level 3 of the system. At this level, the slab requirements of hot rolling mill corresponding to client demand vary in terms A heating trajectory of slab temperature is then constructed offline for each slab according also to the available heating power of each furnace zone, see Figure 4.1. Based on this temperature trajectory, target temperatures of slab at the end of controllable zones: preheating, heating, and soaking zones are predetermined. Data of slab scheduling, slab target temperature at the end of each zone is commanded to level 2 where setpoints of furnace zone temperatures are calculated. At this level, the steady-state model described in Chapter 3 is used to predict temperature of a slab at the end of its current zone by imposing a value of temperature set-point for the zone. Similarly, temperature set-point of the zone is calculated for all other slabs currently being inside the zone. Each slab inside a zone will demand a different zone temperature set-point. Therefore, a weighted averaging technique is applied. Usually, for a time instant k, slabs that situate at the end or entry of a zone will In furnace operation, it is generally better to have rather slightly overheated slabs than under-heated slabs because under-heating of slabs causes incidents in hot rolling steps, which costs a lot in terms of maintenance and production delays. However, overheating has also to be restricted because it can lead to slab melting which also stops the heating process. Moreover, overheating generates more scale formation on slab surfaces than usual, therefore more material is lost. Waste of energy is obviously undesirable effect caused by overheating of slabs.

As a result, slab overheating significantly contributes to the price of final steel products.

In the current situation of the considered furnace, it happens that slabs have a large variety of dimensions, target temperature, and material due to diversity of client demands. Therefore, between two consecutive slabs, there can be considerable difference between desired final temperature. In these cases, the furnace control will impose furnace temperatures that will successfully heat the slab that has higher desired final temperature. Hence, in average, slabs are overheated. The value representing this overheating tendency is a performance criteria of the furnace which reflects also the accuracy of furnace control. A new control strategy which gives higher accuracy will decrease the gap between discharge and desired final temperature and therefore reduce energy consumption of the furnace.

Level 1 of the furnace control system tracks furnace zone temperatures according to their temperature set-points which are calculated by the level 2. Power set-points for each zone will be computed and be sent to lower level of the control system. In the considered furnace, PI controllers are used to control zone temperatures as shown in Figure 4.2. The fact that furnace zones strongly interact with each other influences considerably the performance of PI controllers. Therefore, new controllers that take into account the zone interactions will improve tracking performance of zone temperatures and also reduce energy consumption of the furnace.

The proposed solution for furnace control bases on the same structure as for the current control system as shown in Figure 4.4. At level 2, a nonlinear model predictive controller is designed. The dynamical furnace model presented in Chapter 3 is expected to cope with non-steady-state furnace operations and gives more accurate estimations of slab temperature profiles. This model is also utilized to predict furnace quantities from a sequence As described in Chapter 3, the dynamical model of the furnace is a numerical model which uses finite volume method to calculate slab temperatures. Therefore, objective function cannot be written in analytic form of control variables which are temperature set-points of furnace zones. Especially, partial derivatives of slab temperature will not be calculated directly. Even estimation of the derivatives is expensive in term of computational resources. Therefore, optimization algorithm will have to satisfy time constraint of the problem. Here, direct optimization algorithms are investigated to find a suitable one for the reheating furnace control problem.

For the level 1, a distributed model predictive control method is proposed to replace PID controllers as shown in Figure 4.4. This method benefits from knowledge about the dynamic and coupling of furnace zones to calculate optimal input power to each zone. An iterative algorithm is constructed based on a output coupling model of furnace zones. The optimization problem is solved by quadratic programming.

All these control architecture of levels 1 and 2 will be further described in the next chapters.

Chapter 5

Zone temperature control

Control system of the considered furnace has a cascaded hierarchical structure as described in chapter 4. At level 1, the tracking of zone temperature set-points is currently ensured by PID controllers implemented in decentralized scheme. In reheating furnaces, however, furnace zones are significantly coupled with each other. The decentralized PID controllers cannot therefore explicitly tackle these interactions. In this chapter, distributed model predictive controllers (distributed MPCs) are designed to replace decentralized PID controllers, based on furnace zone models. First, an overview on distributed MPC is given. Second, the modeling of furnace zones and their interactions will be described. Identification procedures are detailed for the considered furnace zones.

Overview on distributed MPC

MPC is an advanced control method granted with a great interest as much in theoretical studies as in industrial application practice. The basics of MPC is utilization of a process model to predict output variables over a finite horizon. Especially, MPC is able to handle process constraints, to deal with multivariable processes. The control problem is transformed to an optimization one. And, optimal control inputs over a finite horizon are derived by solving an optimization problem which usually represents an economical criterion. Principles of MPC applied to the furnace zone temperature control are illustrated in Figure 5.1. At each time instant k, an optimization problem is formulated accounting, e.g. in our specific case, for zone temperature tracking errors and furnace energy consumption over a finite horizon. Optimal input trajectory for time k is derived by solving this optimization problem. Only the first action of the optimal trajectory is employed. The procedure is repeated for the next sample time k + 1. Nowadays, with the fast growth of computational power, process modeling becomes more and more affordable for real-time applications. As a result, MPC method is widely applied to industrial processes replacing the traditional PID controllers, especially in economically important units [START_REF] Rawlings | Coordinating multiple optimization-based controllers: New opportunities and challenges[END_REF].

Reheating furnaces are networked large-scale engineering systems that consist of a number of subsystems interacting with each other through flows of material, energy, and information. The subsystems are represented by furnace zones, flow of material is the steel slab charge moving through the furnace. Information from thermocouples, pressure sensors, etc., allows the monitoring and control of the furnace. Energy flow inside furnace is mainly the heat supplied from fuel combustion. local feedback action provided by local controllers may be sufficient to compensate the coupling effects. Therefore, decentralized MPC strategy is expected to work adequately. For many plants where interactions among subsystems are strong, ignoring these interactions leads to significant loss of control performance. In reheating furnaces, steel slabs are heated in high temperature zones. The radiation is a global and dominant heat transfer phenomenon as reasoning in chapter 3. In addition, circulation of high temperature gas also makes couplings between furnace zones considerable. For these reasons, decentralized MPC will not be suitable for the control of zone temperatures.

In contrast to decentralized MPC, centralized MPC structure uses a MIMO model for large-scale systems, therefore takes into account every coupling effect between sub-systems, see Figure 5.3. However, centralized controllers, are viewed by practitioners as inflexible and monolithic solution. For most large-scale processes, and so reheating furnaces, the primary burdens to centralized control are not computational but organizational [START_REF] Venkat | Distributed Model Predictive Control: Theory and Applications[END_REF].

On the one hand, operators are usually unwilling to deal with the substantial data collection and data handling effort required to design and maintain a valid centralized control system for a large plant. To the best of our knowledge, no such centralized control systems are operational today for any large, networked system. In fact, zones of furnace were individually controlled by PID controllers. And so, a decentralized structure was already in place and plant personnel does not wish to completely redesign the control system to implement centralized MPC. On the other hand, different zones of the furnace are subject to modifications through time, which makes the model development and maintenance effort required for centralized control impractical. Unless these organizational hurdles change in the future, centralized control of networked, large-scale systems stays useful mainly as a benchmark against which other control strategies can be compared and evaluated. On the other hand, a MPC agent can cooperate with others by optimizing the performance criterion of a set of subsystems [START_REF] Rawlings | Coordinating multiple optimization-based controllers: New opportunities and challenges[END_REF][START_REF] Venkat | Distributed Model Predictive Control: Theory and Applications[END_REF], or optimizes independently the performance criterion of one individual subsystem.

In term of optimality of control solution, in a fully cooperative algorithm where each local controller will optimize a global performance criteria, the overall system will converge to a so called Pareto equilibrium [CJKT02, AGBE11, MN14]. On the contrary, in a non-cooperative algorithm where each local controller tries to maximize only the performance of its corresponding subsystem, the overall system will come to the solution called Nash equilibrium [Nas51, ZL07, MN14]. To understand notions of Nash and Pareto equilibriums in an easy way, suppose that there are two controllers which will take action to optimize performance of two coupled subsystems: controller 1 will do action a 1 , controller 2 will do action a 2 , in order to minimize objective functions: J 1 (a 1 , a 2 )

and J 2 (a 1 , a 2 ). A set of actions ( a N 1 , a N 2 ) constitutes a Nash equilibrium iff: J 1 (a N 1 , a N 2 ) ≤ J 1 (a 1 , a N 2 ) for any a 1 and J 2 (a N 1 , a N 2 ) ≤ J 2 (a N 1 , a 2 )
for any a 2 . In other words, Nash equilibrium is where no controllers can act unilaterally to improve their cost knowing action of the other. While a set of actions ( a P 1 , a P

2

) is said to be Pareto optimal if there does not exist another set (a 1 , a 2 ) such as J 1 (a 1 , a 2 ) ≤ J 1 (a P 1 , a P 2 ) and J 2 (a 1 , a 2 ) ≤ J 2 (a P 1 , a P 2 ). Therefore, Pareto solution is the best solution that two controllers can obtain acting together. In term of communication among controllers, in case of highly interactive system-wide, it requires a huge resource for communication. For a less interactive system-wide the degree of communication is reduced. As a result, communication structure of distributed MPC varies according to coupling topology of the system-wide.

Therefore, two communication topologies can be defined for the design of distributed MPCs: 1) MPC agents are fully connected with each other and 2) MPC agents are partially connected in a neighborhood [START_REF] Rawlings | Coordinating multiple optimization-based controllers: New opportunities and challenges[END_REF].

In term of communication protocol, in platforms where hardware resource is abundant, the optimization and communication procedures can be executed in an iterative way to guarantee the optimal solution. Otherwise, non-iterative solution may still be sufficient for control requirements. To summarize different structures that distributed MPC could have, a classification of distributed MPC is given in Table 5.1.

In metallurgical processes, highly interactive system-wide is not the common case. Material, energy and information generally flow sequentially from sub-system to sub-system, which characterizes the major coupling effects.

In fact, the interactions inside a reheating furnace are limited in a neighborhood of a zone. Each furnace zone interacts mainly with its two neighboring zones. By taking advantage of the coupling topology, we can reduce the required communication between the MPC agents. In the following, modeling of furnace zones and their interactions will be presented.

Furnace zone modeling

Modeling of furnace zones concerns identifying dynamic of zone temperatures with respect to variation of input heat powers. Therefore, the modeling involves only furnace zones equipped with fuel burners which are: preheating zone, heating zones and soaking zone. For the considered furnace, heating powers are mainly installed at preheating and heating zones where slabs absorb most of the necessary heat. Two halves of heating zones are separately controlled as shown in Figure 5.5. Soaking zone provides the additional heat for slabs to homogenize Also by observation on real furnace, the dynamic of zone temperature with respect to variations of input power is similar to a first order system with delays and can be expressed as follows: where y j = ∆T j = T j -T j,nom the variation of zone temperature T j from its nominal value T j,nom , u j = ∆P j = P j -P j,nom the variation of zone input power P j from its nominal value P j,nom , L represents the operator of Laplace transformation, and s the Laplace variable.

H j (s) = L[y j (t)] L[u j (t)] = G jj e -τ jj s 1 + λ jj s , ( 5
Couplings among furnace zones are mainly due to convective heat transfer of flue gas and radiative heat transfer inside the furnace. The flue gas flows from exit to entry of the furnace in opposite direction of slab movement as depicted in Figure 5.5. As mentioned in furnace modeling chapter, convective heat transfer depends on difference between participating object's temperature, while radiation heat transfer is proportional to 4 th power of temperature difference of participating surfaces. Therefore, coupling between two neighboring zones j and i = j ± 1 can be modeled with respect to temperatures. This topology of interaction refers to an output coupling model.

The dynamic of a zone's temperature according to variation of its neighboring zone's temperature is observed to be similar to first order transfer function with delays and can be expressed as follows:

H ij (s) = L[y i (t)] L[y j (t)] = G ij e -τ ij s 1 + λ ij s . (5.2)
Interaction of two neighboring zones can also be modeled as input coupling. By observation, dynamic of a zone's temperature with respect to variation of its neighboring zone's input power is also similar to first order transfer function with as shown in Figures 5.8 and 5.7. The system can be represented by

H ij (s) = L[y i (t)] L[u j (t)] = G ij e -τ ij s 1 + λ ij s . ( 5.3)
In conclusion, around a nominal operating condition, furnace zones can be simplified to linear time invariant delayed first order systems. In the following, two models of the considered zones with two distinct coupling structures defined by Equations 5.2 and 5.3 are detailed for the design of distributed MPC strategy.

Input coupling model

Zone model based on input coupling assumption is depicted as shown in Figure 5.9. Three subsystems correspond to preheating, heating 1 and heating 2 zone, respectively. The overall model of three subsystems is written as follows:

Y 1 (s) = H 11 (s)U 1 (s) + H 12 (s)U 2 (s) + H 13 (s)U 3 (s), (5.4a
)

Y 2 (s) = H 21 (s)U 1 (s) + H 22 (s)U 2 (s) + H 23 (s)U 3 (s), (5.4b
)

Y 3 (s) = H 31 (s)U 1 (s) + H 32 (s)U 2 (s) + H 33 (s)U 3 (s), (5.4c)
with

H ii = H i , for i = 1, 2,
Y 1 (s) = H 11 (s)U 1 (s) + H 12 (s)U 2 (s), (5.5a
)

Y 2 (s) = H 21 (s)U 1 (s) + H 22 (s)U 2 (s) + H 23 (s)U 3 (s), (5.5b
)

Y 3 (s) = H 32 (s)U 2 (s) + H 33 (s)U 3 (s).
(5.5c) The input coupling model of the considered furnace has been identified in previous work reported in [DP, JR] as

Y 1 (s) = 0.0140e -59s 1 + 513.1s U 1 (s) + 0.0048e -60s 1 + 460s U 2 (s) (5.6a) Y 2 (s) = 0.0031a -70s 1 + 520s U 1 (s) + 0.0108e -60s 1 + 460s U 2 (s) + 0.005e -65s 1 + 533.3s U 3 (s) (5.6b) Y 3 (s) = 0.0038e -90s 1 + 438s U 2 (s) + 0.0125e -60s 1 + 500s U 3 (s) (5.6c)
In this coupling model of furnace zones, we see that temperature of a zone is more rapidly affected by the input power of that zone than by input power of other zones. Specifically, for the temperature of zone preheating (y 1 ), delay from u 1 to y 1 is 59 s whereas delay from u 2 to y 2 is 60 s. For temperature of the zone heating 1 (y 2 ), delay from u 2 to y 2 is 60 s while delays from u 1 and u 3 to y 2 are 70 s and 65 s, respectively. And for the zone heating 2, delay from u 3 to y 3 is 60 s when delay from u 2 to y 3 is 90 s which is significantly longer. In addition, the static gain from input power of a zone to temperature of that zone is more significant than that to temperature of neighboring zones. For example, in zone heating 1, static gain from u 2 to y 2 is 0.0108 while the gains from u 1 and u 3 to y 2 are only 0.0031 and 0.005, respectively.

Output coupling model

Zone model based on output coupling assumption is depicted as shown in Figure 5.10. Three subsystems correspond to preheating, heating 1 and heating 2 zone, respectively. The overall model of three subsystems is written as follows:

Y 1 (s) = H 11 (s)U 1 (s) + H 12 (s)Y 2 (s) + H 13 (s)Y 3 (s), (5.7a
)

Y 2 (s) = H 21 (s)Y 1 (s) + H 22 (s)U 2 (s) + H 23 (s)Y 3 (s), (5.7b
)

Y 3 (s) = H 31 (s)Y 1 (s) + H 32 (s)Y 2 (s) + H 33 (s)U 3 (s), (5.7c)
with H ii = H i , for i = 1 ÷ 3, is given in Equation 5.1, and H ij , j ̸ = i, is described in Equation 5.2. In identifi-Figure 5.10: Output coupling of furnace zones cation of output coupling model of furnace zones, some coupling terms are obsolete. For reheating furnaces, the flue gas circulates from furnace exit to furnace entry meaning from high temperature zones to low temperature zones. By that, the supplied thermal energy is redistributed to heat cold slabs. In addition, zone temperatures are affected directly by temperature of the flue gas. Therefore, the coupling effects happen mainly in the direction of circulation of the flue gas. As a result, the influence of zone with lower temperature to zone with higher temperature is neglected. Thus, in the Equation 5. 4 (5.8c)

The obtained output coupling model is given as:

Y 1 (s) = 0.0143e -60s 1 + 605s U 1 (s) + 0.210e -20s 1 + 500s Y 2 (s) (5.9a) Y 2 (s) = 0.012e -40s 1 + 742s U 2 (s) + 0.16e -20s 1 + 416s Y 3 (s) (5.9b) Y 3 (s) = 0.0066e -70s 1 + 415s U 3 (s) (5.9c)
These two models are validated against measurement with a different set of data of zone temperatures as shown in Figures 5.12 and 5.11. Generally, it can be seen that both identified models can capture the temperature dynamic of three considered zones. However, the output coupling model fits more the measurement for temperature of zone: heating 1 and heating 2. This can be explained by the fact that temperature of a zone is directly influenced by the temperatures of neighboring zones. For the next sections, output coupling model will be used for controller design.

Distributed MPC design

After examination of furnace zone models, we now generalize the furnace zones as subsystems which are coupled through outputs and distributed MPC strategies are constructed accordingly. First, the predictions of outputs is formulated based on a state space representation of system. Second, the control problem is formulated for two distinct strategies of distributed MPC: cooperative and non-cooperative. Additionally, iterative algorithms are designed assuring solutions with a low level of suboptimality. Also, convergence properties of iterative algorithms and nominal stability of closed-loop system are analyzed for an unconstrained case. Finally, the centralized MPC strategy is constructed and employed as a benchmark for two distributed MPC algorithms.

Networked system modeling

Assuming that the system-wide is composed of M subsystems coupled with each other by outputs, for the subsystem i, the state space representation can be written as

{ x i (k + 1) = A i x i (k) + B i u i (k) + ∑ M j=1,j̸ =i B ij y j (k) y i (k) = C i x i (k) (5.10)
This representation is based on a linearization around the operating point of the furnace. Vectors u i ∈ ℜ n i and y i ∈ ℜ m i are respectively the variations of input power, and temperature of zone i. The term B ij is the coupling coefficient of zone i with zone j. It should be noted that a zone i interacts only with its neighboring zones.

Consequently, B i,j = 0 if j ̸ = (i ± 1). The couple {C i , A i } is assumed to be observable for i = 1 ÷ M . The prediction of output y i (k + s) at instant k is denoted by ŷi (k + s | k). From the state space representation,

ŷi (k + s | k) is calculated as ŷi (k +s | k) = C i A s i x i (k)+ min(s,Nu) ∑ p=1 C i A s-p i B i u i (k +p-1 | k)+ ∑ j̸ =i s ∑ p=1 C i A s-p i B ij ŷj (k +p-1 | k), (5.11)
N u is the control horizon. In a vector form, the predictions of output over prediction horizon N p are gathered as follows:

Ŷi (k | k) = [ ŷ′ i (k + 1 | k), . . . , ŷ′ i (k + N p | k)
] ′ (5.12)

Similarly, the control input sequence over N u steps is put into a vector as

U i (k | k) = [ u ′ i (k | k), . . . , u ′ i (k + N u -1 | k) ] ′ . (5.13)
Output prediction of zone i in Equation 5.11 is therefore rewritten as

Ŷi (k | k) = Ψ i x i (k) + Φ i,i U i (k | k) + ∑ j̸ =i Φ i,j Ŷj (k | k), ( 5.14a)
where

Ψ i = [ (C i A i ) ′ . . . (C i A Np i ) ′ ] ′ , and
Φ i,i =           C i B i • • • 0 . . . . . . . . . C i A Nu-1 i B i • • • C i B i . . . . . . . . . C i A Np-1 i B i • • • C i A Np-Nu i B i           , Φ i,j =        C i B i,j 0 • • • 0 C i A i B i,j C i B i,j . . . . . . . . . . . . . . . 0
C i A Np-1 i B i,j C i A Np-2 i B i,j • • • C i B i,j        .

Observer design for output feedback control

In the considered reheating furnaces, available measurements are zone temperatures which are the controlled outputs. In order to design a controller, an observer is needed to estimate the states of subsystem from temperature measurements. The structure of observer is shown in Figure 5.13. In this case, pole-placement technique is employed by introducing state and output estimations: xi (k), ȳi (k), respectively. The structure of the observer

xi (k + 1) = A i xi (k) + B i u i (k) + L i [y i (k) -ȳi (k)] + ∑ j̸ =i B i,j y j (k) (5.15a) ȳi (k) = C i xi (k) (5.15b)
with L i is a feedback vector. State estimation error is defined as xi (k) = xi (k) -x i (k). From Equations 5.10 and 5.15 the dynamic of estimation error is given as

xi (k + 1) = xi (k + 1) -x i (k + 1) (5.16a) = A i xi (k) -A i x i (k) + L i (y i (k) -ȳi (k)) (5.16b) = A i xi (k) + L i (C i x i (k) -C i xi (k)) (5.16c) = (A i -L i C i )x i (k). (5.16d)
The estimate state xi (k) will converge to the real state x i (k) with L i that satisfies

eig (A i -L i C i ) < 1.
(5.17)

That is, the eigenvalues of (A i -L i C i ) are all inside the unit circle. The dynamic of the observer has to be faster than the dynamic of the subsystem to provide the information for the controller. In general, eigenvalues of

(A i -L i C i )
are chosen and imposed according to eigenvalues of A i , L i is then calculated subsequently.

Non-cooperative distributed MPC

Temperature tracking performance is defined with repects to specific criteria, in particular zone temperature accuracy and energy consumption. In a non-cooperative strategy, a local MPC agent optimizes only the performance of the corresponding zone while taking into account influences of other zones. Therefore, the local controller has to gather output predictions of other zones in order to make an accurate predictions of output of the corresponding zone. In return, the controller has to send out its output predictions so that other controllers can make their own predictions.

In this section, first, the formulation of objective function will be given, then followed by the description of an iterative algorithm. Afterwards, the convergence of the iterative algorithm and closed-loop stability for an unconstrained case are examined.

Problem formulation

Zone temperature control is subject to limitations of zone power and variation speed of zone power, and to upper/lower bounds of zone temperatures. Thus, optimization problem for local MPC agent i is defined as

min u i (k|k),...,u i (k+Nu-1|k) J i (k) = Np ∑ s=1 q i,s [ŷ i (k + s | k) -y i,ref (k + s)] 2 + Nu-1 ∑ p=0 r i,p+1 u 2 i (k + p | k), (5.18a) subject to u i,min ≤ u i (k + p) ≤ u i,max ; p = 0 ÷ N u -1, (5.18b
)

∆u i,min ≤ ∆u i (k + p) ≤ ∆u i,max ; p = 0 ÷ N u -1, (5.18c
)

y i,min ≤ ŷi (k + s | k) ≤ y i,max ; s = 1 ÷ N p , ( 5.18d)
where y i,ref (k + s) is the set-point value of the output of zone i at sample time (k + s) and ∆u i (k + p) is the speed of change of the control input defined by ∆u i (k + p) = u i (k + p) -u i (k + p -1). The coefficients q i,s and r i,p+1 are positive weighting terms on tracking errors and input powers, respectively. From Equations 5.13 and 5.12, the local objective function in Equation 5.18 can be rewritten in a vector form as min

U i (k|k) J i (k) = ∥ Ŷi (k | k) -Y i,ref (k)∥ 2 Q i + ∥U i (k | k)∥ 2 R i (5.19a) subject to a i U i (k) ≤ b i (k) (5.19b) with Q i = diag(q i,1 , ..., q i,Np ), R i = diag(r i,1 , ..., r i,Nu ), and Y i,ref (k) = [ y ′ i,ref (k + 1), ..., y ′ i,ref (k + N p )
] ′ .

a i and b i (k) are determined by the following

b i (k) =             U i,max -U i,min Y i,max -Ψ i x i (k) - ∑ j̸ =i Φ i,j Ŷj (k | k) -Y i,min + Ψ i x i (k) + ∑ j̸ =i Φ i,j Ŷj (k | k) ∆U i,max -∆U i,min             , a i =             I Nu×n i -I Nu×n i Φ i,i -Φ i,i Ξ i -Ξ i             , with U i,max = [ u ′ i,max , ..., u ′ i,max ] ′ , U i,min = [ u ′ i,min , ..., u ′ i,min ] ′ , Y i,max = [ y ′ i,max , ..., y ′ i,max ] ′ , Y i,min = [ y ′ i,min , ..., y ′ i,min ] ′ , ∆U i,max =        ∆u i,max + u i (k -1) ∆u i,max . . . ∆u i,max        , ∆U i,min =        ∆u i,min + u i (k -1) ∆u i,min . . . ∆u i,min       
, and

Ξ i =           I n i 0 • • • • • • 0 -I n i I n i . . . 0
• • • 0 -I n i I n i           Nu.n i ×Nu.n i .
Consequently, a quadratic programming problem is derived as min

U i (k|k) 1 2 U ′ i (k | k)H i U i (k | k) + f ′ i U i (k | k) (5.20a)
subject to (5.20b)

a i U i (k | k) ≤ b i (k), (5.20c)
where

H i = ( Φ ′ i,i Q i Φ i,i + R i ) , f i = Φ ′ i,i Q i ( Ψ i x i (k) + ∑ j̸ =i Φ i,j Ŷj (k -1 | k) -Y i,ref (k)
)

. At each sample time, local MPC controllers communicate the output predictions made at the previous sample. Based on this information, individual quadratic problems are solved providing the optimal input power sequence to each subsystem.

Iterative algorithm

To ensure optimal solution, iterative procedure can be used when communication resource is sufficient. At The local MPC controller i initializes its control input sequence U i (k | k) by the optimal one found at the previous sample like

U (l=0) i (k | k) = U i,opt (k -1 | k -1)
, and output predictions by those of previous sample

Ŷ (l=0) i (k | k) = Ŷi (k -1 | k -1)
, l is the iteration index within a sampling period

Ŷ (l) j (k | k) for j = 1 ÷ M , j ̸ = i. • Step 3: Optimization.
MPC controller i solves its local optimal problem and obtains the optimal control input sequence for the iteration l + 1.

U (l+1) i (k | k) = arg min U i (k|k) 1 2 U ′ i (k | k)H i U i (k | k) + (f (l) i ) ′ U i (k | k) (5.21a)
subject to (5.21b)

a i U i (k | k) ≤ b (l) i (k), (5.21c)
where f

(l) i = Φ ′ i,i Q i ( Ψ i x i (k) + ∑ j̸ =i Φ i,j Ŷ (l) j (k | k) -Y i,ref (k) ) , and
b (l) i (k) =             U i,max -U i,min Y i,max -Ψ i x i (k) - ∑ j̸ =i Φ i,j Ŷ (l) j (k | k) -Y i,min + Ψ i x i (k) + ∑ j̸ =i Φ i,j Ŷ (l) j (k | k) ∆U i,max -∆U i,min             .

• Step 4: Checking and Updating

The predicted outputs are updated with the obtained control input sequence as follows:

Ŷ (l+1) i (k | k) = Ψ i x i (k) + Φ i,i U (l+1) i (k | k) + ∑ j̸ =i Φ i,j Ŷ (l) j (k | k). (5.22) If ∥U (l+1) i (k | k) -U (l) i (k | k)∥ ≤ ϵ i for i = 1 ÷ M
, where ϵ i is the convergence precision, or the iteration index l = l max then MPC controllers end the iterations, otherwise l = l + 1, and the procedure goes back to Step 2.

• Step 5: Application of optimized control input

Once the iteration is terminated, from the optimized control input sequence U i,opt (k), local MPC controller i applies the first action

u i,opt (k) = N i U i,opt (k), with N i = [ I n i 0 n i ×(Nu-1)n i ]
, to the corresponding zone. The steps of the proposed iterative algorithm are recapitulated in the Figure 5.14. The communication between MPC controllers is assured by the a supervisor at high level as shown in Figure 5.15. We shall note that the communications between local controllers are basically predictions of outputs and no control action of one is acknowledged to others.

Convergence analysis

For a unconstrained case, optimization problem in Equation 5.21 becomes with

U (l) i (k | k) = arg min U i (k|k) 1 2 U ′ i (k | k)H i U i (k | k) + (f (l-1) i ) ′ U i (k | k) (5.23)
H i = ( Φ ′ i,i Q i Φ i,i + R i ) , and
f (l-1) i = Φ ′ i,i Q i ( Ψ i x i (k) + ∑ j̸ =i Φ i,j Ŷ (l-1) j (k | k) -Y i,ref (k)
)

. The result of unconstrained quadratic problem is given as

U (l) i (k | k) = ( Φ ′ i,i Q i Φ i,i + R i ) -1 Φ ′ i,i Q i × [ -Ψ i x i (k) -Φ i,j Ŷ (l-1) j (k | k) + Y i,ref (k)
] .

(5.24) For the overall system, we use notations

U (l) (k | k) = [ (U (l) 1 ) ′ (k | k), ..., (U (l) M ) ′ (k | k) ] ′ and Ŷ (l-1) (k | k) = [ ( Ŷ (l-1) 1 ) ′ (k | k), ..., ( Ŷ (l-1) M ) ′ (k | k) ] ′
, as a result the following relation is obtained

U (l) (k | k) = Π(k) + Γ Ŷ (l-1) (k | k), ( 5.25)
where

Γ =        0 -γ 1 Φ 1,2 • • • -γ 1 Φ 1,M -γ 2 Φ 2,1 0 • • • -γ 2 Φ 2,M . . . . . . . . . -γ M Φ M,1 • • • -γ M Φ M,M -1 0        , γ i = ( Φ ′ i,i Q i Φ i,i + R i ) -1 Φ ′ i,i Q i , with i = 1÷M, (5.26)
and

Π(k) =     γ 1 (-Ψ 1 x 1 (k) + Y 1,ref (k)) . . . γ M (-Ψ M x M (k) + Y M,ref (k))     = F x(k) + GY ref (k) , (5.27) with F = diag(-γ 1 Ψ 1 , ..., -γ M Ψ M ), G = diag(γ 1 , ..., γ M ), x(k) = [x ′ 1 (k), ..., x ′ M (k)] ′ , and Y ref (k) = [Y ′ 1,ref (k), ..., Y ′ M,ref (k)] ′ .
In Equation 5.22, by replacing l = l -1, we have

Ŷ (l) i (k | k) = Ψ i x i (k) + Φ i,i U (l) i (k | k) + ∑ j̸ =i Φ i,j Ŷ (l-1) j (k | k).
(5.28)

And for the whole system, it is written like

Ŷ (l) (k | k) = Ω(k) + DU (l) (k | k) + E Ŷ (l-1) j (k | k), ( 5.29)
with

Ω(k) = diag(Ψ 1 , ..., Ψ M )x(k) = Jx(k), D = diag(Φ 1,1 , ..., Φ M,M )
, and

E =        0 Φ 1,2 • • • Φ 1,M Φ 2,1 0 • • • Φ 2,M . . . . . . . . . Φ M,1 • • • Φ M,M -1 0        .
From Equations 5.25 and 5.29, the recursive equation of system outputs are given as:

Ŷ (l) (k | k) = Ω(k) + D ( Π(k) + Γ Ŷ (l-1) (k | k) ) + E Ŷ (l-1) (k | k) = Ω(k) + DΠ(k) + (DΓ + E) Ŷ (l-1) (k | k).
(5.30)

The convergence of Ŷ (l) (k | k) determines the convergence of the algorithm. In the recursive Equation 5.30, Ω(k) and Π(k) do not depend on the iterative index l, therefore the convergence condition of

Ŷ (l) (k | k) is given as ∥ρ (DΓ + E) ∥ < 1 (5.31)
where the spectral radius ρ (DΓ + E) is the maximum eigenvalues of the matrix (DΓ + E).

Nominal stability analysis

From the state space representation of subsystem i in Equation 5.10, the overall system is constructed by gathering the states, inputs, outputs and matrices of subsystems as

x(k) = [x ′ 1 (k) • • • x ′ M (k)] ′ , u(k) = [u ′ 1 (k) • • • u ′ M (k)] ′ , y(k) = [y ′ 1 (k) • • • y ′ M (k)] ′ , A = diag(A 1 , ..., A M ), B u = diag(B 1 , ..., B M ), C = diag(C 1 , ..., C M ).
Then we have the state space representation of the overall system described as

x(k + 1) = (A + B y C)x(k) + B u u(k) y(k) = Cx(k) (5.32)
where

B y =              0 B 1
B M -1,1 . . . . . . B M -1,M -2 0 B M -1,M B M,1 B M,2 . . . . . . B M -1,M 0              .
If the convergence condition in 5.31 is satisfied, Ŷ (l) (k | k) converges to a final value Ŷ * (k | k) and Equation 5.30 becomes

Ŷ * (k | k) = Ω(k) + DΠ(k) + (DΓ + E) Ŷ * (k | k) ⇒ Ŷ * (k | k) = [I -DΓ -E] -1 Υ [Ω(k) + DΠ(k)] .
(5.33)

From Equations 5.33, 5.27 and 5.25, the optimal control input sequence of the whole system U opt (k) is given as:

U opt (k | k) = Π(k) + Γ Ŷ * (k | k) = Π(k) + ΓΥ [Ω(k) + DΠ(k)] = ΓΥΩ(k) + (I + ΓΥD) Π(k) = ΓΥJx(k) + (I + ΓΥD) [F x(k) + GY ref (k)] = [ΓΥJ + (I + ΓΥD) F ] V x(k) + (I + ΓΥD) G P Y ref (k) = V x(k) + P Y ref (k).
(5.34)

The control input to be applied at sample time k

u(k | k) = N U opt (k | k) = N V x(k) + N P Y ref (k) (5.35)
where

N = diag(N 1 , ..., N M ), N i = [ I n i 0 n i ×(Nu-1)n i ]
. In fact, x(k) is only obtained through state observer, therefore, optimal control input in Equation 5.35 is actually calculated by

u(k | k) = N V x(k) + N P Y ref (k) = N V (x(k) + x(k)) + N P Y ref (k) (5.36)
As a result, the dynamic of closed-loop system-wide is defined by

x(k + 1) = (A + B y C)x(k) + B u u(k) = (A + B y C)x(k) + B u [N V x(k) + N V x(k) + N P Y ref (k)] = (A + B y C + B u N V )x(k) + B u N V x(k) + B u N P Y ref (k) . ( 5

.37)

To analyze the stability of a subsystem and its state observer, we have to add an extended state which is the state estimation error as provided in Equation 5.16 as xi (k + 1) = (A i -L i C i ) xi (k). For the overall system, the extended state is x(k + 1) = [x ′ 1 (k + 1), ..., x′ M (k + 1)] ′ which can be expressed as

x(k + 1) = (A -LC) x(k) (5.38)
where L = diag(L 1 , ..., L M ). The dynamic of closed-loop system including distributed MPC controllers and state observers are given as follows

[ x(k + 1) x(k + 1) ] = [ A + B y C + B u N V B u N V 0 A -LC ] [ x(k) x(k) ] + [ B u N P 0 ] Y ref (k).
(5.39) Therefore, the stability conditions in this case is eig

([ A + B y C + B u N V B u N V 0 A -LC ]) < 1 (5.40)
That is, the eigenvalues must be all inside the unit circle to guarantee the nominal stability of the system.

Serial cooperative distributed MPC

As a furnace zone i interacts only with two neighboring zones (i ± 1), and adding to that, for the considered furnace, H 2,1 , H 3,2 = 0, the local MPC controller i can optimize the performance of a neighborhood of three zones (i, i + 1, i -1). That is, local controllers not only communicate with each other to take into account the coupling effects, but also compromise between individual performance of each zone to achieve a best solution for the entire neighborhood.

First, the formulation of objective function is given. Second, an iterative algorithm is presented. However, due to high complexity, the corresponding convergence analysis and nominal stability of closed-loop will not be examined as in non-cooperative strategy.

Problem formulation

Optimization problem for a local MPC controller i is defined based on performance of its corresponding subsystem and two neighboring subsystems i ± 1 as min

U i (k|k) J i± (k) = i+1 ∑ j=i-1 ∥ Ŷj (k | k) -Y j,ref (k)∥ 2 Q j + ∥U i (k | k)∥ 2 R i (5.41a) subject to u i,min ≤ u i (k + p) ≤ u i,max ; p = 0, ..., N u -1, (5.41b)
∆u i,min ≤ ∆u i (k + p) ≤ ∆u i,max ; p = 0, ..., N u -1, (5.41c)

y j,min ≤ ŷj (k + s | k) ≤ y j,max ; s = 1, ..., N p , j = i, i + 1, i -1.
(5.41d) with Q j = diag(q j,1 , ..., q j,Np ), and R i = diag(r i,1 , ..., r i,Nu ). Now, we want to formulate J i± (k) as function of U i (k | k). Let's consider Ŷi+1 (k | k) by substituting Equation 5.14 as

Ŷi+1 (k | k) =Ψ i+1 x i+1 (k) + Φ i+1,i+1 U i+1 (k | k) + ∑ j̸ =i+1 Φ i+1,j Ŷj (k | k) (5.42a) = Θ i+1 (k) Ψ i+1 x i+1 (k) + Φ i+1,i+1 U i+1 (k | k) + ∑ j̸ =i,i+1 Φ i+1,j Ŷj (k | k) + + Φ i+1,i Ŷi (k | k) (5.42b) =Θ i+1 (k) + Φ i+1,i Ŷi (k | k) (5.42c) =Θ i+1 (k) + Φ i+1,i [Ψ i x i (k) + Φ i,i U i (k | k) + ∑ j̸ =i+1 Φ i,j Ŷj (k | k)] (5.42d) =Φ i+1,i Φ i,i U i (k | k) + Θ i+1 (k) + Φ i+1,i [Ψ i x i (k) + ∑ j̸ =i+1 Φ i,j Ŷj (k | k)] Λ i+1 (k)
.

(5.42e)

In a similar approach, Ŷi-1 (k | k) can be formulated as

Ŷi-1 (k | k) = Φ i-1,i Φ i,i U i (k | k) + Λ i-1 (k) (5.43)
where

Λ i-1 (k) = Θ i-1 (k) + Φ i-1,i [Ψ i x i (k) + ∑ j̸ =i-1 Φ i,j Ŷj (k | k)].
The objective function is now expressed by

J i± (k) =∥Φ i-1,i Φ i,i U i (k | k) + α i-1 (k) Λ i-1 (k) -Y i-1,ref (k) ∥ 2 Q i-1 + ∥U i-1 (k | k)∥ 2 R i-1
(5.44a)

+ ∥Φ i+1,i Φ i,i U i (k | k) + α i+1 (k) Λ i+1 (k) -Y i+1,ref (k) ∥ 2 Q i+1 + ∥U i+1 (k | k)∥ 2 R i+1 + ∥Φ i,i U i (k | k) + Ψ i x i (k) + ∑ j̸ =i Φ i,j Ŷj (k | k) -Y i,ref (k) β i (k) ∥ 2 Q i + ∥U i (k | k)∥ 2 R i (5.44b) =U ′ i (k | k)H i± U i (k | k) + 2f ′ i± U i (k | k) + g i± , ( 5.44c)
where

H i± = [Φ ′ i,i (Φ ′ i-1,i Q i-1 Φ i-1,i + Φ ′ i+1,i Q i+1 Φ i+1,i + Q i )Φ i,i + R i ], f i± = Φ ′ i,i (Φ ′ i-1,i Q i-1 α i-1 (k) + Φ ′ i+1,i Q i+1 α i+1 (k) + Q i β i (k)), and
g i± = α ′ i-1 (k)Q i-1 α i-1 (k) + α ′ i+1 (k)Q i+1 α i+1 (k) + β ′ i (k)Q i β i (k) + ∥U i+1 (k | k)∥ 2 R i+1 +∥U i-1 (k | k)∥ 2 R i-1 .
As a result, at sample time k, the optimization problem of Equation 5.41 becomes a quadratic programming problem as min

U i (k|k) 1 2 U ′ i (k | k)H i± U i (k | k) + f ′ i± U i (k | k) (5.45a) subject to a i± U i (k | k) ≤ b i± (k), (5.45b)
where

b i± (k) =                        U i,max -U i,min Y i,max -Ψ i x i (k) - ∑ j̸ =i Φ i,j Ŷj (k | k) -Y i,min -Ψ i x i (k) + ∑ j̸ =i Φ i,j Ŷj (k | k) Y i+1,max -Λ i+1 (k) -Y i+1,min + Λ i+1 (k) Y i-1,max -Λ i-1 (k) -Y i-1,min + Λ i-1 (k) ∆U i,max -∆U i,min                        , a i± =                        I Nu×nu i -I Nu×nu i Φ i,i -Φ i,i Φ i+1,i Φ i,i -Φ i+1,i Φ i,i Φ i-1,i Φ i,i -Φ i-1,i Φ i,i Ξ i -Ξ i                        , with Ξ i , ∆U i,max , ∆U i,min , Y i+1,max , Y i+1,
min defined in Equation 5.19.

Iterative algorithm

Assume that the communication resource is sufficient, then an iterative algorithm should be designed to ensure the optimal solution of the problem 5.45. At each sample time k, each local MPC controller undergoes the procedures of initialization, communication with other controllers, solving the local optimal problem, checking the convergence and application of optimized control input to the corresponding subsystem, i.e.,

• Step 1: Initialization

The local MPC controller i initializes its control input sequence U i (k | k) by the optimal one found at the previous sample like

U (l=0) i (k | k) = U i,opt (k -1), the output prediction Ŷi (k -1 | k) by the prediction of previous sample as Ŷ (l=0) i (k | k) = Ŷi (k -1 | k -1)
, l is the iteration index.

• Step 2: Communication

For iteration l, local MPC controllers exchange their predicted output

Ŷ (l) m (k | k), current state x m (k), and control input sequence U (l) m (k | k) within a neighborhood. Therefore, controller i receives Ŷ (l) j (k | k), U (l) j (k | k),
(l+1) i (k | k) = arg min U i (k|k) 1 2 U ′ i (k | k)H i± U i (k | k) + (f (l) i± ) ′ U i (k | k) (5.46a) subject to a i± U i (k | k) ≤ b (l) i± (k), (5.46b)
• Step 4 :Updating and Checking

The predicted outputs are updated with the new control input sequence as

Ŷ (l+1) i (k | k) = Ψ i x i (k) + Φ i,i U (l+1) i (k | k) + ∑ j̸ =i Φ i,j Ŷ (l) j (k | k). (5.47) If |U (l+1) i (k | k) -U (l) i (k | k)| ≤ ϵ i for i = 1 ÷ M
, where ϵ i is the convergence precision, or the iteration index l = l max then MPC controllers end the iterations and

U i,opt (k) = U (l+1) i (k | k). Otherwise l = l+1,
and the procedure goes back to Step 2.

• Step 5: Application of optimized control input

Once the iteration is terminated, from the optimized control input sequence U i,opt (k), local MPC controller i applies the first action u i,opt (k) = N i U i,opt (k) to the corresponding zone.

The communication between local controllers is illustrated in Figure 5.16. In contrast to non-cooperative distributed MPC, local controllers communicate not only their output predictions but also system states, and control input sequences in order to calculate objective functions associated to control performance of other zones. Nevertheless, the communication is restricted between neighbor zones, which helps reducing the resource requirements.

Centralized MPC as a benchmark

From the state space representation of the overall system given in Equation 5.32, output can be predicted with the following formula

ŷ(k + s | k) = C(A + B y C) s x(k) + min(s,Nu) ∑ p=1 C(A + B y C) s-p B u u(k + p -1 | k), ( 5
J(k) = Np ∑ s=1 q s [ŷ(k + s | k) -y ref (k + s)] 2 + Nu-1 ∑ p=0 r p+1 u 2 (k + p | k), ( 5.49a)
subject to

u min ≤ u(k + p) ≤ u max ; p = 0 ÷ N u -1, (5.49b
)

∆u min ≤ ∆u(k + p) ≤ ∆u max ; p = 0 ÷ N u -1, (5.49c
)

y min ≤ ŷ(k + s | k) ≤ y max ; s = 1 ÷ N p , ( 5.49d)
In a vector form, output prediction of the system-wide is written as

Ŷ (k | k) = [ ŷ′ (k + 1 | k), ..., ŷ′ (k + N p | k) ] ′ (5.50a) = Ψx(k) + ΦU (k | k), ( 5.50b)
where

U (k | k) = [u ′ (k | k), ..., u ′ (k + N u -1 | k)] ′ is the control input sequence at time instant k, Ψ = [ (C(A + B y C)) ′ , ..., (C(A + B y C) Np ) ′ ] ′ , and
Φ =           CB u • • • 0 . . . . . . . . . C(A + B y C) Nu-1 B u • • • CB u . . . . . . . . . C(A + B y C) Np-1 B u • • • C(A + B y C) Np-Nu B u           .
The optimization problem in Equation 5.49 is rewritten in vector form as min

U (k|k) J(k) =∥ Ŷ (k | k) -Y ref (k) ∥ 2 Q + ∥ U (k | k) ∥ 2 R (5.51a) subject to U min ≤ U (k | k) ≤ U max , (5.51b) ∆U min ≤ ∆U (k | k) ≤ ∆U max , (5.51c) Y min ≤ Ŷ (k | k) ≤ Y max .
(5.51d) From Equation5.50, the objective function can be expressed in form of quadratic function as

J(k) = ∥ Ŷ (k | k) -Y ref (k) ∥ 2 Q + ∥ U (k | k) ∥ 2 R (5.52a) = ∥ Ψx(k) + ΦU (k | k) -Y ref (k) ∥ 2 Q + ∥ U (k | k) ∥ 2 R (5.52b) =U ′ (k | k)((Φ ′ QΦ + R H )U (k | k) + 2[Φ ′ Q(Ψx(k) -Y ref (k) f ] ′ U (k | k)+ (5.52c) + (Ψx(k) -Y ref (k)) ′ Q(Ψx(k) -Y ref (k)) (5.52d) =U ′ (k | k)HU (k | k) + 2f ′ U (k | k) + (Ψx(k) -Y ref (k)) ′ Q(Ψx(k) -Y ref (k)).
(5.52e)

As a result, the optimization problem of Equation 5.53 can be solved with a quadratic programming optimization min

U (k|k) J(k) = U ′ (k | k)HU (k | k) + 2f ′ U (k | k) (5.53a) subject to aU (k | k) ≤ b(k), (5.53b) where b(k) =             U max -U min Y max -Ψx(k) -Y min -Ψx(k) ∆U i,max -∆U i,min             , a =             I Nu×n -I Nu×n Φ -Φ Ξ -Ξ            
, and

Ξ =           I n 0 • • • • • • 0 -I n I n . . . 0 . . . 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 0 • • • 0 -I n I n           Nu.n×Nu.n
.

For an unconstrained case, solution of optimization problem 5.53 is calculated by

U opt (k) = H -1 f (5.54)
and the control input applied at sample time k is u opt (k) = N U opt (k) with N given in 5.35. For output feedback control u opt (k) is actually calculated as

u opt (k) = N U opt (k) = N H -1 f x(k) (5.55a) = N H -1 f [x(k) + x(k)] (5.55b)
The whole closed-loop system with centralized MPC controller and state observers is determined by the state space representation of overall system in Equation 5.32 and the dynamic of state estimation error x(k + 1) = (A -LC)x(k) given in Equation 5.38. With the control solution in Equation 5.55, the closed-loop can be written as

[ x(k + 1) x(k + 1) ] = [ A + B y C + B u N H -1 f B u N H -1 f 0 A -LC ] [ x(k) x(k)
] .

(5.56)

The closed-loop stability condition is therefore eig

([ A + B y C + B u N H -1 f B u N H -1 f 0 A -LC ]) < 1.
(5.57)

Simulation results

We want to examine the performance of the proposed non-cooperative and serial cooperative algorithms. Because in a networked system where subsystems interacts weakly the difference between the non-cooperative method and the cooperative method will not be significant. Therefore, a networked system where three subsystems strongly interact with each other through outputs is constructed in order to clarify difference between two methods. The proposed non-cooperative and cooperative distributed MPC algorithms are applied to solve tracking problem of step signals. The results are compared with the centralized MPC to verify the optimality of distributed MPC methods. Afterwards, the algorithms are applied to zone temperature control of the considered walking-beam reheating furnace. The PID controllers are also simulated for zone temperature control so that the results can be appropriately compared.

Illustrative example

Consider a system composed of three subsystems given in the following subsystem 1:

             x 1 (k + 1) =     0.15 0 1 0 0.23 0.6 0 0 0.4     x 1 (k) +     0.41 0.53 0.62     u 1 (k) +     0.03 0.01 0.1 0.03 0.05 0.2     [ y 2 (k) y 3 (k) ] y 1 (k) = [ 0.2 0.16 0.25 ]x 1 (k) (5.58) subsystem 2:        x 2 (k + 1) = [ 0.4 0.16 0 0.25 ] x 2 (k) + [ 0.6 0.2 ] u 2 (k) + [ 0.4 0.4 0.1 0.3 ] [ y 1 (k) y 3 (k) ] y 2 (k) = [ 0.08 0.09 ]x 2 (k)
(5.59) subsystem 3:

                   x 3 (k + 1) =        0.
J i (k) = Np ∑ s=1 q i [ŷ i (k + s | k) -y i,ref (k + s)] 2 + Nu ∑ p=1 r i u 2 i (k + p -1 | k), i = 1 ÷ 3.
(5.61)

With regards to centralized MPC, the performance index for overall system will be

J(k) = J 1 (k) + J 2 (k) + J 3 (k).
(5.62)

The tuning parameters for each subsystem are set with N p = 10, N u = 5, q i = 1, r 1 = 10 -3 , r 2 = 10 -4 , r 3 = 4.10 -3 with (i = 1 ÷ 3), using sampling time T s = 0.1 s, and convergence precision ϵ i = 0.01 (i = 1 ÷ 3).

Simulation on Matlab is carried out to evaluate the proposed distributed MPC algorithms through comparison with centralized MPC. As shown in Figures 5.17

Control of walking beam reheating furnace

The proposed distributed MPC algorithms for a networked system with subsystems coupled with outputs have been proved to be efficient with the illustrative example. Because the furnace zones are not fully coupled we will utilize only non-cooperative distributed MPC to compare with the PID controllers which is currently used. Here, the output coupling model is provided as:

Y 1 (s) = 0.0143e -60s 1 + 605s U 1 (s) + 0.210e -20s 1 + 500s Y 2 (s), (5.63a) Y 2 (s) = 0.012e -40s 1 + 742s U 2 (s) + 0.16e -20s 1 + 416s Y 3 (s), (5.63b
)

Y 3 (s) = 0.0066e -70s 1 + 415s U 3 (s).
(5.63c)

It shall be noted that this model of furnace zone does not have full couplings between subsystems. The interactions are mainly in the direction of gas flow inside the furnace and considerable only between two neighbor zones. Therefore, coupling effects between zones preheating (y 1 ) and heating 2 (y 3 ) are neglected. Similarly, influences of y 2 on y 3 , y 1 on y 2 are obsoleted. In order to use this model for design of MPC controllers, the delays are approximated by P ade method as e -τ s = 1 1 + τ s .

(5.64)

The interaction model of furnace zone is constituted from second order transfer functions, i.e.,

Y 1 (s) = 1 1 + 60s 0.0143 1 + 605s U 1 (s) + 1 1 + 20s 0.210 1 + 500s Y 2 (s) (5.65a) Y 2 (s) = 1 1 + 40s 0.012 1 + 742s U 2 (s) + 1 1 + 20s 0.16 1 + 416s Y 3 (s) (5.65b) Y 3 (s) = 1 1 + 70s 0.0066 1 + 415s U 3 (s) (5.65c)
This model can be transformed into state space representation then discretized with sample time T s = 20. The dimension of the obtained model is 6. The model of Equation 5.63 corresponds with operating point defined In general, a trade-off between tracking performance and power consumption must be fulfilled with a proper choice of weighting matrices. For an energy reduction mode, the control efforts are minimized while sustaining temperature tracking satisfaction. On the contrary, for a tracking mode, the tracking errors are minimized as much as possible which will cost more control effort. In Figures And the second period is the first semester 2014 when the furnace was operating with distributed MPC. In Figure 5.30, we have data of the daily average specific energy consumption of the furnace in function of daily furnace production. It can be noticed that the more the furnace produces in a day, the less energy it consumes to heat one ton of slab. And, the distributed MPC has reduced furnace consumption compared to PID controllers since red dots are concentrated at a lower level than the blue dots. The results can be seen more clearly in the histogram of daily specific consumption shown in Figure 5.31. The general energy reduction of replacing PID controllers with distributed MPC is 2.5% which correspond to the energy consumption 2500 houses.

In addition to reduction of energy consumption, the oxygen ratio inside the furnace is also reduced and less fluctuated by distributed MPC in comparison to result of PID controllers as shown in Figure 5.32. Since lower and more stable oxygen ratio will reduce the decaburization process, the quality of slabs is therefore improved by distributed MPC.

Conclusions

To conclude, in this chapter the interaction effects between furnace zones are examined. Coupling between zones can be modeled as input coupling through zone input power or output coupling through zone temperature. The output coupling mode is validated and adopted for the control of furnace. According to the adopted zone model, two distributed MPC algorithms are developed for system, where subsystems are coupled by outputs. The first algorithm is a non-cooperative strategy in which local MPC controllers do not account for control performance of other controllers. The second algorithm is constructed so that local MPC controllers account for performance of neighboring controllers. This is due to the fact that a furnace zone interacts mainly with its neighbors. Numerical simulation is carried out for an illustrative case using centralized MPC and distributed MPC strategies, which allows verifying the optimality of proposed non-cooperative distributed MPC algorithm. Simulation for a realistic scenario of furnace operation shows improvement of proposed control strategy in comparison with PID decentralized controllers. The proposed non-cooperative distributed MPC has been implemented in industrial conditions for the level 1 of the considered furnace and provides significant reduction of energy consumption for the furnace. In next chapter, MPC strategy is examined for level 2 of furnace control system to optimize zone temperature set-points.

Problem statement

The problem of slab temperature control in this chapter is limited to a fixed scheduling data which comes from level 3. That is, loading and discharging events of a slab are predefined. As a result, residence time of a slab, component of load inside furnace at a given time is known in advance.

The closed-loop of furnace control system is clarified in Figure 6.2. The desired final temperature of slab j is denoted Tj,end ; it is provided by the level 3 (see Figure 6.1). The final slab temperature is estimated as T j,end via a dynamic simulator of furnace, this simulator is designed based on the non-steady-state model presented in chapter 3 and utilizes the measurements of furnace zone power P i and zone temperature T z,i , with i = 1 ÷ M .

The distributed MPC assures the tracking of zone temperature set-points Tz,i by determining the optimal zone power Pi which minimize zone temperature error and energy consumption criterion. In general, the control of

Qualitative control objectives and constraints

A temperature control system of a slab reheating furnace should satisfy the following control objectives:

• Minimum deviation between the desired and the realized final slab temperature profiles.

• Minimum specific energy consumption which is defined by the ratio Energy supplied by fuel combustion Mass of heated slabs

• Minimum loss of material through scale formation

• Minimum decarburization which has an influence on the material quality Chapter 6. Slab temperature control 101

In general, these control objectives are motivated by economical reasons, in particular operating costs and increasing demands in terms of quality and diversity of products. Obviously, the main objective of slab reheating furnace is that sufficient heat is absorbed and distributed homogeneously into slabs. This objective is defined by a desired value which is referred to as desired final slab temperature profile. Typically, it is imposed by downstream process which is hot rolling mills. The second objective is minimum specific energy consumption which is measure by amount of energy consumed to sufficiently heat one ton of steel slab. This implies minimizing fuel consumption and at the same time minimizing material losses. Scale formation is the process of oxidation of steel which creates iron oxides F e x O y on slab surfaces and therefore leads to loss of material. The decarburization is the process of reduction of carbon content in slabs, this can undesirably modify the quality of steel.

As a result, both scale formation and decarburization are phenomenons to be kept at a minimum. Unfortunately, many mathematical models of reheating furnace, including the one suggested in Chapter 3, are not sufficiently detailed to provide quantitative estimation of the two later objectives. However, a good approach to the first two objectives usually leads to good results of the two latter ones.

The reheating furnace is subject to various constraints, specially temperature constraints due to the extreme environment of furnace interior. A slab temperature control scheme should account the following essential constraints:

• Bounds on the temperature of the furnace interior and the wall in order to protect them from heat damage, this is usually expressed as bounds on zone temperature set-points.

• Limitations of the control inputs: fuel and air feeds.

• Constraint on the slab temperature trajectories.

Temperature control is a complicated task due to following reasons:

• The temperatures inside slabs can not be measured, they are only estimated by the furnace model.

• A slab reheating furnace is a nonlinear MIMO distributed parameter system and exhibits a switched dynamic behavior. That is because the furnace steel load varies over time due to changes of number of slabs inside furnace, and characteristic of these slabs. Furnace components span a wide range of dynamics, that is for example: response times range from fractions of a second for fuel/air valves up to several hours associated with the thermal inertia of furnace walls.

• Slab movement is generally discontinuous and sometimes unforeseeable due to scheduled and unscheduled halts or delays of upstream and downstream process in the work-flow.

Furthermore, the references [Soc, Tot04, HBK09, RSCG86, Sta04] report some recent trends in furnace operation which should be considered in furnace temperature control design. These trends are listed as follows: • An increasing diversification of products with more exigent requirements entails large variations of the desired final temperatures as shown in Figure 6.3, the slab geometry as illustrated in Figure 6.4, the material properties, the residence time, the monetary values, etc. Indeed, there are more and more operating scenarios where each slab inside the furnace is the only one of its kind.

• For higher mechanical properties of product, the requirement concerning temperature gradient inside slabs may be tighter

• Utilization of available furnace capacity should be maximized to avoid situation where the reheating furnace may constitute a bottleneck of the production process. However, this problem is out of the scope of this chapter and will be examined in the next one.

• Non-standard operating phases like furnace start-up or shut-down require particular control strategies. An automatic furnace temperature control system should cope with such situations and keep manual interventions at a minimum level.

• Nowadays, furnace operators have to confront with official emission limits. The minimization of exhaust emissions directly links to maximum efficiency and minimum fuel consumption. Therefore, furnace control design striving for the latter challenges is not only an economical but also an environmental solution.

These trends cause an increased deviation from an ideal steady-state operation of the furnace. Here, the steadystate operation is referred to furnace conditions where all processed slabs have the similar properties, dimensions, desired final temperature profiles, residence time, etc., thus the production rate of the furnace is constant, and loading and discharging events occur at regular time intervals. The greater deviation of an operating scenario from the ideal design point of a furnace and the greater diversification of slab properties, the more difficult is the slab temperature control task. In fact, non-steady-state operation of the furnace has become the normal rather than the exceptional operating scenario. In practice, if the difference between the desired temperature of neighboring slabs is considerable, the slab with highest desired final temperature will have the priority.

Finally, it shall be noted that reheating process is sensitive with the thickness of slab because thicker slabs have slower dynamics behavior compared to that of thin slabs. The slab volume grows linearly with the slab thickness, while the heat exchange area, depends only marginally on the thickness because slab width and length are much bigger than slab thickness. As a result, slab thickness has a direct impact on the necessary heating time for obtaining the desired final temperature profile. A furnace is usually designed for specific slabs with a certain value of thickness to assure that residence times do not deviate too much. As an example, the considered furnace heats slab with thickness of 235 mm.

Quantitative control objectives and constraints

The qualitative control objectives and constraints are now restated in a mathematical form. As temperatures are relevant process variables, it seems reasonable to express control constraints in terms of temperature or their time derivatives, even if these constraints correspond to other physical quantities. There are constraints on zone temperatures and slab temperatures. To facilitate the formulation, the constraints will be presented in the continuoustime domain. The transformation to the discrete-time domain is straightforward.

Constraint of zone temperature set-points

The fuel flow rates to the burners are subject to certain limits. Some burner types cannot be automatically switched off, which leads to a lower bound on the fuel flow rate that is considerably above zero. Moreover, furnace operation needs to maintain the zone temperatures at some minimum values. Generally, these bounds on zone Chapter 6. Slab temperature control 104 temperatures vary over time. For these reasons, temperature set-points of controllable zones and their rate of change are limited as Tz,i,min (t) ≤ Tz,i (t) ≤ Tz,i,max (t) (6.1a)

Ṫz,i,min (t) ≤ Ṫz,i (t) ≤ Ṫz,i,max (t), (6.1b)
where i = 1 ÷ M is the zone number and M is the number of controllable zones. For the considered furnace, M = 4 which corresponds to four controllable zones: preheating, heating 1, heating 2, soaking. The constraints Tz,i,min (t), Tz,i,max (t), Ṫz,i,min (t), and Ṫz,i,max (t) are required for preventing the furnace zone from thermal or thermomechanical damage.

The inequalities 6.1 imply that any satisfactory temperature set-point trajectory Tz,i (t) can be followed by the real furnace system. This is the substantial condition because Tz,i (t) (i = 1, ..., M ) are used as control inputs.

In a vector form, zone temperature set-point vector is denoted by In fact, a time-varying design of zone temperature constraints facilitates occasional interventions of furnace operator in some particular situations. For example, the furnace operator may want to temporarily keep the zone temperature at some constant values or within some narrow range. This happens usually for manual shut-down, start-up, or production halts of the furnace. In these conditions, furnace is deviated from nominal operating point, Chapter 6. Slab temperature control 105 which affects the reheating quality of the slabs. The next section presents how the reheating quality of a slab is defined.

Tz (t) = [

Constraints of slab temperatures

We will consider the heating constraints of a single slab j. The main control objective is that the slab has the final temperature profile close to the desired one denoted Tj,end (y j ) ∀ y j ∈ [-D j /2, D j /2] (D j is thickness of slab j) when it is withdrawn from the furnace at the time t j,1 . Let t j,0 be the time at which the slab is loaded into the furnace, T j (y j , t) ∀ y j ∈ [-D j /2, D j /2] is the temperature profile of slab at time t ∈ [t j,0 , t j,1 ]. As detailed in chapter 3 where the heat conduction problem inside the slabs are presented, the slab temperature profile is determined at limited points inside the slab. The more estimated points the more slab temperature profile is accurate with the price of more required computational resource. For control purpose, a considerably small number of points can be sufficient. For instance, temperatures of 10 points inside a slab are estimated to constitute the temperature profile of the slab. Then, the control would achieve

T j (y j , t j,1 ) = Tj,end (y j ), ∀ y j ∈ [-D j /2, D j /2]. (6.3)
Slab temperatures are subject to hard and soft constraints. Violation of soft constraints is only allowed if oth-Figure 6.5: Constraints on slab temperature trajectory erwise a solution satisfying all hard constraints is infeasible. Nevertheless, because the distinction between hard and soft constraints is a qualitative one, hard constraints may also be broken as a last alternative. The restrictions are usually destined to satisfy metallurgical requirements, and requirements of downstream processes, e.g., hot rolling mills. Constraints may be imposed on absolute temperature or relative temperature. Figure 6.5 shows a typical slab temperature trajectory inside the furnace. Slab temperature profile is hard-constrained through the residence period of the slab with the upper limit T j,max preventing waste of energy as well as damage and loss of material, i.e.,

T j (y j , t) ≤ T j,max ∀ t ∈ [t j,0 , t j,1 ], ∀ y j ∈ [-D j /2, D j /2]. (6.4)
In practice, although this constraint is categorized as hard constraint; it is generally easier to satisfy than other constraints which are presented in the following. The final temperature profile of the slab is constrained with soft upper (lower) limit T j,end,max (T j,end,min ) as shown in Figure 6.6. These bounds allow a certain level of tolerance of final slab temperature profile error which comes from the control error. The formal expressions are

T j,end,min ≤ T j (y j , t j,1 ) ≤ T j,end,max ∀ t ∈ [t j,0 , t j,1 ], ∀ y j ∈ [-D j /2, D j /2]. (6.5)
In spite of being a soft constraint, this constraint is more difficult to be satisfied than the hard constraint on slab temperature provided in Equation 6.4 because it demands more precision of the furnace control. For downstream process, e.g., hot rolling mills, a homogeneous temperature distribution is usually preferred. However, in practice, a certain level of temperature inhomogeneity can be allowed. Let ∆T j (t) denote the temperature gradient of the slab which is defined as

∆T j (t) = max y j
{T j (y j , t)} -min

y j {T j (y j , t)}, y j ∈ [-D j /2, D j /2], ∀ t ∈ [t j,0 , t j,1 ]. (6.6)
Large temperature differences within a slab can cause undesirable global deformation like bending of slab.

In addition, excessive inhomogeneities of final temperature profile may lead to unwanted incident during the transport of hot slab, and rejection at the rolling process because average temperature of slab at entrance of rolling mill may be not sufficiently high. Therefore, the inhomogeneity ∆T j (t) of the slab temperature is subject to an upper bound hard constraint and a tighter final bound but considered as a soft constraint as shown in Figure 6.7. The formal expressions are written as:

∆T j (t) ≤ ∆T j,max ∀ t ∈ [t j,0 , t j,1 ], (6.7a)
∆T j (t j,1 ) ≤ ∆T j,end,max , (

where ∆T j,max is a constant, and generally ∆T j,end,max ≤ T j,end,max -T j,end,min . In addition to forgoing constraints and objectives, to reduce energy consumption, one has to take into account the following objectives:

• The slab temperature trajectories should be kept as low as possible, which may also reduce material losses through scale formation and allow lower zone temperatures • The zone temperature is kept at minimum sufficient level, which reduces directly fuel supply into the furnace In Table 6.2, the constraints of a slab temperature profile are recapitulated. In the next section, the presented objectives and constraints will be accounted for a set of slabs in the design of MPC strategy for slab temperature control.

Variables

MPC strategy applied to slab temperature control

This section describes the MPC approach for the control at level 2 of furnace control system. The discussion will begin with how the system can be seen from level 2 point of view by defining its inputs and outputs. Next, the formulation of the problem is given for the design of MPC method. And finally, some important issues on computational requirement and stability are discussed.

Principle

As previously mentioned in chapter 5, MPC is a model-based feedback control concept. For level 2 of furnace control system, control input denoted u(t) are defined by the zone temperatures as

u(t) = [u 1 (t), ..., u M (t)] ′ (6.8a)
= [ Tz,1 (t), ..., Tz,M (t)] ′ . (6.8b) The control problem is clearly the tracking of the desired final slab temperature profile. Preferably, the desired final slab temperature profile would be homogeneous. However, in real process, tolerance of a certain level of inhomogeneity is adopted. As a result, the control problem can be defined by the tracking of an average desired final slab temperatures and the minimization of the slab temperature gradients. The outputs to be controlled are therefore average final slab temperatures, and slab temperature gradients, respectively defined as follows

T s (t) = { mean y j {T j (y j , t j,1 )}, if ∃ j ∈ S, t = t j,1 0, otherwise (6.9a)
g(t) = { ∆T j (t j,1 ), if ∃ j ∈ S, t = t j,1 0, otherwise (6.9b)
where S denotes the set of slabs to be processed and the subscript s stands for slab. It can be noticed that outputs of the system can also be considered as discontinuous variables where discontinuous points correspond with discharging times of slabs. Consequently, set-point for output T s (t), and its bounds can be defined as

Ts (t) = { Tj,end , if ∃ j ∈ S, t = t j,1 0, otherwise (6.10)
Ts,min (t) = { Tj,end,min , if ∃ j ∈ S, t = t j,1 0, otherwise (6.11) Ts,max (t) = { Tj,end,max , if ∃ j ∈ S, t = t j,1 0, otherwise (6.12)
And the constraint on maximum allowed gradient temperature can be defined as

g max (t) = { ∆T j,end,max , if ∃ j ∈ S, t = t j,1 0, otherwise (6.13)
Since slabs are loaded and discharged gradually, the furnace system has actually a switched dynamic. Clearly, the variables of interest of the system at each time t are temperature profiles of all slabs that are currently inside the furnace. However, we see that slab temperature constraints can be expressed in terms of the temperature gradient ∆T j (t) and maximum temperature inside slabs: max y j T j (y j , t). Because most of the constraints concerns the upper bounds of slab temperature profile and the temperature gradient, we chose the vector x(t) = [x ′ 1 (t), x ′ 2 (t)] ′ as the variable of interest, with x 1 (t) and x 2 (t) defined as follows:

x 1 (t) = [max y j+1 T j+1 (y j+1 , t), max y j+2 T j+2 (y j+2 , t), ..., max y j+N in T j in (y j+N in , t)] ′ , (6.14a) x 2 (t) = [∆T j+1 (t), ∆T j+2 (t), ..., ∆T j+N in (t)] ′ , y j ∈ [-D j /2, D j /2], (6.14b)
with T j+1 (y j+1 , t) is temperature profile of slab situating at entrance of the furnace, and T j+N in (y j+N in , t) is the temperature profile of slab being at the exit door of the furnace. N in is the number of slab residing inside the furnace at time instant t and can change over time. S in (t) defined by [j + 1, j + 2, ..., j + N in ] denotes the index set of slabs that are inside the furnace at time t. Model predictive control principles applied to slab temperature control can be now explained by the following steps.

• At the present time τ the variable x(τ ) is estimated by the non-steady state model. Based on this model the system outputs T s (t), g(t) are predicted over the prediction horizon [τ, τ + ∆t P ] with ∆t P the prediction horizon. The basic data necessary for the forecast is composed of: future system control input u(t) ∀ t ∈ [τ, τ + ∆t C ] (∆t C is the control horizon) and future knowledge about furnace configuration, disturbances, etc.

• The control input u(t) ∀ t ∈ [τ, τ + ∆t C ] is calculated so that the predicted system outputs T s (t), g(t)

∀ t ∈ [τ, τ + ∆t P ] match the desired values as much as possible. This matching criterion is specified by an objective function. In general, the control input u(t) is obtained as a solution of an open-loop optimal control problem defined in the interval [τ, τ + ∆t P ).

• The optimized control input u(t) is applied to the real system over the interval [τ, τ + ∆τ s ), with ∆τ s represents the time interval between consecutive optimization events.

• The procedure is repeated for the optimization at τ + ∆τ s to redesign the optimal control input u(t

) ∀ t ∈ [τ + ∆τ s , τ + ∆τ s + ∆t C ].
The whole procedure is presented in Figure 6.9. In the following, formulation of the optimal control problem is clarified with mathematical equations accounting for furnace control objectives and constraints.

.2.1 Problem formulation

In order to understand more the furnace control task, we will first approach the problem in continuous-time domain. To formulate the optimization problem, some terminological definitions are provided for time instant τ :

• Control input sequence U (τ ) = {u(t), ∀t ∈ [τ, τ + ∆t P }. It shall be noted that u(t) = u(τ + ∆t C ), ∀t ∈ [τ + ∆t C , τ + ∆t P ].
• Performance indexes on error of final average slab temperature, on inhomogeneity of final slab temperature, on energy consumption of the furnace are denoted by J 1 (U (τ )), J 2 (U (τ )), J 3 (U (τ )), respectively.

General control objectives of low overall furnace temperature and low slab temperature trajectories are accounted by the term J 4 (U (τ )) and J 5 (U (τ )), respectively. These performance indexes will be described in the next section.

• Figure 6.8 illustrates the sets S in (τ ): index set of slabs that are inside the furnace at time τ , and S 1 (τ ): the index set of slabs that are discharged during the prediction horizon period [τ, τ + ∆t P ].

• Upper bounds on slab temperature trajectories and temperature gradients are denoted T max (t) and ∆T max (t)

respectively. These variable are switched accordingly with loading and discharging events (changes of Chapter 6. Slab temperature control 111 S in (t)).

T max (t) = [T j+1,max , T j+2,max , ..., T j+N in ,max ], (6.15a)

∆T max (t) = [∆T j+1,max , ∆T j+2,max , ..., ∆T j+N in ,max ], (6.15b)
with j ∈ S in (t) = [j 0 , j 0 + 1, ..., j in ] and t ∈ [τ, τ + ∆t P ]. The optimization problem is then defined as minimization of the overall objective function J o (U (τ )) which is the sum of weighted objectives concerning slab temperature, zone temperature, and energy consumption. This minimization problem is subjected to the constraints of zone temperature set-points and slab temperature profiles. The formal expression is given as: min

U (τ ) J o (U (τ )) = α 1 J 1 (U (τ )) + α 2 J 2 (U (τ )) + α 3 J 3 (U (τ )) + α 4 J 4 (U (τ )) + α 5 J 5 (U (τ )) (6.16a) subject to T z,i,min (t) ≤ u i (t) ≤ T z,i,max (t), ∀i = 1 ÷ M, t ∈ [τ, τ + ∆t P ] (6.16b) Ṫz,i,min (t) ≤ ui (t) ≤ Ṫz,i,max (t), ∀i = 1 ÷ M, t ∈ [τ, τ + ∆t P ] (6.16c) T j,end,min ≤ T s (t) ≤ T j,end,max , if j ∈ S 1 (τ ), t = t j,1 , t ∈ [τ, τ + ∆t P ] (6.16d) g(t) ≤ ∆T j,end,max , if j ∈ S 1 (τ ), t = t j,1 , t ∈ [τ, τ + ∆t P ] (6.16e
)

x 1 (t) ≤ T max (t), ∀t ∈ [τ, τ + ∆t P ] (6.16f ) x 2 (t) ≤ ∆T max (t), ∀t ∈ [τ, τ + ∆t P ] (6.16g)
where α 1 , α 2 , α 3 , α 4 , α 5 are positive weighting coefficients. The choice of these coefficients will determine the priority order of different control objectives. In general, the objectives of high homogeneity and accurate final slab temperatures are the most important, then comes the objective of energy efficiency that we want to improve.

Performance indexes

The criterion J 1 (U (τ )) accounts for the deviation of final average slab temperatures from the desired ones. Let j be the index of a slab that is discharged from the furnace during [τ, τ + ∆t P ], i.e., j ∈ S 1 (τ ), then J 1 (U (τ )) is defined by the sum of weighted error of final temperature of slabs discharged during the prediction horizon, i.e.,

J 1 (U (τ )) = ∑ j∈S 1 (τ ) w 1 (z j (τ )) × |T s (t j,1 ) -Ts (t j,1 )|, t j,1 ∈ [τ, τ + ∆t P ], (6.17)
where w 1 (z j (τ )) is a positive weighting coefficient corresponding with slab j, and z j (τ ) is the position of slab j at time τ inside the furnace which refers to the distance from the slab to the furnace exit. The adjustment of w 1 (z j (τ )) is given in Section 6.2.4.

Chapter 6. Slab temperature control 112

The objective of homogeneous slab temperature profile is expressed as

J 2 (U (τ )) = ∑ j∈S 1 (τ ) w 2 (z j (τ )) × g(t j,1 ), t j,1 ∈ [τ, τ + ∆t P ], (6.18)
where w 2 (z j (τ )) is a positive weighting coefficient corresponding to slab j.

The criterion on consumed energy can be formulated directly with the power of furnace zone as:

J 3 (U (τ )) = ∫ τ +∆t P τ w 3 (t) M ∑ i=1 P i (t)dt (6.19)
where P i (t) is the heat power supplied to zone i at time t, w 3 (t) is a weighting coefficient dependent on time instant t. In fact, P i (t) is predicted by the furnace model, therefore, prediction errors may constitute a problem.

As a result, instead of directly calculate total consumption of the the furnace over the prediction horizon, we introduce the time-dependent weighting coefficient w 3 (t) to reduce the effects of prediction errors.

Nevertheless, the zone temperatures are linked proportionally to the energy consumption of the furnace. Thus, an optimization strategy will have to maintain the overall furnace temperature as low as possible. The corresponding objective function J 4 (U (τ )) accounting for zone temperature is given as

J 4 (U (τ )) = ∫ τ +∆t P τ w 4 (t) M ∑ i=1 u i (t)dt (6.20)
where w 4 (t) is a weighting coefficient dependent on time instant t.

The slab temperatures should be kept as low as possible to prevent scale formation and decarburization. Therefore, the corresponding objective function J 5 (x(t)) will take into account the maximum temperatures inside slabs which are represented by the vector x 1 (t) defined in Equation 6.14. The formulation of J 5 (x(t)) is given by

J 5 (U (τ )) = ∫ τ +∆t P τ [x ′ 1 (t) × w 5 (t) × x 1 (t)]dt (6.21)
where w 5 (t) is a weighting positive definite square matrix having the same dimension as length of x 1 (t).

We have introduced different coefficients in the description of objective functions. In contrast to weighting coefficients α i (i = 1 ÷ 5) (c.f Equation 6.16) which determine the priority order of the control objectives in the overall objective function, ω 1 (z j (τ )) and ω 2 (z j (τ )), associating with slab j, define the priority order of the discharged slabs in the corresponding objective functions J 1 (U (τ )) and J 2 (U (τ )). That is because the latter objective functions involve with discontinuous variables: the final temperature gradient g(t) and the final slab temperature T s (t). For the objective functions that concern continuous variables such as: the heat power P i (t), the zone temperature set-points (control inputs) u i (t), and the maximum slab temperature vector x 1 (t), the corresponding weighting coefficient ω 3 (t), ω 4 (t), and ω 5 (t) are time-dependent in order to reduce the effect of the prediction errors.

In the following, constraints in the minimization problem 6.16 are reformulated in forms of penalty terms in order to transform it into an unconstrained minimization problem.

Approximate constraint criterion

In the problem statement, we have identified the constraints of control task. Hard constraints of a slab j are the upper limit of slab temperature T j,max and maximum values of temperature gradient ∆T j,max . Besides, soft constraints of the slab j are the upper (lower) bound of final temperature T j,end,max (T j,end,min ), and the upper bound of temperature gradient ∆T j,end,max . In practice, hard constraints are crucial to prevent incidents such as furnace stoppages due to slab bending or deformation. Fortunately, these hard constraints can be satisfied more easily than soft constraints which needs high accuracy of temperature control. As a results, we will corporate not only the soft constraints as penalty terms but also the hard constraints. In this way, the constrained optimization problem is transformed into an unconstrained one. Firstly, limits of zone temperatures and their rate of change are taken into account as follows:

J c,1 (τ ) = ∫ τ +∆t C τ M ∑ i=1 [max(0, u i (t) -T z,i,max (t)) + | min(0, u i (t) -T z,i,min (t))|]dt (6.22a) J c,2 (τ ) = ∫ τ +∆t C τ M ∑ i=1 [max(0, ui (t) -Ṫz,i,max (t)) + | min(0, ui (t) -Ṫz,i,min (t))|]dt (6.22b)
The idea in these formulations is that if the control inputs or their rate of change get out of their upper and lower limits the penalty terms J c,1 (τ ) and J c,2 (τ ) will take a positive values . Otherwise, they take values of zero.

As defined in Equation 6.9, the final temperature gradient variable g(t) is a discontinuous output of the system.

Therefore, constraints on g(t) can be expressed either in integral continuous form or as sum of discontinuous values, i.e.,

J c,3 (τ ) = ∫ τ +∆t P τ max(0, g(t) -g max (t))dt (6.23a) = ∑ j∈S 1 (τ )
max (0, g(t j,1 ) -∆T j,end,max ) . (

If the final temperature gradient of one slab or more exceeds its upper bound, J c,3 (τ ) will take a positive values, otherwise, it becomes zero.

Similarly, soft constraint on final temperature of slabs is formulated as:

J c,4 (τ ) = ∫ τ +∆t P τ [max(0, T s (t) -Ts,max (t)) + | min(0, T s (t) -Ts,min (t))|]dt (6.24a) = ∑ j∈S 1 (τ )
max(0, T s (t j,1 ) -Ts,max (t j,1 )) + | min(0, T s (t j,1 ) -Ts,min (t j,1 ))|] (

Hard constraints on slab temperatures and temperature gradients can be formulated as

J c,5 (τ ) = ∫ τ +∆t P τ [max(0, [x 1 (t) -T max (t)] ′ [x 1 (t) -T max (t)])]dt (6.25) J c,6 (τ ) = ∫ τ +∆t P τ [max(0, [x 2 (t) -∆T max (t)] ′ [x 2 (t) -∆T max (t)])]dt (6.26)
Global objective function takes into account the control objectives and the temperature constraints:

J o (U (τ )) = 5 ∑ i=1 α i × J i (U (τ )) + µ 6 ∑ i=1 J c,i (U (τ )) (6.27)
where µ is a positive weighting coefficient. It can be noticed that the penalty terms do not contain weighting coefficient in contrast to objective functions (previous section) in which weighting coefficients are introduced to impose a priority order between slabs or to cope with prediction errors. µ is often imposed with huge values compared to α i (i = 1 ÷ 5). The main idea is to be able to penalize the overall objective function

J o (U (τ ))
whenever a single constraint is not satisfied.

The optimization becomes an unconstrained minimization problem.

min

U (τ ) J o (U (τ )) = 5 ∑ i=1 α i × J i (U (τ )) + µ 6 ∑ i=1 J c,i (U (τ )). (6.28)

Discrete-time problem

In order to implement the control strategy online, the problem has to be discritized in discrete-time domain.

Sampling time for level 2 is generally higher than for level 1. Some terminological definitions are added or redefined for discrete time k:

• Sampling time of level 2: ∆t s

• Prediction horizon [k, k + N P ], with N P = int(∆t P /∆t s ) • Control horizon [k, k + N C -1], with N C = int(∆t C /∆t s ) ≤ N P • Control input sequence U (k) = {u(k), ∀t ∈ [k, k + N C -1]}. It shall be noted that u(k) = u(N C ), ∀k ∈ (N C , N P ].
• Discrete times of loading and discharging events of a slab j are denoted: k j,0 and k j,1 , respectively. Their relations with continuous times of loading and discharging events is defined by the follows:

k j,0 = int(t j,0 /∆t s ) (6.29a) k j,1 = int(t j,1 /∆t s ). (6.29b)
A recapitulation of continuous parameters and their corresponding discrete forms is given in Table 6 The unconstrained minimization problem can be written as min

U (k) J o (U (k)) = 5 ∑ i=1 α i × J i (U (k)) + µ 6 ∑ i=1 J c,i (U (k)) (6.30)
where J i (U (k)), i = 1 ÷ 5 and J c,i (k), i = 1 ÷ 6 are straightforwardly obtained from continuous forms:

J 1 (U (k)) = ∑ j∈S 1 (τ ) w 1 (z j (k)) × |T s (k j,1 ) -Ts (k j,1 )|, k j,1 ∈ [k, k + N P ]. (6.31)
J 2 (U (k)) = ∑ j∈S 1 (τ ) w 2 (z j (k)) × g(k j,1 ), k j,1 ∈ [k, k + N P ]. (6.32)
J 3 (U (k)) = N P ∑ i=1 w 3 (k + i) M ∑ j=1 P j (k + i) (6.33) J 4 (U (k)) = N C ∑ i=1 w 4 (k + i) M ∑ j=1 u j (k + i) (6.34) J 5 (U (k)) = N P ∑ i=1 [x ′ 1 (k + i) × w 5 (k + i) × x 1 (k + i)] (6.35) J c,1 (k) = N P ∑ j=1 M ∑ i=1 [max(0, u i (k + j) -T z,i,max (k + j)) + | min(0, u i (k + j) -T z,i,min (k + j))|] (6.36a) J c,2 (k) = N P ∑ j=1 M ∑ i=1 [max(0, ui (k + j) -Ṫz,i,max (k + j)) + | min(0, ui (k + j) -Ṫz,i,min (k + j))|] (6.36b) Chapter 6. Slab temperature control 116 J c,3 (k) = N P ∑ i=1 max(0, g(k + i) -g max (k + i)) (6.37a) = ∑ j∈S 1 (τ )
max (0, g(k j,1 ) -∆T j,end,max ) . (6.37b)

J c,4 (k) = N P ∑ i=1 [max(0, T s (k + i) -T s,max (k + i)) + | min(0, T s (t) -T s,min (k + i))|] (6.38a) = ∑ j∈S 1 (τ ) max(0, T s (k j,1 ) -T s,max (k j,1 )) + | min(0, T s (k j,1 ) -T s,min (k j,1 ))|] (6.38b) J c,5 (k) = N P ∑ i=1 [max(0, [x 1 (k + i) -T max (k + i)] ′ [x 1 (k + i) -T max (k + i)])] (6.39) J c,6 (k) = N P ∑ i=1 [max(0, [x 2 (k + i) -∆T max (k + i)] ′ [x 2 (k + i) -∆T max (k + i)])] (6.40)

Parameters adjustment

Based on the dynamic optimization problem (6.30), optimal control input u(k) = [ Tz,1 (k), ..., Tz,M (k)] is computed for the system-wide. Therefore, the design of cost function is crucial for the performance of the control.

The cost function contains objective functions and penalty terms of furnace constraints that ensure the accurate heating of slabs. The fuel consumption of furnace predicted by non-steady-state model is accounted in the cost function, which assures the energy efficiency of the heating. General objectives that contribute to reduction of energy, loss of material are the minimization of furnace zone temperatures and the slab temperatures. In the formulation of objective functions, we use weighting coefficients associated with slabs as: w 1 (z j (k)), w 2 (z j (k)) with j ∈ S 1 (τ ) the index set of slabs that are discharged during the prediction horizon. The relative ranking of these coefficients is the information of interest rather than their absolute values. That is, these coefficients will be defined by how a slab is considered more important than another. On one hand, a slab that will be discharged from the furnace before another must have a higher heating priority. This is to ensure that desired average final slab temperature ( Ts (t)) and required slab temperature homogeneity (imposed with g max (t)) are achieved. On the other hand, in a normal furnace operation, it happens that two consecutive slabs may have very different desired final temperatures. In this case, due to the slow dynamics of the furnace, it is difficult or even impossible for the furnace control to correctly realize desired final temperature of both. The control may either under-heat one slab and correctly heat the other, or overheat one slab and correctly heat the other as illustrated in Figure 6.10. Under-heated slab risks to be rejected by hot rolling mills, which represents a loss of energy and material. Overheating of both of slabs will need more energy, but generally less expensive than the handling of slab rejections. For this reason, overheating is preferred to under-heating. Therefore, between two slabs, one with higher desired final temperature will be more prioritized. In addition, some slabs may have extra requirements of quality, this could be also taken into account into the weighting coefficients. To recapitulate w 1 (z j (k)), w 2 (z j (k) can have the following forms accounting for criteria of: slab position, desired final slab temperature, quality requirement, i.e.,

w 1 (z j (k)) = γ(z j (k)) × γ T,j (6.41a
)

w 2 (z j (k)) = β(z j (k)), (6.41b)
provided that γ(z j (k)), β(z j (k)) are functions of slab position at time k and γ T,j associates with desired final slab temperature. Let slabs l, j be discharged with desired final temperature Tl,end , and Tj,end , respectively. If the slab l is discharged before slab j meaning z l (k) < z j (k) then the corresponding weighting coefficients of slab l have bigger values, i.e., γ(z l (k)) > γ(z j (k)) and β(z l (k)) > β(z j (k)). This implies also that the more a slab approaches the furnace exit the heavier it is weighted in the objective functions as illustrated in Figures 6.11 and 6.13. γ(z j (k)) can take the following form:

γ(z j (k)) = L -z j (k) L (6.42)
where L is the furnace length. Similarly, if slab j has a higher desired final temperature, i.e., Tl,end > Tj,end then it is weighted more, i.e.,γ T,j > γ T,i . For example, γ T,j can take the following form accounting the desired final slab temperature

γ T,j = Tj,end Tend,nom (6.43)
where Tend,nom is a nominal value of desired final slab temperature.

Furthermore, each furnace zone is dedicated to a specific stage of heating process, e.g., zone soaking heats slabs from 1200 • C to 1250 • C, zone heating 2 from 1000 • C to 1200 • C, zone heating 1 from 600 • C to 1000 • C, zone preheating from 300 • C to 600 • C as shown in Figure 6.12. The zone soaking is dedicated mostly to equalize the

β(z j (k)) = { L-z j (k) L , if j is located in soaking zone η L-z j (k) L , η ∈ [0, 1], otherwise (6.44)
where η is a positive constant rescaling the weighting coefficient of slabs located in heating zones, preheating and recuperation zones. then : w 3 (k + j) < w 3 (k + i), w 4 (k + j) < w 4 (k + i), and w 5 (k + j) < w 5 (k + i).

Simple linear form of w 3 (k), w 4 (k) and w 5 (k) can be defined as

w 4 (k + i) = (N C -i)/N C , ∀i ∈ [1, 2, ..., N C ], (6.46
)

w 3 (k + i) = (N P -i)/N P , ∀i ∈ [1, 2, ..., N P ], w 5 (k + i) = (N P -i)/N P , ∀i ∈ [1, 2, ..., N P ].
A recapitulation of different weighting coefficient is given in Figure 6.14. The weighting coefficients α i (i = 1 ÷ 5) determine the priority order between control objectives. In general, low error of final slab temperatures and low temperature gradients are the most important control objectives. And, reduction of energy consumption comes as the less important control objective, however it is the main objective of the work. Penalty terms are weighted by the coefficient µ which is usually imposed with a big value in order to penalize the overall objective function J o (U (k)) whenever a single constraint is not satisfied. To recap, the relative rank of the coefficients µ and α i (i = 1 ÷ 5) can be described as

µ ≫ α 1 , α 2 > α 3 , α 4 , α 5 . ( 6.47)
In order to determine the prediction horizon, we should note that, if a long optimization horizon N P is chosen, the optimal solution of time instant k, i. e., U (k) may be completely different from optimal solution of time instant k + N P , i. e., U (k + N P ). Especially, in reheating furnaces which have switched dynamics due to the continuous flow of slabs passing through the furnace interior, the control input u(k) applied at time k has significant influence in a near future, but this influence is weakened as time goes on and practically negligible after some time. In the furnace control problem, the prediction horizon should be greater than the time that a slab resides in controllable zones (preheating, heating 1, heating 2, soaking). In average, this period corresponds to 2/3 the residence time of a slab inside the furnace. Therefore we have the conditions for prediction horizon:

∆t P ≥ 2 3 (t j,1 -t j,0 ), (6.48)
with j a representative slab. For the considered furnace the average residence time of slabs is 3 h. Thus, ∆t P = 120 min is an appropriate choice.

The control horizon N C should be chosen with consideration of the available computational resource of the control system. N C × ∆t s is lower bounded by the delay of the furnace which has at least the value of sampling period ∆t s . In fact, ∆t s should be chosen according to the changes of system dynamics. These changes are associated with loading and discharging events. For the considered furnace, average of time interval between two consecutive loading (discharging) events is about 5 min, thus we might chose ∆t s = 5min. As a result, the prediction horizon N P takes a value of ∆t P /∆t s = 24.

Optimization methods

In order to solve the unconstrained minimization problem in 6.30, the objective function is calculated based on variables predicted by a simulation of the furnace. This simulation is constructed based on the non-steadystate model presented in chapter 3. The numerical simulation involves heat exchange between slabs and their surrounding environment, switching of furnace load, movement of slabs inside furnace, energy losses into door openings or the cooling system, etc. Therefore, an analytic formula of the objective function cannot be obtained.

As a result, derivatives of objective function are not explicitly formulated.

In term of optimization methods, there are two major classes: gradient-based methods, and gradient-free methods.

Gradient-based methods demand the function to be continuously differentiable. If only first order gradients are needed, the method is often referred to as gradient search or first order gradient search. Examples of first order methods are method of steepest descent and line search method. If first and second order gradients are demanded, the optimization methods are referred to as second order methods such as Newton-Raphson and quasi-Newton methods.

Gradient-free methods are methods that need only values of objective function or the ranking of these values. with different values of the control input U (k). Because real-time implementation needs the optimization algorithm to provide optimal control input in period of one sample time, the repetition of furnace simulation is limited within a certain number l max . Besides, provided non-steady-state model is more elaborated than the steady-state model and its computation time is also considerably longer. But we also want to benefit the advantage of non-steady-state model in the optimization. Therefore, gradient-based optimization methods will not be studied further. Utilization of gradient-free methods will be more appropriate for optimization problem of reheating furnace.

In the following, a more detailed of furnace control system is described. The existing level 2 of the furnace is presented showing how furnace model is used in furnace observer and set point optimizer. Afterwards, the basic ideas of different direct search methods are presented. Considering the characteristic of reheating furnace optimization problem, an appropriate method will then be chosen and adapted for the considered furnace. are lower than their desired values, the zone temperature set points are incremented with a step ∆U 0 . Inversely, if the simulated end-zone slab temperatures are too much higher than their desired values (overheating), the zone temperature set points are decremented with the step ∆U 0 . When the simulated end-zone temperatures of all slabs located inside a furnace match their requirements, the zone temperature set points are accepted and applied to the furnace. Zone temperatures realized by the temperature control of level 1 are measured and utilized in the observer to estimates the system state (slab temperatures). The system state is then utilized by simulation program to predict end-zone slab temperature. These procedures restart at each sampling time.

Simulation-based optimization

It can be noticed that the method used in the existing level 2 is simply calculating adequate zone temperature set points to achieve the desired final slab temperature. The method is not optimizing the furnace heating process in terms of energy consumption and final slab temperature at the same time. In the following section, the solution for optimization problem of the furnace previously defined in Equation 6.30 is studied. Different optimization methods are described. The appropriate method is then chosen for our problem.

Direct search methods

Direct search methods are reviewed in [START_REF] Lewis | Direct search methods: then and now[END_REF][START_REF] Kolda | Optimization by direct search: New perspectives on some classical and modern methods[END_REF]. The authors show the distinguishing characteristic of direct search methods compared to other methods where either explicit gradients or their approximated values are used. From practical point of view, three classical direct search methods concerning unconstrained minimization are classified as: pattern search method, simplex method (different from simplex method for linear programming), and method with adaptive sets of search directions.

Pattern search methods

A generalized formulation of pattern search methods is given in [START_REF] Torczon | On the convergence of pattern search algorithms[END_REF]. We consider a problem of minimizing a continuously differentiable function

min x∈ℜ n f (x) (6.49)
where f (x) : ℜ n → R. (6.50)

We denote by ℜ, Q, Z, and N the sets of real, rational, integer, and natural numbers, respectively. If we use direct search methods for this problem then neither analytic nor approximated derivatives of f (x) are needed.

Basic idea of pattern search methods is utilization of a series of exploratory moves that consider the behavior of the function f at a pattern of points (see Figure 6.17 ) to find the minimum.

Firstly, to define a pattern, we need two components called a basic matrix and a generating matrix. The basic matrix can be given by any nonsingular matrix B ∈ ℜ n×n . And the generating matrix is denoted C k ∈ Z n×p , Figure 6.17: Example of search pattern in ℜ 2 with a given step length parameter ∆ k with k is the iterate index and p > 2n, that is C k may be actualized every iterate k. And C k is partitioned into components as

C k = [M k | -M k | L k ] = [Γ k L k ] (6.51)
where M k ∈ M ⊂ Z n×n , with M is a finite set of nonsingular matrices, and L k ∈ Z n× (p-2n) and has at least one column of zeros. Note that both B and C k have rank n.

Then, we define pattern P k by

P k = BC k = [BM k | -BM k | BL k ] = [BΓ k | BL k ] (6.52)
Because B and C k have rank n, the column vectors of P k span ℜ n . Provided ∆ k ∈ ℜ, and ∆ k > 0 as a step length parameter, we define a trial step s i k as a vector calculated by

s i k = ∆ k Bc i k , ( 6.53)
where c i k denotes a column of

C k = [c 1 k ...c p k ]
. We can see that Bc i k determines the direction of the trial step. Then, a trial point is defined as any point having the form x i k = x k + s i k , with x k is the current iterate. In Figure 6.17, an example of search pattern in ℜ 2 is illustrated. From he current iterate x k , eight trial points are defined with the step length ∆ k . In this example, the basic matrix B is the identity matrix, i.e., B = [1 0; 0 1]. The generating matrix C k is fixed for all k. In this case, C k = C 0 , where

C 0 = [ 1 0 -1 0 1 1 -1 -1 0 0 1 0 -1 1 -1 -1 1 0 ] .
These parameters are defined for the coordinate search algorithm provided in the Appendix B. The pattern P k is used to guide a series of search from the current iterate x k for a new iterate x k+1 = x k +s k giving a lower function value, i.e., f (x k+1 ) ≤ f (x k ). These series of search is referred to the explanatory moves mentioned previously.

At one extreme, an effective explanatory moves algorithm would be one that somehow guesses which of the and c i k ≤ C, ∀k, ∀i = 1, ..., p, and lim k→+∞ ∆ k = 0. Then, the convergence of x k produced by the algorithm of generalized pattern search method satisfying the strong hypotheses of explanatory moves is given as lim k→+∞ ∥▽f (x k )∥ = 0. (6.58) One can be sure that when a reflected vertex produces a strict decrease on objective function compared to the current best vertex, a new candidate for function minimizer is obtained. For a constrained problem, a variant of Nelder-mead simplex algorithm which copes with bounds of variables is mentioned in [START_REF] Luersen | Globalized nelder-mead method for engineering optimization[END_REF]. The bounded variables can leave the domain after either the reflection or the expansion operation. The problem is handled by projection Due to its heuristic characteristic, Nelder-Mead simplex method lacks an explicit mathematical convergence proof. In [START_REF] Lagarias | Convergence properties of the nelder--mead simplex method in low dimensions[END_REF], convergence results are obtained for strictly convex functions in dimensions 1 and 2. The Nelder-Mead simplex method is proved to converge to a minimizer in 1-dimensional problem. And various limited convergence results is obtained for dimension 2. In [START_REF] Mckinnon | Convergence of the nelder-mead simplex method to a nonstationary point[END_REF], the method failed to converge to the minimizer Chapter 6. Slab temperature control 130 from a specific starting point for a family of strictly convex functions in two dimensions. Some efforts of improving the original Nelder-Mead method are found in [LLR04, Tse99, PCB02, NT02]. In these studies, strategies for solving the defects of original Nelder-Mead simplex method include using a sufficient decrease condition for acceptance of a new vertex rather than a simple decrease and restarting if the simplex becomes ill-conditioned, e.g., degenerated into subspace.

Simplex methods

{ if (x k < x min k ), then x k = x min k if (x k > x max k ), then x k = x max k ( 6
On the other hand, with its simple implementation and effectiveness in most of practical applications, Nelder-Mead simplex became the most known direct search method. In december 2014, Google Scholar displayed more than 1200 papers published in 2014 that referred to the Nelder-Mead method for a wide variety of applications.

The algorithm is often found in optimization numerical software packages [START_REF] Lewis | Direct search methods: then and now[END_REF] Rosenbrock's method conducts a series of stages. For each stage, a set of search directions is defined. These search directions are updated from stage to stage to make use of the information acquired about the objective function.

For the initial stage of Rosenbrock's method, the search directions are defined with the coordinate directions. In Figure 6.23 [VB], an example of Rosenbock's algorithm is illustrated for 2-dimensional minimization problem.

The algorithm search for improvement of the objective function in axis x 1 , x 2 until it finds at least one successful (points 5, 7, 9, 11, 13 for direction of axis x 1 and points 2, 4, 6, 8, 10 for direction of axis x 2 ) and one unsuccessful step in each of the search directions (the points 3 for axis x 1 and 12 for axis x 2 ). Then, the stage terminates. If there exists new point that gives the lowest of function value (point 13), then it is used as initial point for the next stage.

For the next stage, instead of using the same search directions, new search directions is generated by rotating the current ones to adapt with the function curvature. For that, the connection between the initial point of the previous stage (point 1) and the new best point (point 13) is taken as a new search direction, with its orthogonal direction, they constitute new set of search directions (defined by points: 1, 13, 18). It should be noted that search directions are orthogonal and therefore they are linearly independent. The search continues and converges to the minimizer of the objective function.

Powell's method uses function modeling as the basic approach. The main intention is to be able to use conjugate directions to accelerate the convergence if the method is applied to a convex quadratic function. Like Rosenbock's The first n searches are executed along each of linearly independent directions (from P 0 to P 1 and from P 1 to P 2 as illustrated in Figure 6.24 [START_REF] Flannery | Numerical recipes in c[END_REF]). The last search (n + 1) is conducted along the direction connecting the point obtained at the end of the first n searches (point P 2 in the example) and the starting point (P 0 ). The search is accelerated finding new minimum point of the stage (point P 3 ). At the end of the stage, one of the first n directions is replaced by the last one (connection P 0 P 2 ). In a new stage, the newly found minimum point of the previous stage is used as the starting point, and the whole steps are repeated to search for the minimizer of the objective function. As illustrated in Figure 6.25, the optimization by MPC strategy for the reheating furnace is also a simulationbased method as for the existing method. The part of setpoint optimizer in the system includes mainly furnace simulation block, and optimization algorithm block. The furnace simulation calculates the objective function J o (k) via the non-steady-state furnace model. As formulated in the forgoing sections, J o contains objective and penalty terms concerning accuracy of final slab temperature, energy consumption and control constraints. To calculate these terms, for sample time k, the simulation of furnace predicts the final slab temperature profiles and furnace energy consumption over the prediction horizon [k, k + N p ]. This simulation is executed repetitively from given control input sequence U (k) and the estimations of system state (current temperature profiles of slabs, zone temperatures, etc.). Control input sequence is varied according to the chosen optimization algorithm. Once the control input sequence is optimized according to a given tolerance, it is applied to the real furnace.

Adopted optimization strategy

On the other hand, a furnace observer takes into account the measured heating power supplied to the furnace zones in order to estimate current temperature profiles of all slabs inside the furnace. Also, this observer estimates the system outputs which are the final temperature profile of the slabs.

J o (U (k)) is a nonlinear function with the following properties:

• Calculation of J o (U (k)) is very time-consuming.

• Exact first partial derivatives of J o (U (k)) cannot be calculated.

• Numerical approximation of the gradient of J o (U (k)) is impractical and also time-consuming since it requires evaluations of J o (U (k)) at different points.

• Values of J o (U (k)) are noisy due to numerical calculation of heat radiation which can be slightly different from one simulation to another.

In application for considered reheating furnace, a highly accurate solution is not practically desired because of high cost simulation. On the other hand, it may be impossible because of uncertainties and errors in the model or data. Actually, what we want is improvement of reheating furnace rather than optimization. Therefore, the crucial issue is finding a better answer quickly for each iteration; asymptotic convergence properties are not relevant.

Among the presented direct search methods, except simplex search methods, the convergence property of algorithms is assured. However, while simplex search methods only need no more than n + 1 evaluations of objective function to complete the search for improvement of an iterate, the other methods require from 2n to 2 n evaluations to do the same thing [START_REF] Lewis | Direct search methods: then and now[END_REF].

Consider the requirement of the application and characteristic of the available optimization algorithms, we choose Nelder-Mead simplex method to solve the minimization problem in the MPC strategy for level 2 of the reheating furnace. For that, we expect to obtain good control input sequence in a limited number of function evaluations.

The design of Nelder-Mead simplex method for reheating furnace case consists of determining the initial simplex and executing a series of reflection, expansion, contraction, shrink operation. For sample time k, the initial simplex can be constructed from the solution of time k -1, i.e., U 0 (k) = U (k -1). As in [START_REF] Luersen | Globalized nelder-mead method for engineering optimization[END_REF], a simplex of size a can be initialized at P 0 based on the following rule

U i = U 0 + ρe i + n ∑ j=1,j̸ =i (qe j ) for i = 1, .., n, (6.61)
where e i is a unit base vector and ρ, q determined by

ρ = a n √ 2 (√ n + 1 + n -1 ) q = a n √ 2 (√ n + 1 -1 ) . (6.62)
The size of the simplex a is taken between 2% and 10% of the smallest domain dimension. This latter parameter is in fact determined by the difference between upper and lower limits of control inputs. For the considered furnace, the differences between upper and lower limits of zone temperature set-points are 90

U i = ϑ[U i + U l ] for i = 1, n (6.66)
where shrink coefficient ϑ = 1/2. If the convergence tolerance is satisfied then stop, otherwise go to step 1 and the algorithm restarts

The convergence criterion used to halt the procedure is calculated as √ 1/(n + 1) ∑ n i=1 (y i -y) 2 . That is, the algorithm terminates the iteration when the simplex become sufficiently small.

Numerical result

In order to validate the adopted MPC strategy on the considered reheating furnace, we simulate the furnace operation by constructing a furnace simulator from the non-steady-state model presented in Chapter 3. The data used to validate the simulator is provided from the database of the considered furnace. From the scheduling table, the information of the heated slabs, fuel composition, etc., and the zone power set-points calculated by the level 1 controllers, the furnace simulator calculates instantaneous temperature profile of slabs, exact position of each slab, the zone temperatures, and also the final slab temperature profiles as depicted in Figure 6.26. At each sample time, the set-points optimizer calculates the optimal zone temperature set-points based on the current position, temperature of slabs, and the current zone temperatures. In addition, to compare MPC strategy and the existing This section presents firstly the simulation environment, operational scenario of the furnace, and the validation of the simulator against the industrial data. Specific tuned parameters of MPC strategy will then be given. Consequently, the results of MPC strategy will be analyzed and compared to the results of the existing method.

Simulation conditions

The whole simulation is coded in Intel Fortran 90 to facilitate the programming task because the non-steadystate model is originally developed in Fortran 90. The structure of simulation program is outlined in Figure 6.26. One may ask why the level 1 controllers are simulated as PID instead of distributed MPC which has been studied and implemented. It is because the available data used in the simulation was recuperated when the furnace was still running with PID controllers. In addition, from the point of view of level 2, control of level 1 is considered transparent. Thus, the choice of controller for level 1 will not affect controller level 2. The simulation Chapter 6. Slab temperature control 136 environment captures the whole cascaded control system described in Section 6.1 (cf. Figure 6.2). The main difference is that the observer is not required because all the system variables are available from the furnace simulator itself. The simulation environment has been installed on a standard PC with a 2.5GHz CPU and 3 GB of RAM. On this machine, to simulate 1 h of furnace operation the non-steady-state model needs 5 s. This is considerably more time-consuming than steady-state mode which is around 10 5 times faster. We have chosen the prediction horizon of ∆t P =2 h and sample time ∆t s =5 min according to the average reheating time of slabs and time interval of loading and discharging events. Thus, predictions of slab temperatures over the prediction horizon will need 10 s.

On the other hand, we want the simulation environment to be implementable for real-time operation. Therefore, the optimization of MPC controller has to deliver the optimized zone temperature set-points in less than a period of one sample. As a result, the maximum number of function evaluations we have is l max = ∆t s × 60/10 = 30.

The heating scenario uses database of three-day furnace operation in which around 900 slabs are heated and The reheating times of slabs t j,1 -t j,0 , or in other word, the residence times of the slabs inside the furnace are given in Figures 6.29 and 6.30. The average reheating time is around 2.5 h. Besides, we can see slabs such as the slabs j = 210 to j = 247 of which residence times are increased to around 3 h. This is because the furnace has been stopped for 30 mn when it contained slabs j = 210 to j = 247. In a longer stop, reheating time can be dramatically increased to 4 h and 5 h as the case of slabs j = 600 to j = 641. These long residence times may lead to overheating and unwanted damage of the material quality. The desired final slab temperatures are presented in Figures 6.33 and 6.34. For all the slabs, the desired temperature profile Tj,end (y j ) is homogeneous with respect to y j . And, the bounds on the final temperature profile are determined by T j,end,max = Tj,end + 20 • C and T j,end,min = Tj,end -20 • C. In the Figures 6.33 and 6.34, these bounds are indicated by two limits of the grey bars. In practice, these bounds could be relaxed according to operating scenario, T j,end,max in particular because over-heating is less restricted than under-heating which causes rejections of slabs in hot rolling mills. Generally, we can see that desired temperatures of slabs are highly irregular with exception of slabs j = 618 to j = 687. The irregularity between desired final slab temperatures will make the control task more difficult to achieve high temperature accuracy. Temperature range is from 1170 • C, as of slabs j = 30, l = 31, to 1230 • C as of slabs j = 57, j = 58. Due to slow thermal inertia of the furnace, it would be desirable that consecutive slabs have as much overlap between the allowed final temperature ranges [T j,end,min , T j,end,max ] as possible. We see that most of slabs have some required temperature overlaps with their neighboring ones, with exception of some slabs like j = 31 and j = 32 or slab j = 295. In fact, the arrangement of slabs according to their desired final temperature is one of important task of level 3, which, however, is not in the scope of this dissertation. The upper bound of slab temperature is T j,max = 1350 For the simulation, all slabs have initial temperature close to ambient temperature which is fixed at 20 • C. At time t=0 h, the furnace is assumed to have inside it 31 (average number) fictive slabs with standard dimensions of (700 x1860 x235) mm and standard time interval between discharging events of 5 min. These fictive slabs are given with a certain temperature profile according to its position inside the furnace. The slabs of real system are subsequently entered the simulator for the rest of the simulation. On the other hand, we should also look at the consumed power calculated by the simulator in comparison to the measurement of the real system. In Figure 6.37, we can see that the simulator has been able to reproduce the dynamic of the measurement although some static errors are present. As given in Table 6.5, the average consumed power calculated by the simulator is 97.5 MW which is close to average measurement value of 99 MW.

Simulation results

Zone temperature set-points calculated by the MPC controller are shown in Figures 6.38 The black solid curve is the simulated zone temperatures following set-points of MPC controller given by the PID controllers of level 1. We see generally the set-points optimized by MPC controller respect their bounds.

The temperature of the soaking zone is given in Figure 6.38. We see that set-points calculated by MPC are lower than those given by the existing controller. The PID controllers are, in most of time, able to follow set-points provided by MPC controllers. Due to upper limits of zone powers PID may sometimes fail to attain the set-points like at the intervals [2 h, 3 h] and [25 h, 26 h].

For zone heating 2, temperature set-points calculated by MPC are generally higher than values given by the existing controller as shown in Figure 6.39. The level 1 controllers succeed to follow the temperature set-points.

In zone heating 1, the set-points given by MPC are mostly lower than that of the existing controller, especially for the time intervals Unlike other zones, we see clearly that PID controllers of preheating zone fail to follow the set-points given by the MPC controller which are a lot lower than values calculated by the existing controller. In fact, the temperature of this zone is not only affected by its fuel flow which has lower bound greater than zero, but also temperature of other zones through the flow of exhaust gas. Therefore, low temperature set-points are not followed by PID The average zone temperature set-points are given in The average final slab temperatures mean y j T j,end (y j ) are shown in Figures 6.42 On the other hand, existing controller often overheats the slabs like for the slabs j = 102 to j = 182 in Figure 6.42 or the slabs j = 630 to j = 695 in Figure 6.45. This tendency is the preventive strategy adopted in the existing controller to avoid slab under-heating and slab rejections in the hot rolling mill. However, it may become excessive and conservative as for slabs j = 755 to j = 766 in Figure 6.45, in these cases energy is wasted. MPC controller, on the contrary, gives less slab overheating and less deviations from the desired temperatures, as example of slabs j = 75 to j = 95 in Figure 6.42 or slabs j = 752 to j = 782 in Figure 6.45.

The slabs that have high residence time compared to the average due to a long furnace stop, like the slabs j = 600 to j = 700 when a stop of 2 h occurred (c.f. Figure 6.30), are often largely overheated with the existing controller as shown in Figure 6.45. The MPC controller, in contrast, is able to reduce the overheating caused by long stop.

More specified statistics are given in Figure 6.46, we have the histograms of final slab temperature errors (difference between calculated and desired values) obtained with MPC strategy and with the existing controller. On average, we see that both methods slightly overheated slabs. And the overheating is reduced from 13.7 • C with the existing controller to 7.4 • C with MPC controller as specified in Table 6.7. This is a crucial improvement which may contribute to the reduction of overall energy consumption of the furnace. In addition, standard deviation of temperature errors are reduced from 10 • C to 8 • C. This indicates that MPC strategy is more accurate that the existing controller. To measure the success of the heating, we consider the bounds of temperature errors which is the range of -20 • C to 20 • C. For MPC strategy, 92% of slabs is successfully heated while the existing controller only obtained the results of 75%. More success on final slab temperature leads to less slab rejections at the downstream process, e.g., the hot rolling mills. In conclusion, compared to the existing controller, the MPC Homogeneity of slabs is reflected with the final temperature gradients inside slabs as shown in Figure 6.47.

The red curve represented final slab temperature gradients simulated with MPC controller, the blue one are those simulated with the existing controller. In general, MPC controller gives lower temperature gradients with exceptions as of slabs j = 700 to j = 747. Some slabs have particularly more homogeneous temperature with MPC controller than with existing controller such as slabs j = 40 to j = 95 and slabs j = 200 to j = 238. The slabs j = 600 to j = 700 are affected by a long stop of 2 h, therefore, their final temperature gradients are much smaller than other slabs in the results of both of controllers. The improvement of MPC controller is confirmed Finally, a general observation for the whole simulation of 900 slabs is given in Figure 6.61 and Figure 6.62.

These figures show the average temperature trajectory and the temperature gradient of 900 slabs. With the given information of these figures, we can confirm that the MPC controller heat the slabs more gradually than the existing controller. The average temperature of slabs and the corresponding gradient is reduced in general and therefore energy can be saved.

The energy consumption of the furnace zones simulated with two controllers are shown in Figure 6.63. The red columns represent results of MPC controller and the blue ones are results of the existing controller. For two heating zones the MPC controller requires more fuel power than the existing controller. On the other hand, power of other zons, especially preheating zone, are higher when simulated with existing controller than with the MPC controller. MPC controller requires nearly no supply of fuel for the preheating zone which means that preheating zone is heated mainly by the exhaust gas. This can also explain why the MPC controller consumes for in the zone heating 1, since the exhaust gas heating the preheating zone come directly from the zone heating 1. Total fuel consumption of the furnace is shown in Figure 6.64 for the representative time interval [33h 42h]. We see that with MPC controller the furnace fuel consumption is more smooth, and in general lower than that simulated with the existing controller. The specific energy consumption, i.e, consumed energy per ton of steel heated is given in Table 6.9. The provided numerical values are normalized with a nominal value of specific energy consumption relating to the considered furnace. For the whole simulation of three-day furnace operation, the MPC controller

Conclusions

The MPC strategy is designed for the level 2 of reheating furnace. The problem is the control of slab temperature so that it attains the desired final value. Firstly, the solution consists in formulating an optimization problem which is subject to furnace constraints. The cost function of the optimization problem takes into account the error of final slab temperature, temperature gradient inside slabs, energy consumption of the furnace, and general objectives such as reduction of furnace temperature level and slab temperature level. The furnace constraints concern the limits on zone temperatures and their rate of changes. Constraints on temperature of slabs are maximum temperature inside slabs, maximum temperature gradient of slabs, bounds on final slab temperature and final temperature gradients. These constraints can be accounted as penalty terms.

Secondly, the dynamic optimization required by MPC strategy is based on the simulation of furnace by a computer program in which the non-steady-state model is implemented. This model demands considerable computational resources. Therefore, gradient-based optimization methods are not those of interest because they need either explicit or approximation of derivatives of the cost functions. On the contrary, direct search methods need only the evaluations of cost function. The Nelder-Mead simplex method is well-known for its simple implementation, flexibility with constraints, here proved to be efficient for the design of MPC level 2 of reheating furnace. The simulation results on industrial scenarios show that energy consumption can be reduced up to about 5%, with higher homogeneities and accuracy of final slab temperatures.

In the next chapter, we reconsider the scheduling of slabs in order to boost the furnace productivity while assuring the heating requirement of slabs. The MPC controller proposed in the present chapter will be the basic of the controller constructed for scheduling task. Beside the criterions on accurate final slab temperatures, low temperature gradient and energy efficiency, the maximum throughput rate of furnace is the new criterion which is equally expressed by finding the minimum time interval between slab discharging events. A new control input added to the optimization is the average time interval between slab discharging events which is lower bounded according to the requirement of the hot rolling mill.

Part III

Improvement and Perspective

two zones, the nonlinear model predictive strategy utilizing zone temperatures and forwarding speed as control input is presented. The constraints of upstream and downstream are not taken into account.

Another nonlinear predictive control scheme for furnace control is presented in [START_REF] Balbis | Nonlinear predictive control of steel slab reheating furnace[END_REF]. The control inputs are the fuel flows of furnace zones and the conveying speed. The developed method firstly finds the steadystate conditions that maximize the production profit (productivity, energy consumption, heating quality). Then a linearized model around the steady-state point is used to construct a dynamic optimization. For this later problem, the objective function accounts for the furnace constraints, and the deviation of the system from the steady-state optimum (optimized speed and fuel flows) found by the static optimization. The obtained simulation results

showed that the fuel consumption can be reduced and throughput rate of furnace is increased. This is because the system can respond more quickly to changes in material flows and delays.

Recently, a method for slab scheduling is developed in [START_REF] Steinboeck | Optimized pacing of continuous reheating furnaces[END_REF] to find the minimum reheating time of slabs. A switched analytic model of reheating furnace is able to compute enthalpy absorbed by slabs. Based on this model, a quadratic programming problem is formulated taking into account furnace constraints. The optimal solution is the smallest reheating times needed to attribute sufficient heat for slabs. An example scenario in which three furnaces operate in parallel obtained significant results like: higher accuracy of final slab temperature, lower temperature gradient and increased throughput rate.

In summary, most of the published works utilized either furnace temperature or slab velocity as control inputs.

Only few combine these two strategies. We want to contribute to the filling in of this gap by utilizing the MPC method. In this chapter, based on the numerical non-steady-state model presented in chapter 3, we design MPC strategy to improve at the same time furnace productivity, heating quality and energy consumption. The control inputs are time gaps between discharging and loading events of slabs, and the furnace zone temperature set-points.

In the following, constraints and objective of scheduling task in reheating furnace are clarified.

Problem statement

Boosting productivity of the furnace and optimizing furnace temperature control cannot be done as separated tasks. This is because maximum reheating power of a furnace is in fact a dynamic quantity, which depends, at each sampling time, on furnace load and the current furnace thermal state (slab, walls and gas temperature, etc.). As a result, the optimization of furnace scheduling must take into account the heating requirement, furnace constraints as well as productivity objective and constraints of hot rolling mill, i.e.,

• The temperature trajectories of slabs are upper bounded by T j,max .

• Final slab temperature profiles are limited with T j,end,max and T j,end,min .

• Trajectory of temperature gradient of slabs are upper bounded by ∆T j,max so that slab deformations are avoided.

Chapter 7. Boosting of productivity 166

• Final slab temperature gradient should be as small as possible and stay below a certain required value ∆T j,end,max .

• Maximizing the throughput rate of processed steel slabs. This implies that forwarding velocity should be as high as possible, or time intervals between consecutive discharging events should be as small as possible.

• Furnace is limited with certain productivity due to constraint of hot strip mill. Therefore, the time interval between discharging events are limited with a lower bound δt sch,min .

• Slabs are often arranged into montage according to its material properties, with different montage the cylinders of hot rolling mill need to be changed; so a space between these sets of slabs is imposed giving time for necessary replacements. Therefore, the time interval between these montage will also be imposed meaning not dependent on δt sch .

• Stoppages of furnace can be scheduled or unscheduled events, and they are unavoidable and need to be included in the simulation.

Problem formulation

The sequencing of slabs through the furnace is determined with a table as example shown in Figure 7.1, in which time intervals between consecutive loading and discharging events are predefined. Therefore, the smaller these time intervals are the more throughput rate of the furnace is increased. Because the furnace can only support certain level of slab load, the discharging and loading are usually synchronized, meaning if one slab is discharged, another slab at the furnace entry would be loaded. As a result, the average time interval between discharging and loading of consecutive slabs are in general equal. Consequently, we can use either the time intervals of slab discharging events or those of slab loading events as extra control inputs along with zone temperature set-points.

Considering the complexity of the non-steady-state model that we want to use for the optimization of MPC controller, it would be preferable that the number of input variable has a low dimension to reduce number of furnace simulation needed. Therefore, we should use only an average value of the time intervals as the extra control input. More specifically, for each furnace simulation, we can redefine the slab sequencing by replacing the current values of time interval in table of Figure 7.1 with an value provided by the optimization algorithm.

Because the slab scheduling is reconsidered each sampling time, reheating time of a slab can change constantly.

As a result, if we use a fixed prediction horizon for the optimization problem, the number of slabs discharged of the furnace can change from one furnace simulation to another provided with the different time intervals.

For example, if the time interval of discharging is small then more slab will be discharged during the prediction horizon period and inversely, if the time intervals are long then there will be less slabs discharged in the same period of time. In the formulation objective functions, number of slabs discharged is a crucial quantity which affects the calculated results. We shall want the number of discharged slabs during the prediction horizon to be constant. Thus, MPC strategy will be applied with a variable prediction horizon. The following variables are considered for scheduling optimization tasks:

• Average time interval of discharging events δt sch (k). And obviously, δt sch (k) ≥ 0

• Lower bound of time interval of discharging events δt sch,min (k)

• Sampling time of level 2: ∆t s

• Number of slabs considered in the prediction horizon is a constant number N s , whereas the prediction horizon [k, k+N P (k)] defined with N P (k) = ⌈N s ×δt sch (k)/∆t s ⌉ is time-varying dependent of δt sch (k).

• Fixed control horizon [k, k + N C -1] with N C ≤ N P (k)
• Control input u(k) is now defined by the zone temperature set-points and the average time interval of discharging events, i.e.,

u(k) = [δt sch (k), T z,1 (k), T z,2 (k), ..., T z,M (k)] (7.1a) = [u 1 (k), u 2 (k), ..., u M +1 (k)]. (7.1b) • Control input sequence U (k) = {u(k), ∀t ∈ [k, k + N C -1]}. It shall be noted that u(k) = u(N C ), ∀k ∈ [N C , N P ].
• Loading and discharging events of a slab j correspond to: k j,0 and k j,1 , respectively, with

k j,0 = int(t j,0 /∆t s ) (7.2a) k j,1 = int(t j,1 /∆t s ). (7.2b)
In comparison to the MPC method proposed in Chapter 6, the difference of MPC constructed in this chapter dedicated to the boosting of productivity is that we use one more variable concerning slab sequencing δt sch (k) and the prediction horizon N P is a time-varying parameter.

The objective and constraint concerning throughput rate are added along with the heating objectives and constraints that have been formulated in Chapter 6. Especially, the weighting coefficients used in the objective functions are chosen by the same approach defining the priority of each slab according to its position inside the furnace.

The requirement of final slab temperatures are formulated as the sum of weighted final temperature errors of slabs that are discharged during the prediction horizon, i.e.,

J 1 (U (k)) = N P (k) ∑ i=1 w 1 (k + i) × |T s (k + i) -Ts (k + i)|, (7.3a) = ∑ j∈S 1 (k) w 1 (k j,1 ) × |T s (k j,1 ) -Ts (k j,1 )|, k j,1 ∈ [k, k + N P (k)]. (7.3b)
And, objective of homogeneity of slab is the sum of weighted temperature gradients as

J 2 (U (k)) = N P (k) ∑ i=1 w 2 (k + i) × g(k + i) (7.4a) = ∑ j∈S 1 (k) w 2 (k j,1 ) × g(k j,1 ), k j,1 ∈ [k, k + N P (k)]. (7.4b)
The criterions on energy consumption is expressed as

J 3 (U (k)) = N P (k) ∑ i=1 w 3 (k + i) M ∑ j=1 P j (k + i). (7.5)
And, general objective function of reducing furnace temperature is formulated as

J 4 (U (k)) = N C ∑ i=1 w 4 (k + i) M +1 ∑ j=2 u j (k + i). ( 7.6)
The general objective to reduce slab temperature trajectory is written as

J 5 (U (k)) = N P (k) ∑ i=1 [x ′ 1 (k + i) × w 5 (k + i) × x 1 (k + i)]. ( 7.7)
In particular, objective of maximum throughput rate is translated into minimum average time interval between discharging events. That is, the criterion on scheduling objective can be given as

J 6 (U (k)) = N C ∑ i=1 w 6 (k + i) × u 1 (k + i) (7.8)
where w 6 (k) is a time-dependent coefficient defined with

w 6 (k + i) = N C -i + 1 N C (7.9)
The constraints of zone temperature is now expressed as This is the extra constraint added to the heating problem formulated in chapter 6. If the time interval is lower than the limit, the furnace discharge too fast and the downstream process, e.g., will not be able to tolerate it.

Therefore, the penalty term J c,7 (k) will have a positive value every time the upper limit of time interval is not respected.

The scheduling optimization task is constructed as an unconstrained minimization problem given by min

U (k) J o (U (k)) = 6 ∑ i=1 α i × J i (U (k)) + µ 7 ∑ i=1 J c,i (U (k)) (7.16)
The weighting coefficient µ usually has high value to penalize the cost function when constraints are violated.

Other coefficients α i (i ∈ [1,6]) are chosen according to the importance of the objectives. Most important objectives are success of desired final slab temperatures and low temperature gradients. In second places, there are objectives of maximizing production rate and reducing energy consumption (reducing zone temperature and fuel supply). Final slab temperature gradient is shown in Figure 7.9. We see that in general the final slab temperature gradients simulated with MPC controller are more regular between slabs compared to the existing controller. The average temperature gradient of slab simulated with MPC controller is reduced of 7 • C with a standard deviation decreased from 22 • C to 7 • C as shown in Figure 7.10 and Table 7.1. In term of throughput rate, different results are presented in Figures 7.11,7.12 and Table 7.2. It should be noted that the estimation of throughput rate is synchronized with discharging events and not with time because, at each optimization, MPC controller modifies the scheduling of the slabs. The throughput rate can be increased by 15 t/h with MPC in comparison to the existing controller. The total heating time and average reheating time are also reduced. The price of higher production rate is higher energy consumption which is increased from 1.25 Gj/t to In term of reheating time of slabs, the results obtained with MPC strategy is compared with that of the real data as shown in Figure 7.15. In general, with MPC strategy, reheating times of slab are reduced. For both results of MPC and the real data,when the reheating time of slab is increased suddenly, it means that there is a stop. A short stop of furnace can be about 20 min. And, a long stop is more than 50 min. We can see that when a short stop happens MPC does not reduce reheating time of slabs. However, when a long stop happens, MPC strategy can reduce significantly the reheating time of slabs. This means MPC reacts faster to the stops than the current method. For this reason, the production rate of the furnace is increased as shown in Table 7.2.

Conclusions

This chapter has widened the furnace control from a predefined scheduling strategy to a variable scheduling strategy where MPC method is utilized to optimize the zone temperature set-points, and at the same time to find the maximum possible throughput rate of the furnace. The MPC controller designed for boosting of productivity is basically developed from the solution constructed in the chapter 6. Heating objectives and constraints are the same as those formulated of predefined scheduling.

Beside the zone temperature set-points, the time interval of discharging events is added as a control input. The latter variable actually is the direct index of furnace throughput rate since the faster slabs are discharged at the required temperatures the higher furnace productivity is. Therefore, an new objective of control is to minimize the time interval between discharging events. On the other hand, the hot rolling mill is limited with a upper bound. As a result, an additional constraint is represented by a lower bound of the time interval.

When the new control input, the time interval between discharging events, is varied from one furnace simulation to another, the number of discharged slabs over a fixed prediction horizon changes. Therefore, the calculation of objective functions will be affected since they account for final temperature and temperature gradient of every individual discharged slab. As a result, instead of utilizing a fixed prediction horizon, we used a variable horizon changing according to the time interval while the number of discharged slabs over the prediction horizon is constant.

The design is tested with the industrial data of three-day furnace operation. Simulation results show that the furnace throughput rate can be increased up to 15 ton per hour while furnace consumption is also increased of 0.02 Gj/ton. Considering the economical benefit of increased throughput rate compared to the price of an energy unit, it is usually preferable to push the furnace to its maximum capacity. In term of heating requirement, the MPC controller shows to be more accurate over the final slab temperature than the existing controller and especially the slab temperature gradients are reduced significantly.

Chapter 8 Extension to other metallurgical furnaces

Besides slab walking-beam reheating furnaces, other types of process in which steel products are heat treated are also big energy consumers of a steel plant. These processes are often divided into zones which are controlled at different level of temperatures. The application of proposed distributed MPC algorithm for control of zone temperatures may lead to reduction of energy consumption and improvement of product quality. In this chapter, we present the furnaces for continuous annealing, continuous galvanizing and the hot blast stoves in which we want to use the MPC strategy to improve heating quality and reduce energy consumption. Continuous annealing is the process of heating rolled strip in a sequence of zones to a predefined temperature according to steel grade and strip dimensions. The final result of the process is an increased of ductility and elimination of strains in the strip that lead to failure in subsequent processes. Figure 8.1 shows the steps of continuous annealing process. Because the total processed strip can be several kilometers long, the furnace is vertical, withe the strip making several traverses within each zone. Firstly, the strip goes through cleaning section which is dedicated to removal of all dirt, grease, mill scale, and other oxides form on the surface of steel. Accumulators act as buffers between static steel coils (at entry and exit) and the continuous strip running through the furnace. As the entry coil is emptied and is replaced with a new full coil and both strips welded together, the entry accumulator unwinds to provide continuous strip. Similarly, the exit accumulator winds up while the exit coil is full and is unloaded at the exit end.

The annealing furnace is divided in different stages of heating and cooling as shown in Figure 8.1. These stages are necessary to give steel the desired properties by heating to particular temperatures and profiles thereby determining the grain structure within the metal, and to prepare for further processing such as galvanizing. In the heating zone, the cold strip is heated to the highest temperature of the annealing profile. Because of the danger of oxidation and flame length hot spots to the strip, the heat is radiated to the strip using ceramic radiant tubes, themselves fired by gas.

The soaking zone is dedicated to maintain strip temperature at 700 • C to 800 • C. and is electrically heated. The first of the cooling section is the slow cooling chamber, where two variable speed blowers circulate the atmosphere gas through two water-cooled heat exchangers. The second of the primary cooling sections is the fast cooling chamber, where the strip is cooled by a variable speed fan and the cooling length is determined by the positions of three dampers along its vertical length. This chamber is also split vertically for temperature profile control across the width of the strip using a scanning pyrometer. It is desirable for productivity to maintain strip movement through the furnaces at all times. However, it is inevitable that faults occur and the strip stops. At this stage, there is a risk of the strip being overheated and breaking, causing major disruption. The control system, therefore, looks ahead and takes action at low line speed, line deceleration and line stop. These essential actions are to minimise the heating furnace firing and maintain atmosphere gas pressure due to the cooling strip.

To prevent oxidation of the strip, the atmosphere around the strip inside the furnaces is a controlled mixture of H2 and N2 (cracked ammonia NH3 providing 5% H2 and 95% N2) although hydrogen only is sometimes used.

The furnace themselves are sealed gas tight and each of them is maintained under slight positive pressure.

Continuous galvanizing

Galvanizing is the practice of immersing clean, oxide-free iron or steel into molten zinc in order to form a zinc coating that is metallurgically bonded to the iron or steel's surface. The zinc coating protects the surface against corrosion by : 1) shielding the base metal from the atmosphere; 2) providing a cathodic or sacrificial protection since zinc is more electropositive than iron or steel.

In fact, the main difference of galvanizing process with annealing process is that it have a step of galvanizing with zinc pot as shown in Figure 8. The contiuous galvanizing processes a different type of product compared to that of continuous annealing process.

The final products have more value-added than that of continuous annealing process. The strip is heated in the furnace and then immersed in zinc bath. The final step consists of quenching, removing excess zinc and inspection. The major control involves Strip temperature, bath immersion time, bath withdrawal rate and steel cooling rate, which is more exigent than the control in furnace of continuous annealing process.

Furnace control

The control system of furnaces in the continuous galvanizing and annealing is quite similar to that of the slab reheating furnace. There are two main levels. The level 1 controls the temperature of the zones by varying the fuel flow rates, and the level 2 calculate the temperature set-points of furnace zones in order to achieve desired temperature of strip. The furnace zone temperatures are usually controlled by PID controller as depicted in Figure 8.4. Temperature of furnace zones are coupled with each other. In furnaces of annealing and galvanizing process, number of controllable zones are greater than that of reheating furnaces. Another important difference is that in continuous galvanizing and annealing radiant tubes are used as indirect heat sources while reheating furnaces use gas burners which is direct heat sources. In addition, the gaseous environment inside furnace of annealing and galvanizing is highly controlled and monitored to avoid surface reaction of steel strip and there is no circulation of gas. Therefore, the coupling effects between zones are not due to gas circulation but mainly to the flow of hot strip.

The control challenge of furnaces annealing and galvanizing is that it requires high control precision of strip temperature since temperature error constraints can be very tight (5 • C). On the other hand, the strip is transported through the furnace at high speed (500 m/s). As a result, the dynamic of these furnaces is much faster than that of the reheating furnaces.

In order to successfully apply the proposed hierarchical MPC design in the galvanization furnaces, one must developed highly accurate coupling model of the furnace zones. The sampling period must be smaller than that utilized in reheating furnaces. we can see that coke and iron ore are grinded and mixed together. The mixture of coal and iron ore is then transfered and loaded into blast furnace from the upper part. Hot combustion air (enriched with oxygen), also called blast air, is supplied to the blast furnace from the lower part at high pressure and at an averaged temperature of 1200 • C. To do the task the hot blast stoves firstly burn a gas mixture to generate the necessary heat. Then, the heat generated from combustion process is stocked in a solid structure called the checkerwork by heat transfer as shown in Figure 8.6. The latter solid structure is constructed from different layers of refractory brick of different types. One typical setup of different type of refractory bricks from top to bottom is mullite and silica, super-duty and high-allumia. These layers contain several flues so that air or gas can flow through. After the heat is stocked in the solid mass, the cold combustion air needed for blast furnace is now blown into the hot blast stoves and then heated by the hot solid mass. Therefore, a hot blast stove functions through two alternate cycles of heating and cooling, or also called gas phase and blast phase referring to periods when the stove is filled with gas and blast air, respectively. Specifically, these phases can be described as follows:

• Gas phase: the flow of gas in this phase is shown in Figure 8.6. The stove is heated by combustion of fuel gas in the combustion chamber of the stove. The wast gas which is the combustion product passes through the dome (top) of the stove and then goes down through the refractory bricks. At the end of the gas phase, most energy of combustion has been absorbed by the refractory layers. Therefore, temperature of waste gas is kept at low level which is around 300 • C.

• Blast phase: in this phase, cold blast air is blown into the stove, the path of blast air is in the opposite direction of waste gas in gas phase. Cold blast air firstly is heated by the hot refractory layers at the bottom of the stove. Subsequently, it goes up to the dome of the stove then delivered to the blast furnace as hot blast air of 1200 • C. After gas phase, there is a changeover phase which is about five minutes long. The main task in this phase is to clean the impurities and dust produced from the waste gas. For that, a small flow of cold wind is blown into the combustion chamber. After the blast phase, other changeover phases are used to rebalance the pressure inside the stove, especially inside the refractory brick structure. Temperatures of the refractory layers have to be restricted The difference between temperature of fluid and solid is reduced from the bottom to the top of the stove. One hot blast stove can only blow hot blast during its blast phase whereas the blast furnace needs a constant supply of hot blast. So several stoves must run in dephased conditions to guarantee a constant hot blast supply to the blast furnace. The minimum number of stove to ensure it is two but generally the 3 stoves configuration is the most utilized. For this configuration there is always one stove on blast phase while the two others are on gas or on changeover phase as shown in Figure 8.8.

Control task

The control task of hot blast stove is to determine the minimum amount of fuel necessary during gas phase to achieve the blast air flow rate and temperature requirements for the next blast phase. It is wanted that at the end of the blast phase, cold blast air flow is zero. However, in practice, the later flow can not precisely reach zero value, so a small amount of additional heat is put into the stove to ensure that the temperature and flow rate requirements of the hot blast are met. The additional heat can be measured by the by-pass valve position at the end of the blast phase.

Control constraints will be firstly the upper and lower bounds of fuel flow rates at the gas phase. Then, upper and lower temperature limits of the dome, the interface and the bottom of the stoves. These latter values are given in As the fuel consumption of the hot stoves represents from 10% to 15% of the total energy consumption of a steel plant, the economical and environmental challenges for control are considerable. The MPC strategy can be considered in order to optimize the fuel flow rate profile of the hot stoves during the gas phase. The control is located at high level of the system and therefore is similar to the level 2 of the slab reheating furnace, because what are needed to be optimized are the fuel flow rate set-points given to subordinate levels. Currently, elaborated numerical models of hot stoves have been developed by ArcelorMittal R&D. These models can be used for the predictions of different process variables. Based on these predictions, an objective function can be formulated and computed taking into account the fuel consumption, constraint of fuel flow rates, the temperature of the delivered hot blast, temperature and constraints of different part of the stoves. Then optimal profiles of fuel flow rates in gas phases of the stoves are calculated via an appropriate optimization algorithm. One particular point is that the stoves function in cycles, the ratio between gas phase and blast phase can also be an important variable affecting the performance of the stoves which can be considered as input to the optimization problem.

Chapter 9

Conclusions and perspectives

Summary

The objectives of this work were to develop the MPC strategy under a hierarchical distributed architecture of a slab walking-beam reheating furnace. The function of slab reheating furnaces is to bring temperature of steel slabs from an initial temperature to around 1250 • C with a high homogeneity before entering hot rolling mills.

The modeling of the furnace is crucial for observation of process variables such as: slab temperature profile, wall temperature, etc., and for the temperature control of the furnace. A non-steady-state model replaces the steady-state model of the furnace. The new model can cope with the increasing diversification of slabs in terms of desired final temperature, dimensions, material properties, etc. The non-steady-state model is validated against industrial data. This model is then utilized in the design of furnace control.

The control task of reheating furnaces concerns mainly two levels. The level 1 controls furnace zone temperatures by varying the fuel flows supplied to the furnace zones. And the level 2 controls the final slab temperature by varying the zone temperature set-points.

The level 1 should accurately track the temperature set-points of furnace zones, and it should take into account other important control objectives as minimum energy consumption, and minimum operating costs. To solve the stated control problem, the distributed MPC strategy is applied based on a model of furnace zones. This later model takes into account the interaction between furnace zones. Two approaches are carried out to capture coupling effect between furnace zones. The first approach considers that temperature of one zone is a function of not only the heat input into that zone but also that of its neighboring zones. This leads to a so-called input coupling model. The second approach is more direct; it considers that the temperature of one zone is a function of the heat input into that zone and the temperature of its neighboring zones. This leads to a so-called output coupling model.

By observations in the real furnace, the temperature of furnace zones behaves as a first order transfer function with respect to input powers and zone temperatures. Based on the structure of first order transfer function and the furnace data, least square method is utilized to identify the parameters of two models. These models are then validated against measurement of zone temperatures and zone heat powers. The results of the validations showed that output coupling model is better than input coupling model in capturing the dynamic of furnace temperatures.

Therefore, the output coupling has been used for the subsequent design of distributed MPC.

A non-cooperative distributed MPC algorithm is constructed for networked system which has sub-systems interacting with each other through outputs. Local MPC controllers communicate with each other about their actions.

Each local MPC controller optimizes only the performance of its corresponding sub-system knowing how others sub-systems will behave and their influences. In order to cooperate actions of local MPC controllers, a cooperative algorithm is designed for serially connected networked system. In this later algorithm, each local MPC controller try to optimize not only the performance of its own sub-system but also two other neighbors of it. The communication between controllers is increased compared to non-cooperative strategy. However, cooperative distributed strategy brings the overall system to a better optimum point. Two algorithms are simulated for a illustrative example showing their stability and convergence property.

Using a simulation on Matlab environment, the proposed distributed MPC algorithms have been tested in an example scenario taken from real furnace system. The simulation results obtained with a validated model of furnace zones allow a comparison with the performance of the PI controllers. In general, better tracking performance and lower energy consumption is observed for distributed MPC compared to PI controllers. Moreover, the proposed method is implemented in the level 1 of a real furnace. The results obtained from this application are compared with the results of the level 1 that has been previously used. The comparison showed that distributed MPC can save up to 3 % of energy consumption while assuring good tracking of zone temperature set-points. In addition, the oxygen ratio in the furnace is reduced, which leads to less scale formation and decarburization and therefore less material losses. The global gas demand fluctuation is reduced, which leads to more stable pressure of gas supply station. Especially, automatic mode usage is increased by 10 % in soaking zone, which means that furnace operator needs less manual intervention in furnace operation.

At level 2 of the furnace, zone temperature set-points are calculated so that slabs can reach their desired final temperatures. Computation of zone temperature set-points is crucial for the success of heating task as well as for improvement of energy efficiency of the furnace. This task is subject to heating requirement of slabs, slab temperature and zone temperature constraints, and fuel supply constraints. Provided with a predefined schedule of slabs, MPC strategy utilizes the presented non-steady-state model to predict slab temperatures and furnace consumption from a given control input sequence of zone temperature set-points. The heating objectives are formulated as minimization criterions, and the limits and bounds on temperatures are accounted as penalty terms.

Therefore, the control problem of level 2 is constructed as a unconstrained minimization problem.

Particular attention should be drawn to the formulation of the objective functions. Slabs are taken into account with different weighting coefficients based on their priority. In general, slabs are accounted by their discharging order, or in other word according to their position inside the furnace, the more a slab is close to the furnace exit, the more it is important because there is less time to heat it. In addition, a slab with higher desired final temperature is more prioritized than another neighboring slab because that would avoid under-heating and slab rejection at

Fig. 4 :

 4 Fig. 4 : Coordonnées de la brame i

  concernent le profil de température de la brame au moment de son enfournement t i

  et temporelle est nécessaire pour résoudre (8) et (9). Par exemple, si nous calculons la température en 2n+1 points dans l'épaisseur de la brame, l'échantillonnage spatial est défini par

Fig. 5 :

 5 Fig. 5 : Discrétisation spatiale dans l'épaisseur de la brame

Fig. 6 :

 6 Fig. 6 : Rayonnement thermique dans un four de réchauffage Avec ces hypothèses, le problème de rayonnement thermique dans une zone est le même que celui dans un espace fermé, comme illustré Fig. 7.

Fig. 7 :

 7 Fig. 7 : Espace fermé équivalent d'une zone utilisé par la méthode stationnaire Les débits thermiques absorbés par les surfaces d'une brame dans une zone sont constitués des éléments suivants :

Fig. 9 :

 9 Fig. 9 : Structure distribuée de la commande MPC L'approche MPC distribuée nécessite tout d'abord une modélisation de la dynamique de la température des zones, ainsi que des couplages entre zones. Les algorithmes permettant la communication et la coopération entre les régulateurs MPC locaux sont ensuite élaborés.

Fig. 10 :

 10 Fig. 10 : Réponse du premier ordre de la température des zones où y 1 , y 2 , y 3 et u 1, u 2 , u 3 représentent les variations des températures et puissances des zones par rapport au point de fonctionnement du four. Ce modèle traduit un couplage par les sorties. Pour construire la loi de commande MPC distribuée, il convient de transformer (17) sous forme de représentation d'état comme suit :

Fig. 11 :

 11 Fig. 11 : Algorithme non-coopératif

Fig. 13 ;

 13 Fig. 13 ; Communication au sein d'un MPC distribué coopératif

Fig. 15 :

 15 Fig. 15 : Comparaison des comportements MPC distribué et régulateurs PID

  Fig. 16 : Consommation énergétique MPC et PID D'autres mesures du four sont comparées Fig. 17 et Fig. 18. Le ratio d'oxygène issu de la commande MPC est plus faible et surtout en général moins fluctuant que dans le cas PID. Cela peut permettre de réduire la formation de calamine sur les brames. En termes d'utilisation du mode automatique, les opérateurs du four ont moins besoin d'intervenir en manuel avec la commande MPC, particulièrement dans la zone d'égalisation.

Fig. 17 :

 17 Fig. 17: Ratio d'oxygène dans les fours

Fig. 20 :

 20 Fig. 20 : Contraintes sur la température de brames

Fig. 24 :

 24 Fig. 24 : Schéma structurel de la simulation (regroupant les trois niveaux hiérarchisés) Les Fig. 25 et Fig. 26 illustrent les épaisseurs, les largeurs et les longueurs des brames sur ces trois jours. On constate que l'épaisseur des brames ne varie pas significativement (229 mm à 235 mm), le four étant conçu spécialement pour des brames avec une certaine épaisseur. En revanche, la longueur et la largeur des brames changent de façon non négligeable d'une brame à l'autre, par exemple, de 1 m jusqu'à 2 m en largeur, et de 6 m jusqu'à 10 m en longueur.

Fig. 25 :

 25 Fig. 25 : Epaisseur de brames

Fig. 34 : 5 :

 345 Fig. 34 : Température de défournement des brames

Fig. 35 : 6 :

 356 Fig. 35 : Histogramme de l'erreur de la température de défournement des brames En termes d'homogénéité de la température de défournement des brames, nous pouvons constater sur la Fig. 36 que la structure MPC induit des gradients plus faibles. En moyenne, les gradients de température des brames sont réduits de 5°C comme le montre le Tab. 6.

Fig. 37 :

 37 Fig.37 : Trajectoire de la température de brame j=400

Fig. 42 :

 42 Fig. 42 : Trajectoire moyenne du gradient de température de brames En termes de consommation d'énergie, nous pouvons constater sur la Fig. 43 que la consommation d'énergie des zones de préchauffage et d'égalisation est diminuée alors que celle de deux zones de chauffage principales est augmentée. La consommation totale du four est montrée Fig. 44. Nous voyons que la structure MPC consomme en général moins que la commande actuelle. Le Tab. 7. récapitule les réductions obtenues via le MPC par rapport à la commande existante. La structure MPC induit en moyenne une diminution de 5% de l'énergie consommée.

Fig. 43 :Tab. 7 :

 437 Fig. 43 : Comparaison de la consommation d'énergie des zones

  Cadencement fixe des brames Par exemple, considérons deux brames consécutives labellisées 14P144 et 14P145. Le cadencement fixe impose que la brame 14P145 soit chargée dans le four 2 minutes après la brame 14P144. De même, la brame 14P145 sera évacuée du four 4 minutes après la brame 14P144.

60 )Fig. 46 :

 6046 Fig. 46 : Simulation globale avec optimisation de la cadence La Fig. 46 montre le schéma de la simulation incluant la structure MPC hiérarchisée et l'optimisation de la cadence du four. Le régulateur MPC niveau 2 optimise non seulement les consignes de température mais aussi l'intervalle de temps entre les défournements des brames. Cette dernière variable est ensuite envoyée au simulateur du four pour changer le cadencement des brames. La Fig. 47 montre que le régulateur MPC donne alors généralement des intervalles de temps entre défournements plus petits, ce qui correspond à une augmentation de la productivité du four comme le montre la Fig. 48.

Figure 1 . 1 :

 11 Figure 1.1: Hot rolling of steel slabs

Chapter 1 .

 1 Introduction 16

Figure 1 . 2 :

 12 Figure 1.2: Hierarchical cascade structure of control system

Figure 2 . 1 :

 21 Figure 2.1: Lateral view of a walking beam slab reheating furnace

Figure 2 . 2 :

 22 Figure 2.2: Top view of a walking-beam slab reheating furnace

Figure 2 . 3 :

 23 Figure 2.3: Example of path-time diagram of slabs

Figure 2 . 4 :

 24 Figure 2.4: Heating curve of a slab with respect to its position

Figure 2 . 5 :

 25 Figure 2.5: Instrumented bloom with thermal box before charging into the furnace

Figure 2 . 6 :

 26 Figure 2.6: Energy flows in a walking-beam reheating furnace

HVS04,Chapter 3 .

 3 JLKK10] [JK07]. Third, 2-dimensional or 3-dimensional temperature profile inside furnace is calculated with techniques of fluid dynamics [HCK10b, HBK09, JKZ05, Kim07]. Temperature profile of slabs is usually assumed as 1-dimensional along the thickness direction of the slabs [vDSK02, Led86a, RSCG86, SWKK11, WK99, YNY91, ZCX + 02]. In [SWKK10], some conditions that justify this assumption are given. Other models in [CCL05, ESYS97, GSR + 98b, JK07, Kim07, Li86] provide 2-dimensional temperature distribution in the slabs. These models allow an analysis of the influence of skids on the in-homogeneity of the slab temperature profile. In [HCK10b, JKZ05], a 3-dimensional temperature profile is computed. Frequently, heat conduction equation is solved by means of the finite difference method [CCL05, ESYS97, GSR + 98b, HBK09, JLKK10, Led86a, WK99, YL88, ZCX + 02]. The finite element method is used in [JK07, Modeling of the furnace 29 ZIS97

[ ESYS97 ,

 ESYS97 GSR + 98b, KKY + 00, YNY91]. Accurate formulation of energy balance also requires a model of heat losses through the furnace walls and the cooled skids [CCL05, ESYS97, ZIS97, GSR + 98b, HHLW08]. The convective heat transfer happens between the gas flow and the surfaces inside the furnace as presented in [WK99, HCK10a, KMP94, HHLW08, HVS04]. According to [HCK10a, HHLW08, KHK00, ZIS97], 5% of heat input to slabs comes from convection.

Furnace

  models for steady-state operation are usually computationally low-cost, which makes them suitable for control design [Li86, WCGS99, Kim07]. Steady-state operation implies constant slab speed, regular discharge events, similarity of product dimensions and material [SGW + 11]. In that case, the slab temperature inside furnace is a function of slab position [SGW + 11].

Figure 3 . 1 :

 31 Figure 3.1: Simplified geometry of walking beam reheating furnace

Figure 3 . 2 :

 32 Figure 3.2: Geometry and position of slab i

Figure 3 . 3 :

 33 Figure 3.3: Temperature profile of 2n + 1 points inside a slab

Figure 3 . 4 :

 34 Figure 3.4: Radiative heat exchange between slabs and environment

Figure 3 . 5 :

 35 Figure 3.5: Assumption of local radiation for a slab i residing in zone j

Figure 3 . 6 :

 36 Figure 3.6: Overheating of slab due to stops; waste of energy

Figure 3 . 7 :

 37 Figure 3.7: Grid-point discretization for one-dimensional heat conduction problem

Figure 3 . 8 :

 38 Figure 3.8: Variation of temperature with time for three different schemes

Figure 3 .

 3 Figure 3.10: Averaged slab temperature from validation experiment of the dynamical furnace model

•

  Control of furnace zone temperatures and slab temperatures • Control of temperature gradients inside slabs • Throughput rate and movements of slabs • Air and fuel ratio control which assures the proper ratio between combustion air and fuel in the air/fuel mixture • Scheduled or unscheduled downtime or delays, which may be caused by upstream or downstream processes • Control of heat recuperation units • Control of opening and closing of furnace doors • Control of temperature, pressure, and flow rates of combustion air and fuel

  presented control system directly defines reference signals for fuel flow rates, which are applied to the furnace by open-loop control. The furnace zone temperatures are measured only for safety reasons, e.g., to prevent slabs from overheating. Therefore, open-loop control of the slab temperatures is reasonable as long as the furnace zone temperatures are feedback controlled or at least monitored. Feedback control strategies control directly slab temperatures by feedback as in [ESYS97, FPSF90, KMP94, Led86a, PW98, WWC04]. The current slab temperatures can be obtained either by measurement, e.g., radiation pyrometers, from estimation of an online simulator or a state observer. It shall be noted that only slab surface temperatures are accessible by means of pyrometers. From slab surface temperatures, the slab temperature profile can be further estimated by heat conduction models. Chapter 4. Control system of the furnace 53 The control schemes presented in [WYY09] use temperature feedback from radiation pyrometers installed inside the furnace. The radiation pyrometers are also used in feedback control application by measuring slab surface temperatures after discharge from furnace or after the first roughing mill in order to compensate for remaining model uncertainties [She94, Sta04, Wö, WYY09, WWC04]. This strategy of tracking and feeding back slab surface temperatures is particularly suitable for steady-state operation. For feedback control, reference trajectories of slab temperature profiles may be provided online or preplanned and stored in the database [CXZ + 03, Bar, She94, SII + 99]. Each slab may be reheated according to its individual reference trajectory depending on parameters like the slab thickness, the desired final temperature, and the scheduled residence time. The challenging task of this strategy is the development of a MIMO controller that transforms the heat input requested by each individual slab into optimal values of control variables, e.g., zone temperatures or fuel flows. Other feedback control strategies have been developed for steady-state furnace operation or based on steady-state furnace models [BBG08, Bar, CA08, DCM + 92, ESYS97, FBF + 83, HZ82, Led86a, RB92, SII + 99, VJKM08, WCGS99, WWC04, Wö]. Obviously, if the furnace has to operate in non-steady-state conditions, the basic assumption is no longer satisfied and therefore furnace performance may significantly decline. If open-loop control structures are utilized in these cases, the problem can be particularly severe because the heating history of the furnace is totally ignored.

  Cascaded hierarchical control structure is widely utilized in the majority of published furnace control methods [BBG08, DBS88, CXZ + 03, vDSK02, FRLD06b, HZ82, IZS04, Led86a, MPS04, Pri80, She94, SS04, Sta04, VGB05, VJKM08, WYY09, WWC04, ZCX + 02, ZXWS01]. This control structure may facilitate the design of linear controllers for nonlinear systems. Particularly, cascaded hierarchical structure is useful for the control of Chapter 4. Control system of the furnace 54 interconnected subsystems who exhibit significantly different dynamics such as, in case of reheating furnaces, furnace zones, steel slabs, fuel burners [GGS01, Lev96]. For reheating furnace control, the inner loops usually control the furnace zone temperatures, for instance by PI or PID controllers [BDM99, vDSK02, DCM + 92, FRLD06b, GSR + 98a, HZ82, IZS04, Sta04, VGB05, VJKM08, WYY09, WWC04, ZCX + 02]. Furnace zones can be modeled as linear first order systems as argued in [NRAD + 12a, VJKM08]. The design approach for inner loops presented in [NRAD + 12a] accounts for the dynamic interaction between furnace zones capturing the dynamic relation between fuel inputs and zone temperatures. In some studies, set-point values of furnace zone temperatures are stored in databases as lookup tables or predefined curves [Bar, CA08, HZ82, Led86a, Pri80, She94, VGB05]. The curves or lookup tables data can be determined off-line using steady-state planning and optimization or empirically from operator's experiences. The averaging method used in these cases accounts for slab geometry, scheduled reheating time, slab material, slab desired final temperature [Bar, FPSF90, HZ82, Pri80, WWC04]. The advantage of this approach is that it can save computer power and memory. With some penalties on slab reheating quality, the furnace control task can be simplified into individual controls for each furnace zone. Disadvantage of this strategy is that the dynamic interactions between furnace zones are neglected [Bar, CA08, FBF + 83, HZ82, PW98, Sta04, SII + 99, WWC04]. Zone-based feedback control methods impose that the slab should reach predefined (optimal) temperature set-point values at the end of furnace zones [CXZ + 03, vDSK02, DCM + 92, FBF + 83, MPS04]. Therefore, the error used for feedback control is usually formulated in terms of estimated slab temperatures (at the end of furnace zones). The feedback controller either defines set-point values for the zone temperature [CXZ + 03, vDSK02, Pri80, She94, Sta04, WYY09, WWC04] or directly for the heat input power to the corresponding zone [DCM + 92, FBF + 83, Wes88] .

Chapter 4 .

 4 Control system of the furnace 56 Also based on hierarchical cascaded structure,[NvWS] presented a nonlinear model predictive controller used in the primary loop where furnace zone temperature set-points are optimized. A transient furnace model is utilized in dynamic optimization of both furnace zone temperatures and production rate. The quadratic cost function accounts for control errors in terms of final slab temperatures and for furnace zone temperatures to reduce the energy consumption.

  objective function. In the context of furnace control, one can find the utilization of linear programming (simplex algorithm) [CXZ + 03, ESYS97, WCGS99, WWC04], in [FUK + 10] with mixed integer linear programming, in [BBG08, YNY91, YUU94] with quadratic programming, or in [IZS04, NvWS] by nonlinear programming.Real-time implementations of solving a nonlinear optimization problem are described in[START_REF] Yoshitani | Optimal slab heating control for reheating furnaces[END_REF][START_REF] Yoshitani | Optimal slab heating control with temperature trajectory optimization[END_REF] with method of steepest descent and in[START_REF] Steinböck | Model-based control and optimization of a continuous slab reheating furnace[END_REF][START_REF] Steinboeck | Nonlinear model predictive control of a continuous slab reheating furnace[END_REF] with quasi-Newton method.

Figure 4 . 1 :Figure 4 . 2 :

 4142 Figure 4.1: Heating trajectory of slab inside reheating furnace

Figure 4 . 3 :

 43 Figure 4.3: Averaging technique of current furnace control

Figure 4 . 4 :

 44 Figure 4.4: MPC structure for walking-beam reheating furnace

  Afterwards, distributed MPCs are constructed with a non-cooperative strategy in which a local MPC controller only optimizes the performance of one zone, or a cooperative strategy in which a local MPC controller not only optimizes the performance of one zone but also that of the neighboring zones. Iterative algorithms are also proposed for each strategy to increase the optimality of the solutions. These algorithms are accompanied with closed-loop stability and convergence analysis. Simulation on an illustrative example is carried out showing the efficiency of the proposed distributed MPCs algorithms. Experimental results show the improvements obtained with distributed MPCs in comparison to decentralized PID controllers.

Figure 5 . 1 :

 51 Figure 5.1: Model predictive control principles

Figure 5 . 2 :

 52 Figure 5.2: Structure of decentralized MPC

Figure 5 . 3 :

 53 Figure 5.3: Structure of centralized MPC

Figure 5 . 4 :

 54 Figure 5.4: Structure of distributed MPC

Chapter 5 .

 5 Tracking of zone temperatures by distributed model predictive control 65 the temperature gradient. In fact, the soaking zone is isolated from the other zones to ensure the homogeneity of temperature. Therefore, as we want to utilize distributed MPC strategy to reduce energy consumption of the furnace, only preheating zone and heating zones are considered in this part of the dissertation.

Figure 5 . 5 :

 55 Figure 5.5: Controllable zones of walking-beam reheating furnace

Figure 5 . 6 :

 56 Figure 5.6: Imposed space between slabs

. 1 ) 5 .Figure 5 . 7 :Figure 5 . 8 :

 155758 Figure 5.7: Step signal of zone preheating

Figure 5 . 9 :

 59 Figure 5.9: Input coupling of furnace zones

Figure 5 .

 5 Figure 5.11: Validation of output coupling model

Figure 5 .

 5 Figure 5.13: Structure of observers

  each sample time k, each local MPC controller executes five main actions: initializing, communicating with other controllers, solving the local optimization problem, checking the convergence, and applying the optimized control input to the corresponding subsystem [NRAD + 12b, ZL07]. These steps are explained as follows • Step 1: Initialization

Figure 5 .

 5 Figure 5.14: Diagram of the non-cooperative iterative algorithm

Figure 5 .

 5 Figure 5.15: Communication between MPC controllers

  and x j (k) with j = i ± 1 to calculate necessary terms of quadratic programming problem in Equation5.45.• Step 3: Optimization.Each MPC controller solves its local quadratic problem expressed in Equation5.45 and obtains the new control input sequence for the iteration.

U

Figure 5 .

 5 Figure 5.16: Diagram of cooperative distributed MPC algorithm

  Figure 5.17: Results of non-cooperative distributed MPC

  , 5.18 and 5.19, distributed MPC gives good performance in

Figure 5 .

 5 Figure 5.18: Results of cooperative distributed MPC

Figure 5 .

 5 Figure 5.19: Output and control signals under centralized MPC

Figure 5 .

 5 Figure 5.20: Number of iteration

Figure 5 .

 5 Figure 5.22: Input power of furnace zones

Figure 5 . 5 . 5 Chapter 5 .Figure 5 .Figure 5 .

 555555 Figure 5.23: Zone temperatures

Figure 5 .

 5 Figure 5.28: Utilization of automatic mode

Chapter 5 .Figure 5 .Figure 5 .

 555 Figure 5.30: Specific energy consumption

  level 1 is considered transparent to level 2 meaning the tracking of zone temperature set-points has high accuracy and fast response in comparison to calculation frequency of level 2. MPC strategy in level 2 computes the optimal zone temperature set-points which minimize error of final slab temperature and energy consumption. This optimization of temperature set-points is based on a simplified version of the non-steady-state model. In the following, different control objectives and constraints are specified.

Figure 6 . 2 :

 62 Figure 6.2: Cascaded closed-loop of the furnace control system with MPC strategy

Figure 6 . 3 :

 63 Figure 6.3: Variation of desired final slab temperatures in an operation of the considered furnace

Figure 6 . 4 :

 64 Figure 6.4: Variation of slab width and length

Figure 6 . 6 :

 66 Figure 6.6: Constraints on final slab temperature profile

Figure 6 . 7 :

 67 Figure 6.7: Constraints on inhomogeneity of slab temperature

Figure 6 . 8 :

 68 Figure 6.8: Distinction between S in (τ ) and S 1 (τ ) for furnace operation during prediction horizon τ, τ + ∆t P

Chapter 6 .Figure 6

 66 Figure 6.9: MPC applied for slab temperature control: a) at time τ , b) at time τ + ∆τ s

Figure 6 .

 6 Figure 6.10: Heating of consecutive slabs that have big difference on desired final temperature

Figure 6 .

 6 Figure 6.11: Variation of weighting coefficient on final slab temperature error according to slab position

Figure 6 .

 6 Figure 6.13: Variation of weighting coefficient on temperature gradient according to slab position

Figure 6 .

 6 Figure 6.14: Weighting coefficients according to slab position and time

Chapter 6 .

 6 These methods are also referred to as direct search methods[START_REF] Lewis | Direct search methods: then and now[END_REF][START_REF] Torczon | Multi-directional search: a direct search algorithm for parallel machines[END_REF][START_REF] Hooke | direct search''solution of numerical and statistical problems[END_REF][START_REF] Kolda | Optimization by direct search: New perspectives on some classical and modern methods[END_REF]. Slab temperature control 121 For reheating furnaces, utilization of gradient-based methods will need finite-differencing of the objective function. This choice will require repeating furnace simulation in order to evaluate the objective function J o (U (k))

Figure 6 . 6 . 6 .

 666 Figure 6.15: Architecture of the furnace control system

Figure 6 .

 6 Figure 6.18: Example of a simplex in ℜ 2 and ℜ 3

Chapter 6 .Figure 6 . 3 Figure 6 .

 6636 Figure 6.19: Reflection of worst vertex on the centroid of the opposite face in ℜ 2 and ℜ 3

Figure 6 .

 6 Figure 6.21: Reflection, expansion, and contraction moves of Nedlder-Mead simplex method

  upper bounds of the function variable. Authors have also pointed out that by projecting on bounds, the simplex may collapse into subspace of the saturated variables. To tackle the problem, a convergence test based on re-initialization at the point of convergence is then required.

Figure 6 .Figure 6 .

 66 Figure 6.23: Rosenbock's method in problem of dimension 2

Figure 6 .

 6 Figure 6.25: Simulation-based optimization of the MPC strategy

Figure 6 .

 6 Figure 6.26: Software-based simulation environment

Figure 6 .

 6 Figure 6.27: Thickness of slabs

Figure 6 .

 6 Figure 6.28: Width and length of slabs

Figure 6 .

 6 Figure 6.27 shows the thicknesses D j of the slabs. The thicknesses vary between 229 mm and 235 mm and most common value is 235 mm, this variation is not considerable because, as mentioned previously, the considered furnace is dedicated to the heating of slab with thickness of 235 mm. The width and length of the slabs are spread over the range of 0.99 m to 1.93 m and 5.90 m to 10.75 m, respectively as shown in Figure 6.28. Despite highly regular thicknesses, the slab width and length vary greatly. Thus, the transient effects on the heating process are unavoidable.

Figure 6 .

 6 Figure 6.29: Residence time of slabs 1 to 399

  The usual position-time diagram of representative slabs are shown in Figure6.31. The lines represents the position of each slab in function of time. The position of a slab is the distance from its center to the entry of the furnace.Due to the long stop at time 46.3 h, the slabs j = 600 to j = 606 remained at the same position for about 2 h.For the considered reheating scenario, the temperature set-point of zone preheating, heating 1, heating 2 and soaking are limited respectively between Tz,min =[830, 1085, 1200, 1210] • C and Tz,max =[1100 , 1250, 1355, 1300]

Figure 6 .Figure 6 .

 66 Figure 6.30: Residence time of slabs 400 to 783

Figure 6 .

 6 Figure 6.32: Position-time diagram of slabs 506 to 606

  • C, and the temperature gradient of slabs should remain below ∆T j,max = 350 • C. Maximum inhomogeneity of final slab temperature profile ∆ j,max,end is dependent of slab thickness. Since the considered furnace heats most slabs that are 235 mm thick, ∆ j,max,end is fixed at value 100 • C. The Nelder-Mead simplex algorithm is used as optimization method. For initialization of simplex algorithm, we use a simplex with size a = 10 • C. The initial control inputs U 0 of the simulation are imposed with the averaged values of upper and lower limits, i.e., U 0 = ( Tz,max + Tz,min )/2=[965, 1167.5, 1277.5, 1255] • C. Prediction horizon ∆t P = 120 min which represents average residence time of a slab inside burner zones. The sampling period is ∆t s = 5 min which represents the average discharging time (or loading time) interval between 2 successive slabs. The control horizon is ∆t C = 5 min, which means that N C = 1. About the parameter tuning of the objective function, weighting factors on error of final slab temperatures and final slab temperature gradients are α 1 = 70, α 2 = 10. Others weighting factors concerning energy consumption are α 2 = α 3 = α 5 = 10. Weighting factor on penalty terms µ is imposed with great value, here, µ = 1000.

Figure 6 .Chapter 6 .

 66 Figure 6.33: Desired final temperature of slabs 1 to 399

Figure 6 .

 6 Figure 6.35: Averaged final slab temperatures

Chapter 6 .Figure 6 .

 66 Figure 6.37: Consumed power of the furnace

  , 6.39, 6.40 and 6.41 in comparison with values designed by the existing controller for the time interval [0 h, 64 h]. The dashed lines show the upper and lower bounds of the temperature set-points. The red solid curve shows the set-points calculated by MPC controller, whereas the blue solid curve shows the set-points given by the existing controller.

  [2 h, 5 h] and [50 h, 60 h]. The PID controllers have no problem tracking set-points provided by MPC controllers.

Figure 6 .

 6 Figure 6.38: Temperature set-point of soaking zone calculated by MPC controller compared to that of the existing controller

Figure 6 .

 6 Figure 6.41: Temperature set-point of zone preheating calculated by MPC controller compared to that of the existing controller

  , 6.43, 6.44 and 6.45. The final slab temperatures simulated with MPC controller are represented by the red dotted curve and the blue one represents the temperatures simulated with the exiting controller. Most of the temperatures provided by MPC controllers stay inside of the expected bounds with exceptions of slabs of which desired final temperatures are particularly different with those of neighboring slabs such as slabs j = 269 to j = 275 or slabs j = 384 to j = 386 in Figure 6.43.

Figure 6 .

 6 Figure 6.43: Final temperature of slabs 203 to 404

Chapter 6 .Figure 6 .

 66 Figure 6.44: Final temperature of slabs 405 to 606

Figure 6 .Figure 6 .Figure 6 .

 666 Figure 6.45: Final temperature of slabs 607 to 783

Figure 6 .Figure 6 .

 66 Figure 6.52: Temperature gradient of slab j = 200

Figure 6 .

 6 Figure 6.54: Temperature gradient of slab j = 400

Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .

 6666666 Figure 6.55: Temperature trajectory of slab j = 600

J c, 1 (Chapter 7 .

 17 u i (k + j) -T z,i,max (k + j)) + | min(0, u i (k + j) -T z,i,min (k + j))|] (7.10a) ui (k + j) -Ṫz,i,max (k + j)) + | min(0, ui (k + j) -Ṫz,i,min (k + j))|]. (7.10b) Boosting of productivity 169The constraints on final slab temperature gradient is given asJ c,3 (k) = N P (k) ∑ i=1 max(0, g(k + i) -g max (k + i)) (7.11a) = ∑ j∈S 1 (k) max (0, g(k j,1 ) -∆T j,end,max ) . (7.11b)The upper and lower bound on final slab temperatures are formulated asJ c,4 (k) = N P (k) ∑ i=1 [max(0, T s (k + i) -T s,max (k + i)) + | min(0, T s (t) -T s,min (k + i))|] (7.12a) = ∑ j∈S 1 (k) max(0, T s (k j,1 ) -T s,max (k j,1 )) + | min(0, T s (k j,1 ) -T s,min (k j,1 ))|].(7.12b)The slab temperature trajectory and temperature gradient are limited by upper bounds T max , ∆T max , respectively.J c,5 (k) = N P (k) ∑ i=1 [max(0, [x 1 (k + i) -T max (k + i)] ′ [x 1 (k + i) -T max (k + i)])] (7.13) J c,6 (k) = N P (k) ∑ i=1 [max(0, [x 2 (k + i) -∆T max (k + i)] ′ [x 2 (k + i) -∆T max (k + i)])] (7.14) Important constraint of scheduling task is the minimum time interval between discharging events δt sch,min , i.e., J c,7 (k) = N C ∑ i=1 | min(0, δt sch,min (k + i) -δt sch,min )|. (7.15)

Figure 7 . 3 :Figure 7 . 4 :

 7374 Figure 7.3: Implementation diagram of MPC strategy

Figure 7 . 5 :Figure 7 . 6 :

 7576 Figure 7.5: Discharging time interval of events j = 451 to j = 900

Figure 7 . 7 :Figure 7 . 8 :

 7778 Figure 7.7: Average final temperature of slabs j = 451 to j = 900

Figure 7 .

 7 Figure 7.9: Final slab temperature gradient

Chapter 7 .Figure 7 .Figure 7 .

 777 Figure 7.11: Furnace throughput rate measured at discharging events j = 1 to j = 450

Figure 7 .

 7 Figure 7.13: Average temperature of slab number 570

Figure 7 .

 7 Figure 7.14: Average temperature gradient of slab number 570

Figure 8 . 1 :

 81 Figure 8.1: Continuous annealing process

Figure 8 . 2 :

 82 Figure 8.2: Temperature trajectory of the strip through different zones

Chapter 8 .

 8 Extension to other metallurgical furnaces 181 The over-aging chamber provides electrical heating to maintain the strip at an intermediate temperature at around 400 • C. The temperature trajectory of the strip is shown in Figure 8.2 [CKC00].

  3.

Figure 8 . 3 :

 83 Figure 8.3: Continuous galvanizing process

Figure 8 . 4 :

 84 Figure 8.4: Temperature control of furnace with PID controllers

Figure 8 . 5 :

 85 Figure 8.5: Hot blast stoves with blast furnace

Figure 8 . 6 :

 86 Figure 8.6: Internal structure of hot blast stoves

Chapter 8 .Figure 8 . 7 :Figure 8 . 8 :

 88788 Figure 8.7: Temperature of the fluids (blast air, waste gas) and the solid (refractory layers) at different levels of a hot blast stove

  Le problème d'optimisation MPC sur l'horizon de prédiction ∆𝑡 𝑃 , l'horizon de contrôle ∆𝑡 𝐶 est défini par la minimisation du critère suivant :

	1 x t ( )	1 12 1 ( , ) , max max j j y y j j T y t      2 2 ( , ) ,..., ( , ) max j j j N j N in in y j N in T y t T y t        	T	(32)
	2 x t () Δ ( 1 j T y t T y t 1 2 2 , ),Δ ( , ),...,Δ j j j T      ( , ) , y T j N j N in in y t    	j	2,  2 , j j D D   
	Où [j+1, j+2,…, j+N

in ] est l'ensemble des indices des brames dans le four à l'instant t.

  S 1 représente l'ensemble des brames défournées pendant la période [𝜏, 𝜏 + ∆𝑡 𝑃 ]. Les termes J 1 , J 2 de l'équation (33) représentent les objectifs sur la température de défournement et son homogénéité. Ces termes sont exprimés par les relations :

	(34)
	p   t
	1 x t T t t τ τ max ( ) ( ), , Δ t  p     
	2 x t () Δ T t t τ τ max ( ), , Δ t  p     

où

  3.

	Nom de paramètre	Continu	Discret
	Instant de temps d'enfournement	t	j	,0	k	j	,0
	Instant de temps de défournement	t	j	,1	k	j	,1
	Instant de temps	τ			k
	Période d'échantillonnage	Δ τ	s	Δ t	s
	Horizon de prédiction	Δ t	p	N	p
	Horizon de commande	Δ t	c	N	c
	Tab. 3 : Variables discrètes et continues correspondantes
	On en déduit alors la forme discrète des termes de la fonction de coût comme suit :

  sont nettement inférieures à celles données par la commande actuelle. Les valeurs moyennes des consignes sont récapitulées dans le Tab. 4. Consignes de température de zone préchauffe Les températures de défournement des brames sont comparées Fig.34. En général, la structure MPC donne des températures de brame plus précises. Ceci est confirmé par le Tab. 5, pour lequel la loi MPC, obtient pour 92% des brames la température dans la marge des limites, tandis que le régulateur existant n'en obtient que 75%.

							Tab. 4 : Moyennes des consignes de température de zone
							Fig. 32 : Consignes de température de la zone de chauffe 1
							Fig. 33 :
						Set-points MPC			Upper bound			Lower bound			set-point database
		1320
	(°C)	1260 1290
		1230
		1200
			0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21
													Temps (h)
						Set-points MPC				Upper bound			Lower bound			set-point database
			1400
			1350
	(°C)	1250 1300
			1200
			1150
			0	1	2	3	4	5	6	7	8	9	10 11 12 13 14 15 16 17 18 19 20 21
													Temps (h)

  3.1 System of equations for radiative heat transfer inside reheating furnace . . . . . . . . . . . . . . 5.1 Classification of distributed MPC structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Comparison between PID controllers and distributed MPC . . . . . . . . . . . . . . . . . . . . 6.1 Constraints on zone temperature set-points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Constraints of slab temperature profile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Correspondence of continuous and discrete parameters . . . . . . . . . . . . . . . . . . . . . . 6.4 Difference between steady-state model used in furnace observer and set-point optimizer . . . . . 6.5 Average consumed power of the furnace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Average zone temperature set-points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.7 Deviation of final slab temperature from the desired value . . . . . . . . . . . . . . . . . . . . 6.8 Inhomogeneity of final slab temperatures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.9 Normalized specific energy consumption . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Temperature constraints of hot blast stove . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	7.1 Final slab temperature and temperature gradient with MPC strategy and the existing controller .
	7.2 Production rate and energy consumption indicators . . . . . . . . . . . . . . . . . . . . . . . .
	8.1

  The rolling mill needs a diversity of products in terms of steel grade, dimensions, desired rolling temperature. The steady-state model cannot provide accurate temperature estimations for furnace control. Consequently, safety margins of slab temperatures are used to avoid under-heated slabs. The strategy is rather conservative, and it leads to general over-heating of slabs, especially during stops of the furnace as shown in Figure3.6. Energy cost reduction of reheating furnace

	Discharged temperature (°C)	1180 1200 1220 1240 1260	Overheated slabs due to stops		Plant estimated temperatures Desired temperatures
		1160
		0	20	40	60	80	100	120	140	160	180	200
						Slab number

4 g Table 3 .

 g3

1: System of equations for radiative heat transfer inside reheating furnace

Table 5 .

 5

1: Classification of distributed MPC structure

  3, given in Equation 5.1, and H ij , j ̸ = i described in Equation 5.3. Because interactions are mainly between neighboring zones, H 13 and H 31 are neglected. Therefore, in Equation 5.4, coupling terms H 13 and H 31 are neglected. The input coupling model of zones is simplified to

  Tracking of zone temperatures by distributed model predictive control

							91
	Zone	Average temperature error( • C) Power consumption(Watt)
		PID DMPC Improvement	PID	DMPC	Gain
	Preaheating	8	6	2	9572	9163	4.87%
	Heating 1	6	4	2	26026 25768	1%
	Heating 2	4	3	1	25678 25224 1.78%

5.22

, it is noticed that distributed MPC does not have delay on control action as seen with PID controllers, this can be explained due to anticipation feature of MPC strategy. Moreover, tracking performance with distributed MPC can be visually observed in Figure

5

.23 to Chapter 5.

Table 5 . 2 :

 52 Comparison between PID controllers and distributed MPC be improved from that of PID controllers. Especially in transient phase, distributed MPC reduces the overshoot effects which is produced with PID controllers. More precisely, as shown in Table5.2, although only slight improvement on tracking temperature is obtained, energy consumption can be saved from 1% to 5% according to zones.

Table 6 .

 6

	1: Constraints on zone temperature set-points
	Other modifications of Tz,i,min (t), Tz,i,max (t), Ṫz,i,max (t), and Ṫz,i,min (t) may come from furnace operators.

Table 6 . 2

 62

Description Nominal values T j,max . . . . . . . . . . . upper bound on slab temperature . . . . . . . . . . . . . . . . . . . . . . 1350 • C Tj,end,max -Tj,end upper margin of final slab temperature . . . . . . . . . . . . . . . . . 20 • C Tj,end -Tj,end,min lower margin of final slab temperature . . . . . . . . . . . . . . . . . 20 • C ∆T j,max . . . . . . . . . maximum inhomogeneity of slab temperature profile . . . . 350 • C ∆T j,end,max . . . . . . maximum inhomogeneity of final slab temperature profile 100 • C : Constraints of slab temperature profile

Table 6 .

 6 .3 Parameter Names Continuous notation Discrete notation Loading time of slab j . . . t j,0 . . . . . . . . . . . . . . . k j,0 Discharging time of slab j t j,1 . . . . . . . . . . . . . . . k j,1 Time instant . . . . . . . . . . . . τ . . . . . . . . . . . . . . . . . k Sampling time . . . . . . . . . . ∆τ s . . . . . . . . . . . . . . ∆t s Prediction horizon . . . . . . . ∆t P . . . . . . . . . . . . . . N P Control horizon . . . . . . . . . ∆t C . . . . . . . . . . . . . . N C 3: Correspondence of continuous and discrete parameters

  , e.g., in function fminsearch in Matlab optimization toolbox. In general, what is noted about the method is that it tends to be able to satisfy optimization purpose by producing a rapid initial declines of function value in comparison to other direct search The algorithms of this class accelerate the search by adapting search directions to the curvature of the objective function obtained during the course of the search. Underlying algorithms are Rosenbrock's and Powell's methods

	methods [Wri10].
	6.3.2.3 Adaptive methods
	[Pow64, Ros60].

  • C, 155 • C, 175 • Cand 270 • C for zone soaking, heating 2, heating 1, and preheating, respectively. As a result, we take an average value a = 10 • C. Different operations are defined with the coefficients of reflection, expansion, contraction and shrink. In general, reflection coefficient has value 1, expansion coefficient takes value 2, contraction and shrink operations have coefficient of 1/2. Specific steps of the algorithm is given in the following• Step 1-Initialization of simplexThe initial simplex in n -dimensional space is defined by (n + 1) points {U 0 , U 1 , ..., U n } which are calculated in Equations 6.62. The objective function is evaluated at each point of the simplex:y i = J o (U i ) for i = 0, n. y h = J o (U h )and y l = J o (U l ) represent respectively highest and lowest value of cost is the reflection coefficient. If the function value at the new point is lower than the current highest point, i.e. U l < y * = J o (U * ) < U h then U h is replaced by U * and the algorithm restarts with new simplex. If y * ≤ y l , go to step 3, if y * ≥ y h , go to step 4. If y * * < y l , U h is replaced by U * * and the algorithm restarts. If y * * > y l , U h is replaced by U * before restarting.Where β is contraction coefficient. If y * * ≤ y h then P h is replaced by U * * and the algorithm restarts. On the contrary, if the highest point U h cannot be replaced with any of reflection, expansion, or contraction operations. Shrink operation is executed in order to generate a new simplex, i.e.

	function. U is centroid point of the simplex, i.e., U = 1/n	∑ n j=0 U j .
	• Step 2-Reflection operation
	The highest point U h is reflected through the centroid point U to a point U * which is determined as
	U * = (1 + α)U -αU h	(6.63)
	where α = 1 • Step 3-Expansion operation
	U * is expanded to U * * by a relation
	U * * = γU * + (1 -γ)U	(6.64)
	where the expansion coefficient γ = 2. • Step 4-Contraction operation
	U h is replaced by U * . Contracted point U * * is determined by
	U * * = βU h + (1 -β)U .	(6.65)

Table 6 . 5

 65

			Average consumed power (MW)
		Soaking Heating 2 Heating 1 Preheating Total
	Simulator	10	28.91	42.27	16.35	97.53
	Measurement	10.97	28.06	43.12	16.89	99.04

: Average consumed power of the furnace

Table 6

 6

	Chapter 6. Slab temperature control

.6. In average, the MPC controller has redistributed the zone temperatures towards the exit of the furnace, especially at the zone heating 2. By increasing the average temperature set-points of zone heating 2, MPC controller is able to decrease that of other zones, especially in the front zones: preheating and heating 1. Now, we should be interested in the heating results of MPC controller.

Table 6 . 6

 66

: Average zone temperature set-points Figure 6.42: Final temperature of slabs 1 to 202

Table 6 .

 6 7: Deviation of final slab temperature from the desired value

Table 6 .

 6

9: Normalized specific energy consumption

Table 7 .

 7 1: Final slab temperature and temperature gradient with MPC strategy and the existing controller

Table 7 . 2 :

 72 1.27 Gj/t. In Figure7.13, we can see the average temperature of slab number 570 with MPC and with existing controlled in the same time reference. The slab is heated faster with MPC controller and therefore discharged sooner than with Quantity MPC Existing controller Average throughput rate (t/h) 296 281 Average reheating time (min) 169 179 Total heating time (h) 68 72 Specific energy consumption (Gj/t) 1.27 1.25 Production rate and energy consumption indicators

Table 8 .

 8 Table 8.1. Finally, a maximum allowed value of final by-pass valve position at the end of the blast phase is imposed.

	Location Maximum constraints Minimum constraints
	Dome	1500 • C	900 • C
	Interface 1200 • C	800 • C
	Bottom	350 • C	200 • C

1: Temperature constraints of hot blast stove

2.1. La conduction thermiqueSachant que la largeur et la longueur des brames sont largement plus importantes que l'épaisseur, seul le problème de conduction dans l'épaisseur de la brame est considéré. La Fig.4illustre les coordonnées associées à la brame i, de dimensions (W i x L i x D i ). La température en un point de la brame dépend uniquement de sa position y et du temps t, notée donc T i (y, t) (avec -D i /2 ≤ y ≤ D i /2) pour désigner la température d'un point de la brame i.

Au niveau 2 du four, les consignes de température de zones sont calculées de sorte que les brames puissent atteindre leurs températures visées. La satisfaction de cette spécification est cruciale pour le succès de la tâche de chauffage ainsi que pour l'amélioration de l'efficacité énergétique du four. Cette tâche est par ailleurs soumise à des contraintes de température de zones, de température de brames, et de puissance de zones. La stratégie MPC de ce niveau 2 utilise le modèle de transfert thermique transitoire de façon à prédire la température des brames ainsi que la consommation du four. Les objectifs de chauffage sont formulés via des critères qu'il convient de minimiser, et toutes les contraintes mentionnées sont prises en compte comme des termes de pénalité incorporés dans le critère d'optimisation. Par conséquent, le problème de commande de niveau 2 se ramène à un problème d'optimisation sans contrainte.Une attention particulière doit être portée à la formulation des fonctions de coût. Les brames doivent être prises en compte avec des coefficients de pondération en fonction de leur priorité. En effet, en général, les brames sont comptabilisées par leur ordre de défournement, ou en d'autres termes en fonction de leur position à l'intérieur du four, plus une brame est proche de la sortie du four, plus elle est importante, car moins il reste de temps pour la chauffer à sa température visée. De plus, une brame ayant une température visée plus importante est plus prioritaire qu'une autre brame voisine, le but étant de la chauffer suffisamment pour éviter le sous-chauffage et le rejet de la brame au laminoir à chaud. La zone d'égalisation est particulièrement dédiée à l'homogénéisation de la température de brames, par conséquent, dans la formulation du critère de minimisation relatif au gradient de température, des brames situées dans la zone d'égalisation sont plus prioritaires que celles des autres zones.

Part II

Control design

Chapter 6

Slab temperature control

The previous chapter has dealt with zone temperature control, so-called level 1 of Figure 6.1. In this chapter, the tracking problem of zone temperature set-points has been resolved with distributed MPC strategy. Zone temperature set-points are in fact calculated by the level 2 of the control system (see Figure 6.1). This latter level is also mentioned as slab temperature control because it determines heating process of slabs through out the furnace. Therefore, optimization of temperature set-points calculation is crucial to obtain efficient furnace operation. In this chapter dealing with the elaboration of these zone temperature set-points, the problem statement is firstly given describing qualitative and quantitative control objectives of level 2. Secondly, model predictive control strategy is clarified for slab temperature control. Next, mathematical optimization methods are presented, among which Nelder-Mead simplex algorithm is further utilized as the optimization method of the MPC strategy.

A numerical simulation is carried out for an example scenario of the real furnace system. The obtained results showed promising improvements in term of energy consumption and temperature accuracy. steps defined by the columns of the matrix ∆ k P k will produce a simple decrease (to be differed from sufficient decrease as for gradient-base methods like steepest descent method), and subsequently evaluate f (x) at only one point. This will reduce a lots of required computational resource. However, it is not always easy to do such guessing. At the other extreme, an explanatory moves algorithm would evaluate all possible moves defined by

k is a column of P k ), and returns the trial step that produces the smallest function values.

A pattern search method is distinguished partly with other methods by its explanatory moves algorithm. In a generalized description of pattern search methods, we want to allow the maximum possible choice of explanatory moves algorithm and still maintain the required properties for the convergence analysis. These requirements are given in the following hypothesis of explanatory moves.

Hypostheses of explanatory moves

• (1)

The algorithm of generalized pattern search method is given as follows:

Algorithm of generalized pattern search method

Let initial variables x 0 and ∆ 0 be given a priori. Then for k = 0, 1, ...,

• (2) Determine a step s k utilizing an explanatory moves algorithm.

A particular pattern search method specified by: the basic matrix, the generating matrix, the explanatory moves algorithm utilized to determine the step s k , and the algorithms for updating C k and ∆ k .

Two components: the generating matrix C k and the step length ∆ k need to be updated every iterate k. Generally, we want to update ∆ k so that ϵ k > 0 which means the iteration k is successful, otherwise it fails to find new iterate that decreases the function value. It should be recalled that to accept a step, we only need a simple, as different from sufficient, decrease of function value. The idea for updating ∆ k is to reduce it when the iteration Chapter 6. Slab temperature control 126 is unsuccessful and to increase it when the iteration is successful. Therefore, updating algorithm can be given by

Updating of step length ∆

θ and λ k can be determined with given parameters τ ∈ Q and {ω 0 , ω 1 , ..., ω N } ⊂ Z, with N < +∞, τ > 1, ω 0 < 0, and ω i ≥ 0, i = 1, ..., N , as follows θ = τ ω 0 (6.54)

The algorithm for updating C k varies with each particular pattern search method. It is sufficient to choose C k satisfying the form in 6.51 with corresponding conditions placed on matrices

For analysis of general convergence results, we define the set L(x 0 ) as

Assume that L(x 0 ) is compact. And that function f (x) is continuously differentiable on a neighborhood of L(x 0 ). Then for the sequence x k produced by the algorithm of genralizedd pattern search method given above, the following result of convergence is assured:

To obtain a stronger convergence results, the hypotheses of explanatory moves must be strengthened as

Strong hypotheses of explanatory moves

Assume that L(x 0 ) is compact, and function f (x) is continuously differentiable on a neighborhood of L(x 0 ). In addition, if the columns of the generating matrix

with Figure 6.48, and Table 6.8. Average gradient is reduced from 73 • C to 68 • C and 3 more percent of slabs is the theoretical upper limit of final temperature gradient which is fixed at 100 • C. simulated with MPC controller, whereas the blue one is that given by existing controller. The corresponding temperature gradients are given in Figures 6.50, 6.52, and 6.54. The thin lines with markers are the maximum and minimum temperature of slabs. Generally, we see that the temperature trajectories respect the upper limit and the final temperature gradients satisfy the required upper bounds. The temperature gradients hardly exceed the upper bound.

The average slab temperatures given by MPC controller are higher than that simulated with the existing controller.

Similarly, the maximum temperature of slab given by existing controller is also higher than that simulated with the MPC controller. Temperature gradients of the considered slabs given by existing controller increase rapidly since slabs enter the furnace. This is one of the causes leading to energy inefficiency. On the contrary, with MPC controller, the gradients evolve more slowly compared to those simulated with the existing controller. This tendency will lead to reduction of energy consumption.

For the slabs that are affected by the furnace stops as slabs j = 600, j = 601 and j = 641, their representative Chapter 7

Boosting of productivity

Slab temperature control may involve control of furnace temperature or control of slab velocity or both. In part II, we have focused on the problem of optimizing furnace temperature set-point in level 2 and control them in level 1 with a predefined slab scheduling. That is, the sequencing of slabs is known in advance with a table providing the time interval between loading and discharging events of consecutive slabs as shown in Figure 7.1. For example, we consider two consecutive slabs labeled 14P144 and 14P145. Then, the predefined schedule imposes that the slab 14P144 will be loaded into the furnace 2 minutes after the slab 14P145. Similarly, the slab 14P144 will be discharged from the furnace 4 minutes after the slab 14P145.

In practice, when downstream process, e.g., hot rolling mill does not constitute a bottleneck of the production chain, the priority of reheating furnace is maximizing the throughput rate, and obviously, at the same time assuring the requirement of heating quality. Therefore, we will retake the solution previously constructed in chapter 6 and adapt it to the problem of maximizing throughput rate. This chapter will start with overview on published work on scheduling problem of reheating furnace. Afterwards, the problem statement clarifies the objectives and constraints of furnace control concerning slab scheduling. Then, the formulation of the problem is given. And finally, simulation results based on industrial data are presented and discussed.

Literature review

In published works for reheating furnace control, there are generally two classes of control: production scheduling and furnace temperature control. The production scheduling may ignore the effect of furnace temperature control and the non-uniformity of distributed heat power on the productivity of the furnace. Whereas, works on furnace temperature control often consider that production scheduling is predefined.

In [START_REF] Tang | A review of planning and scheduling systems and methods for integrated steel production[END_REF], a review on scheduling method for integrated steel plant is provided. The authors emphasized on the direct path of slab from continuous caster to hot rolling mill, where reheating furnaces are not present or only used to slightly heat a small amount of slabs. This is hot-charging problem because at the furnace entrance, slabs are already at high temperature. Another review is given in [START_REF] Dutta | A survey of mathematical programming applications in integrated steel plants[END_REF] where production planning and scheduling of integrated steel plants are solved by mathematical methods.

In [START_REF] Lopez | The hot strip mill production scheduling problem: A tabu search approach[END_REF], the production scheduling task of two identical reheating furnaces in a hot strip mill is determined by the tabular search heuristic. In this method, each slab is attributed with a desired forwarding velocity, and each furnace is run at the desired velocity of the slowest slab it contains.

In [START_REF] Petersen | Inter-process synchronization in steel production[END_REF], a heuristic scheduling system for a single slab reheating furnace and a rolling mill is presented. The system uses a modified greedy algorithm with small decision trees and total enumeration. An objective function that accounts for the empty space in the furnace, the discharging times of slabs and the deviations from ideal values of slab residence times is used.

In [START_REF] Assaf | Steel scheduling optimization for ipsco's rolling mill and reheat furnace[END_REF], an algorithm for production scheduling of a reheating furnace and a hot rolling mill is designed based on a heuristic model. A cost function taking into account reheating times and fuel consumption and scheduling constraints are considered in the sequencing problem. The problem is solved by enumeration and a modified and pruning procedure. The constraints keep the problem size in a reasonable range. The utilized heuristic model of furnace requires that the residence time of slabs in each furnace zone is lower bounded.

In [START_REF] Broughton | A paradigm for the scheduling of a continuous walking beam reheat furnace using a modified genetic algorithm[END_REF], the scheduling task of a reheating furnace is developed. A furnace model is used to calculate the temperature gradient and mean error of final slab temperature. An objective function tacking into account the later quantities is minimized with a genetic algorithm.

In [START_REF] Jaklič | The influence of the space between the billets on the productivity of a continuous walking-beam furnace[END_REF], a reheating furnace for billets is studied. The influence of predefined space between billets on furnace productivity is analyzed. The study showed that an optimum space between billets exists, i.e., maximal furnace productivity can be achieved while assuring the heating requirements. For the considered furnace, the products to be heated are steel slabs. Spaces between slabs are predefined and do not significantly influence the furnace productivity because these spaces are much smaller than slab's width and length. Therefore, this space between slabs will not be included in the productivity study of the considered furnace.

In the past, the optimum furnace zone temperatures and forwarding velocity for certain kinds of products are often tabulated for steady-state operation as presented in [START_REF] Leden | A control system for fuel optimization of reheating furnaces[END_REF]. In [START_REF] Knoop | Nonlinear pi-controller design for a continuous-flow furnace via continuous gain scheduling[END_REF], a nonlinear PI-type furnace controller is proposed that uses only the slab velocity as a control input. In the control of a small reheating furnace with We test the forgoing design of MPC strategy for the scheduling optimization of the three-day furnace operation presented in chapter 6. The simulation scheme is given in Figure 7.2. In this scheme, in addition to optimized set-points of zone temperature, the optimization also delivers the optimal time interval of the discharging events.

Simulation results

The latter value will be used to update the slab scheduling. In practice, it is a little difference since scheduling optimization by MPC strategy, instead of directly modifying the slab sequencing, will propose to the furnace operator the new time interval that increases productivity and at the same time satisfies the heating requirements, see hot rolling mill. In the formulation of minimization criterion concerning temperature gradient of slabs, soaking zone is especially dedicated to homogenization of slab temperature, therefore, slabs located in soaking zone are more prioritized than those of other zones.

In term of optimization algorithm, direct search algorithms are studied because the optimization problem is based on simulation of the furnace. Numerical non-steady-state model provides more accurate estimation temperature variables but also demands more computational resource than the steady-state model. Using gradient-based optimization methods needs analytic form of cost function gradient or their estimation by finite differencing.

Nelder-Mead algorithm is utilized due to its advantages like: simple implementation, fast decline of function value.

Using a simulation on Intel Fortran environment, the proposed MPC strategy and optimization algorithm for level 2 is tested in a scenario of real furnace system. In this scenario, 900 slabs are processed within 3 days with an average heating time of 180 min. As mentioned before, the sequencing of slabs is predefined with an average throughput rate of 281 ton/h. The simulation obtained with a validated simulator of the furnace allow comparison of MPC strategy and the existing strategy of level 2. The obtained results show that utilization of MPC controller can save 5 % of furnace energy consumption while achieving higher slab final temperature and slab temperature homogeneity.

A study aiming at maximization of furnace productivity is carried out. In this study, the sequencing of slabs and optimization of zone temperature set-points are considered at the same time. In addition to the cost function used in the optimization zone temperature set-points, the average time interval between discharging events is accounted. As a new control input beside from the zone temperature set-points, the average time interval between discharging events is subject to a lower bound due to constraint of downstream process, e.g., hot rolling mill.

Particular point in the design MPC for scheduling task is that the number of slab taken into account in each execution of optimization is constant whereas the prediction horizon is variable from one furnace simulation to another according to the time interval.

The simulation in Intel Fortran environment shows that furnace production rate can be increased with a cost on energy consumption. In addition, accuracy of final slab temperature and its homogeneity can be significantly improved. In practice, these results prove that MPC can be used to inform the furnace operator to adjust furnace throughput rate to instantaneous furnace capacity. This information will help the furnace operator to adjust the scheduling of slabs so that furnace can operate at its maximum rate.

Conclusions

The presented results for level 1 of the furnace demonstrate the feasibility and the capabilities of the proposed distributed MPC method. The output coupling model is able to capture the coupling effects between furnace zones with a system of first order transfer functions. Although expected to give a limited optimality compared Chapter 9. Conclusions and perspectives 190 to cooperative algorithm, the non-cooperative algorithm is able to provide satisfactory results for the considered walking-beam reheating furnace.

The MPC designed for level 2 has been tested successfully with data of real system. The MPC controller shows more accurate and homogeneous final slab temperature profiles compared to the results of the existing controller.

As a result, overheating of slab is reduced significantly and consequently the energy efficiency is enhanced.

Despite highly expensive non-steady-state furnace model used in the optimization, the Nelder-Mead simplex algorithm is able to deliver the satisfactory temperature set-points in period of a sample time meaning the proposed MPC controller for level 2 can be implemented in the real furnace.

The slab scheduling can be improved with the MPC method. Instead of using a fixed prediction horizon, we use a variable one whereas number of slabs taken into account is constant. The choice of an average time interval of discharging events as the additional control input is proved to be relevant. The obtained results show that the throughput rate can be increased up to 15 ton/h while temperature heating requirement are improved, especially the homogeneity of slabs. The propose MPC controller will not actually serve as furnace controller of level 2 but to inform furnace operator to modify the slab scheduling so that the furnace can operate at its maximum capacity.

Perspectives

The proposed distributed MPC strategy has been implemented and currently runs in full time in the furnaces of hot rolling mill of Florange plant ArcelorMittal. With the concluding results of MPC level 2 simulated based on the data of the real system, it is planned to install the MPC controller on the considered furnace for experimental validation. On the same hot rolling mill there is another reheating furnace having a similar design but different dimensions of zones. Adaptation of the solution to this new furnace will only involve the parameter modification of the non-steady-state model of the furnace.

In the proposed solution, the non-steady-state model used for furnace simulation and optimization of temperature setpoints is quite time-consuming. This restricts the optimality of the solution. Therefore, we want to use another model which is less time-consuming but is able to capture the transient of the furnace operation. In the presented heat transfer model, around 90% of calculation time is dedicated to radiative heat transfer. Therefore, the simplification of model will involve mainly the calculation of radiative heat transfer in the furnace. Other type of heat transfer like conduction or convection will remains the same as in the developed non-steady state model.

On the other hand, reducing the calculation time of the model will permit the utilization of gradient-based methods which can lead to faster convergence of the optimization problem. For that, the objective function will need to be differentiated by finite-differencing.

Chapter 9. Conclusions and perspectives

191

In addition, simplified model will give more margin on real-time constraints. The control horizon can be increased, which means that longer period of temperature setpoints can be planned at each execution of optimization.

Robustness of the MPC controller level 2 in response to the uncertainty in the model is also and important direction that should be exploited. In practice, the uncertainty can come from unscheduled stop of the furnace which is due to furnace incidents. The recurrent changing of process data should be also analyzed more in simulation studies.

These process data could involve degradation of furnace wall material, change of cooling systems, etc.

In the development of MPC for level 2, we have used constant limits of zone temperatures as constraints in the optimization problems. In real-time implementation, these limits will actually vary over time. At a specific time instant, the limits of zone temperatures depends not only on power available of furnace zones but also on the slab charge located in each zone. To tackle this issue, the furnace model may effectively give estimations of these limits which might then be utilized in the optimization problems.

The large perspective of this work is the implementation of MPC level 2 for all reheating furnaces of ArcelorMittal steel plants, which is expected to generate significant economical and environmental impact.

Other type of reheating furnaces which are dedicated to different type of products such as blooms, billets, tubes, etc., may be studied in order to apply distributed MPC and hierarchical MPC.

Appendix A

BISRA Table 193 Appendix B

Coordinate search algorithm

Let consider the unconstrained minimization problem of a continuously differentiable function

where f (x) :

Coordinate search is a pattern search method. To define the pattern for coordinate search, the basic matrix B is usually the identity matrix, i.e., B = I n . For example, with n = 2, we have

3)

The generating matrix for each iterate k is fixed as C k = C 0 which contains in its columns all possible combination of -1, 0, 1. Thus, C k has p = 3 n columns and contains both I n and -I n . Matrix M can be defined as M = I n , and L consists of the remaining 3 n -2n columns of C. For example, with n = 2,

and

All possible trial points defined by the pattern matrix P = BC 0 , for a given step length ∆ k are depicted in Figure To guide the search for descent point, the algorithm of explanatory moves is given as following algorithm where e i , i = 1, ..., n, denote the unit coordinate vectors:

Explanatory moves algorithm for coordinate search.

The current iterate is given with x k , ∆ k , f (x k ), and B, set s k = 0, ϵ k = 0, and min = f (x k ).

For iteration k, explanatory moves may compute no less than n trial point as the first case (on the top left), and no more than 2n points as the last case (on the bottom right). The updating of step length ∆ k is determined by

Appendix C

Nelder-Mead simplex search algorithm

The simplex algorithm follows the following steps:

• Step 1-Initialization

The initial simplex in n -dimensional space is defined by (n + 1) points {P 0 , P >7 ? 8@ %7 9 7% 7 8 87 9 9 7%9 7 7 %< " 8 9 9 9 " &89 7% 8 7 7 89 7 ()*+,-7 87 % 9 %9 8 < 8 8 9 A 7 9 7%9 9 8 7 8 < %< 7 >@ " 7 7 9 9 7 " ( 78 " ) A 7 9 9 9 9 " 7 8 78 78 898 < 8 < 89 8 87 %< 9 8 7 7 9 % " 8 0 BC-1 9 7 < % " 9 ; 8 7 9 9 9 % " 8 9 7 " " ( 8 7? &8 7 9 @ 8 9 A 9 % 7 8 7> 87% 9 89 8 9 % >8 &89 > 8