Monsieur A Ridha Mahjoub M'

Keywords: Blending Problem, Robust Optimization, Polyhedra, Infeasibility Analysis, Visualization

de mener ce travail. Je tiens à remercier aussi Stefan Janaqi qui m'a recruté pour cette thèse et qui a partagé avec moi ses idées et sa connaissance.

J'adresse mes vifs remerciements à Mesdames Nadia Maïzi et Mériam Chèbre pour avoir bien voulu juger ce travail.

A mes employeurs de la Banque Centrale du Mexique, en particulier à Manuel Galan et Jorge de la Vega vont mes remerciements les plus sincères pour leur soutien.

Finalement, un grand Merci à DDelphine Moussard pour la révision du texte en français.

Ce travail est dédié à Aline Castro, ma collègue et ma femme qui est à l'origine de cette aventure.

Resumé

Le procédé de fabrication de mélanges consiste à déterminer les proportions optimales à mélanger d'un ensemble de composants pour que le produit obtenu vérie un ensemble de spécications. Le problème des mélanges (PM) apparaît typiquement dans l'industrie pétrochimique où plusieurs types de composants bruts sont mélangés pour produire de l'essence et autres huiles à propriétés spéciées.

Le système de mélange est généralement constitué de trois sous-systèmes fonctionnels : le sous-système de planication, un optimiseur en ligne et le sous-système de contrôle. Le sous-système de planication prévoit la production générale de la ranerie et d'une façon plus spécique, il génère la recette objective initiale pour toute opération de mélange programmée. Plusieurs sources d'incertitude (voir ci-dessous) apparaissent tout au long du procédé. L'optimiseur en ligne est alors nécessaire pour mettre à jour la recette cible qui peut devenir sous-optimale, voire irréalisable, à cause de ces perturbations dans le procédé. Le feedback du système est basé sur des mesures des propriétés du mélange actuel et celles des composants qui sont recueillies par des analyseurs en ligne. Enn, le sous-système de contrôle est chargé de régler le taux de débit des composants an d'atteindre la recette cible courante qui a été produite par l'optimiseur en ligne. Dans la Figure 1, nous montrons l'intégration de ces systèmes et le diagramme qui décrit un procédé de mélange typique.

Dans ce travail, nous nous concentrons sur l'analyse du système d'Optimisation en Temps Réel (RTO) qui est formé par l'optimiseur en ligne et le sous-système de contrôle. Une hypothèse principale du travail est que les propriétés d'un mélange sont déterminées comme une combinaison linéaire des propriétés des composants, c'est-à-dire que nous supposons que les lois de mélange sont linéaires. En général elles ne le sont pas, mais elles peuvent être approchées par des fonctions linéaires d'autres propriétés de mélange ; voir, par exemple, [START_REF] Odeh | Modélisation mathématique des propriétés de mélanges : Bsplines et optimisation avec conditions de forme[END_REF] et [START_REF] Gary | Petroleum Rening: Technology and Economics, Fifth Edition[END_REF]. Nous tenons à mentionner que cette hypothèse est loin d'être une simplication excessive du problème. En fait, l'optimiseur en ligne du logiciel Anamel qui est utilisé par certaines unités de mélange au sein du groupe TOTAL (voir [START_REF] Chèbre | Feedback control and optimization for the production of commercial fuels by blending[END_REF]) utilise ce modèle de mélange linéaire. En outre, nous nous concentrons sur la production d'un seul mélange à partir d'un ensemble de composants et nous considérons le cas des composants stockés dans les "running tanks", c'est-à-dire le cas où les bacs des composants sont reliés directement aux procédés en amont, où les composants sont en train d'être fabriqués. Ceux-ci s'écoulent directement dans ces bacs en même temps qu'ils sont mélangées. Une caractéristique importante dans le procédé de mélange est l'existence d'exigences strictes ("bornes dures") sur certaines propriétés du mélange qui sont contrôlées par des normes environnementales, réglementaires et technologiques. Si le mélange produit ne satisfait pas les exigences imposées par les "bornes soft", il y a un coût que nous encourons pour la sur-qualité (giveaway ) du mélange. Néanmoins, si le mélange ne satisfait pas les bornes dures, alors le mélange n'est pas vendable et on doit le remélanger. Le remélange réduit la capacité de la ranerie et il a plusieurs coûts associés, (par exemple, la location du bac, les coûts logistiques, les coûts d'inventaire,...). Donc, il faut considérer les bornes dures comme des contraintes qui doivent être satisfaites. D'autre part, le mélange de pétrole est un procédé complexe où plusieurs facteurs inconnus et incertains ont une incidence sur les propriétés du mélange. Outre l'erreur produite par la linéarisation des lois de mélange, il existe d'autres sources d'incertitude : la précision des instruments de mesure sur les composants et les propriétés ; une connaissance incertaine des propriétés des composants en raison des variations des procédés en amont et des conditions de mélange non contrôlables tels que la température de l'air, l'humidité, etc. Ces sources d'incertitude sont typiques de tout problème RTO (cf. [START_REF] Zhang | Results analysis for trust constrained real-time optimization[END_REF]).

Les bornes dures sur les propriétés du mélange en combinaison avec ces sources d'incertitude sont les caractéristiques idéales pour appliquer les techniques d'Optimisation Robuste (RO) (voir, [START_REF] Ben-Tal | Robust Optimization[END_REF]). L'optimisation robuste est une méthodologie pour faire de l'optimisation avec des données incertaines. Une caractéristique principale dans les techniques de RO est de considérer un modèle d'incertitude déterministe et basé sur des ensembles. Nous ne faisons aucune hypothèse probabiliste sur l'incertitude et la solution obtenue est op-v timale pour toute réalisation de l'incertitude dans un ensemble donné. Les techniques de RO semblent être ad hoc pour résoudre un problème RTO avec incertitude, où la faisabilité du problème est la principale préoccupation.

Plusieurs travaux portent sur la faisabilité et l'optimisation du PM (voir, par ex., [START_REF] Slaoui | Application des techniques mathématiques à la gestion des mélanges : histosplines et optimisation[END_REF][START_REF] Girard | Zonotopes and Mixtures management. International Series of Numerical Mathematics[END_REF][START_REF] White | Gasoline blending optimization cuts use of expensive components[END_REF][START_REF] Szafran | Zonoèdres : de la géométrie algorithmique à la théorie de la séparation[END_REF]). La faisabilité de plusieurs mélanges à produire à partir du même ensemble de composants est considérée dans [START_REF] Lacolle | Modélisation géométrique de la faisabilité de plusieurs mélanges[END_REF]. Plus tard, dans [START_REF] Daoudi | Zonotopes et zonoïdes : études et applications aux processus de la séparation[END_REF], l'auteur est préoccupé aussi par l'incertitude des propriétés. Il développe une analyse de sensibilité de la recette lorsque les erreurs de mesure dans les propriétés des composants varient dans un ellipsoïde. Dans [START_REF] Serpemen | Blending technology key to making new gasolines[END_REF], une méthode d'optimisation en ligne est utilisée pour mettre à jour la recette et réduire la divergence, entre le modèle et le réel, causée par l'incertitude. Dans [START_REF] Hendrix | Finding robust solutions for product design problems[END_REF], une approche géométrique est présentée pour les problèmes de conception de produits et de mélange, mais en considérant seulement l'incertitude due à des imprécisions de mesure. Pour traiter le problème de l'incertitude des propriétés des composants, [START_REF] Vermeer | Blendcontrol system all but eliminates reblends for Canadian rener[END_REF] et [START_REF] Singh | Model-based real-time optimization of automotive gasoline blending operations[END_REF] ont proposé respectivement des modèles RTO linéaire et non-linéaire basés sur les pronostics des propriétés des composants. D'autre part, [START_REF] Wang | On-line optimization model design of gasoline blending system under parametric uncertainty[END_REF] a présenté un modèle de contraintes en probabilités qui est résolu par une méthode heuristique. Plus récemment, [START_REF] Chèbre | Feedback control and optimization for the production of commercial fuels by blending[END_REF] introduit un algorithme de contrôle qui gère ce type d'incertitudes au moyen d'un estimateur de propriétés des composants. Un modèle non-linéaire basé sur la programmation stochastique et qui couvre diérents types d'incertitude est présenté dans [START_REF] Zhang | Real-time optimization under parametric uncertainty: a probability constrained approach[END_REF].

L'objectif principal de ce travail est d'étudier le procédé de mélange et la production de mélanges robustes d'un point de vue polyédral. Plus précisément, nous nous intéressons au développement des outils informatiques pour optimiser en ligne la production de mélanges robustes qui résistent à l'incertitude qui survient lors du procédé de mélange.

Dans le chapitre 2, nous présentons les Polytopes de Mélange inhérents au problème de mélange et l'équation de base qui sera utilisé pour introduire les données mises à jour dans la méthode RTO développée au chapitre 3. Nous modélisons diérentes sources d'incertitude liées au procédé de mélange et nous proposons une méthode robuste RTO basé sur les techniques de RO. La méthode présentée est destinée à éviter le remélange et on mesure sa performance en termes de sur-qualité du mélange produit et de coût de la recette.

En outre, nous analysons la faisabilité d'un mélange à fabriquer à partir d'un ensemble de composants lorsque le fond du bac d'un mélange précédent est utilisé dans la composition du nouveau mélange. Des informations essentielles pour le système de contrôle sont alors produites. Par exemple, savoir le volume minimal à couler avant d'obtenir un mélange faisable robuste réduit les interventions inutiles du système de contrôle. Celles-ci tentent de maintenir le mélange dans les spécications en modiant la recette produite par l'optimiseur en ligne, mais souvent elles sont fondées sur une vision limitée du vi problème.

Enn, au chapitre 4, nous analysons le cas où la méthode RTO ne parvient pas à produire un mélange réalisable. Nous mesurons et produisons des visualisations de l'infaisabilité du PM et proposons une méthode pour guider le procédé de mélange vers le meilleur mélange robuste.

Polytopes de Mélange

Désignons par n et m le nombre de composants et de propriétés respectivement (pour l'application industrielle, m et n ne dépassent pas 20).

Soit S 1 le simplexe déni comme

S 1 = {u ∈ R n | 0 ≤ u ≤ 1, 1 T • u = 1}.
Nous allons représenter une recette de mélange (recette) par un vecteur u ∈ S 1 tel que u j est le pourcentage du composant j présent dans le mélange alors que les propriétés du mélange (propriétés) seront désignées par un vecteur y ∈ R m . De cette façon, l'ensemble des composants sera décrit par la matrice B de taille (m, n), où B i,j est la i ème caractéristique physique ou chimique du composant j.

Lors de la fabrication d'un mélange, il est habituel d'utiliser le fond du bac d'un mélange précédent dans la production du nouveau mélange. Le procédé commence alors avec un volume V 0 d'un mélange précédent avec des propriétés b 0 ∈ R m et continue en ajoutant progressivement un volume v du nouveau mélange pour obtenir un volume du produit nal V total = V 0 + v.

Comme on l'a annoncé, dans ce travail nous supposons que les propriétés du mélange sont une combinaison linéaire (ou plus précisément, une combinaison convexe) des propriétés des composants. En vertu de ces hypothèses, pour tout v ≥ 0 les propriétés de la recette u sont déterminées par l'équation de base dénie par l'application linéaire y : R n → R m tel que :

y(u) = α v b 0 + β v Bu (0.1) avec α v = V 0 V 0 +v et β v = v V 0 +v
. Pour tout horizon de planication xe et déterminé par le volume (v), on peut obtenir la recette de coût minimum comme la solution du Problème de vii mélange (PM) :

min c T • u s.t. u ≤ u ≤ u (0.2)
y ≤ y(u) ≤ y (0.3)

1 T • u = 1, 0 ≤ u ≤ 1. (0.4)
Ici le vecteur c ∈ R n représente les coûts des composants. Une recette u doit vérier des bornes de disponibilité et hydrauliques (voir contraintes sur les composants (0.2)). De la même manière, les contraintes sur les propriétés (0.3) sont déterminées par des bornes inférieures (y) et supérieures (y). Ces bornes des propriétés sont étiquetées comme des bornes hard et soft. Les bornes hard, sont liées aux spécications commerciales et environnementales devant être satisfaites alors que les bornes soft peuvent être violées, mais en produisant un mélange avec de la sur-qualité. Enn, (0.4) sont les contraintes de pourcentage.

Géométriquement, l'intersection des contraintes sur les composants (0.2) avec S 1 détermine le polytope

U = {u ∈ R n | u ≤ u ≤ u , 1 T • u = 1, 0 ≤ u ≤ 1}
qui représente l'ensemble de toutes les recettes qui peuvent être produites. D'autre part, les contraintes sur les propriétés dénissent le polyèdre (car il peut être non borné) :

C U = {u ∈ R n | y ≤ y(u) ≤ y}.
Tous les deux sont dénis dans l'espace des composants.

Dans l'espace des propriétés, on a les images de U et C U par l'équation de base :

P = {y(u) ∈ R m | u ∈ U } et C Y = {y(u) ∈ R m | u ∈ C U }.
[V min , V M ax ] tel que P (v) ∩ C = ∅, ∀v ∈ [V min , V M ax ].
Le PM est réalisable seulement pour ces volumes (voir la Figure 2).

Le fait d'avoir les valeurs V min et V max permet de sélectionner un hori-b0 P (VMin) zon de planication approprié à chaque boucle du RTO. Habituellement, le volume à produire est xé par le département de planication, mais l'imprécision dans le modèle et l'incertitude tout au long du procédé peuvent requérir une révision constante de ces volumes. Choisir un horizon de planication avec un volume v < V min produit un PM irréalisable et rend nécessaire l'utilisation d'une recette alternative avec une éventuelle détérioration de la performance globale de la méthode RTO. De l'autre côté, V max est le volume maximum à mélanger lorsqu'on cherche un mélange qui répond aux spécications. Après que V max a été coulé, P (v) s'éloigne de C et les propriétés du mélange se détériorent progressivement.

P (VMax) C P (∞)
De la même manière, nous pouvons analyser le rôle de l'utilisation de V 0 dans la faisabilité du PM. En faisant u 0 = V 0 v , on peut réécrire l'équation de base d'une manière utile :

y(u) = u 0 b 0 1 + u 0 + Bu 1 + u 0 .
(0.6)

x

En utilisant cette équation, on peut déterminer l'intervalle [u 0M in , u 0M ax ] pour lequel P ∩ C = ∅ (d'une manière équivalente, U ∩ C U = ∅). Ces valeurs nous donnent toutes les combinaisons de volumes (V 0 , v) produisant un mélange réalisable. Désignons par

C U H = {u ∈ R n | y H ≤ y(u) ≤ y H }
le polyèdre déni par les contraintes dures. À toute boucle RTO on procède comme suit :

Algorithme 1 Méthode RTO I On incorpore la nouvelle information disponible (V 0 , b 0 , B et v) dans l'équation de base. Si V 0 = 0, on résout les problèmes (0.8), (0.10) et (0.11) dans cet ordre jusqu'à obtenir une solution réalisable. Sinon, on passe à l'étape (II).

II On calcule V M in (V M ax ) en résolvant le problème de programmation linéaire :

max (min) u 0 (0.7) s.t. u ∈ U ∩ C U u 0 ≥ 0. Si le problème (0.7) est irréalisable ou si v / ∈ [V M in , V M ax ] alors on passe à l'étape (IV).
III On calcule la recette optimale :

min c T • u (0.8) s.t. u ∈ U ∩ C U .
IV On oublie les coûts des composants et on se concentre exclusivement sur les contraintes dures. Donc, on calcule l'intervalle [V M in , V M ax ] en résol-xi vant le problème :

max (min) u 0 (0.9) s.t. u ∈ U ∩ C U H u 0 ≥ 0. Si ce problème est irréalisable ou si v / ∈ [V M in , V M ax ]
, alors on passe à l'étape (VI).

V On calcule la recette u ∈ U ∩ C U H qui produit le mélange avec sur-qualité minimale.

min y(u) -y T 1 (0.10) s.t. u ∈ U ∩ C U H .
VI On trouve u ∈ U pour produire un volume de mélange v avec des propriétés aussi près que possible de y T :

min y(u) -y T 1 (0.11) s.t. u ∈ U.

Incertitude dans le problème de mélange

En suivant ces étapes, la méthode RTO produit toujours une recette qui guide le processus de contrôle. Pour réduire l'imprécision du modèle produit par les incertitudes, la méthode met à jour les propriétés des composants et du mélange courant via l'équation de base. Cependant, l'actualisation du modèle peut échouer même à garantir une solution réalisable. Une des raisons principales de cet échec est l'hypothèse implicite que les données restent inchangées à l'intérieur de chaque boucle RTO. Lorsque les composants sont stockés dans des running tanks, la valeur de B uctue avec le temps parce que les composants sont issus de diérents procédés présentant aussi des perturbations. De plus, et conformément à [START_REF] Zhang | Results analysis for trust constrained real-time optimization[END_REF], l'incertitude dans tout système RTO peut être de quatre types :

1. L'incertitude du procédé : les propriétés des composants, la température, l'humidité, etc.

2. L'incertitude de mesure. Dans ce travail on fait face à l'incertitude due à la linéarisation de la loi des mélanges (incertitude du type 3) à travers de la mise à jour du modèle. On considère également les incertitudes découlant de la mesure des composants et des propriétés du mélange (incertitude du type 2) et l'incertitude causée par la connaissance imprécise des propriétés des composants (incertitude du type 1). Toutefois, d'autres types d'incertitude (e.g., l'incertitude due à des facteurs non contrôlés extérieurs au procédé comme la température ou l'humidité, l'incertitude sur les prix des composants ou l'incertitude dans les disponibilités des composants) se manifestent géométriquement d'une façon similaire et peuvent donc être traités de la même manière.

Pour xer notre méthode RTO contre l'incertitude de mesure et la connaissance imprécise de B, nous utilisons les techniques d'optimisation robuste (voir [START_REF] Ben-Tal | Robust Optimization[END_REF]). L'idée consiste à calculer pour le polytope C U sa région robuste convexe RC U telle que tout point dans cette région résiste à l'incertitude de B et aux erreurs de mesure, c'est-à-dire que tout point u ∈ RC U est garanti de rester à l'intérieur du polytope réel pour toutes les réalisations possibles de B et toute erreur de mesure, dans des limites raisonnables.

Dans la Figure 3, nous montrons les polytopes nominaux U et C U , le "robustope" RC U et le polytope réel C U . La diérence entre le polytope nominal et le polytope réel C U et C U est le résultat d'une diérence entre la matrice nominale B et la matrice réelle B. Tout point dans RC U (par exemple la recette nominale robuste u r ) se transforme en un point dans son voisinage (la vraie recette robuste u R ) en raison des erreurs de mesure. Le robustope RC U est calculé de telle sorte que tout point de son voisinage (la région en pointillés) reste à l'intérieur du polytope réel C U .

Nous pouvons utiliser diérents ensembles et diérentes normes pour modéliser l'incertitude. Le niveau de conservatisme (combien nous voulons être protégés de l'incertitude) et la complexité du problème dépendent de ces choix. Prenons comme ensemble d'incertitude des intervalles et la norme L ∞ , les régions robustes obtenues sont des polytopes et la complexité du modèle est préservée au détriment d'être probablement trop conservateur (nous sommes protégés des pires déviations de toutes les valeurs de B et des plus grandes erreurs de mesure survenant toutes, en même temps).

En modélisant l'incertitude de cette façon, on peut remplacer les polytopes robustes dans la méthode RTO et l'on obtient directement une nouvelle méthode RTO robuste et de même complexité. De plus, tous les outils développés pour l'analyse de l'infaisabilité du PM (voir plus bas), seront valables aussi xiii En ce qui concerne le niveau de conservatisme, on peut comparer le prix de la robustesse (la diérence des coûts entre la recette optimale robuste et la recette optimale nominale) avec le coût du remélange. Si le dernier est inférieur au premier, on peut ajuster les paramétres du modèle pour diminuer le niveau de conservatisme et vice versa. Le coût du remélange est dicile à estimer en pratique, alors on propose d'utiliser le prix de la robustesse comme une valeur de référence. On peut par exemple, comparer le prix de la robustesse avec une fraction du coût de remélange : le coût de la recette à produire pour corriger le mélange fabriqué au préalable et qui est hors spécications. Il peut arriver, comme on le montre dans un exemple réel, que cette valeur soit déjà supérieure au prix de la robustesse.

RC U U C U C U u R u R

Infaisabilité du Problème de Mélange

Une instance du PM (B, v, b 0 , V 0 , c, u, u, y, y) est déterminée par les composants B, le volume v, le fond du bac b 0 et son volume V 0 , les coûts des composantes c, les bornes hydrauliques et de disponibilité des composants u xiv et u et par les bornes sur les propriétés y et y. A chaque itération, la méthode RTO résout le PM an de produire, si possible, la recette de coût minimum. Si le PM est infaisable, certaines bornes sont assouplies et un mélange faisable avec surqualité minimale est recherché. La méthode tient compte des possibles modication de l'instance du PM pour obtenir de meilleurs résultats. Plus précisément, les volumes V 0 et v sont utilisés comme variables de contrôle dans le modèle an de produire une recette plus économique ou faisable lorsque les valeurs originales de V 0 et v ne le permettent pas. Dans certains cas, même avec ce degré de liberté supplémentaire, le PM reste infaisable pour les volumes requis. Alors, pour assurer la continuité de l'opération, la méthode propose de produire le mélange infaisable le plus proche du mélange cible y t .

Dans ce chapitre, on s'intéresse à mesurer et visualiser l'infaisabilité du PM et à l'utilisation d'autres paramètres de l'instance du PM pour contrôler le procédé de mélange. Outre les volumes v et V 0 , la disponibilité des composants, les bornes hydrauliques ou encore les composants peuvent être utilisés pour piloter le processus de mélange. Ici, on analyse exclusivement la modication des bornes de disponibilité et hydrauliques.

Les niveaux d'infaisabilité du problème de mélange

Dénissons le polytope

C U 1 = C U ∩ S 1
dans l'espace des composants et le polytope

C Y 1 = {y(u) ∈ R m | u ∈ C U 1 } dans l'espace des propriétés.
Si le PM est infaisable pour une paire de volumes xes (V 0 , v) et si l'on suppose que U = ∅, alors on peut avoir 3 diérents niveaux d'infaisabilité qui peuvent être énoncées en termes de polytopes de mélange :

Niveau 1. C U 1 = ∅ et U ∩ C U 1 = ∅. Niveau 2. C U = ∅ et C U 1 = ∅. Niveau 3. C U = ∅.
Pour les trois niveaux d'infaisabilité, on peut modier les propriétés des composants (la matrice B) et forcer le polytope C U à avoir une intersection non vide avec le polytope U en produisant de cette manière un mélange faisable. D'autre part, comme le simplex S 1 contient toutes les recettes faisables, alors C U 1 représente l'ensemble de toutes les recettes ayant leurs propriétés ( C Y 1 ) xv dans les spécications. Un polytope U peut atteindre, par des modications des bornes, un point u ∈ C U \ U ssi u ∈ S 1 , c'est-à-dire si et seulement si le PM a une infaisabilité de niveau 1.

Mesure de l'infaisabilité et directions faisables

Pour chaque niveau d'infaisabilité, on peut dénir l'infaisabilité du PM en termes des polytopes de mélange : l'infaisabilité de niveau 1 est dénie par

D U = D(U, C U 1 ) et D Y = D(P, C Y 1 )
l'infaisabilité de niveau 2 par

D U 2 = D(U, C U ) et D Y 2 = D(P, C Y ) et l'infaisabilité de niveau 3 par D Y 3 = D(P, C).
Dénissons une direction faisable dans l'espace des composants comme une direction de translation des polytopes U ou C U produisant sa collision et par conséquent la faisabilité du PM. Une direction de translation des polytopes P ou C Y produisant sa collision est appelée également une direction faisable (dans l'espace des propriétés), car elle peut être vue comme l'image par l'équation de base d'une direction faisable dans l'espace des composants.

Comme on ne considère que des modications des bornes pour transformer d'infaisable en faisable une instance du PM, dans ce qui suit, nous allons seulement développer l'analyse de l'infaisabilité de niveau 1. Dans l'espace des composants on s'intéresse au calcul de l'écart entre les polytopes U et C U 1 . Les deux polytopes sont décrits dans leurs formes duales. Donc, on peut calculer sa distance par la résolution du problème :

min t -u (0.12) s.t. u ∈ U, t ∈ C U 1 . Si u * ∈ U et t * ∈ C U 1 est une solution du problème (0.12), alors D U = d(u * , t * ) = t -u et w U = t * -u * D U (0.
s.t. u ∈ U, t ∈ C U 1 on obtient p * = y(u * ) ∈ P et q * = y(t * ) ∈ C Y 1 tel que D Y = d(p * , q * ). Alors, w Y = q * -p * D Y (0.
H 2 = {u(s) | u(s) = u M + w T • s} avec s = (s 1 , s 2 ) T ∈ R 2 , w = (w 1 , w 2 ) T et u M le milieu du segment [u * , t * ].
Désignons par U 2 et C 2 les coupes de U et C U 1 avec H 2 . Pour produire la visualisation de l'infaisabilité dans l'espace des composants, on dessine les polygones U 2 et C 2 dans le plan H 2 avec u M pris comme centre de coordonnées et les vecteurs de w 1 et w 2 comme axes (voir Figure 4).

Il est important de remarquer les caractéristiques de cette visualisation. Visualisation de l'infaisabilité dans l'espace des propriétés Considérons les directions faisables dans l'espace des propriétés

w 1 = w Y et w 2 = w U Y et denotons H 2 le 2D-plan déni par H 2 = {y(s) ∈ R 2 | y(s) = y M + w T • s} avec s = (s 1 , s 2 ) T ∈ R 2 , w = (
U N EW = {u ∈ R n | u -∆ ≤ u ≤ u + ∆ , 1 T • u = 1, 0 ≤ u ≤ 1}
en choisissant les variables élastiques appropriées ∆, ∆ ∈ R n+ de façon à avoir On observe que, en utilisant les vecteurs u -, u + ∈ R n+ , on peut restreindre les modications des bornes en cherchant la meilleure recette u * dans U + ∩C U (voir Equation 0.15) au lieu de C U 1 .

U N EW ∩ C U 1 = ∅.
U + = {u ∈ R n | 0 ≤ u -u -≤ u ≤ u + u + ≤ 1 , 1 T • u = 1} (0.15)
En outre, dans la Section (4.7) on construit des repères concernant l'impact des bornes des composants sur le coût et la sur-qualité du mélange et sur le volume nécessaire à couler pour avoir le premier mélange faisable. Cette xix information pourrait être utile pour décider si l'on applique ou non les modications des bornes.

Parmi les critères d'optimisation pour chercher la meilleure recette u * ∈ C U 1 on peut minimiser l'infaisaibilité (la distance entre les polytopes correspondants), minimiser les modications des bornes ou l'on peut utiliser les mêmes critères de la méthode RTO. Une fois la recette u * déterminée, on obtient les modications optimales des bornes

∆ * = max{0, u -u * } et ∆ * = max{0, u * -
u ∈ C U 1 u 0 ≥ 0. pour trouver u 0M ax et u M in1 (u 0M in et u M ax1
). Puis, à partir de ces valeurs et de la relation u 0 = V 0 /v on obtient l'intervalle de volumes

[V M in1 , V M ax1 ] pour lesquels C U 1 = ∅. Les mélanges y Min1 = y(u Min1 ) et y Max1 = y(u Max1 )
sont le premier et le dernier mélange, en termes de volume, que l'on peut produire en modiant les bornes des composants. Enn,

∆ * = max{0, u -u Min1 } et ∆ * = max{0, u Min1 -u}
sont les modications minimales des bornes pour atteindre le mélange y Min1 au volume V M in1 . La relation entre ces volumes et l'impact des bornes des composants sur le procédé de mélange peut être mieux expliqué en renvoyant à la Figure 6. P (V M in ) Dans cette gure, le polytope P va de P (0) = b 0 à P (∞) formant le cône :

P (V M in1 ) P 1 (V M in1 ) C Y 1 CP ∩ C C CP 1 CP
CP = y ∈ R m | y = conv V P (∞) ∪ b 0 .
Les points y Min ∈ P (V M in ) ∩ C et y Max ∈ P (V M ax ) ∩ C sont le premier et le dernier mélange que l'on peut produire. D'autre part, comme U ⊆ S 1 alors P (v) ⊆ P 1 (v) = y(S 1 ), ∀v ≥ 0. Par conséquent, P 1 (v) décrit le cône :

CP 1 = y ∈ R m | y = conv V P 1 (∞) ∪ b 0
qui comprend tous les mélanges productibles quand les bornes hydrauliques et de disponibilité ne sont pas respectées. Le premier et le dernier de ces mélanges étant dans les spécications (dans le polytope C) sont les mélanges y Min1 et y Max1 . L'intersection de CP 1 avec C est le polytope C Y 1 et la visualisation dans l'espace des propriétés presentée dans la Section (4.5), c'est une coupe des polytopes

CP (v) et C Y 1 (v) pour un volume v xe.
Enn, supposons que c * et q * (c * 1 et q * 1 ) sont le coût et la sur-qualité minimaux qui peuvent être atteints par un mélange dans

CP ∩ C (dans C Y 1 ). Alors, comme P (v) ⊆ P 1 (v) pour tout v ≥ 0, il résulte que P ∩ C ⊆ C Y 1 . Par conséquent, V M in1 ≤ V M in , c * 1 ≤ c * et q * 1 ≤ q * .
Ainsi, nous pouvons xxii dire que les bornes des composants sont responsables de l'augmentation du coût par (c *c * 1 ), de l'augmentation de la sur-qualité par (q *q * 1 ) et du volume supplémentaire à verser (an d'avoir le premier mélange faisable) par (V M in -V M in1 ). Donc, on contourne la diculté à déterminer les coûts de modication des bornes, en fournissant les indices

K C = c * -c * 1 c * 1 , K Q = q * -q * Contents Introduction (French version

Introduction

The oil blending process consists in determining the optimal proportions to blend from a set of available feedstocks or components such that the nal product obtained fullls a set of specications on their properties. Blending problems (BP) arise typically in the petrochemical industry where several types of raw components are blended to produce gasoline and other oils with specied properties.

The blending system is typically constituted by three functional subsystems: the scheduling subsystem, the on-line optimizer and the control subsystem. The scheduling subsystem is the one in charge of the general renery production planning and specically, it generates the initial target recipe for any blend operation scheduled. Several sources of uncertainty (see below) pervade the process. The on-line optimizer is then required to update the target recipe which may became sub-optimal, or even infeasible, due to these disturbances in the process. The feedback is based on measures of the blends' and components' properties gathered by on-line analyzers. Finally, the control subsystem is in charge to adjust the components ow rates in order to attain the current target recipe produced by the online optimizer. In Figure 1.1, we show the integration of these systems and the ow chart which describes a typical blending process.

In this work we concentrate on the analysis of the Real Time Optimization (RTO) system formed by the on-line optimizer and the control subsystem. A main assumption here is that the properties of a blend are determined as a linear combination of the components properties, that is, we suppose that properties blend linearly. In general they do not, but they can be replaced by approximated by linear functions of other blend's properties; see, e.g., [START_REF] Odeh | Modélisation mathématique des propriétés de mélanges : Bsplines et optimisation avec conditions de forme[END_REF] and [START_REF] Gary | Petroleum Rening: Technology and Economics, Fifth Edition[END_REF]. We would like to mention that this assumption is by no means an over simplication of the problem. Actually, the linear blending model is on the basis of the on-line optimizer of the Anamel software which is used by some blending units within the TOTAL group (see [START_REF] Chèbre | Feedback control and optimization for the production of commercial fuels by blending[END_REF]). Additionally, we focus on the production of a single blend from a set of components and we consider the case when the components are stocked in "running" tanks. That is, when the components' tanks are connected directly to upstream processes where the components are being produced. These ones are owing into their tanks at the same time that they are blended. An important characteristic in the oil blending process is the strict requirements ("hard bounds") over some blend's properties controlled by environmental, legal and technological standards. If the blend does not satisfy these hard bounds, then it is not sellable and one must reblend it. Reblending reduces the renery's capacity with several involved costs (e.g., the tank rent, the logistic and inventory costs,...). So, one should consider these blends' bounds as hard constraints which must be satised.

On the other hand there are "soft bounds". If the blend produced does not fulll the requirements imposed by these bounds, then it is better than what is required and there is a cost we incur for the blend's quality giveaway.

The oil blending is a complex process where several unknown and uncertain factors aect the blend's properties. In addition to the plant-model error produced by the linearization of blending laws, there are other sources of uncertainty: measure errors on components and properties caused by instruments' precision; uncertain knowledge of the components' properties due to upstream process variations and uncontrolled blending conditions such as air temperature, humidity, etc. All these are typical uncertainties arising in any RTO problem (cf. [START_REF] Zhang | Results analysis for trust constrained real-time optimization[END_REF]).

Hard constraints on the blends' properties combined to these sources of uncertainty are the ideal characteristics to apply the Robust Optimization (RO) techniques (see [START_REF] Ben-Tal | Robust Optimization[END_REF]). Robust Optimization is an approach to optimize under uncertainty. A main issue in the RO techniques is to consider a deterministic and set-based uncertainty model. No probabilistic assumption is made over the uncertainty and the solution obtained is optimal for any realization of the uncertainty in a given set. RO techniques seem to be ad hoc to address an uncertain RTO problem where the feasibility is the primary concern.

Several previous works (see, e.g., [START_REF] Slaoui | Application des techniques mathématiques à la gestion des mélanges : histosplines et optimisation[END_REF][START_REF] Girard | Zonotopes and Mixtures management. International Series of Numerical Mathematics[END_REF][START_REF] Szafran | Zonoèdres : de la géométrie algorithmique à la théorie de la séparation[END_REF][START_REF] White | Gasoline blending optimization cuts use of expensive components[END_REF]) deal with the BP feasibility and optimization. The feasibility of several blends to be produced from the same set of components is considered in [START_REF] Lacolle | Modélisation géométrique de la faisabilité de plusieurs mélanges[END_REF]. Later, in [START_REF] Daoudi | Zonotopes et zonoïdes : études et applications aux processus de la séparation[END_REF] the author is concerned with the properties' uncertainty too. He develops a sensitivity analysis of the recipe when measurement errors in the components properties vary in an ellipsoid. In [START_REF] Serpemen | Blending technology key to making new gasolines[END_REF] an on-line optimization method is used to update the recipe and reduce the plant/model mismatch caused by uncertainty. In [START_REF] Hendrix | Finding robust solutions for product design problems[END_REF], a geometric approach is presented for the product and mixture design problems but only considering the uncertainty due to measurement imprecisions. To deal with components' properties uncertainty, [START_REF] Vermeer | Blendcontrol system all but eliminates reblends for Canadian rener[END_REF] and [START_REF] Singh | Model-based real-time optimization of automotive gasoline blending operations[END_REF] proposed a linear and non-linear blending RTO systems based on predictions of feedstocks properties whereas [START_REF] Wang | On-line optimization model design of gasoline blending system under parametric uncertainty[END_REF] presented a chance constraint model and an hybrid: neural networks -genetic algorithm solution. More recently, [START_REF] Chèbre | Feedback control and optimization for the production of commercial fuels by blending[END_REF] introduced a linear blending control algorithm which handles this type of uncertainties via an estimator of the components' properties. A more general method for the non-linear blending model based on stochastic programming which covers various types of uncertainty is presented in [START_REF] Zhang | Real-time optimization under parametric uncertainty: a probability constrained approach[END_REF].

This work tries to answer a number of real world questions coming from the petroleum blending industry. The answers to these questions are required to optimize and control, in real time, the production of robust blends. This explains the need for developing ecient computational tools. Another important point of this work is the generation of supportive visualizations of the objects under analysis which are intended to help the decision-making.

In Chapter 2 we present the Blending Polytopes inherent to the Blending Problem and the Basic Equation which will be used to introduce the updated data in the RTO method developed in Chapter 3. Here we analyze the feasibility of a blend to be produced from a set of components when the heel of a previous blend is used in the composition of the new blend. Critical information for the control system is then produced. For instance, knowing the minimal volume to pour before getting a robust feasible blend reduces the unnecessary control system interventions. These ones try to maintain the blend within specications by modifying the recipe produced by the on-line optimizer but often based on a limited view of the problem. Later we model dierent sources of uncertainty arising in the blending process and we propose a robust RTO method based on the RO techniques. The method presented is intended to avoid reblending and we measure its performance in terms of the blend quality giveaway and the recipe's cost. We propose to compare the price of robustness (the cost dierence between the robust and the nominal recipes) with the reblending cost in order to adjust the level of conservatism in the model. Finally, some numerical examples are presented.

In Chapter 4 we analyze the case when the RTO method fails to produce a feasible blend. We exploit the Blending Polytopes to measure, classify and visualize the BP infeasibility. A ne analysis of the components bounds modications is conducted to guide the process towards the best robust blend. A set of indices and visualizations provide a helpful support for the decision maker.

Finally, in the Appendix A we present a customized algorithm for computing the vertices of the blending polytopes and in the Appendix B we present some screenshots of the visualization system VISAMEL in which we have implemented some of the ideas presented in this work.

Chapter 2

Blending Polytopes

After a set of preliminary denitions, four intrinsic polytopes of the blending problem are introduced. These polytopes constitute the core of this work because they dene a natural framework for the study of the feasibility and infeasibility of the BP. The distance between the pairs of polytopes gives a measure of the BP infeasibility and determines a set of signicant points and polygons allowing its visualization.

Preliminaries

We briey review some of the basic denitions on polytope theory and convexity used in this work. A complete treatment of these elds can be found in [START_REF] Ziegler | Lectures on Polytopes (Graduate Texts in Mathematics)[END_REF] and [START_REF] Webster | Convexity[END_REF]. For the basic notions on complexity theory we refer to [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF].

Throughout this work we will denote by x p the L p norm of a vector x ∈ R d . Three important cases are the L 1 norm

x 1 = d i=1 |x i |,
the Euclidean norm

x 2 = d i=1 x 2 i and the L ∞ norm x ∞ = max i=1,...,d |x i |.
For S ⊆ R d , the convex hull conv(S) of S is the smallest convex set in R d containing S. A V-polytope is the convex hull of a nite set of points in some R d . A subset P of R d is called an H-polyhedron if it is the intersection of a nite set of halfspaces in some R d . An H-polytope is an H-polyhedron that is bounded. A polytope is a set P ⊆ R d which can be presented either as a V-polytope or as an H-polytope.

The dimension of a polytope is the dimension of its ane hull and a dpolytope is a polytope of dimension d in some R e (e ≥ d).

Chapter 2. Blending Polytopes

Let P ⊆ R d be a polytope. A linear inequality ax ≤ b is valid for P if it is satised for all x ∈ P . A face of P is any set of the form

F = P ∩ {x ∈ R d | ax = b}
where ax ≤ b is a valid inequality for P . The dimension of a face of P is the dimension of its ane hull: dim(F ) = dim(a(F )). The faces of dimensions 0, 1, dim(P ) -2, and dim(P ) -1 are called vertices, edges, ridges and facets respectively. We denote the set of vertices of P by V .

Every polytope has two representations. The primal representation as the convex hull of the nite set of its vertices:

P P = {v ∈ R d | v = s k=1 λ k v k , s k=1 λ k = 1, λ k ≥ 0, k = 1, . . . , s}
and the dual representation as the intersection of a nite set of half-spaces:

P D = {x ∈ R d | Ax ≤ b}.
The question of how a polyhedron is represented has very important practical consequences. As a simple example consider the d-cube:

C d = {x ∈ R d | -1 ≤ x i ≤ 1, i = 1, ..., d}
and the dual representation is described by 2d linear inequalities, while to describe it in its primal form, 2 d points are required (the d-vectors all of whose components are 1 or -1). Thus the size of the two representations diers greatly, for large d. This fact has important consequences in practice when computing time is a crucial issue.

A d-polytope is called simple if each vertex is contained in exactly d facets and is called simplicial if each facet contains exactly d vertices. A vertex lying on more than d facets is known as a degenerate vertex. In Figure 2.1 we sketch the 2-simplex

S 1 = {x ∈ R 3 | 0 ≤ x ≤ 1, x 1 +x 2 +x 3 = 1} = conv{(1, 0, 0), (0, 1, 0), (0, 0, 1)}
and the 3-cube

C 3 = {x ∈ R 3 | -1 ≤ x ≤ 1} = conv (-1,-1,-1), (-1,-1,1), (-1,1,-1), (-1,1,1), (1,-1,-1), (1,-1,1), (1,1,-1), (1,1,1)
in R 3 . Both polytopes are simple but only the simplex is simplicial. Transforming a polytope from its primal form to its dual form and vice versa are known as the facet enumeration (FE) and the vertex enumeration (VE) problems respectively. These problems are equivalent under polar duality: From a polytope P given in its dual form, we obtain directly the primal representation of its dual (P △ ). Then, solving the facet enumeration problem for P △ we obtain its dual representation and from this we get the primal form of P . Since the number of vertices of a d-polytope P with m facets is O(m ⌊d/2⌋ ) [START_REF] Seidel | The upper bound theorem for polytopes: an easy proof of its asymptotic version[END_REF], then it is common to measure the eciency of a vertex enumeration algorithm in terms of both the size of the input and the size of the output (see e.g., [42]). Under this scheme, a VE algorithm is polynomial if its running time is bounded by a polynomial in d, m and the number of vertices v.

(1,0,0) (0,1,0) (0,0,1) (1,1,1) (1,1,-1) (-1,1,1) (-1,1,-1) (1,-1,-1) (-1,-1,1) (-1,-1,-1) (1,-1,1)
The vertex enumeration problem for polyhedra is known to be a NP-Hard problem [START_REF] Khachiyan | Generating All Vertices of a Polyhedron Is Hard[END_REF] and for general polytopes the complexity is an open problem. However, for some particular types of polytopes the vertex enumeration can be made eciently. For instance, if P is a polytope simple or a 0/1-polytope, we can nd its vertices in linear time on v using the reverse search method [START_REF] Avis | A pivoting algorithm for convex hulls and vertex enumeration of arrangements and polyhedra[END_REF] and a backtracking approach [11] respectively. A 0/1-polytope in R d is one having all its vertices in {0, 1} d .

Later we will need to compute the vertex set for the blending polytopes introduced in Section 2.3. These are not simple polytopes but their structure and their size can be exploited by a backtracking algorithm which we present in the Appendix A.

Distance Between Polytopes

The distance between the polytopes A and B is given by

D(A, B) = min{ x -y | x ∈ A, y ∈ B}
with the most suitable norm for the context. Here we will need to compute D(A, B) for the Euclidean norm and the L 1 norm. So, in the following we discuss the tractability of the distance problem for these norms and both representations of the polytopes A and B.

Let's assume that A p , B p , A d and B d are the primal and dual forms of the polytopes A and B. If A and B are given in their dual representation and the distance is measured by the Euclidean norm, then the distance problem is a convex quadratic programming problem:

min{ x -y 2 2 s.t. x ∈ A d , y ∈ B d }.
If the L 1 norm is used, then D(A, B) can be obtained by solving the linear programming problem

min 1 T s (2.1) s.t. s ≥ x -y s ≥ y -x x ∈ A d , y ∈ B d .
On the other hand, if A and B are given in their primal form, we can proceed in a dierent way. By the separating hyperplane theorem (see, e.g., [START_REF] Webster | Convexity[END_REF]), A ∩ B = ∅ if and only if there is an hyperplane

H(w, b) = {x ∈ R d | w T • x -b = 0} strictly separating A from B, i.e., such that w T • x -b ≥ 0 ∀x ∈ A and w T • x -b < 0 ∀x ∈ B.
We dene the gap of H(w, b) as the sum of the shortest distances to H from the nearest point in each polytope. In Figure 2.2, we show two hyperplanes separating the polytopes A and B. The hyperplane of maximum gap is H * and its gap (represented as the distance between the parallel dotted lines to

H * ) is d(x i , H * ) + d(y j , H * ) = D(A, B).
Then, as shown in [START_REF]Computation of certain measures of proximity between convex polytopes: a complexity viewpoint[END_REF] and the development of support vector machines [START_REF] Burges | A Tutorial on Support Vector Machines for Pattern Recognition[END_REF], nding the distance between A and B (using the norm • ) is equivalent to computing the separating hyperplane with the largest gap: where,

(i) A ∩ B = ∅ ⇔ ∃ w ∈ R d , b ∈ R such that (w, b) ∈ H. (ii) A ∩ B = ∅ ⇒ D(A, B) = 2/min{ w ′ | (w, b) ∈ H} x i A y j B H * H
H = {(w, b) ∈ R d+1 | w T • x -b ≥ 1 ∀x ∈ A, w T • x -b ≤ -1 ∀x ∈ B}
is the family of hyperplanes separating strictly A from B and • ′ is the dual norm to • . Using the Euclidean norm, D(A, B) reduces to solve a quadratic problem and if the L 1 norm is used, the resulting problem is min{ w ∞ | (w, b) ∈ H} which can also be stated as a linear program similar to problem (2.1) (see [START_REF]Linear Programming (Series of Books in the Mathematical Sciences)[END_REF]). Finally, knowing that the convex quadratic programming problem can be solved in polynomial time [START_REF] Kozlov | Polynomial Solvability of Convex Quadratic Programming[END_REF], it follows that computing D(A, B) can be done eciently either if both polytopes are given in their primal forms or in their dual forms.

Blending Polytopes

Let's denote by n and m the number of components and properties respectively (usually, both m and n are lower than 20).

Let S 1 be the simplex dened as

S 1 = {u ∈ R n | 0 ≤ u ≤ 1, 1 T • u = 1}.
We will represent a blend's recipe (recipe in short) by a vector u ∈ S 1 such that u j is the percentage of the component j present in the blend whereas the blend's properties (properties) will be denoted by a vector y ∈ R m . Then the set of components is described by the m × n matrix B where, B i,j is the i th When blending it is customary to use the heel of a previous blend to produce a new one. The process starts then with a volume V 0 of a previous blend with properties b 0 ∈ R m and continues adding gradually a volume v of the new blend to obtain a nal product's volume

V total = V 0 + v.
In this work we suppose that blend's properties are linear combination (or more precisely, a convex combination) of components' properties. Under this assumptions, for any v ≥ 0 the blend's properties of the recipe u are determined by the basic equation dened as the linear mapping y : R n → R m such that:

y(u) = α v b 0 + β v Bu (2.2)
where,

α v = V 0 V 0 +v and β v = v V 0 +v
. Having poured a blend's volume v i with properties b i , we can use the basic equation to tune up the production of the remaining volume to blend. To do this, we just need to update b 0 and V 0 and select the new volume to plan:

v ∈ [0, v -v i ]
. This is the kernel of the RTO method presented in Section 3.2. For this planning horizon (v), we can search the optimal cost recipe by solving the Blending Problem (BP):

min c T • u s.t. u ≤ u ≤ u (2.3) y ≤ y(u) ≤ y (2.4) 1 T • u = 1, 0 ≤ u ≤ 1.
(2.5)

Here the vector c ∈ R n represents the components' costs. A recipe u is subject to components availability and hydraulic bounds which impose the components constraints (2.3). In a similar way, the properties constraints (2.4) are determined by quality lower (y) and upper (y) bounds. These properties bounds are labeled as hard (y H , y H ) or soft (y S , y S ). Hard bounds are related to commercial and environmental specications which must be satised whereas soft bounds can be violated but incurring into quality giveaway. Finally, (2.5) are the percentage constraints.

Geometrically, the intersection of the regions determined by the components constraints and S 1 denes the recipe's polytope :

U = {u ∈ R n | u ≤ u ≤ u , 1 T • u = 1, 0 ≤ u ≤ 1}
which represents the set of all the recipes that can be produced. On the other hand, the properties constraints dene the polyhedron (as it may be unbounded):

C U = {u ∈ R n | y ≤ y(u) ≤ y}.
Both are dened in the components space.

In properties space, the images of U and C U under the basic equation are the polytopes

P = {y(u) ∈ R m | u ∈ U } and C Y = {y(u) ∈ R m | u ∈ C U }.
It is important to note here, that the last two polytopes are not given in a primal or dual form. Maybe, this is the reason why they were not studied previously. We show further how the basic equation and the structure of U and C U are exploited in order to treat them.

Finally, the target polytope is dened by

C = {y ∈ R m | y ≤ y ≤ y}.
In this polytope we distinguish the target blend y T which is a blend provided by the scheduling department having the minimal quality giveaway and other desirable characteristics. This blend will be considered as the ideal blend to produce when the BP fails to generate a feasible recipe.

Any point in P represents the properties of a blend issued from one or more recipes in U whereas C depicts the set of properties within specications. The polyhedron C U contains all the points (not necesarily recipes) in the components space that are mapped by y onto properties within specications forming the polytope C Y ⊆ C.

We notice that the polytopes, P , C U and C Y depend on v. For these polytopes as for the basic equation we omit the v when it is considered to be xed.

At this point we have identied the polytopes underlying the BP and we observe that for a xed volume v, a feasible solution of BP is a recipe

u ∈ U ∩ C U with properties y(u) ∈ P ∩ C Y . Equivalently, BP is infeasible ⇐⇒ P ∩ C Y = ∅ ⇐⇒ U ∩ C U = ∅.
(2.6)

Example 2.1 Blending Polytopes

Let's consider the production of (V total = 10 m 3 ) of blend from two components and with only one property to control. The components are determined by the matrix B = (4, 8) and their bounds are u = (0, 0.3) T and u = (0.6, 1) T .

We are not using any previous blend (V 0 = 0) then we need to produce (v = 10 m 3 ) of fresh blend. We suppose that the blend's property must belong to the interval [y, y] = [START_REF] Ben-Tal | Robust Optimization[END_REF]11] and that the blend with minimal quality giveaway is y T = 9.

In Figure 2.3 we present the blending polytopes in both spaces. In compo- nents space (R 2 ),we exhibit the polytope

U = {(u 1 , u 2 ) ∈ R 2 | 0 ≤ u 1 ≤ 0.6, 0.3 ≤ u 2 ≤ 1, u 1 + u 2 = 1}
and the unbounded polyhedron

C U = {(u 1 , u 2 ) ∈ R 2 | 9 ≤ 4u 1 + 8u 2 ≤ 11}.
In properties space (R 1 ), we show the polytope P and the polytope C Y which coincides this time with the target polytope C. Here, the vertices of U are the points (0.6, 0.4) and (0, 1) which transforms under y in the vertices of P : 5.6 and 8. The gap between the polytopes U and C U and between P and C Y show us the infeasibility of this BP instance. Moreover, in this case the polytope C U does not contain any recipe as C U ∩ S 1 = ∅. Finally, we observe that using the recipe (0, 1), that is, producing the blend with 100% of the second component, we obtain the closest blend to y T . This infeasible blend has then the minimal quality giveaway among all the blends that can be produced.

Chapter 3

Robust RTO Method

The blending system is typically constituted by three functional subsystems: the scheduling subsystem, the on-line optimizer and the control subsystem.

The scheduling subsystem purpose is to plan the renery operations for different time horizons. In the long-term, production targets are determined using long-term forecasts of product demands and oil prices. In the mediumterm planning, these targets are revised and adjusted regarding the estimates of components and tanks availabilities. Finally, in the short-term the planning concentrates on the production of target recipes for any blend operation scheduled. The volume to produce, the components to blend and all the other parameters needed to compute the target recipe are determined at this level and according to the general renery planning. On these basis, one should expect the BP instance to have always a feasible solution. However, the target recipe computation is most of the times based on predictions and occasionally on measurements of the components properties. In both cases the accuracy of the model is aected: in the rst case the model inherits the imprecisions of the predictions and in the second case because the components measurements can be also inaccurate and delayed [START_REF] Chèbre | Feedback control and optimization for the production of commercial fuels by blending[END_REF]. Actually, as the components are in production at the same time that they are blended, therefore, the imprecision in the model induced by the defective knowledge of the components properties remains during all the process. For this and other types of uncertainties, the target recipe may became sub-optimal or even infeasible.

In order to reduce the mismatch between the model and the real system, the on-line optimizer updates the target recipe in a closed loop cycle using the on-line measures of the blend's and components' properties. The updated recipe will serve the controller to guide the blending process.

In this chapter we focus on the RTO system formed by the on-line optimizer and the control subsystem. Specically, in Section 3.1 we analyze the blending process when the heel of a previous blend is used in the composition of the new blend. This analysis turns out to be very important for the development of the RTO method presented in Section 3.2. The RTO method is based on the basic equation and the blending polytopes, it is intended to avoid reblending and we measure its performance in terms of the blend quality giveaway and the feedstocks prices. In Section 3.3, we analyze from a geometric point of view several types of uncertainty arising in the blending process and we outline how to x the RTO method against these types of uncertainty. In 16 Chapter 3. Robust RTO Method Section 3.4 we construct the corresponding uncertainty sets for these types of uncertainty using the robust optimization techniques. These sets determine the robust regions for the blending polytopes in which the new robust RTO (RRTO) method is based. Finally, in Section 3.5 we present some numerical examples to illustrate and compare some key aspects of the RTO method and its robust counterpart. The dierence between the robust and the nominal recipes (the price of robustess) is compared with a component of the reblending cost. This iformation can be used to adjust the level of conservatism in the model.

Blending from a previous blend

If the blend is produced from scratch or the volume to pour is too big compared with the heel's volume, that is, if V 0 = 0 or v → ∞ then the basic equation reduces to y(u) = Bu and the polytope P (v) does not depend on V 0 . However, the fact of using the actual blend's properties b 0 in the basic equation to update the model inside the RTO loop makes blending from a previous blend the regular case. Thus, we have V 0 = 0 and P (v) is a polytope homothetic to P (∞) with b 0 as center of homothecy. In properties space the polytope C doesn't depend on v. There exists then an interval (possibly empty) of volumes

b 0 P (V M in ) P (V M ax ) C P (∞)
[V M in , V M ax ] such that P (v) ∩ C = ∅, ∀v ∈ [V M in , V M ax ].
The BP is feasible only for these volumes (see Figures 3.1 and 3.2).

Having the values of V M in and V M ax helps to select the appropriate plan-

b 0 C P (V M in ) P (∞) Figure 3.2: Polytope P (v) with b 0 ∈ P (∞).
ning horizon at each RTO loop. Usually, the volume to produce is xed by the scheduling department but the inaccuracy in the model and the uncertainty pervading the process may require a constant review of these volumes. Choosing a planning horizon by taking a volume v < V M in produces an infeasible BP and makes it necessary to generate an alternative recipe with a possible deterioration of the overall performance of the RTO method. On the other hand, V M ax is the maximum volume to blend when looking for a blend within specications. After V M ax has been poured, P (v) moves away from C and the blend's properties deteriorate progressively.

In practice, when a blend is o-spec, the control system tries to correct it and make it to be within specications. To do this, a small subset of components is used (the more economical or those having the properties to correct the blend faster) and often based on a limited view of the problem. A bonus of knowing V M in and V M ax is that it helps the control system to optimize their interventions.

In a similar way, we can analyze the role of the use of V 0 in the feasibility of the BP. From the basic equation, we notice that letting xed the volume v, the polytope P moves from P (∞) to b 0 as the volume V 0 increases. In this case, there exists also an interval (possibly empty) of volumes [V 0M in , V 0M ax ] of the previous blend allowing us to produce a blend within specications. Making u 0 = V 0 v , we can rewrite the basic equation in an useful way:

y(u) = u 0 b 0 1 + u 0 + Bu 1 + u 0 .
(3.1) Using this equation, we can determine the interval [u 0M in , u 0M ax ] for which P ∩ C = ∅ (and equivalently U ∩ C U = ∅). These values give us all the combinations of volumes (V 0 , v) producing a feasible blend.

Finally, we dedicate some attention to the cost of reblending a blend which failed to be within specications. The reblending cost has several components, some of which are dicult to estimate. First of all, reblending aects the renery's scheduling and particularly, the sales program. The blend o-spec may require to be transfered to a dierent tank and needs to be stocked until its reprocessing. This generates logistic, storage and inventory expenses. Additionally, reblending may demand the use of non standard components with the suitable properties in order to rectify the blend o-spec. Using additives or modifying the upstream processes to give the components the desirable characteristics which they lack, increase their costs. Hence, the cost dierence between the optimal recipes computed when the heel's volume is used and when it is not, is also a component of the reblending cost. Here, we continue the example 2.1 but considering the use of (6m 3 ) of a previous blend with the property b 0 = 15. In this case, the polytopes U and C Y do not change but the polytopes C U and P do.

C U = {(u 1 , u 2 ) ∈ R 2 | 9 ≤ 6 10 × 15 + 4 10 (4u 1 + 8u 2 ) ≤ 11}
while the vertices of P are y(0.6, 0.4) = 11.24 and y(0, 1) = 12.2.

In Figure 3.3 we show these polytopes and we observe that the BP for the new volumes is again infeasible. Moreover, looking the blending polytopes in the properties space we can deduce that we are taking too much volume from the previous blend or that we need to pour some extra volume in order to have a feasible blend. So, let's compute the extreme values of V 0 and v for which the BP is feasible.

If we use the basic equation 3.1 to describe the polytope C U in terms of u 0 and u, we have

9 -Bu 6 ≤ u 0 ≤ 11 -Bu 4 .
Thus, the minimum and maximum values of u 0 such that U ∩ C U = ∅ are obtained computing u M = max u∈U {Bu} and u m = min u∈U {Bu} respectively.

Here, we can nd directly these values by evaluating Bu in the vertices of U . Then, we get u M = B(0, 1) T = 8 and u m = B(0.6, 0.4) T = 5.6. Hence,

u 0M in = max{0, 9-8 6 } = 1
6 and u 0M ax = max{0, 11-5. 6 4

} = 1.35.
Now, by means of the relation u 0 = V 0 v we can select a suitable pair of volumes (V 0 , v) producing a feasible blend. For instance, using all the available previous blend (V 0 = 6 m 3 ), we need to produce at least V M in = 4.44 m 3 and no more than V M ax = 36 m 3 in order to have a feasible blend. Otherwise, producing v = 4 m 3 constraints the use of V 0 as follows: V 0M in = 0.66 m 3 and V 0M ax = 5.4 m 3 .

In Figure 3.4, we depict the evolution of polytope P from b 0 to P (∞). All the polytopes P producing a feasible blend are in the shaded region. The rst (P f ) and last (P l ) touching C, are obtained from any pair of volumes (V 0 , v) satisfying V 0 v = u 0M ax and V 0 v = u 0M in respectively.

RTO Method

Let's denote by

C U H = {u ∈ R n | y H ≤ y(u) ≤ y H }
the polyhedron dened by the hard constraints.

At any RTO loop we proceed as follows: II Compute V M in (V M ax ) by solving the linear programming problem:

max (min) u 0 (3.2) s.t. u ∈ U ∩ C U u 0 ≥ 0. If Problem (3.2) is infeasible or if v / ∈ [V M in , V M ax ] then go to step IV.
III Compute the optimal recipe:

min c T • u (3.3) s.t. u ∈ U ∩ C U .
IV Forget the components costs and focus only on the hard constraints. Com-3.2. RTO Method 21 pute the interval [V M in , V M ax ] by solving the problem:

max (min) u 0 (3.4) s.t. u ∈ U ∩ C U H u 0 ≥ 0.

If this problem is infeasible or if v /

∈ [V M in , V M ax ] then go to step VI.

V Compute the recipe u ∈ U ∩ C U H producing the blend with minimal quality giveaway.

min y(u)y T 1

(3.5) s.t. u ∈ U ∩ C U H .
VI Find u ∈ U to produce a volume v of blend with properties as near as possible to y T :

min y(u) -y T 1 (3.6) s.t. u ∈ U.
Following these steps, the RTO method always produces a recipe that guides the control process (see Figure 3.5). In step I, V 0 = 0 and solving the Problems (3.3), (3.5) and (3.6) we obtain the feasible recipe of minimal cost, the blend with minimal quality giveaway or the closest o-spec blend to y T . If BP is feasible for some volume v ∈ [V M in , V M ax ] then the method generates the feasible blend of minimal cost (step III). Otherwise, the hard constraints become the priority and the method searches a blend satisfying them while minimizing the blend's quality giveaway (step V). Finally, if there is no blend satisfying the hard constraints, the recipe producing the blend with properties as close as possible to y T is proposed (step VI).

Notice that if we can let the blend's volume to be free in some interval

[v, v] ⊆ [V M in , V M ax ]
then solving Problem (3.3) with v free and the additional constraint:

V 0 /v ≤ u 0 ≤ V 0 /v
we get at the same time the optimal blending horizon v * ∈ [v, v] and the optimal recipe for this volume. Similarly, in steps V and VI, we can obtain

I. Feedback (V 0 ,b 0 ,B,v) II. Compute [V M in , V M ax ] III. Solve min c T • u s.t. u ∈ U ∩ C U . IV. Compute [V M inH , V M axH ] V. Solve min ||y(u) -y T || 1 s.t. u ∈ U ∩ C U H . VI. Solve min ||y(u) -y T || 1 s.t. u ∈ U.

Minimal cost recipe u *

Minimal quality giveaway blend y(u * ) Closest blend y(u * ) to y T V 0 = 0 the optimal volume v * ∈ [v, v] and the blend with minimal quality giveaway in that range. To do this, we perform a dichotomic search over v solving at each step (with v xed) the Problems (3.5) and (3.6) respectively. We nish this section by stressing the possible application of the RTO method to determine the appropriate heel's volume to use at the beginning of the blending process. In our analysis, V 0 is considered as xed but we can solve the Problem 3.3 plus the constraints

V 0 == 0 v ∈ [V M in , V M ax ] v / ∈ [V M in , V M ax ] Infeasible Feasible v ∈ [V M inH , V M axH ] v / ∈ [V M inH , V M axH ] Infeasible Feasible Feasible
u 0M in ≤ u 0 ≤ u 0M ax (3.7) V 0 ≥ 0, v ≥ 0 in order to obtain the cheapest recipe with u * 0 ∈ [u 0M in , u 0M ax ].
Then we can nd (by means of the relation u 0 = V 0 v ) a suitable pair of volumes (V 0 , v).

BP Uncertainty

In the previous section we proposed a blend RTO method based on the Blending Polytopes and its evolution with the blended volume. Looking to reduce the model deviations produced by some uncertainties, the method uses blend's and components' properties updates to feedback the model via the basic equa-tion. However, model updating may fail to guarantee even a feasible solution.

A main reason of this failure is the implicit assumption that data remains unchanged inside each RTO loop.

When blending from running tanks, B ′ s values uctuate with time because the components are issued from dierent processes presenting also perturbations. To address this problem, [START_REF] Singh | Model-based real-time optimization of automotive gasoline blending operations[END_REF] proposed a blending RTO method which updates the model with predictions of the components' properties. Although this method improves the model accuracy, it continues to be non robust as it depends on the quality of the predictions. Moreover, the uncertainty in the blending process aects other than B ′ s values as we will see below.

In accordance with [START_REF] Zhang | Results analysis for trust constrained real-time optimization[END_REF], uncertainty in any RTO system may be of four types:

1. Process uncertainty: components properties, external factors (temperature, humidity, . . . ).

2. Measurement uncertainty.

3. Model uncertainty.

4. Market uncertainty: components availabilities and prices, blends demands, etc.

In this work we address the model uncertainty (type 3) through the RTO model update schema. We consider also the uncertainties arising from components measurement and blend's properties measurement (type 2) and the uncertainty caused by the imprecise knowledge of the components properties (type 1).

Measurement and components properties uncertainties manifest geometrically in dierent ways: For the rst type, the real recipe (the real percentages of each component used in the blend) and its properties may dier from the nominal ones. So, the real recipe and its properties are located in neighborhoods of the nominal recipe and its properties respectively. For the second type, the real matrix B (the real components' properties) diers from the nominal matrix B. Hence the real polytope C U is dierent from the nominal polytope C U . In both cases, when a nominal feasible recipe u ∈ U ∩ C U is computed, the real recipe may fall outside the polytopes U and C U whereas the real blend's properties may be outside the polytopes P and C. Then, the real recipe results to be infeasible.

An intuitive idea to x our RTO method against measure and B ′ s uncertainty follows from the previous geometric information. The idea consists in computing for the polytope C U its convex robust region RC U such that any point in this region resists to B ′ s uncertainty and to measurement errors.

Chapter 3. Robust RTO Method

RC U U C U C U u R u R Figure 3.6: Robustope RC U .
That is, any point u ∈ RC U is guaranteed to remain inside the real polytope for all possible realizations of B and any measurement error, within some reasonable limits. In Figure 3.6, we show the nominal polytopes U and C U , the robustope RC U and the real polytope C U . The dierence between the nominal and real polytopes C U and C U results from a dierence between the nominal matrix B and the real matrix B. Any point in RC U (e.g. the nominal robust recipe u R ) transforms in a point in its neighborhood (the real robust recipe u R ) because of measurement errors. The robustope RC U is computed such that any point in its neighborhood (the dotted region) remains inside the real polytope C U . Here we are not computing the robust region of polytope U because the components constraints are considered as soft constraints. The real recipe u R is usable even if it falls slightly outside the polytope U .

In this way, replacing the polytope C U by RC U in the RTO method presented in Section 3.2, we have a new robust RTO method. To develop this idea, rst we need to model and measure the uncertainty we would like to be protected of and then to compute the appropriate robust region of C U . This development follows the robust optimization theory (see [START_REF] Ben-Tal | Robust Optimization[END_REF]).

Robust Regions and Robust RTO

Let u be a nominal recipe and y(u) its properties. To model the components and properties measurement uncertainties, we suppose that the real recipe u lies in the ball S(u, δ u ) of radius δ u and center u whereas the real blend's properties y lies in the ball

S(y(u), δ y ) = {y ∈ R m | y ∈ [y(u) -δ y , y(u) + δ y ]}
with center y(u). Finally, to model matrix B uncertainty, we use interval sets. That is, we suppose that for each i = 1, . . . , m and j = 1, . . . , n, the real value B i,j is comprised in the interval

I i,j = [B i,j -ε - i,j , B i,j + ε + i,j ]
around its nominal value B i,j for some ε - i,j , ε + i,j ∈ R + . Here we can use dierent sets and any norm to model the uncertainty. The level of conservatism (how much we want to be protected of uncertainty) and the problem complexity depend on these selections. Taking interval sets and the norm L ∞ , the robust regions obtained are polytopes and the complexity in the model is preserved at the expense of being probably too conservative (we are protected from the worst deviations of all B ′ s values and from the biggest measurement errors occurring all at the same time). Now we proceed to construct the robust region of polytope C U .

Any point u ∈ C U is robust regarding B uncertainty i

y ≤ y(u) ≤ y.
That is, i for any B such that B i,j ∈ I i,j ,

y i ≤ α v b 0 i + β v B i • u ≤ y i , ∀i = 1, . . . , m (3.8)
holds. Notice that any B i ∈ I i can be expressed parametrically as

B i (z) = B i + z -T Q - i + z +T Q + i (3.9) with Q - i = diag(ε - i,1 , . . . , ε - i,n ), Q + i = diag(ε + i,1 , . . . , ε + i,n ), z -= min{0, z} and z + = max{0, z} for some z ∈ R n such that z ∞ ≤ 1. Hence, replacing (3.9) in (3.8) we have that u ∈ C U is B-robust i y i -β v (z -T Q - i +z +T Q + i )•u ≤ y i (u) ≤ y i -β v (z -T Q - i +z +T Q + i )•u, ∀i = 1, . . . , m is valid for all z ∈ R n such that z ∞ ≤ 1. Now, as min z∈R n : z ∞≤1 {z -T Q - i + z +T Q + i } = -ε - i , ∀i = 1, . . . , m and max z∈R n : z ∞≤1 {z -T Q - i + z +T Q + i } = ε + i , ∀i = 1, . . . , m, then we have that u ∈ C U is B-robust i y i + β v ε - i • u ≤ y i (u) ≤ y i -β v ε + i • u, ∀i = 1, . . . , m. (3.10)
From now on we denote

γ - i (u) = β v ε - i • u and γ + i (u) = β v ε + i • u, ∀i = 1, . . . , m.
In addition to B uncertainty, u ∈ C U resists also to components' measure uncertainty if and only if any point in S(u, δ u ) satises (3.10). That is, i

y i + γ - i (u) ≤ y i (u + z) ≤ y i -γ + i (u), ∀i = 1, . . . ,
{B i • z} = δ u B i 2 , ∀i = 1, . . . , m
we have that (3.4) is true i

y i + γ - i (u) + δ i ≤ y i (u) ≤ y i -γ + i (u) -δ i , ∀i = 1, . . . , m (3.11)
with

δ i = β v δ u B i 2 .
Equivalently, any recipe u ∈ C U is robust regarding B uncertainty and properties' measure uncertainty i

y i + γ - i (u) ≤ y i (u) + Z i ≤ y i -γ + i (u), ∀i = 1, . . . , m (3.12)
is valid for all Z ∈ R m such that |Z i | ≤ δ yi . As previously, computing the minimum and maximum on Z, (3.12) holds i

y i + γ - i (u) + δ yi ≤ y i (u) ≤ y i -γ + i (u) -δ yi , ∀i = 1, . . . , m (3.13)
Finally, letting ∆ i = max{δ i , δ yi } we obtain the robustope RC U of polytope C U :

RC U = {u ∈ R n | y i +γ - i (u)+∆ i ≤ y i (u) ≤ y i -γ + i (u)-∆ i , ∀i = 1, . . . , m}.
Any u ∈ RC U resists to B uncertainty and to components and properties measurement uncertainty.

To summarize, let's consider a nominal feasible recipe u ∈ U ∩ C U . If we suppose that u ∈ S(u, δ u ) and y ∈ S(y(u), δ y ) for some δ u ∈ R + and δ y ∈ R m+ and if B i,j ∈ I i,j for some ε - i,j , ε + i,j ∈ R + , ∀i = 1, . . . , m, ∀j = 1, . . . , n, then the real recipe u will satisfy the real hard constraints: u ∈ C U .

The RTO method proposed in Section 3.2 transforms then into a robust RTO method by a simple substitution of C U with RC U . This robust RTO method depends completely on the robustope RC U and to obtain it we only need to determine the values of ε -, ε + , δ u and δ y . It's worthwhile to mention that while δ u and δ y are xed values independent of the RTO loop's length, the parameters ε -and ε + depend on it. As uctuations on B may accumulate over time, the longest the loop's length is, the biggest these uctuations can be.

We notice that when v and V 0 are considered as free variables in the RTO method, the value of β v is not determined. In this case, we need to use the identity β v = 1 1+u 0 and the Equations (3.11) and (3.13) to describe the polytope RC U .

In the previous discussion, we limit the analysis to measurement and components properties uncertainties. However, other types of uncertainty manifest geometrically in the same way and thus can be treated identically. For instance, when the uncertainty is due to uncontrolled factors like temperature or humidity, the real blend's properties are located in a ball around the nominal blend's properties. We can include these uncertainties by taking δ y as the maximum radius of these balls. To model the uncertainty in the components prices, we can transform the optimization problem to one with certain objective function and such that uncertainty appears as a constraint which may be considered as a row of the matrix B. We can proceed similarly to deal with the uncertainty in the components availabilities. In this case the robust region RU of polytope U needs to be determined.

S 1 U C U RC U u N u R -c Figure 3.7: Computing a robust blend in R 2 .
suppose that we want to empty the container of the previous blend (V 0 = 6 m 3 ) producing no more than v = 6 m 3 of the new blend and let's assume that the cost of the component u 1 is three times the cost of component u 2 (c = (3, 1)). The resulting blending polytopes for these volumes and the optimal nominal recipe u N = (0.25, 0.75) with cost c * N = 1.5 are shown in Figure 3.7. If we suppose that components and properties measurement errors are bounded by δ u = 0.01 and δ y = 0.01 and that B can vary in the interval [B -, B + ] = [[3.9, 7.9], [4.1, 8.1]], then the recipe u N risks to be infeasible. To generate a robust blend which resists to these uncertainties, we compute

γ - 1 (u) = γ + 1 (u) = 0.05u 1 + 0.05u 2 and ∆ 1 = δ 1 = 0.
0447 in order to obtain the robustope

RC U = {(u 1 , u 2 ) ∈ R 2 | 0.7922 ≤ u 1 + 2.02u 2 , u 1 + 1.97u 2 ≤ 1.6855}.
The optimal robust recipe is u R = (0.2933, 0.7067) and has a cost of c * R = 1.59. In this case, producing the robust recipe is 5.7% more expensive than producing the nominal one.

Example

In this section we present an example to illustrate and compare some key aspects of the RTO method and its robust counterpart. The BP (1) consists in producing a xed volume V total = 5000 m 3 of blend from 8 components and V 0 = 2000 m 3 of the heel's volume from a previous blend. Each component and the previous blend has 7 properties to be controlled during the process and they are represented by the 7 × 8 matrix B and vector b 0 respectively. Vectors y, y and u, u stand for the properties and components bounds respectively while vector c denotes the components costs. From these results we are interested in comparing the price of robustness with a component of the reblending cost (the cost dierence between the recipes obtained when the heel's volume is used and when it is not). We notice that this comparison is biased against the price of robustness as the components for the blending process are selected, obviously, according to the properties of the previous blend. Thus, trying to produce a blend from the same components without using the previous blend, may result more expensive or even infeasible. In a fair comparison, we should select the best available components to produce the blend from scratch. We do not have the elements to do this, but in order to provide a comparison as fair as possible, the price of robustness is obtained by taking V 0 = 0 (no reblending involved) and the reblending cost from the nominal recipes (no robustness involved).

B =              
              [b 0 , y, y] =           30
To compute the price of robustness we conduct a blending simulation over 36 RTO loops of 2-hours length. In Table (3.1) we show the average recipe's cost over the 36 periods and the percentage cost increase (∆ c ) from the nominal recipe to the robust recipe for dierent θ values and taking V 0 = 0.

Then we generate the optimal nominal recipes when V 0 = 168.44 for θ = 0.1) and V 0 = 0 are used in the production of 5000 m 3 of blend. The costs of these recipes are 88.29 and 87.04 respectively. We observe that for this case, the nominal recipe's cost with reblending is (88.29 -87.04)/87.04 = 1.43% greater than the nominal recipe's cost without reblending whereas the relative price of robustness is of only 0.35% for a signicant θ value of 0.1.

(V 0M ax θ c * R c * N ∆ c (%)
In Tables (3.2) to (3.5), we show the relative price of robustness ∆ c for a set of BP instances considering all of them the use of a previous blend. In the following, we comment these results.

First, we notice that the price of robustness for the BP (4) is zero for all the θ values considered. This is explained by recalling that the feasible region of the robust BP is U ∩ RC U ⊆ U ∩ C U . Thus, if the optimal robust solution to the BP is a vertex of U , then it is also the optimal solution for the nominal BP. For BP (2) and BP (4), the price of robustness is less than 1% whereas for BP (3) it reaches 4.53% for the highest level of conservatism.

On the other hand, xing V 0 = 0, the only feasible instance is BP (5). For the BP instances 2 to 4, it is impossible to produce the blend from scratch using the same components that were used to reblend. The optimal recipe's cost for BP (5) is 86.9480 which is 0.14% more expensive than the optimal nominal blend using the heel's volume.

As we have already noticed, using the same set of components to produce the blend from scratch and from a previous blend is not the appropiate way to compute the reblending cost. However, considering that the comparison scheme is biased against the price of robustness and that we are only comparing a component of the reblending cost, we think that the results obtained from these comparisons comfort the idea that reblending might be much more expensive than producing a robust blend even if the uncertainty is modeled by interval sets, which may be too conservative.

Finally, we highlight the importance of computing the reblending cost using the best available components to produce independently the blend from scratch and using heel's volume. This value can be used to adjust the level of conservatism in the model via the parameters δ u , δ y and θ. In the absence of this value, we can use the price of robustness as a benchmark for the reblending cost. The price of robustness is the highest value to pay for reblending.

After these considerations, we can conclude that for BP (1) and BP (4), the cost of producing the robust blend is lower than the reblending cost. For the BP instances 2, 3 and 5 we should produce the robust recipe unless its reblending costs are lower than 0.77%, 4.53% and 0.1% respectively.

Blending Problem Infeasibility

A BP instance (B, v, b 0 , V 0 , c, u, u, y, y) is determined by the components B, the volume v, the heel b 0 and its volume V 0 , the components cost c, the components' availability and hydraulic bounds u and u and by the properties bounds y and y. At each loop, the RTO method presented in Section 3.2 solves the BP in order to produce, if possible, the minimal cost recipe. If the BP is infeasible, some bounds are relaxed and a feasible blend with minimal quality giveaway is searched. The method considers the possible modication of the BP instance to achieve better results. Specically, the volumes V 0 and v are used as control variables in the model in order to produce a more economical recipe or a feasible one when the original values of V 0 and v do not allow it. In some cases, even with this extra degree of freedom, the BP remains infeasible for the required volumes V 0 and v. Then, in order to assure the continuity of the operation, the method proposes to produce the closest infeasible blend to the ideal target blend y T .

In this chapter, we are interested in measuring and visualizing the BP infeasibility and in using other parameters of the BP instance to control the blending process. Besides the volumes v and V 0 , the components' availability, the hydraulic capacities or even the components itself can be used to drive the Blending Process. Here, we consider only components' availability and hydraulic bounds modications.

The organization of this chapter is as follows: First, we review the main features of the infeasibility analysis of linear systems in Section 4.1. Then, in Section 4.2 we classify the BP infeasibility and we state the necessary conditions under which bounds modications can produce a feasible BP instance. In Section 4.3 we measure the BP infeasibility in terms of the blending polytopes and we highlight the existence of two vectors related to the infeasible BP. These vectors will be exploited in Sections 4.4 and 4.5 to construct 2D-visualizations of the BP infeasibility. Later, in Section 4.6 we design customized procedures to modify the components' bounds and transform an infeasible BP instance into a feasible one according with the optimization criteria used in the RTO method. We observe also that replacing the Blending Polytopes by their corresponding robustopes these procedures can be directly applied to produce a robust blend. In Section 4.7 we use the bounds' modications procedure to drive the blending process towards the blend within specications requiring to pour the minimum blends' volume. This analysis 36 Chapter 4. Blending Problem Infeasibility leads us to investigate the impact of the components bounds on the blending process. We propose then a set of indices to measure the impact of the components bounds on the blends' qualities and we accompany these measures with a supportive visualization. An example is given in Section 4.8.

Measuring and Analyzing the BP Infeasibility

In the industrial applications, when a linear programming problem results to be infeasible, the decision maker faces the problem of nding a way to transform the infeasible model into a feasible one. The infeasibility of the problem can be due to a bad formulation (e.g., a model not reecting the real problem or a model not well-posed for its numerical solution by a certain method).

Most linear solvers scale the constraints (see below) in a preprocessing stage to improve the characteristics of the model and make it best suited for its solution. However, the infeasibility can be also the result of a set of incompatible requirements in the model which, if seen individually, are completely pertinent. The detection of these incompatibilities and the study of the required modications to repair the problem instance is known as Infeasibility Analysis (see ( [START_REF] Chinneck | Feasibility and Infeasibility in Optimization:: Algorithms and Computational Methods[END_REF]) for a recent review of this subject). The potential risk of inconsistencies and the complexity to detect them increases with the problem size. The linear programming models arising in industrial applications today are increasingly complex and can have up to millions of constraints. For this reason, it is becoming more necessary to have automated methods to detect and handle the infeasibility.

A dierent reason for developing these methods, which particularly concern us, is the widespread use of linear solvers in the kernel of a RTO system.

In this situation, the optimizer should propose, in real time, a solution that guarantees the continuity of the operation even if the problem is infeasible.

A method to identify the causes of infeasibility and to propose how to repair the model will be very useful to drive the process towards a tailored feasible solution.

The most usual method to measure the infeasibility of a point x violating an individual constraint A i x ≤ b i is the dierence A i xb i . This method extends to measure the infeasibility of a point over a set of constraints in two main ways:

• The sum of the infeasibilities over all the constraints.

• The number of violated constraints.

Measuring and Analyzing the BP Infeasibility 37

In practice, a constraint is considered as violated if the violation exceeds a tolerance τ . That is, if A i xb i > τ . This test is very easy to compute but it is very sensitive to scaling factors as we can see in the following example.

Let's consider the constraint x ≤ 1 evaluated at x = 1 + τ 2 . Numerically, this constraint is considered as satised. However, the same constraint scaled by 10 (10x ≤ 10) and evaluated at the same point produces a violation of the same order than the scaling factor: 10 + 10 × τ 2 -10τ = 4τ .

To avoid this problem, one can scale the constraints to ensure all the coecients to have approximately the same magnitude before measuring the infeasibility. A dierent approach, exempt from scaling issues, consist in measuring the infeasibility as the Euclidean distance between x and the closest point in the feasible region dened by the set of constraints. In Section 4.2

we will use this idea to measure the infeasibility of the whole BP.

Among the tools developed to detect and repair the infeasibility there are three main approaches ( [START_REF] Chinneck | Feasibility and Infeasibility in Optimization:: Algorithms and Computational Methods[END_REF]):

• Identication of an Irreducible Infeasible Subset (IIS) of constraints.

• Identication of a Maximum Feasible Subset (Max FS) of constraints.

• Best Repair proposal of the model.

An IIS has the property of being infeasible but becomes feasible if any single constraint is removed. Detection of IIS allows the decision maker to concentrate in small subsets of inconsistencies which are easier to understand. Notice that an infeasible model can have several IIS's, then repairing an IIS may not repair the whole system. A brute-force method to repair the whole system consists in detecting and repairing iteratively an IIS until feasibility. The identication of a Max FS of constraints can also be viewed as the detection of the minimum number of constraints whose removal leaves a feasible system. This is the minimum unsatised linear relation problem (MIN ULR) which is known to be NP-Hard (see [12]). Both approaches, IIS and Max FS detection have a similar aim, to reduce the scope of the problem and let the decision maker to use its expertise to repair the model. This intervention turn them in semi-automated methods. An empirical analysis of three methods (including IIS) to identify and isolate the infeasibility in a Blending Problem was conducted in 1992 ( [START_REF] Greenberg | An empirical analysis of infeasibility diagnosis for instances of linear programming blending models[END_REF]). The conclusions of the study points out the superiority of the IIS method over the others, but it also highlights the diculty to release useful information from the isolation of the infeasibility.

The best repair proposal approach diers, from the later two, in this sense.

It focuses on determining algorithmically the best modication to the problem' constraints with respect to a predened criterion. Among the methods of this type there are those that only modify the right hand side vector (RHS Methods) and those looking the feasibility of the problem by modifying all the coecients of the constraints. The last problem, much more dicult than the former one is usually addressed by minimizing the norm of a correction matrix H (see [START_REF] Amaral | A reformulation-linearizationconvexication algorithm for optimal correction of an inconsistent system of linear constraints[END_REF]). Somme common choices for the matrix norm are L 1 and L ∞ and the ∞-norm :

H ∞ = max i j |h ij |. In ([4]
) the Frobenius norm is used. The optimal matrix correction (according to some criterion) gives a relatively similar feasible system of constraints, although it is not always clear how to use this information. Here, we will focus on the RHS Methods.

RHS Methods

In [START_REF] Murty | Infeasibility Analysis for Linear Systems, a survey[END_REF], the authors distinguish four dierent criteria to transform the infeasible system Ax ≤ b onto a feasible one using only RHS's modications:

The Smallest Changes Model. This model seeks to modify the minimum number of RHS terms (equivalent to the MIN ULR problem).

The Smallest Penalty Model. If the penalty cost f i of changing the RHS b i is xed for all i (does not depend on the modication's magnitude), this method seeks to change the set of b i 's making the problem feasible while minimizing the associated sum of penalties.

The Smallest Variable Cost Model. If the penalty cost depends on the amount of each change, this model seeks to minimize the total variable cost of all the changes to the b i 's which makes the problem feasible. This can be written as the following linear program:

min m i=1 c + i e + i + c - i e - i + m+p i=m+1 c - i e - i s.t. A i x + e - i -e + i = b i , i = 1, . . . , m A i x + e - i ≥ b i , i = m + 1, . . . , m + p e - i , e + i ≥ 0, i = 1, . . . , m + p where c + i (c - i
) denote the non-negative cost per unit increase (decrease) in the value of b i . The variables e + i and e - i are called elastic variables as they let the constraints to stretch to make the feasible region nonempty. Notice that for the b i 's associated to the inequality constraints, we only need the elastic variables increasing the feasible region. So, we let e + i = 0 for i = m + 1, . . . , m + p. If (x * , e - * i , e + * i ) is an optimal solution of this problen, then b ′ i = b ie - i + e + i for all i is the optimal modication of b under this model and x * is a feasible solution of the modied problem.

Infeasibility Levels 39

The Smallest Variable Cost Model with Bounds. This is the same as the previous model, except that bounds are imposed on the changes. Then, letting

u + i ≥ 0 (u - i ≥ 0)
be the maximum possible increase (decrease) on b i for all i, we have:

min m i=1 c + i e + i + c - i e - i + m+p i=m+1 c - i e - i s.t. A i x + e - i -e + i = b i , i = 1, . . . , m A i x + e - i ≥ b i , i = m + 1, . . . , m + p 0 ≤ e - i ≤ u - i , 0 ≤ e + i ≤ u + i , i = 1, . . . , m + p.
If some b i cannot be increased, decreased or modied at all, we x u

+ i = 0, u - i = 0 or u + i = u - i = 0.
Other dierent measures of the required adjustment to the RHS can be considered, e.g., min (b-Ax)

+ 1 or min m i=1 [(b i -A i x) + ] 2 with (•) + = max{0, (•)}.
All these methods have dierent desirable characteristics that distinguish them (robustness, ease of calculation, etc.). However, the best way to repair an infeasible problem depends on the particular situation.

Infeasibility Levels

Let's dene the polytope

C U 1 = C U ∩ S 1
in components space and the polytope

C Y 1 = {y(u) ∈ R m | u ∈ C U 1 } in properties space.
If the BP results to be infeasible for a xed pair of volumes (V 0 , v) and assuming that U = ∅, then one may have 3 dierent levels of infeasibility which can be stated in terms of the blending polytopes: Level 1.

C U 1 = ∅ and U ∩ C U 1 = ∅. Level 2. C U = ∅ and C U 1 = ∅. Level 3. C U = ∅.
For the three levels of infeasibility, we can modify the components properties (the matrix B) and force the polytope C U to have a non empty intersection with the polytope U producing in this way a feasible blend. On the other hand, as the simplex S 1 contains all the possible recipes, then C U 1 represents the set of all the recipes having their properties (C Y 1 ) within specications. A polytope U can reach, by bounds modications, a point u ∈ C U \ U i u ∈ S 1 . That is, i the BP has infeasibility of level 1.

In Figure 4.1 we show the three levels of infeasibility for a BP instance having two components and two properties. When the infeasibility is of level 1, we observe that modifying u 1 and u 2 the polytope U expands and we can reach any recipe u ∈ C U 1 . In contrast, having infeasibility of levels 2 and 3 we need to change the components (B) in order to get a non empty C U touching the polytope U . In this case, it is the polytope C U which is modied while U stays unchanged.

1 1 S 1 Level 3 u 1 u 1 u 2 u 2 U S 1 Level 2 U C U S 1 Level 1 U C U C U 1
For each infeasibility level, the infeasibility of a recipe u ∈ U can be measured in both spaces as the Euclidean distance between u and the closest point in the corresponding polytope. In the components space, for levels 1 and 2, as the distance between u and the closest point in C U 1 and C U respectively. In the properties space, for the three levels of infeasibility, as the distance between y(u) and the closest point in C Y 1 , C Y and C.

In a natural way, we can dene then the BP infeasibility in terms of the blending polytopes: the BP infeasibility of level 1 is dened by

D U = D(U, C U 1 ) and D Y = D(P, C Y 1 )
is a feasible direction in components space which expresses the minimal eort moving direction of U in order to meet C U 1 .

We observe that p * = y(u * ) ∈ P and q * = y(t * ) ∈ C Y 1 . Thus,

w U Y = q * -p * q * -p * 2
is a feasible direction in properties space.

To compute the distance between polytopes P and C Y 1 , we should notice that for these polytopes we don't have neither the primal nor the dual representations. However, we can compute their distance indirectly from polytopes U and C U 1 and the basic equation. Solving the problem min y(t)y(u) 1

s.t. u ∈ U, t ∈ C U 1
we obtain p * = y(u * ) ∈ P and q * = y(t

* ) ∈ C Y 1 such that D Y = d(p * , q * ). Then w Y = q * -p * D Y (4.3)
is the shortest moving direction of P in order to meet C Y 1 and

w Y U = t * -u * t * -u * 2
is the corresponding feasible direction in components space.

The distances D U 2 and D Y 2 can be obtained in the same way: The polytope C U 2 is given in its dual form and C Y 2 is the image of C U 2 under the basic equation. In contrast, the polytope C has no preimage when the infeasibility is of level 3. Hence, D Y 3 should be computed directly in the properties space. The vertices of C are easy to obtain and in Section 4.5, we will show how to get those ones of P . Then, computing the separating hyperplane of largest gap between the polytopes P and C in their primal form (see Section 2.2), we have D Y 3 .

In this section we highlight the existence of the BP inherent feasible directions:

(i) w U and w Y U in components space and (ii) w Y and w U Y in properties space.

Next we are going to exploit these directions to obtain 2D-visualizations of the BP infeasability in each space.

Infeasibility Visualization in Components Space

Let w 1 = w U and w 2 = w Y U . These vectors dene a 2D-plane in the components space whose parametric equation is

H 2 = {u(s) | u(s) = u M + w T • s} with s = (s 1 , s 2 ) T ∈ R 2 , w = (w 1 , w 2 ) T and u M the midpoint of the segment [u * , t * ] (see Equation (4.2)).
Let's denote by U 2 and C 2 the cuts of U and C U 1 with H 2 . The polygone U 2 is then described by the set of points u(s) ∈ H 2 verifying:

1 T • u(s) = 1 (4.4) u ≤ u(s) ≤ u.
(4.5) Equation (4.4) is satised by all u(s) ∈ H 2 since u M ∈ S 1 and 1 T • w = (0, 0). Thus the vertex set of U 2 can be calculated explicitly because it is the intersection of the set of lines in H 2 described by (4.5). In the same manner, we can nd the vertices of the polygone

C 2 = {u(s) ∈ H 2 | y ≤ y(u(s)) ≤ y, u(s) ∈ S 1 }
as the intersection of a set of lines in H 2 . We know that

1 T • u(s) = 1, ∀u(s) ∈ H 2 ,
therefore C 2 is completely characterized by the set of points u(s) ∈ H 2 such that:

0 ≤ u(s) ≤ 1 (4.6) y ≤ α v b 0 + β v Bu(s) ≤ y. (4.7)
Finally, replacing u(s) in (4.6) and (4.7) we get the set of lines in H 2 whose intersection determines the vertex set of C 2 :

-Bu M ≤ w T • s ≤ 1 -Bu M y -α v b 0 β v -Bu M ≤ (Bw) T • s ≤ y -α v b 0 β v -Bu M .
To obtain the infeasibility visualization in components space (see Figure 4.2), we simply draw these polygones in the plane H 2 with u M taken as the center of coordinates and vectors w 1 and w 2 as the axis. Notice that these axis are not orthogonal in general.

It is important to remark the characteristics of this visualization. As we obtain it by cutting polytopes U and C U 1 with a plane dened by w U and w Y U which passes through the nearest points in both polytopes, the distance between the polygones U 2 and C 2 is D U . The distance between the polytopes U and C U 1 in the original n-space. By the same reason, the angle between directions w U and w Y U is preserved.

U 2 C 2 D U w 1 w 2

Infeasibility Visualization in Properties Space

Following the same lines of the preceding section, let's consider the feasible directions in the properties space w 1 = w Y and w 2 = w U Y and let's denote by H 2 the 2D-plane dened by

H 2 = {y(s) ∈ R 2 | y(s) = y M + w T • s} where s = (s 1 , s 2 ) T ∈ R 2 , w = (w 1 , w 2
) and y M is the midpoint of the segment [p * , q * ] (see Equation (4.3)). We can reuse the same development of the previous section to obtain the polygons P 2 and C 2 by cutting P and C Y 1 with H 2 on condition to have the dual description of P and C Y 1 . We don't have neither the primal nor the dual description of these polytopes but they can be obtained (at least theoretically) as follows.

Let V U and V P be the vertex sets of polytopes U and P . First, we observe that V P is contained in the image of V U under the basic equation as stated in the next Property 4.5.1. If V P and V U are the vertex sets of polytopes P and U respectively then V P ⊆ {y(u) | u ∈ V U }.

Proof. Let y(u) ∈ V P and suppose u / ∈ V U . Then u is a convex combination of V U with at least two coecients strictly positive:

u = u i ∈V U λ i u i with u i ∈V U λ i = 1, λ i ≥ 0 and λ 1 , λ 2 > 0.
Applying the basic equation to u we have

y(u) = α v b 0 + β v Bu = u i ∈V U λ i (α v b 0 + β v Bu i ) with u i ∈V U λ i = 1, λ i ≥ 0 and λ 1 , λ 2 > 0.
So, y(u) is a convex combination of vectors in P which contradicts our assumption.

Then, we obtain the primal and dual representations of P from the dual of U in three steps:

(I) Compute V U = {u 1 , u 2 , . . . , u N } and then V = {v i | v i = y(u i ) , i = 1, . . . , N }.
(II) Delete from V iteratively each v i which can be expressed as a convex combination of V \ v i to obtain V P .

(III) Solve the facet enumeration problem for the set V P .

We can proceed in the same way to construct the primal and dual forms of

C Y 1 .
In Section 2.1 we state that vertex and facet enumeration are dicult problems for general polytopes. If we want to reuse the last section procedure to produce the infeasibility visualization in properties space, we need to solve both problems. In practice, the vertex enumeration for polytopes U and C U 1 Chapter 4. Blending Problem Infeasibility can be done eciently. We can partially explain this situation by noticing that polytope U has a particularly simple structure that can be exploited by a customized algorithm. It is the intersection of a parallelepiped with the simplex S 1 and all its vertices are either vertices of the parallepiped or they lie in the interior of its edges (see Appendix A).

On the contrary, the facets of P and C Y 1 can not always be obtained eciently. Occasionally, we encounter polytopes P and C Y 1 for which the facet enumeration is very time consuming (several minutes). For the practical application this computing times are unacceptable, thus we need to modify the procedure to generate the desired visualization in the properties space avoiding the dual construction. Specically, in the following algorithm we produce the polygons P 2 and C 2 by cutting P and C Y 1 with H 2 but having P and C Y 1 in their primal form. 4.3a). Let P 2 = {p 2 (1), . . . , p 2 (k)} be the polygon obtained with at least one vertex (p * is the optimal solution of Problem (4.8) for d T = [1, 0]). If P 2 has more than one vertex, then we set i = 1 and we try to expand P 2 iteratively in the next step. Else, we go to End.

y M + w T • s = p i ∈V P λ i p i p i ∈V P λ i = 1 λ i ≥ 0, i = 1, . . . , |V P | (4.8) for d T = [1, 0], d T = [0, -1], d T = [-1, 0] and d T = [0, 1] (see Figure

Expand

We look for a new vertex in the outer normal direction d out to the edge (p 2 (i), p 2 (j)), where j = (i + 1) mod (|V P 2 |). This is done by solving Problem (4.8) with d = d out . Let q * be the optimal solution for this problem. If d T • p 2 (i) < d T • q * then q * is a new vertex in P 2 (see Figure 4.3b). We insert q * in the vertex list at position i + 1. We proceed similarly to obtain the polygon C 2 and nally produce the infeasible visualization in the properties space (see Figure 4.4). This visualization also preserves the same important characteristics as in components space: The distance D Y between polygones P 2 and C 2 which is the same as the distance between P and C Y 1 and the angle between directions w Y and w U Y . We nish this section remarking that we can use the algorithm (4.1) to produce the visualizations for the infeasibility of levels 2 and 3. For the level 2, D U 2 and D Y 2 determine 4 feasible directions (2 in each space) that dene the 2-D plans with which to cut the pair of polytopes (U, C U ) in the components space and (P, C Y ) in the properties space. For the infeasibility of level 3, we have only the feasible direction generated by the closest points in P and C (the polytope C U 1 is empty). Then, a second direction is choosen conveniently in order to get an infeasibility visualization in the properties space.

Otherwise, we set

i = i + 1. If i < |V P 2 |,

Becoming Feasible and Robust

In the previous sections we showed how to measure and visualize the BP infeasibility. Let's now consider how to modify the hydraulic and availability bounds in order to improve the blend's quality and to produce, if possible, a feasible recipe. Specically, we are going to transform the polytope U into the polytope:

U N EW = {u ∈ R n | u -∆ ≤ u ≤ u + ∆ , 1 T • u = 1, 0 ≤ u ≤ 1}
by choosing the appropriate elastic variables ∆, ∆ ∈ R n+ and looking to have

U N EW ∩ C U 1 = ∅.
We notice that hydraulic and availability bounds can be adjusted individually and can reach their lower and upper limits: [0, 1]. In other words, the polytope U N EW can attain any recipe in C U 1 as C U 1 ⊆ S 1 . On the other hand, these bounds are settled in accordance with the global renery planning and by operational requirements (e.g., to promote the use of plentiful or inexpensive components, or in order to empty one component's tank). The impact of the bounds violation on the overall renery's performance is dicult to measure. Even harder is to estimate the individual impact of each bound violation. For this reason, we decide to proceed as follows: First we compute the best recipe u * ∈ C U 1 to produce according to some optimization criterion and then we obtain the bounds' modications to have u * ∈ U N EW (see Figure 4.5). We observe that choosing the vectors u -, u + ∈ R n+ , we can constrain the bounds modications by searching the best recipe u * in U + ∩C U (see Equation 4.9) instead of C U 1 .

U + = {u ∈ R n | 0 ≤ u -u -≤ u ≤ u + u + ≤ 1 , 1 T • u = 1} (4.9)
In particular, xing u - i = 0 or u + i = 0 we forbid the modication of the lower or upper bound for the i th component. Aditionally, in Section 4.7 we construct some benchmarks of the impact of the components bounds on the blends' cost, on the quality giveaway and on the blending volume. These benchmarks might be useful to decide whether or not to apply the bounds' modications.

Some particular optimization criteria to search the best recipe are the minimization of the infeasibilities of level 1 in both spaces and the minimization of the bounds' modications. This last one is the smallest variable cost model (see Section 4.1) with all the costs having the value 1.

Another natural choice is to select u * ∈ C U 1 using the same optimization criteria applied in the RTO method: we start searching a feasible recipe of minimal cost by solving the problem

min c T • u (4.10) s.t. u ∈ C U 1 .
If the problem 4.10 is infeasible (C U 1 = ∅), then we search the blend with minimal quality giveaway but considering exclusively the hard constraints:

min y(u) -y T 1 (4.11) s.t. u ∈ C U 1H .
Finally, if the problem 4.11 is infeasible (C U 1H = ∅), then we can not transform U into U N EW to reach a feasible recipe but solving min y(u)y T 1 (4.12) s.t.

u ∈ S 1 we obtain the recipe with properties as close as possible to y T .

Once the recipe u * is determined, we get the optimal bounds' modications ∆ * = max{0, u -u * } and ∆ * = max{0, u * -u} allowing us to produce u * . We denote by

E = n i=1 (∆ * i + ∆ * i ) (4.13)
the required eort to perform the BP instance transformation. The approach described in this section does not consider the BP uncertainty which may cause the feasible blend obtained for the nominal data to be in reality infeasible. This situation is reected by the fact that polytopes U N EW and C U 1 are just touching (see Figure 4.5). So, even a little dierence between the nominal and the real data may result in a false intersection. In order to overcome this fragility and construct a robust blend, we only need to substitute in the previous procedure the polytopes by their corresponding robustopes as determined in Section 3.4. In this way, the recipe u * will be selected deeply inside the polytope C U 1 and will resist to any variation in the data, within the pre-established limits. In Figure 4.6 we illustrate the result obtained from this procedure. In this case, the components' bounds are also considered as hard constraints which must be satised. So, after selecting the robust recipe u * ∈ RC U 1 we transform the robustope RU onto RU N EW in order to have u

* ∈ RU N EW ∩ RC U 1 . U RU C U 1 RC U 1 RU N EW u *

Overseeing the Blending Process

The bound's modication procedure of the previous section was developed to x the BP infeasibility in accordance to the optimization criteria of the RTO method. Here, we present other meaningful objective toward which to drive the Blending Process.

When blending, one may be interested in getting from the current infeasible blend (b 0 ) a feasible one while pouring the least possible blend's volume. This criterion turns out to be very practical in case of unforeseen events, like a component shortage or a pump halt, requiring to interrupt prematurely the blending process. Moreover, it allows to sell the blend as soon as possible thus reducing the inventory costs. Next we describe the procedure to compute the rst reachable blend (in terms of volume) within specications. First, solving From the basic equation (Equation 2.2), we have that the initial blend b 0 dilutes in the new blend as the poured volume increases. Therefore, the possible modications to b 0 provided by the new blend reduce with time. This traduces in a reduction of the movement speed of polytopes P and P 1 as they move away from b 0 . For the same reason, the only variation in this gure when we use a dierent initial volume V 0 , is the movement speed of P and P 1 : reducing V 0 speeds up the movement of both polytopes as we need to pour a smaller amount of blend's volume in order to modify b 0 .

b 0 y M ax y M ax1 y M in y M in1 P (V M in ) P (V M in1 ) P 1 (V M in1 ) C Y 1 CP ∩ C C CP 1 CP
Finally, let's suppose that c * and q * (c * 1 and q * 1 ) are the minimal cost and minimal quality giveaway that can be achieved by any blend in CP ∩ C (in C Y 1 ). Hence, as

P (v) ⊆ P 1 (v) for all v ≥ 0, it follows that P ∩ C ⊆ C Y 1 . Consequently, V M in1 ≤ V M in , c *
1 ≤ c * and q * 1 ≤ q * . Thus, we can say that the components bounds have an impact of (c *c * 1 ) in the blend's cost, of (q *q * 1 ) in the blend's quality giveaway and of (V M in -V M in1 ) in the required volume to pour in order to have the rst feasible blend.

It is worth to mention that we can generate this supportive visualization at any time in the blending process. All that we need to do is to compute the vertex sets V C, V P (∞) and V P 1 (∞). Then, adding the point b 0 to V P and to V P 1 , we have the primal description of CP and CP 1 . Next taking polytope C. Let's suppose that b 0 is the actual blend and y * ∈ P ∩ C is the optimal blend to produce. Then, cutting CP , CP 1 and C with the plane dened by vectors w1 = y *b 0 and w2 = e i (the i th canonical vector in R m ) and passing through b 0 , we obtain the desired cut (see Figure 4.9). In this visualization, we distinguish some new information elements: The dotted line between b 0 and y * represents the nominal path of the i th blend's property i.e., the path to be followed by the i th blend's property in absence of uncertainty. The cone CP shows the limitations imposed by the components bounds on the range of attainable values for the property i and the cone CP 1 shows all the possible values of property i if the components bounds are disregarded. The relative location of the target blend with respect to the blending polytopes (here, we observe that y T is out of reach even with bounds modications) and nally, the location of y * in C and CP 1 give us some information about the non-robustness of y * and the capacity to correct it. Here, the blend y * is located on the boundary of C and is prone to become infeasible, however, with bounds modications a more robust blend can be produced.

Now we explain how to produce some benchmarks of the impact of the produce a volume V total = 5000 m 3 of blend and using V 0 = 2000 m 3 of the previous blend. The infeasible blend y 1 producing the minimal quality giveaway was then computed. The cost, the quality giveaway and the required eort (see Equation 4.13) to reach y 1 are grouped in its quality vector :

Q(y 1 ) = [90.7162 422.9910 0] T .

Then we compute the interval [u 0M in , u 0M ax ] for which a feasible blend exist. From these values and the relation u 0 = V 0 /v, we can compute all the couples (V 0 , v) producing a blend within specications. In particular, a minimum volume v = 34578 m 3 needs to be produced if we want to use completely the heel's volume. All this is reected in Figure 4.10, where the blend y Min ∈ P (34578) ∩ C is the rst feasible blend on the path followed by the polytope P .

From the material presented in this chapter, we have more information about the BP (1) infeasibility. The problem has infeasibility of level 3 as the polytope C U (3000) is empty. This means that for v = 3000 m 3 and V 0 = 2000 m 3 , bounds modications do not suce to attain a feasible blend. For these volumes, the best we can do is to extend the polytope U to reach the recipe u 2 ∈ S 1 with minimal quality giveaway but remaining infeasible. To do this, we solve the Problem (4.12) getting the blend y 2 with quality vector Otherwise, solving the Problem (4.14) we get (u 0M ax1 , u 0M in1 ) and from them, taking V 0 = 2000 m 3 we compute the range of volumes

[V M in1 , V M ax1 ] = [5100, ∞)
from which bounds modications, including the null modication, can produce a feasible blend. In Figure 4.10, the blend y Min1 ∈ P 1(5100) ∩ C is the rst reachable blend by bounds modications. For any volume less than 5100, the cone CP 1 does not intersect C.

For v = 5100, the problem continues to be infeasible but of level 1. The polytopes U and C U 1 are at distance D U = 0.5404 whereas the distance between P and C Y 1 is D Y = 234. Taking w 1 = w U and w 2 = w Y U we cut the polytopes U and C U 1 to get the infeasibility visualization in components space (see Figure 4.11). The distance between the polytopes U and C U 1 living in R 8 is exactly the gap depicted between the cuts of U and C. The last one consists only of the point t * .

In properties space, the visualization (see Figure 4.12) is generated using the vectors w 1 = w Y and w 2 = w U Y . Here again, the cut of C Y 1 consists only of the point q * . The gap separating the doted lines is the distance between Having infeasibility of level 1, we can compute the minimal cost recipe u 3 = 0 0 0 0 0 0.3389 0 0.6611 T in C U 1 (5100). This recipe can be attained by the polytope U applying the following bounds' modications: and the visualization of this transformation is presented in Figure 4.13. Here we are cutting the polytopes U , U N EW and C U 1 using w 1 = w U and w 2 = u 3 -u M . In this way, we get in the same cut the cheapest recipe u 3 and the closest one t * . This explains why the cuts of U and C U 1 are dierent from the ones of the Figure 4.11. Here, the cut of U consists only of the point u * whereas the cut of C is the segment [t * , u 3 ].

(∆, ∆) =             0.
Now we proceed to measure the impact of the components' bounds on the blending process. In Figure 4.10, we can see at a glance how the components bounds impose limitations not only on the size of the feasible region CP ∩ C but also to the required time (that is the volume to pour) to reach a feasible blend. Numerically, the impact of the components bounds on the optimal blend's cost, on the quality giveaway and on the required volume to pour before getting the rst feasible blend are:

K C = 9.71% K Q = 234.38% K V = 578.10%
These indices tell us that components' bounds are responsible of an increase of almost 10% in the blend's cost, of 234% in the blend's quality giveaway and of 578% in the blending volume required to produce the rst blend within specications. Now, we measure the impact of each bound on the blend's qualities. That is, we compare the optimal blend's qualities for the BP (1) when all the components bounds are disregarded versus the BP (1) when only one bound at a time is active. These indices are showed (in percentage) in Table 4.1. There we can read, for instance, that the upper bound 8 has no incidence on the blend's cost, and that it increases the quality giveaway in 107% and the minimal blending time in 26%. This does not mean that allowing the modication of u 8 will reduce the quality giveaway and the blending time by these quantities. What this really means is that u 8 prevents to reduce by more than half the quality giveaway and by 26% the blending time, even if all the other bounds are free to be modied.

Regarding the uncertainty in the process, let's suppose that we want to be protected from measurement and B's uncertainty within the following limits: δ y = 0.1200 0.0003 4.5000 0.0826 0.0280 0.000049 0.2000 δ u = 0.01 and θ = 0.01 as we did at the beginning of the example in Section 3.5. Then we compute the robust volumes intervals:

T Component LB UB K C K Q K T K C K Q K T
[V M inR , V M axR ] = [55304, ∞)
and

[V M in1R , V M ax1R ] = [5784, ∞).
This means that we need to pour a surplus of 784 m 3 of blend in order to be in range of getting by bounds modications the rst robust blend. Let's produce 5784 m 3 of blend and let 93.8729 44.3127 1.017 From these vectors we observe that the cost of the robust blend y 6 is 5.01 % higher than the cost of the nominal blend y 7 . In contrast, we realize that the cost of the nominal blend y 3 (for v = 5100) is 10.79 % higher than the cost of y 7 . Thus, producing the blend's surplus of 784 m 3 we can reach a robust and cheaper blend than y 3 .

In RC U 1 taking the vectors w 1 = w U and w 2 = u 6u M . The cut of RC U 1 being the line segment [t * , u 6 ]. Notice that the robust recipe u 6 is located on the boundary of the polytope C U 1 . Contrary to the intuition, this can happen because the boundary where u 6 is located is a boundary of S 1 : recall that C U 1 = C U ∩ S 1 and observe that u 6 is the strict convex combination of two vertices of S 1 . Thus, not being u 6 on the boundary of C U , it can be robust to properties measurement and to B's uncertainty. Moreover, as it is impossible to produce a recipe outside S 1 , then u 6 is robust also to components measurement uncertainty. Now, in Figures 4.15 T which is the rst feasible blend that can be produced. We notice that at the beginning of the process b 0 has 5 of 7 properties within specications and that the property 5 will be the last property to reach its specications. In this example, the cuts of CP are very thin and are confused with the nominal path of y * . We observe also that the target blend is out of reach even with bounds modications (y T / ∈ CP 1 ). Finally, regarding the robustness of y * , it seems that the properties at most risk are the properties 5 and 6. Chapter 5

Conclusions

In this work, we introduced the Blending Polytopes inherent to the Blending Problem. Based on these polytopes we proposed a real time optimization (RTO) method for the linear oil blending process. This method is intended to avoid reblending while minimizing the blend's cost and the quality giveaway.

A main characteristic of our RTO method is the integration of the case when the blend has to be produced using the heel of a previous blend. This feature provides meaningful information for the control system, for instance, to determine the appropriate heel's volume to use in the blend or to have an estimate of the required volume to pour before getting the rst blend within specications.

Several sources of uncertainty pervade the blending process. We identied and modeled some dierent types of them following the RO techniques. Using interval sets, the robust regions are polytopes and a straightforward robust version of the RTO method is obtained.

Reblending a previous blend which failed to be within specications has several costs associated: The tank rent, the logistic and inventory costs . . . Another one is the recipe's cost increase for using the heel's volume of the previous blend. Here we compared this component of the reblending cost with the price of robustness to stress the convenience of producing robust recipes and prevent reblending. The example presented and the limited set of other simulations that we conducted are not conclusive and further studies are needed before considering the implementation of our method in the renery blending circuit. However, we are only comparing a component of the reblending cost (the recipe's cost increase) and the results obtained comfort the idea that reblending might be much more expensive than producing a robust blend. Then, using a more conservative technique like RO may improve the global performance of the blending process. The recipe's cost rises by taking the robust recipe in place of the nominal one but the reblending expenses cancel.

In the second part of this work, we exploited the Blending Polytopes to develop a ne analysis of the BP infeasibility. We use them to classify, measure and visualize the BP infeasibility and to generate a component's bounds modications procedure to drive the blending process towards the best robust blend. A set of indices measuring the impact of the component's bounds For the general case, there are two main classes of algorithms for solving VE: the pivoting methods and the insertion methods. Apparently, the rst published method to solve VE is the insertion method [START_REF] Motzkin | The double description method[END_REF], now known as the double description method. Te DD method involves the insertion of constraints iterativelly. At each step, a vertex list is updated: the old vertices which become infeasibles in the new system are deleted whereas the vertices generated by the new constraint are inserted. To detect the new vertices at each step, the DD method and its variants must determine whether two former vertices x i and x j are adjacent. The complexity of this test depends strongly on the number of present vertices at each step and thus on the insertion order of the constraints. In [START_REF] Dyer | The Complexity of Vertex Enumeration Methods[END_REF], an example is presented to show that the worst-case complexity of the insertion methods is not bounded by a polynomial function on m, n and v. However, assuming that the number of present vertices at intermediate steps remains in O(v), the whole complexity of the DD method is in O(m 2 v 3 ) (see [START_REF] Fukuda | Double description method revisited[END_REF] for ecient implementations).

The pivoting methods, inspired by the simplex method, exploit the relation between the polytope's vertices and the Basic Feasible Solutions (see, e.g., [START_REF]Linear Programming (Series of Books in the Mathematical Sciences)[END_REF]). For simple polytopes this relation is one to one whereas for the degenerate case, a degenerate vertex is related to several BFS. These methods start producing an arbitrary vertex and continue iteratively generating all the neighbor BFS produced by a pivoting step. Each neighbor BFS determines a vertex which might have been discovered at a previous step. Thus, the method must include a test to prevent duplicated vertices to be output. In the worst case a polytope may have a unique vertex and up to m d BFS. Thus, the pivoting methods are not ecient for degenerated polytopes.

For some special classes of polytopes, ecient algorithms have been proposed. For example, if P is a polytope simple, [START_REF] Dyer | The Complexity of Vertex Enumeration Methods[END_REF] presented the rst linear- time algorithm (on v) to solve the VE problem. The algorithm performs a depth rst search on the graph of P and uses a balanced tree (see, e.g., [3]) to store the polytope vertices. Later, the reverse search method was presented in [START_REF] Avis | A pivoting algorithm for convex hulls and vertex enumeration of arrangements and polyhedra[END_REF]. This method enumerates the vertices of P in time O(ndv) and has the advantage of using a memory space which is bounded by a polynomial in the input size. For polytopes arising from networks, [START_REF] Provan | Ecient enumeration of the vertices of polyhedra associated with network LP's[END_REF] gives an algorithm in O(v 2 ) and for polythedra associated with dual LI(2) systems (systems having at most two nonzero coecients per constraint), [START_REF] Abdullahi | Listing Vertices of Simple Polyhedra Associated with Dual LI(2) Systems[END_REF] presents a pivoting algoritm running in O(nv). For polytopes having all its vertices in {0, 1} d , a backtracking approach is presented in [11] which runs in polynomial time on the input size and in linear time on the number of vertices.

Available codes implementing some of these methods can be found in [START_REF] Fukuda | cdd/cdd+ Reference Manual[END_REF]5,[START_REF] Christof | PORTA: Polyhedron representation transformation algorithm[END_REF]30].

A.2 ZerOne Algorithm

Let P be a polytope described in its dual form and let's denote by P the class of polytopes which are combinatorially equivalent to a d-cube of appropriate dimension. That is, it exists a bijection between the faces of the d-cube and those ones of P ∈ P that preserves the inclusion relation. Finally, let Π the class of polytopes which can be obtained taking the convex hull of a set of vertices from a polytope P ∈ P . In Figure A.1 we depict the construction in dimension 3 of a polytope P in this class.

Figure 1 :

 1 Figure 1: Système de mélange.

et le polytope cible C

 cible = {y ∈ R m | y ≤ y ≤ y}. Dans ce dernier polytope on distingue le mélange cible y t qui est un mélange fourni par le département de planication ayant la sur-qualité minimale et d'autres caractéristiques souhaitables. Ce mélange sera considéré comme le mélange idéal à produire quand on ne parvient pas à générer une recette viii réalisable pour le PM. À ce point nous avons identié les polytopes sous-jacents au PM et nous observons que, pour un volume xe v, une solution réalisable du PM est une recette u ∈ C U ∩ U avec des propriétés y(u) ∈ P ∩ C Y . De façon équivalente, PM est infaisable ⇐⇒ P ∩ C Y = ∅ ⇐⇒ U ∩ C U = ∅.(0.5) Une Méthode RTO Robuste Le système de mélange est généralement constitué de trois sous-systèmes fonctionnels : le sous-système de planication, l'optimiseur en ligne et le soussystème de contrôle. L'objectif du sous-système de planication est de planier les opérations de ranage pour diérents horizons temporels. Dans le long terme, les objectifs de production sont déterminés en utilisant des prévisions de la demande des produits et des prix du pétrole. Dans la planication à moyen terme, ces objectifs sont révisés et ajustés par rapport aux estimations de disponibilité des composants, des réservoirs,. . . Enn, dans le court terme la planication se concentre sur la production de recettes cibles pour toute opération de mélange prévue. Le volume de production, les composants à mélanger et tous les autres paramètres nécessaires au calcul de la recette cible sont déterminés à ce niveau et selon le planning général dans les raneries. Dans ces conditions, on devrait s'attendre à avoir toujours une solution réalisable pour l'instance du PM. Toutefois, le calcul de la recette cible est la plupart du temps basé sur des prévisions et, occasionnellement, sur des mesures des propriétés des composants. Dans les deux cas, la précision du modèle est aectée : dans le premier cas, le modèle hérite des imprécisions des prédictions et dans le second cas les mesures des composants peuvent également être inexactes et retardées [13]. En fait, comme les composants sont en production en même temps qu'ils sont mélangés, l'imprécision dans le modèle induit par la connaissance défectueuse des propriétés des composants reste tout au long du procédé. Pour cela et d'autres types d'incertitudes, la recette cible peut devenir sous-optimale, voire irréalisable. An de réduire la disparité entre le modèle et le système réel, l'optimiseur en ligne met à jour la recette cible, dans un cycle en boucle fermée et au moyen de mesures en ligne des propriétés des composants et du mélange actuel. La nouvelle recette servira au contrôleur à guider le procédé de mélange. Avant de présenter le modèle RTO, on discute de la production d'un mélange à partir d'un mélange fabriqué préalablement. Si le mélange est produit à partir de zéro ou si le volume à couler est trop grand par rapport au volume du fond du bac, c'est-à-dire, si V 0 = 0 ou v → ∞, alors l'équation de base se réduit à y(u) = Bu et le polytope P (v) ne dépend pas de V 0 . D'autre ix part, si V 0 = 0 alors P (v) est un polytope homothétique à P (∞) avec b 0 comme centre d'homothétie. Dans l'espace des propriétés, le polytope C ne dépend pas de v. Il existe alors un intervalle (éventuellement vide) de volumes

Figure 2 :

 2 Figure 2: Evolution du polytope P (v) avec le volume coulé.

xii 3 .

 3 L'incertitude du modèle.4. L'incertitude des marchés : les disponibilités et les prix des composants, la demande des diérents mélanges, etc

Figure 3 :

 3 Figure 3: Robustope RC U .

2 Figure 4 :

 24 Figure 4: Visualisation de l'infaisabilité dans l'espace des composants.

Figure 6 :

 6 Figure 6: Impact de la modication des bornes des composants sur le procédé de mélange.

3 . 3

 33 Blending from a previous blend in R 2 . . . . . . . . . . . . . . . 18 3.4 Evolution of polytope P with the blending volume. . . . . . . 19 3.5 RTO method. . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.6 Robustope RC U . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.7 Computing a robust blend in R 2 . . . . . . . . . . . . . . . . . 28 4.1 Infeasibility levels. . . . . . . . . . . . . . . . . . . . . . . . . 40 4.2 Infeasibility visualization in components space. . . . . . . . . . 44 4.3 Algorithm 2D-Cut from primal form. . . . . . . . . . . . . . . 47 4.4 Infeasibility visualization in properties space. . . . . . . . . . . 48 4.5 Transforming U onto U N EW to reach the feasible blend u * . . . 49 4.6 Transforming RU onto RU N EW to reach the robust blend u * . . 51 4.7 Impact of the components' bounds on the blending process. . . 53 4.8 Monitoring the blend's properties. . . . . . . . . . . . . . . . . 54 4.9 Monitoring the blend's properties (new information elements). 55 4.10 BP (1): Impact of the components' bounds. . . . . . . . . . . 57 4.11 BP (1): Infeasibility visualization in components space. . . . . 59 4.12 BP (1): Infeasibility visualization in properties space. . . . . . 60 4.13 BP (1): Reach the cheapest recipe. . . . . . . . . . . . . . . . 61 4.14 BP (1): Reaching a robust recipe. . . . . . . . . . . . . . . . . 63 4.15 BP (1): Monitoring the property 1. . . . . . . . . . . . . . . . 64 4.16 BP (1): Monitoring the property 2. . . . . . . . . . . . . . . . 64 4.17 BP (1): Monitoring the property 3. . . . . . . . . . . . . . . . 65 4.18 BP (1): Monitoring the property 4. . . . . . . . . . . . . . . . 65 4.19 BP (1): Monitoring the property 5. . . . . . . . . . . . . . . . 66 4.20 BP (1): Monitoring the property 6. . . . . . . . . . . . . . . . 66 4.21 BP (1): Monitoring the property 7. . . . . . . . . . . . . . . . 67 A.1 Obtaining a polytope in class Π from the 3cube . . . . . . . 73 A.2 Recipe's polytope and its boundary polytope P . . . . . . . . 75 A.3 Search tree generated by algorithm 1. . . . . . . . . . . . . . . 75 A.4 Recipe's polytope U without a known boundary polytope. . . 76 xxv Chapter 1

Figure 1 . 1 :

 11 Figure 1.1: Blending system.

Figure 2 . 1 :

 21 Figure 2.1: The 2-simplex and 3-cube in R 3 .

Chapter 2 .

 2 Blending Polytopes

Figure 2 . 2 :

 22 Figure 2.2: Maximum gap dened by the hyperplane H*.

Chapter 2 .

 2 Blending Polytopes physical and/or chemical characteristic of the component j.

Figure 2 . 3 :

 23 Figure 2.3: Blending polytopes in R 2 .

Figure 3 . 1 :

 31 Figure 3.1: Polytope P (v) with b 0 / ∈ P (∞).

Chapter 3 .

 3 Robust RTO Method

Example 3 . 1

 31 Blending from a previous blend

Figure 3 . 3 :

 33 Figure 3.3: Blending from a previous blend in R 2 .

Figure 3 . 4 :

 34 Figure 3.4: Evolution of polytope P with the blending volume.

Algorithm 3 . 1

 31 Blending RTO Method I Incorporate the new available information (V 0 , b 0 , B and v) in the basic equation. If V 0 = 0, solve the problems (3.3), (3.5) and (3.6) in this order until having a feasible solution. Otherwise, go to step II.

Figure 3 .

 3 Figure 3.5: RTO method.

Example 3 . 2

 32 Computing a robust blendIn the example 3.1 we determine the interval [u 0M in , u 0M ax ] = [1 6 , 1.35] from which we can compute all the couples (V 0 , v) producing a feasible blend. Let's Chapter 3. Robust RTO Method

Tu

  = 0.1428 0.0819 0.0352 0.0881 0.0338 0.0100 0.0352 0.0100 T u = 0.4225 0.1930 0.0704 0.1690 0.2000 0.2000 0.3592 0.2000 T In order to produce a robust recipe, we assume that components and properties measurement errors are bounded by δ u = 0.01 and δ y = 0.1200 0.0003 4.5000 0.0826 0.0280 0.000049 0.2000 T respectively. Regarding B uncertainty, we dispose of B -and B + the absolute lower and upper limits of matrix B. So, we dene τ -= B -B -, τ + = B + -B and T = B + -B -. We suppose that B ′ s real values are comprised in the intervals [Bε -, B + ε + ] with ε -= min{τ -, θT } and ε + = min{τ + , θT } for some θ ∈ [0, 1]. As we stated previously, the values of δ u and δ y are xed costs have decreased. The dierence c * N -c * N 0 is a component of the reblending cost.

Figure 4 . 1 :

 41 Figure 4.1: Infeasibility levels.

Figure 4 . 2 :

 42 Figure 4.2: Infeasibility visualization in components space.

  Algorithm 4.1 2D-Cut Primal Form Start Compute the extreme points of P ∩ H 2 in the directions w 1 and w 2 by solving the problem: max d T • s s.t.

Figure 4 . 3 :

 43 Figure 4.3: Algorithm 2D-Cut from primal form.

2 Figure 4 . 4 :

 244 Figure 4.4: Infeasibility visualization in properties space.

Figure 4 .

 4 Figure 4.5: Transforming U onto U N EW to reach the feasible blend u * .

Figure 4 . 6 :

 46 Figure 4.6: Transforming RU onto RU N EW to reach the robust blend u * .

Figure 4 . 7 :

 47 Figure 4.7: Impact of the components' bounds on the blending process.

Figure 4 . 9 :

 49 Figure 4.9: Monitoring the blend's properties with new information elements.

Q(y 2

 2 ) = [111.0754 156.1711 1.0614] T .

Figure 4 .

 4 Figure 4.12: BP (1): Infeasibility visualization in properties space.

Figure 4 . 13 :

 413 Figure 4.13: Bounds modications of polytope U to attain the cheapest recipe u 3 ∈ C U 1 .

u 6 =T

 6 be the minimal cost of robust and nominal recipes that can be reached by bounds modications. The quality vectors of these blends are:Blend Cost Quality Giveaway Eort

Figure 4 .

 4 [START_REF] Chinneck | Feasibility and Infeasibility in Optimization:: Algorithms and Computational Methods[END_REF] we show the infeasibility visualization in components space for v = 5784. This time we are cutting the polytopes U , U N EW , C U 1 and 4

Figure 4 . 14 :

 414 Figure 4.14: Transforming polytope U to reach the robust recipe u 6 .

to 4 .

 4 21 we show the visualizations to monitor the blend's properties of BP (1). Here, we are tracking the blend y * = 33.2649 0.2305 579.2266 24.3628 432.0940 1.4379 41.8209

Figure 4 . 15 :

 415 Figure 4.15: BP (1): Monitoring the property 1.

Figure 4 . 16 :

 416 Figure 4.16: BP (1): Monitoring the property 2.

Figure 4 . 17 :

 417 Figure 4.17: BP (1): Monitoring the property 3.

Figure 4 . 18 :

 418 Figure 4.18: BP (1): Monitoring the property 4.

Figure 4 . 19 :

 419 Figure 4.19: BP (1): Monitoring the property 5.

Figure 4 . 20 :

 420 Figure 4.20: BP (1): Monitoring the property 6.

Figure 4 . 21 :

 421 Figure 4.21: BP (1): Monitoring the property 7.

Figure A. 1 :

 1 Figure A.1: Obtaining a polytope in class Π from the 3cube .

  w 1 , w 2 ) et y M le milieu du segment [p * , q * ].

	xviii
	résoudre). Cependant, la structure des polytopes U et C U 1 et leurs tailles
	(m, n ≤ 20), permettent l'obtention de leurs sommes d'une manière ecace.
	Au contraire, les facettes de P et C Y 1 ne peuvent pas toujours être obtenues
	de manière ecace. Parfois, nous rencontrons des polytopes P et C Y 1 pour
	lesquels l'énumération de facettes est très coûteuse en temps de calcul. Ainsi,
	on modie la procédure pour obtenir les polygones P 2 et C 2 en coupant les
	polytopes P et C Y 1 à partir de leur forme primale.
	Devenir faisable et robuste
	Dans la dernière section on a montré comment mesurer et visualiser l'infaisabi-
	lité du PM. On va voir maintenant les modications des bornes hydrauliques
	et de disponibilité an d'améliorer la qualité du mélange et de produire, si
	possible, une recette faisable. Plus précisément, on va transformer le polytope
	U dans le polytope
	Pour produire la visualisation de l'infaisabilité dans l'espace des propriétés, on
	peut réutiliser la procédure suivie dans l'espace des composants, à condition
	d'avoir la description duale des polytopes P et C Y 1 . Nous n'avons, ni la des-
	cription primale ni la description duale de ces polytopes, mais on peut les avoir
	(au moins théoriquement) à partir des polytopes U et C U 1 et de l'équation
	de base. Le problème d'énumération des sommets d'un polytope quelconque
	est un problème ouvert (on ne connaît pas d'algorithme polynomial pour le

  u}

	xx
	U
	U N EW
				C U 1
				u *
	qui nous permettent de produire u * . On désigne par
	E =	n	(∆ * i + ∆	* i )
		i=1
	l'eort requis pour eectuer la transformation de l'instance du PM.
	L'approche décrite précédemment ne considère pas l'incertitude du PM qui Figure 5: Transformation du polytope U pour atteindre le mélange faisable peut provoquer que le mélange faisable obtenu pour les données nominales soit en réalité infaisable. Cette situation s'explique par le fait que les polytopes u * ∈ C U 1 .
	U N EW et C U 1 se touchent à peine (voir Figure 5). Donc, toute petite dié-
	rence entre la valeur nominale et les données réelles peut entraîner une fausse
	intersection. Pour remédier à cette fragilité et construire un mélange robuste, comme suit. Tout d'abord, on résout les problèmes
	il sut de substituer dans la procédure précédente, chaque polytopes par son robustope correspondant tels qu'ils ont été déterminés dans la Section (3.4). max (min) u 0 (0.16)
	De cette façon, la recette u * sera choisi bien à l'interieur du polytope C U 1 et s.t.
	résistera à toute variation des données, dans les limites pré-établies.
	Guidage du procédé de mélange
	Lors du procédé de mélange, on peut vouloir obtenir un mélange faisable à par-
	tir d'un mélange infaisable (b 0 ) tout en versant le moins de volume possible.
	Ce critère se révèle être très pratique en cas d'événements imprévus, comme
	le manque d'un ou plusieurs composants ou l'arrêt d'une pompe, obligeant
	à interrompre prématurément le procédé de mélange. En outre, il permet de
	vendre le mélange le plus tôt possible en réduisant ainsi les coûts d'inventaire.
	Pour trouver le premier mélange faisable (en termes de volume), on procède

  Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Distance Between Polytopes . . . . . . . . . . . . . . . . . . . 2.3 Blending Polytopes . . . . . . . . . . . . . . . . . . . . . . . . The 2-simplex and 3-cube in R 3 . . . . . . . . . . . . . . . . . . 7 2.2 Maximum gap dened by the hyperplane H*. . . . . . . . . . 9 2.3 Blending polytopes in R 2 . . . . . . . . . . . . . . . . . . . . . 12 3.1 Polytope P (v) with b 0 / ∈ P (∞). . . . . . . . . . . . . . . . . . 16 3.2 Polytope P (v) with b 0 ∈ P (∞). . . . . . . . . . . . . . . . . . 17

		List of Figures
	)	iii
	1 Introduction
	2 Blending Polytopes
	2.1 5 Conclusions
	A Vertex Enumeration
	B VISAMEL
	Bibliography
	xxiii

3 Robust RTO Method 3.1 Blending from a previous blend . . . . . . . . . . . . . . . . . 3.2 RTO Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 BP Uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Robust Regions and Robust RTO . . . . . . . . . . . . . . . . 3.5 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Blending Problem Infeasibility

4.1 Measuring and Analyzing the BP Infeasibility . . . . . . . . . 4.2 Infeasibility Levels . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Gap Computation and Feasible Directions . . . . . . . . . . . 4.4 Infeasibility Visualization in Components Space . . . . . . . . 4.5 Infeasibility Visualization in Properties Space . . . . . . . . . 4.6 Becoming Feasible and Robust . . . . . . . . . . . . . . . . . . 4.7 Overseeing the Blending Process . . . . . . . . . . . . . . . . . 4.8 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 ZerOne Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Vertex Enumeration of Blending Polytopes . . . . . . . . . . . 1.1 Blending system. . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.1

Table 3 .

 3

	0.01 87.3851 87.1097	0.32
	0.02 87.3542 87.0763	0.32
	0.03 87.3847 87.1036	0.32
	0.04 87.3632 87.0795	0.32
	0.05 87.4058 87.1194	0.33
	0.06 87.4134 87.1249	0.33
	0.07 87.4232 87.1309	0.33
	0.08 87.4258 87.1306	0.34
	0.09 87.4297 87.1307	0.34
	0.10 87.4077 87.1062	0.35

1: BP (1): Relative price of robustness.

Table 3 . 2 :

 32 BP (2): Relative price of robustness.

	θ	c * R	c * N	∆ c (%)
	0.01 86.0184 85.8682	0.17
	0.02 86.0567 85.8682	0.22
	0.03 86.1158 85.8682	0.29
	0.04 86.1749 85.8682	0.36
	0.05 86.2342 85.8682	0.43
	0.06 86.2935 85.8682	0.50
	0.07 86.3529 85.8682	0.56
	0.08 86.4124 85.8682	0.63
	0.09 86.4720 85.8682	0.70
	0.10 86.5317 85.8682	0.77
	θ	c * R	c * N	∆ c (%)
	0.01 84.8437 82.2032	3.21
	0.02 84.9591 82.2032	3.35
	0.03 85.0749 82.2032	3.49
	0.04 85.1911 82.2032	3.63
	0.05 85.3079 82.2032	3.78
	0.06 85.4251 82.2032	3.92
	0.07 85.5427 82.2032	4.06
	0.08 85.6609 82.2032	4.21
	0.09 85.7916 82.2032	4.37
	0.1 85.9239 82.2032	4.53

Table 3 . 3

 33

	θ	c * R	c * N	∆ c (%)
	0.01 87.2632 87.2632	0
	0.02 87.2632 87.2632	0
	0.03 87.2632 87.2632	0
	0.04 87.2632 87.2632	0
	0.05 87.2632 87.2632	0
	0.06 87.2632 87.2632	0
	0.07 87.2632 87.2632	0
	0.08 87.2632 87.2632	0
	0.09 87.2632 87.2632	0
	0.1 87.2632 87.2632	0

: BP (3): Relative price of robustness.

Table 3 .

 3

			Chapter 3. Robust RTO Method
	θ	c * R	c * N	∆ c (%)
	0.01 86.8841 86.8267	0.07
	0.02 86.8872 86.8267	0.07
	0.03 86.8903 86.8267	0.07
	0.04 86.8933 86.8267	0.08
	0.05 86.8964 86.8267	0.08
	0.06 86.8995 86.8267	0.08
	0.07 86.9025 86.8267	0.09
	0.08 86.9067 86.8267	0.09
	0.09 86.9113 86.8267	0.1
	0.1 86.9159 86.8267	0.1

4: BP (4): Relative price of robustness.

Table 3 .

 3

5: BP (5): Relative price of robustness. Chapter 4

Table 4 . 1 :

 41 BP (1): Individual bounds impact on the blend's qualities.

Remerciements

Figure B.4: Analysis of bounds modications. Bibliography

Chapter 3. Robust RTO Method during the process whereas θ is directly related to the RTO loops' length. We take θ = 0.01 for the 2-hours length RTO loops.

Solving Problem (3.2) we get

The fact of having u 0M inR = 0 has a special meaning that we are going to exploit in this example. As u 0 = V 0 V , then taking V 0 = 0 (producing the blend from scratch) there exists a robust recipe for any volume v > 0. For the moment, let's consider V 0 = 2000 m 3 . This gives the robust feasible volumes interval:

The corresponding interval for the nominal case is:

This means that we need to produce at least 34578 m 3 (55304 m 3 ) in order to get a (robust) blend within specications which uses completely V 0 . For the problem with only the hard constraints, we obtain similar intervals.

If we decide to produce V total = 5000 m 3 using V 0 = 2000 m 3 , we solve Problem (3.6) and we obtain a recipe with cost 90.72 per m 3 but producing a blend out of specications. Instead of this, we can compute the biggest heel's volume allowing us to produce 5000 m 3 of robust blend: From relations u 0M axR = V 0 V and V total = V 0 + v we obtain

Next, xing V 0 = 174.5 m 3 and v = 4825.5 m 3 we solve Problem (3.3) to obtain the optimal robust recipe: Therefore, producing a robust recipe induces a cost increase of 2.68% over the nominal recipe's cost. This is the price of robustness. However, taking V 0 = 0 (no reblending), the robust recipe cost is c * R0 = 87.70 and the nominal recipe cost is c * N 0 = 87.04. Hence, the price of robustness is only 0.76%. Moreover, we observe that both, the nominal and the robust the infeasibility of level 2 by

and the infeasibility of level 3 by D Y 3 = D(P, C).

These distances give us relevant information for the analysis of infeasibility, although of dierent kind: In the components space, D U and D U 2 are estimates of the required eort to obtain a feasible BP instance by hydraulic and availability bounds modications and by components modications respectively. These distances are only an estimate because, as we will see later, the required eort to transform the instance depends on the type of transformation applied and on the way to measure the eort. In the properties space, D Y , D Y 2 and D Y 3 measure the blend's quality losse that can be saved by components bounds modications and/or by components modications.

Gap Computation and Feasible Directions

Let's dene a feasible direction in component's space as a translation direction of polytopes U or C U producing its collision and therefore the BP feasibility. A translation direction of polytopes P or C Y producing their collision is called also a feasible direction (in properties space) as it can be seen as the image under the basic equation of a feasible direction in components space.

As we have stated early, here we are considering only bounds modications to transform an infeasible BP instance into a feasible one. In the following, we will only develop the analysis for the infeasibility of level 1. However, we notice that we can measure and construct the visualizations (when they exist) for the infeasibility of levels 2 and 3. We will highlight the dierences in the procedure explained below.

In components space we are interested in computing the gap between polytopes U and C U 1 . Both polytopes are described in their dual form. So, as discussed in Section 2.2, we can compute this distance by solving:

we obtain u 0M ax and u M in1 (u 0M in and u M ax1 ). Then, from these values and the relation u 0 = V 0 /v we get the volumes interval

for which C U 1 = ∅. The blends y Min1 = y(u Min1 ) and y Max1 = y(u Max1 ) are the rst and last blends, in terms of volume, that can be produced by means of "some" bounds' modications. Finally,

are the minimal bounds modications to reach the blend y Min1 at volume

Regarding robustness, we only need to replace the polytope C U 1 by its robustope RC U 1 in Problem (4.14) in order to recover the rst and last pro- ducible robust blends when no bounds are considered.

The procedure to obtain [V M in1 , V M ax1 ] is the same we used in the RTO method to compute the feasible volumes interval [V M in , V M ax ] but taking U ∩S 1 instead of U ∩ C U . The relation between these volumes and the impact of the components bounds on the blending process can be better explained by referring to Figure 4.7. In this gure, the polytope P runs from P (0) = b 0 to P (∞) forming the cone

The points y Min ∈ P (V M in ) ∩ C and y Max ∈ P (V M ax ) ∩ C are the rst and last feasible blends that can be produced.

On the other hand, as U ⊆ S 1 then P (v) ⊆ P 1 (v) = y(S 1 ), ∀v ≥ 0. Therefore, P 1 (v) describes the cone

which consists of all the producible blends when the hydraulic and availability bounds are disregarded. The rst and last of these blends being within specications (in polytope C) are y Min1 and y Max1 . The intersection of CP 1 with C is the polytope C Y 1 and the visualization in properties space presented in the directions w 1 = y Minb 0 and w 2 = y Min1b 0 and by means of the algorithm (4.1), we obtain the cuts of C, CP and CP 1 with the 2D-plane H 2 dened by

where s = (s 1 , s 2 ) T ∈ R 2 and w = (w 1 , w 2 ).

In the blending practice, it is common to monitor graphically the properties of the blend during its production (see Figure 4.8). These graphs oer at a glance useful information for the control department. For instance, one can In this paragraph, we propose an improvement for these graphs inspired by Figure 4.7. The enhancement consists in adding at the end of the observed data, a suitable cut of the blending cones CP and CP 1 and of the target Chapter 4. Blending Problem Infeasibility components bounds on the blends' cost, on the quality giveaway and on the blending volume. These benchmarks might be useful to decide whether or not to apply the bounds' modications.

First we construct a cost-related benchmark K C . To do this, let's consider the Problem (3.3) with the additional constraints (3.7):

Let u * be the optimal solution and c * its cost. This is the cheapest recipe we can produce while satisfying the components bounds. Then, replacing the constraints 4.16 and 4.17 by

we ignore the components' bounds and we get the optimal recipe u *

measures the relative cost increase produced by the components bounds.

Similarly, by means of a dichotomic search over v and solving the problem

xed at each step, we get the blend in P ∩ C with minimal quality giveaway q * . Then, after replacing U by S 1 , and

] we obtain the blend in C Y 1 with minimal quality giveaway q * 1 . In this case, the ratio

measures the relative quality giveaway increase produced by the components bounds.

Finally, the quotient

)/V M in1 measures the impact of hydraulic and availability bounds on the blending volume required to produce the rst feasible blend.

We highlight the possibility to compute these benchmarks in order to know 4.8. Example 57 the impact of a single bound or a set of bounds on the three blend's qualities: cost, giveaway and volume. To do this, we simply x the corresponding variables u - i and/or u + i to zero for all the components bounds under analysis and we let the other bounds to be variable.

Example

Let's consider the BP instance BP (1) = (B, v, b 0 , V 0 , u, u, y, y, c) presented in the example of Section 3.5.

The next discussion can be followed in Figure 4.10 which was produced by cutting the blending cones and the target polytope with the plane determined by the vectors y Minb 0 and y Min1b 0 passing through b 0 , (see 4.7). In this visualization, we observe how polytopes P and P 1 evolve in time as the blending volume increases. From b 0 = P (0) = P 1 (0) until P (∞) and P 1 (∞) they form the cones CP and CP 1 respectively.

As we have seen in Section 3.5, the BP (1) is infeasible if we want to Properties Space on the recipe's cost, on the blend's quality giveaway and on the blending time required to produce the rst feasible recipe have been proposed. These indices might be useful to decide whether or not to apply the bounds' modications.

We accompany these measures with supportive visualizations for the decision maker.

Appendix A

Vertex Enumeration

In this section we are interested in computing the vertex set for the blending polytopes U and C U 1 . For completeness, in Section A. Every polytope has then two representations. The primal representation as the convex hull of the nite set of its vertices:

and the dual representation as the intersection of a nite set of half-spaces:

The transformation A.1 to A.2 is known as the vertex enumeration (VE) problem and from A.2 to A.1 is known as the facet enumeration (FE) problem.

The dimension of a polytope is the dimension of its ane hull and a dpolytope is a polytope of dimension d in some R e (e ≥ d).

A d-polytope is called simple if each vertex is contained in exactly d facets 71 the complexity of a feasibility check.

We don't know if the blending polytopes belong to the class Π (see below).

However, we continue to present this algorithm as it is the basis for the vertex enumeration algorithm that we introduce in Section A.3.

Similar to the dual description of the d-cube, we can describe any polytope P ∈ P as

where A L , A U ∈ R n×n , and each vertex of P is uniquely determined by n inequalities satised with equality.

Geometrically, we have n pairs of inequalities that represent opposite facets of P and every vertex v of P is determined by a partition L ∪ U of the indices {1, . . . , n} for which

Given a polytope P ∈ Π, we add the redundant description of its boundary polytope P . The algorithm 1 generates recursively all the partitions L ∪ U of the rows' indices using the constraints

each vertex of P whereas its feasibility is tested using Ax ≤ b.

For a set of rows' indices L , U , let's dene the polytope

Then, the algorithm can be described as follows:

Algorithm 1 Vertex enumeration for polytopes in class Π end if 10: end if 11: end {enum} The complexity of the algorithm 1 (O(nvT LP )) follows from the fact that each vertex v ∈ V is generated only once and requires at most n feasibility tests. In order to exemplify the algorithm, let's consider the recipe's polytope U belongs to the class Π as it can be obtained from the polytope

in class P (showed in dotted lines). Then, using both representations U and ∅ u 3 (0.2,0,0.8)

(1,0,0) V U = {(0.2, 0, 0.8), (0.2, 0.8, 0), (1, 0, 0)}.

For this example, we were able to obtain the boundary polytope U ∈ P of the recipe's polytope U by shrinking all its redundant bounds. Adding U to the descrption of U , all the vertices in V U become degenerate. Each vertex is determined by three bounds and all of them satisfy the equality u 1 + u 2 + u 3 = 1. Thus, we can ignore the equality constraint for the location of the vertices and use it only for the feasibility test.

This procedure does not work in general and we don't know if the blending polytopes belong or not to the class Π. In the Figure A.4 we show the recipe's polytope

In this case, the constraints 0 ≤ u 1 ≤ 0.8, 0 ≤ u 2 ≤ 1, 0 ≤ u 3 ≤ 0.8 form a polytope in P but U 2 can not be obtained as a convex combination of a subset of its vertices. The polytope U 2 has four non-degenerate vertices, (0.2, 0, 0.8), (0, 0.2, 0.8), (0.8, 0, 0.2) and (0.8, 0.2, 0) which require the use of the equality constraint for their location.

In conclusion, as we don't have a method to produce U (if this exists), then we can not apply the zerOne algorithm to the recipe's polytope. In the next section we present a very fast VE algorithm for the recipe's polytope which explores the rows partitions and exploits the fact that the feasibility test involves only one addition. Before presenting the algorithm, we discuss some issues regarding the number of degenerate vertices of the recipe's polytope.

So, let's consider the parallepiped L n,k dened by u = s/k and u = (1s)/(nk) for some 0.5 < s < 1 and

The number of vertices in L n,k determined by k upper bounds and nk lower bounds is n k and by construction all of them are in S 1 . Then, as n k reaches its maximum when k = ⌊n/2⌋, the recipe's polytope U can have up to n ⌊n/2⌋ degenerate vertices. These pathological polytopes have n ⌊n/2⌋ vertices all of them degenerate and each vertex is determined by exactly n n-1 = n partitions L ∪ U of n -1 rows' indices.

A.3 Vertex Enumeration of Blending Polytopes

In this section we introduce a customized VE algorithm for the recipe's polytope U . The algorithm is inspired from the zerOne method and exploits the structure of the recipe's polytope. As each vertex u ∈ U is determined by n bounds or by n -1 bounds and the percentage constraint (1 T • u = 1), then u is either a vertex of the parallepiped L or it lies in the interior of one of its edges. Thus, we need to search the vertices of U in the edges of L in an ecient manner. The second property that we use is the simplicity of the feasibility test: a partition p ∈ L ∪ U is feasible i the coordinates of the vertex v(p) determined by p sum one.

Let's denote by s = u -u the slack between the components bounds and by e i the i-th canonical vector in R n . For convenience, we denote a partition p ∈ L ∪ U as a natural number whose binary representation has the i-th bit set to 1 i i ∈ U . The set of 2 n partitions encode the vertices of L and the function enumU presented in the Algorithm 2 conducts a lexicographic depth rst search on them in order to detect those which are also vertices of U .

The function enumU is called with three parameters. The partition p to be analyzed, the most signicant active bit of p (which corresponds with the biggest index of the upper bounds used in p) and the sum of the bounds determining the partition p: sum p = n i=1 p i u i + (1 -p i )u i . Initially we call enumU with p = 0, last = 0 and sum 0 = n i=1 u i . That is, the method starts Appendix A. Vertex Enumeration the dfs analyzing the vertex of L determined by all the lower bounds. We suppose that sum 0 < 1 otherwise L has no intersection with S 1 and U = ∅ or u is the only vertex of U . In general, each partition p being analyzed satises sum p < 1. The analysis of p consists in scanning all the neighbor partitions q such that p < Lex q. If sum q < 1 then we proceed to analyze q. If sum q = 1 then q determines a (degenerate) vertex of U and if sum q > 1 then the edge (p, q) of L contains the vertex v(p) + (1sum p ) e k of U , where k is the only dierent coordinate between p and q.

Algorithm 2 Vertex enumeration for the recipe's polytopes 1: function enumU (p, last, sum p ) 2: for i = last + 1 to n do 3:

enumU (q, i, sum q ) 7: end if 14: end for 15: for i = 1 to last -1 do 16:

if (p(i) == 0) and (sum q > 1) then To show the validity of the method, we notice that any vertex of U is either the extreme point q or an interior point of an edge (p, q) of L where sum p < 1 ≤ sum q . The algorithm analyzes all the partitions p such that sum p < 1 and scans exactly once all the edges (p, q) with sum q ≥ 1. Therefore, any vertex of U is enumerated by the algorithm. Finally, a partition q such that sum q = 1 is detected only once (line 11 of the algorithm) because of the use of the lexicographic order to scan the partitions. In this way we avoid the duplication of any degenerate vertex.

A.3. Vertex Enumeration of Blending Polytopes 79

In the worst case, the complexity of the algorithm is exponential in the input size. For example, the algorithm requires to analyze all the edges of L (n2 n-1 ) in order to nd the only vertex of the polytope

In spite of this, the algorithm has a very good behavior in comparison with two of the best general methods. For n = 20, which is the maximum number of components used in practice, the algorithm nds the 184, 756 vertices of the highly degenerate polytope L 20,10 described at the end of the previous section in 58 seconds. In contrast, the reverse search method ([5]) nds these vertices in more than 29 minutes and the double description method ( [START_REF] Fukuda | cdd/cdd+ Reference Manual[END_REF]) takes more than 19 minutes to enumerate the 48, 620 vertices of L 18,9 .

All the previous discussion and the algorithm 2 concern only the recipe's polytope U . However, we need to solve also the VE for the polytope C U 1 .

Unfortunately, the structure of this polytope depends completely on the components matrix B from which we don't have any information.

In order to enumerate the vertices of C U 1 , we use a backtracking algorithm similar to the algorithm 1. This algorithm generates all the partitions of the indices of n -1 rows and it tests its feasibility solving the simplex phase 1 in order to establish if the partition denes a vertex or not. If feasible, the test provides a vertex which is determined by the partition and by the equality constraint n i u i = 1. However, by using this approach, we need to deal with the redundant partitions producing the same degenerate vertex. Theoretically, the algorithm can behave badly, nevertheless it performs well in practice. (It solves the VE for any polytope C U 1 with n ≤ 20 in less than 1 minute).

VISAMEL

In this section we present some screenshots of the decision support system VISAMEL (Blends Visualization) developed for the TOTAL group. Some of the ideas presented in this work have been already implemented in this application while others are being evaluated for their future inclusion.

Robustness and visualization of blends production

Abstract: The oil blending process (BP) consists in determining the optimal proportions to blend from a set of available components such that the nal product fullls a set of specications on their properties. Two important characteristics of the blending problem are the hard bounds on the blend's properties and the uncertainty pervading the process. In this work, a real-time optimization method is proposed for producing robust blends while minimizing the blend quality giveaway and the recipe's cost. The method is based on the Robust Optimization techniques and under the assumption that the components properties blend linearly. The blending intrinsic polytopes are exploited in order to measure, visualize and characterize the infeasibility of the BP. A ne analysis of the components bounds modications is conducted to guide the process towards the best robust blend. A set of indices and visualizations provide a helpful support for the decision maker.