Modélisation en programmation par contraintes
Modélisation en programmation par contraintes
Année 2010
Modélisation et résolution en
programmation par contraintes de
problèmes mixtes continu/discret
de satisfaction de contraintes et
d’optimisation
THÈSE
pour obtenir le grade de
DOCTEUR DE L’UNIVERSITÉ DE NANTES
Discipline : INFORMATIQUE
Nicolas B ERGER
le Jeudi 7 Octobre 2010
au LINA
Nicolas B
Université de Nantes
Nicolas B
Modélisation et résolution en programmation par contraintes de problèmes mixtes
continu/discret de satisfaction de contraintes et d’optimisation
xv+146 p.
Ce document a été préparé avec LATEX2ε et la classe these-IRIN version 0.92 de l’association
de jeunes chercheurs en informatique LOGIN, Université de Nantes. La classe these-IRIN est
disponible à l’adresse :
http ://login.irin.sciences.univ-nantes.fr/
Révision pour la classe : $ Id : these-IRIN.cls,v 1.3 2000/11/19 18 :30 :42 fred Exp
Tous sous le ciel, connaissant le beau comme le beau : voici le laid !
Tous connaissant le bien comme le bien : voici le mal !
C’est ainsi que l’être et le non-être naissent l’un de l’autre,
Que le difficile et le facile s’accomplissent l’un par l’autre,
Que mutuellement le long et le court se délimitent,
Le haut et le bas se règlent,
Le ton et le son s’accordent,
L’avant et l’après s’enchaînent.
7
un grand nombre de sujets, donnant toujours plus de sens aux moments que nous passons (et passerons
encore !) ensemble. Dans tous les cas, un grand merci à tous et à toutes.
Comment ne pas penser maintenant à mes parents Brigitte et Jean-Pierre, mon frère Antoine, ainsi
que ma soeur Marie, son bien-aimé Jean-François et leur adorable petit Ruben ? Mais aussi tous les
membres de l’accueillante et très nombreuse famille à laquelle j’ai la chance d’appartenir. Sans vous, je
n’y serais sûrement jamais arrivé aussi facilement, même si j’ai eu fort à faire ces (bientôt) dix dernières
années. Je sais très bien que tout le monde n’a pas cette chance en ce monde et qu’en plus de toutes les
épreuves que j’ai eu à traverser, d’autres moins chanceux que moi étaient seuls et ont dû en plus lutter
pour se vêtir, se nourrir et avoir un toit. Rien de cela pour moi qui ai atterri en thèse sans jamais avoir eu à
me préoccuper de ce genre de problèmes qui auraient facilement pu empoisonner davantage un quotidien
déjà difficilement respirable... Mais c’est également le moment d’adresser un énorme merci à ce cher
dragon d’Erwan, que j’ai rencontré chemin faisant et qui, bien que nous ne soyons pas d’un même sang,
m’a lui aussi apporté le soutien et l’affection indéfectibles qui m’ont aidé à garder la tête hors de l’eau
quand je manquais de couler sous les assauts des vagues, toujours plus profondes et puissantes, de cette
quête de sens, vitale, qui m’anime au quotidien et, pourquoi ne pas l’espérer, pour encore très longtemps.
Pour finir, un énorme merci à L-T, ainsi qu’à ses nombreux traducteurs. Pour sa rigueur et sa
discipline, mais aussi pour l’admiration et la profondeur que m’inspirent de longue date ses écrits, qui
m’ont plus que jamais été d’une aide précieuse pour garder le cap tout au long de la rédaction de ce
manuscrit, ultime étape de mon doctorat. Oui, merci beaucoup, Vieux-Maître1 ! !
1
Traduction française communément admise pour Lao-Tseu.
S OMMAIRE
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
2 État de l’art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3 Modélisation et résolution d’un problème de conception en robotique . . . . . . . . . . . . . . . . . . . . . . 31
4 Collaboration de deux solveurs discret et continu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5 Utilisation de la contrainte AllDifferent dans un solveur continu . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
6 Spécialisation aux domaines d’entiers du filtrage basé sur l’arithmétique des intervalles . . . 87
7 Conclusion et perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
vii
I
I NTRODUCTION
Les contraintes apparaissent spontanément dans nombre de situations auxquelles sont confrontés les
êtres humains. Elles formalisent d’une manière transparente et naturelle les dépendances qui existent
dans le monde physique ou ses représentations mathématiques abstraites [Bar99]. Ce sont des restric-
tions de l’ensemble des possibles, des éléments de connaissance par lequel le champ des possibles se
réduit : les contraintes sont un moyen très générique de représenter les régularités qui gouvernent notre
monde, que ce soit au niveau informationnel, physique, biologique ou social [Dec03]. Elles matéria-
lisent une information qui existe indépendamment d’elles, et fixent les règles de ce qui est possible ou
non. Les contraintes peuvent ainsi porter sur des objets, des personnes, des grandeurs physiques ou des
abstractions mathématiques :
• la somme des angles d’un triangle vaut 180 degrés ;
• une personne ne peut pas se trouver dans plusieurs lieux simultanément ;
• Suzanne ne peut pas être mariée à la fois avec Jacques et avec Bernard ;
• chaque ligne d’une grille de sudoku contient exactement une fois chacun des chiffres de 1 à 9 ;
• etc.
Les contraintes ont quelques intéressantes propriétés qui leur sont caractéristiques [Bar]. Elles sont
non-directionnelles, c’est-à-dire que les variables sur lesquelles une contrainte porte sont tantôt l’origine,
tantôt la destination de l’information que celle-ci véhicule. Une autre de leur propriété est qu’elles sont
déclaratives : elles spécifient une relation qui doit être vérifiée, sans pour autant définir la manière
d’atteindre ce but. Elles sont également additives, i.e. l’ordre de définition des contraintes ne compte pas,
seule la conjonction finale des contraintes posées est importante.
Étant donné un ensemble de contraintes, portant chacune sur un sous-ensemble donné d’un certain
nombre de variables auxquelles on cherche à affecter une valeur, une question centrale est de savoir s’il
existe une possibilité de satisfaire simultanément toutes ces contraintes et, le cas échéant, quelles valeurs
on peut affecter à chaque variable pour être sûr de satisfaire toutes les contraintes. L’exemple 1.1 montre
un problème bien connu1 qui se base sur la notion de contraintes.
Exemple 1.1. Cinq maisons sont alignées, de couleurs différentes. Dans chaque maison vit un homme
de nationalité différente. Chaque personne boit une boisson différente, chaque personne fume un type
de cigarettes différent, chaque personne élève un animal différent. Pour chaque maison, on cherche à
connaître sa couleur, la nationalité de la personne qui y vit ainsi que sa boisson préférée, son type de
cigarettes et son animal. Les contraintes sont les suivantes :
1. L’Anglais vit dans la maison rouge.
2. L’Espagnol possède un chien.
3. On boit du café dans la maison verte.
1
Ce problème est parfois attribué à Einstein ou à Lewis Carroll, et est aussi connu sous le nom de problème du zèbre.
1
2 I
au type de problème considéré et on trouve autant de sous-disciplines que de familles de problèmes : par
exemple, la programmation en nombre entiers (IP), la programmation linéaire (LP), etc.
La programmation par contraintes est également une discipline dont le but est, comme son nom
l’indique, de résoudre des problèmes posés sous la forme de contraintes. C’est cette branche de l’infor-
matique, dédiée à la résolution de tels problèmes, qui nous intéresse ici particulèrement. Son objectif est
la résolution exacte de problèmes sous contraintes dans lesquels la notion de qualité des solutions n’est
pas obligatoirement présente.
de résolution ou plus largement les interactions avec l’utilisateur, ou bien encore la mise au point de
techniques d’apprentissage automatique.
• Nous avons eu l’opportunité de faire l’application des techniques de programmation par contraintes
pour la modélisation et la résolution d’un problème de conception en robotique pour lequel la PPC
n’avait jamais été mise à contribution jusqu’ici. Une fois que nous sommes parvenus à modéliser
le problème, nous avons eu à modifier le cadre classique de résolution d’un solveur de contraintes
continues destiné à la satisfaction de contraintes uniquement, pour y intégrer la notion de recherche
rigoureuse d’optimum. Des expériences nous ont permis de mettre en évidence que notre approche
était compétitive par rapport à d’autres outils, à la pointe des techniques d’aujourd’hui en matière
d’optimisation.
• Nous avons mis en oeuvre une collaboration de deux solveurs, l’un spécialement dédié aux pro-
blèmes discrets, l’autre conçu en particulier pour les problèmes continus. L’idée directrice qui
nous a poussé à entreprendre cette réalisation est extrêmement naturelle : il semble évidemment
très prometteur de pouvoir profiter des avantages de chacun des outils dans leur champ de prédi-
lection respectif. Les techniques orientées continu sont relativement peu efficaces pour traiter la
partie discrète des problèmes mixtes et, réciproquement, les techniques orientées discret ne sont
pas efficaces pour traiter la partie continue de ces problèmes. En pratique, cette collaboration nous
permet d’abord de faciliter la phase de modélisation en offrant la possibilité de s’abstraire d’un
mécanisme de modélisation rigide et difficile à utiliser. En outre, cette collaboration nous permet
de résoudre plus efficacement certains types de problèmes mixtes.
• Nous avons ensuite cherché à faire un usage optimal d’une contrainte globale entière dans un sol-
veur continu. Pour ce faire, nous avons effectué une comparaison des différentes modélisations
possibles de la contrainte globale discrète alldifferent, ainsi que des techniques de filtrage qui
sont dédiées à chacune d’entre elles. Cette étude expérimentale nous a conduit à résoudre des pro-
blèmes d’optimisation issus du monde réel et à implémenter une variété d’algorithmes de filtrage
dans un solveur de contraintes continues. Il ressort de nos expériences qu’en ce qui concerne la
résolution de ces problèmes difficiles du monde réel, même s’il reste avantageux, surtout en termes
de robustesse, d’opter pour une modélisation globale de la contrainte et non son éclatement en une
clique de contraintes binaires, il n’est pas intéressant d’opter pour un algorithme de filtrage trop
complexe comme c’est le cas de l’algorithme de filtrage de type consistance de bornes le plus
répandu, même dans sa version la plus optimisée.
• Nous avons également réfléchi à une spécialisation des techniques de filtrage basées sur l’arithmé-
tique des intervalles, dédiée au cas des contraintes arithmétiques discrètes et mixtes. Réécrivant
dans un premier temps de nouvelles définitions à certaines opérations de l’arithmétique des inter-
valles, en les adaptant au contexte particulier d’intervalles d’entiers, et réfléchissant à une prise
en compte plus tôt au cours du calcul sur les intervalles des contraintes d’intégralité sur les va-
riables entières, nous avons pu dans un deuxième temps implémenter ces nouvelles idées dans
un solveur de contraintes continues. Cela nous a ainsi permis d’observer l’impact significatif que
peuvent avoir ces optimisations naturelles du calcul, d’une part sur l’efficacité théorique des mé-
thodes de filtrage basées sur cette arithmétique des intervalles, et d’autre part en pratique sur les
temps effectifs de résolution des problèmes.
6 I
II – État de l’art
(a) présentation des éléments essentiels de modélisation des problèmes qui nous intéressent en
programmation par contraintes et en programmation mathématique également,
(b) présentation des techniques majeures de résolution de ces problèmes par la programmation
par contraintes ;
(e) mise en évidence expérimentale des effets sur la résolution de problèmes mixtes des nouvelles
techniques mises en oeuvre.
La programmation par contraintes est un paradigme générique de résolution de problèmes dont l’une
des principales caractéristiques est de distinguer deux grandes tâches principales :
• modéliser le problème à résoudre : cette première tâche consiste à formaliser le problème à l’aide
de concepts dédiés, dans le but de le reformuler de manière plus abstraite ;
• résoudre le modèle ainsi obtenu : cette seconde étape vise à appliquer des techniques de résolution
qui vont permettre de trouver une (des) solution(s) s’il en existe.
L’objet du présent chapitre est de présenter en détails chacune de ces deux étapes. Dans la première
section, on s’intéressera donc à la modélisation des problèmes en programmation par contraintes, ainsi
qu’en programmation mathématique. La seconde section sera quant à elle consacrée aux techniques de
résolution de ces problèmes par la programmation par contraintes, qu’il s’agisse de méthodes dédiées
aux problèmes discrets, continus ou mixtes.
9
10 É ’
Définition 2.1 (CSP). Un problème de satisfaction de contraintes (CSP) est un triplet (X,D,C) avec :
• X un ensemble {X1 , X2 , . . . , Xn } de variables ;
• D un ensemble {D1 , D2 , . . . , Dn } de domaines : Xi ∈ Di , i = 1, . . . , n ;
• C un ensemble {C1 , C2 , . . . , Ct } de contraintes où chaque Ci définit un sous-ensemble du produit
cartésien des domaines des variables sur lesquelles elle porte : Ci (Xi1 , . . . , Xik ) ⊆ Di1 × · · · × Dik .
La notion de domaine désigne l’ensemble de valeurs que peut potentiellement prendre une variable.
Cette notion générique fait référence à des ensembles de natures potentiellement très différentes :
• un ensemble de valeurs symboliques : on peut vouloir représenter des couleurs, e.g. D = {rouge,
bleu, vert}, ou encore des jours de la semaine, par exemple D = {lundi, mercredi, samedi} ;
• un ensemble d’entiers non contigus : il est possible d’utiliser un ensemble d’entiers quelconques,
e.g. D = {5, 8, 12} ;
• un intervalle d’entiers : on peut également définir un ensemble d’entiers contigus, en compré-
hension, ne spécifiant que les bornes inférieure et supérieure de l’intervalle, par exemple D =
{1, . . . , 10} ;
• un intervalle de réels1 : de la même façon, on peut définir en compréhension un ensemble de réels,
en ne déclarant que ses bornes inférieure et supérieure, e.g. D = [−π, π].
Chaque variable possède ainsi un domaine initial qui contraint la valeur qu’on peut lui affecter. Si
sa valeur est à chercher dans N, dans Z ou dans tout ensemble discret de valeurs, on parlera de variable
discrète. Si sa valeur est à trouver dans R ou dans tout sous-ensemble continu de valeurs, on parlera
de variable continue. En fonction du type des variables du problème, on pourra donc distinguer trois
grandes catégories de problèmes de satisfaction de contraintes : selon qu’il ne contient que des variables
discrètes, que des variables continues, ou bien à la fois des variables discrètes et continues, on dira d’un
CSP qu’il est discret, continu ou mixte.
Une fois qu’on a défini pour chaque variable son domaine initial, on peut poser des contraintes les
reliant les unes aux autres, restreignant les combinaisons de valeurs autorisées. Là aussi il existe des
contraintes de différentes natures :
1
On reviendra dans la section suivante de ce chapitre sur l’utilisation en machine de cette notion d’intervalle de réels.
2.1. Modélisation des problèmes 11
Définition 2.2 (Solution). Une solution d’un CSP est un n-uplet S = (s1 , s2 , . . . , sn ), avec si ∈ Di ∀i,
pour lequel toute contrainte C j ∈ C est satisfaite : aucune contrainte C j (X j1 , . . . , X jk ) n’est violée quand
chacune des variables X ji sur lesquelles elle porte se voit affecter sa valeur correspondante s ji dans la
solution.
Selon les situations, résoudre un CSP peut consister à trouver une, toutes, ou la meilleure des solu-
tions. Dans ce dernier cas, on doit d’ailleurs fournir un moyen d’évaluer la qualité d’une solution, moyen
qui prend très souvent la forme d’une expression arithmétique définie sur un sous-ensemble des variables
du problème et qu’on va utiliser pour quantifier la qualité d’une solution donnée.
2.1.1.2 Exemples
À présent, nous allons examiner quelques exemples de CSP, purement théoriques ou au contraire
extrêmement concrets, pour comprendre en pratique comment on peut se servir de cette notion pour
formaliser des problèmes.
Exemple 2.2. [Mac77] Soit le CSP discret défini par :
• X = {X1 , X2 , X3 , X4 , X5 },
• D = {D1 , D2 , D3 , D4 , D5 },
avec D1 = D2 = {a, b, c} et D3 = D4 = D5 = {a, b},
• C = {C13 , C23 , C34 , C35 , C45 },
avec C13 = C23 = {(b, a), (c, a), (c, b)}, C34 = C35 = {(a, b)} et C45 = {(b, a)}.
Dans l’exemple ci-dessus, une contrainte Ci j entre deux variables Xi et X j est définie en extension par
l’ensemble des couples (u, v) qui la satisfont, le premier (resp. deuxième) élément du couple étant une va-
leur du domaine de Xi (resp. X j ). Remarquons qu’une autre façon d’exprimer en extension une contrainte
est d’expliciter des couples non pas de valeurs, mais de couples variable-valeur e.g. {(Xi , u), (X j , v)}.
Ce problème utilise des variables à domaine symbolique et des contraintes en extension définissant
explicitement les combinaisons de valeurs autorisées. En outre, ce problème n’a aucune solution. En
effet, on constate que les paires de valeurs autorisées par les contraintes ne permettent pas de trouver
une affectation de l’ensemble des variables du problème telle que toutes les contraintes soient satisfaites
2
Ce type de contrainte est également appelée contrainte conditionnelle.
3
Il existe plusieurs centaines de contraintes globales correspondant à des usages et des applications diverses et variées. Elles
sont exhaustivement répertoriées dans [BCDP07] par exemple.
12 É ’
simultanément, e.g. il n’y a aucune valeur pour X5 qui soit simultanément compatible avec les valeurs
autorisées pour X3 et X4 par les contraintes C35 et C45 . Cet état de fait est d’ailleurs plus facilement
visible si l’on considère une représentation graphique du problème :
Exemple 2.3. [Mac77] Deux représentations possibles pour le CSP de l’exemple précédent :
En effet, considérant par exemple le schéma de droite, où les sommets représentent les variables du
problème et les combinaisons de valeurs autorisées sont exprimées par un arc reliant ces valeurs dans
le domaine des variables concernées, on s’aperçoit qu’aucune valeur de X5 n’est reliée à la fois à une
valeur de X3 et à une valeur de X4 . Le schéma de gauche, où les tables que portent les arcs sont les paires
de valeurs autorisées pour les variables reliées, représente quant à lui une autre forme de visualisation
possible pour le CSP, qui pourra s’avérer plus adéquate dans certaines circonstances, en particulier dans
le cas où les domaines des variables du CSP comporteraient un grand nombre de valeurs. Avant de
poursuivre, notons à propos de ce problème purement discret que les contraintes qui le définissent, ici
données en extension sous la forme d’une liste de tuples valides, peuvent également dans le cas présent
être exprimées en intension sous la forme C xy : x > y, en utilisant l’ordre lexicographique.
L’exemple suivant est, quant à lui, purement continu et montre une autre représentation graphique
possible pour un CSP :
Exemple 2.4. [Gou00] Soit le CSP continu défini par :
• X = {x, y, z, k},
• D = {D x , Dy , Dz , Dk }
avec D x = [3, 4], Dy = [4.5, 8], Dz = [−7, 9]
et Dk = [0, 12],
x2 = y
x + y = z
.
• C=
z=k+y
t = k − 5
Les contraintes de ce CSP continu sont les sommets de son graphe représentatif, et un arc étiqueté
par un ensemble de variables relie deux contraintes partageant ces variables. En outre, ce CSP n’admet
aucune solution. Sans entrer dès maintenant dans les détails à propos du filtrage des domaines, qui sera
abordé en profondeur dans la suite de ce chapitre, on peut d’ores-et-déjà constater que le domaine D x =
[3, 4] de x ne contient aucun réel tel qu’il existe un y compatible au sens de la contrainte x2 = y :
intuitivement, on comprend bien que les valeurs compatibles de y se trouvent dans le domaine [9, 16],
dont l’intersection avec Dy = [4.5, 8] est vide.
2.1. Modélisation des problèmes 13
Pour finir, voici un problème avec une signification plus concrète, qui de plus s’avère être un CSP
mixte i.e. faisant intervenir à la fois des variables discrètes et des variables continues :
Exemple 2.5. [HSG01] Soit un circuit électrique composé d’une résistance R1 de 0.1 MΩ connectée en
parallèle à une résistance variable R2 de valeur comprise entre 0.1 MΩ et 0.4 MΩ. Un condensateur
K est connecté en série à ces deux résistances. On dispose d’un kit de composants électroniques ne
contenant que les condensateurs de valeurs 1 µF, 2.5 µF, 5 µF, 10 µF, 20 µF et 50 µF. On cherche à ce
que le temps de charge du condensateur soit compris entre 0.5 s et 1 s, i.e. que la tension aux bornes du
condensateur atteigne 99% de sa valeur entre 0.5 s et 1 s. Le CSP correspondant est :
• X = {t, R2 , K},
• D = {Dt , DR2 , DK }
avec Dt = [0.5, 1], DR2 = [0.1 × 106 , 0.4 × 106 ],
et DK = {1×10−6 , 2.5×10−6 , 5×10−6 , 10×10−6 , 20×10−6 , 50×10−6 },
• C = {Cohm , Ccharge }
avec Cohm : R1 = R11 + R12 la loi d’Ohm pour des résistances en
parallèle,
−t
et Ccharge : 1 − exp( R·K ) = 0.99 la loi de charge d’un condensateur
en série.
La seule possibilité pour le condensateur cherché ici est de K = 10 µF. Il existe cependant une
infinité de couples de valeurs réelles pour (t, R2 ) qui satisfont toutes les contraintes du problème pour
cette valeur de K. Dans ces conditions, on peut imaginer pertinent de formuler différemment le problème,
par exemple en ajoutant l’objectif de minimiser R. Nous allons maintenant voir, dans la suite de cette
section, que cette notion d’objectif est prépondérante en programmation mathématique.
Définition 2.3 (Problème d’optimisation sous contraintes). Un problème d’optimisation sous contraintes
a la forme générique suivante :
Un problème d’optimisation consiste à trouver les valeurs des variables qui minimisent la valeur
correspondante de la fonction objectif. En outre, cette minimisation peut faire intervenir ou non des
contraintes à satisfaire. Certaines variables peuvent être entières, d’autres réelles, et à l’instar d’un CSP,
14 É ’
un problème d’optimisation peut être discret, continu ou mixte selon le type des variables qui le consti-
tuent.
Le langage de modélisation en usage pour exprimer des problèmes d’optimisation se veut le plus
simple possible et utilise une syntaxe basée essentiellement sur celle des mathématiques. Notons quand
même l’utilisation de primitives de modélisation plus complexes issues de la programmation par contraintes,
i.e. les contraintes globales, comme une tendance apparue récemment [Hoo07].
Les lecteurs intéressés par une description plus approfondie de ce modèle et en particulier la signifi-
cation physique des différentes contraintes utilisées pourront se reporter à [KK94]. Quant à la solution
optimale pour ce problème, on obtient x1 = 10, x2 = 0.283 et x3 = 1.180, correspondant à une valeur de
l’objectif de 2.798.
2.1.3 Généralisation
Nous pouvons à présent généraliser la notion de problème de satisfaction de contraintes et celles
de problème d’optimisation sous contraintes en une notion plus générale dont l’un et l’autre sont des
instances particulières :
Notion de consistance
Les pionniers de cette technique sont W, M, M, et F dans [Wal72],
[Mon74], [Mac77], et [Fre78]. Si c’est à M qu’on doit la notion de consistance pour désigner
5
Cela sera détaillé au début de la sous-section suivante à propos de l’exploration de l’espace de recherche.
16 É ’
la propriété qui sert à distinguer les valeurs qu’on garde et celles qu’on élimine, il ne fait que poser un
mot sur une notion qu’ont déjà définie précédemment W et M. Ce sont les mêmes F
et M qui reviennent en 1985 dans [MF85] sur la complexité de leurs propres algorithmes et
de leurs évolutions. Dans [MH86], on trouvera ensuite une nouvelle version de leur algorithme, puis
notamment dans [HDmT92] et dans [Bes94].
En effet, depuis son apparition cette notion de consistance a été utilisée sous de nombreuses décli-
naisons. En fonction de leur efficacité à réduire les domaines et de leur rapidité à le faire, ces différentes
évolutions du constat initial de M se sont transformées, ont été améliorées, revues et corrigées
les unes après les autres et sont plus que jamais en usage aujourd’hui.
Les techniques de filtrage varient selon qu’il s’agit d’un CSP discret, continu ou mixte. Dans cette
sous-section, nous passons en revue les différentes méthodes mises en application pour chaque type de
CSP.
Les valeurs a des variables X1 , X2 et X4 ne sont compatibles avec aucune valeur de X3 ni de X5 , et ont
donc été supprimées.
L’algorithme de filtrage par consistance d’arc dans sa version optimale AC7 a dans le pire des cas
une complexité temporelle en O(ed2 ) et spatiale en O(ed), avec e le nombre de contraintes et d la taille
du plus grand domaine [BFR99].
Parallèlement, M a défini dans [Mon74] une propriété très forte des solutions des réseaux
de relations. C’est cette propriété, que M a rebaptisé ensuite consistance de chemin. Elle est
bien plus forte que la consistance d’arc car elle permet cette fois de retirer des domaines des n-uplets de
valeurs incompatibles, qui n’apparaissent dans aucune solution. Mais elle est en revanche plus difficile à
implémenter et coûteuse à calculer.
Définition 2.6 (Consistance de chemin). Etant donné un CSP P = (X, D, C), la paire de domaines
(Dα , Dω ) des variables Vα et Vω est chemin-consistante par rapport à un ensemble ordonné E = {C1 , . . . ,
Cm } de contraintes Ci portant chacune sur un ensemble Vi de variables tel que Vα ∈ V1 , Vω ∈ Vm et
Vi ∩ Vi+1 , ∅ ∀i < m ssi pour toute paire de valeurs de (Dα , Dω ) qui satisfait toutes les contraintes de E
portant sur Vα ∪ Vω , il existe au moins une combinaison de valeurs des variables des Vi telle que toute
contrainte de E est satisfaite.
À l’instar de la consistance d’arc (AC), la consistance de chemin (PC) permet de filtrer les domaines
pour en éliminer les valeurs inconsistantes qui ne peuvent participer à aucune solution. Ce mécanisme
est toutefois encore plus coûteux que le précédent. Son algorithme initial a lui aussi été revu et corrigé de
nombreuses fois. Les deux algorithmes ont évolué de concert et, aujourd’hui, l’algorithme PC5 a dans le
pire des cas une complexité temporelle en O(n3 d3 ) et spatiale en O(n3 d2 ), avec d la taille du plus grand
domaine et n le nombre de variables [Sin96].
En outre, M a démontré l’équivalence entre la consistance de chemin et la 2-consistance
de chemin i.e. le cas particulier de la consistance de chemin quand les contraintes de F ne portent que
sur une variable hormis Vα et Vω . Cela signifie par exemple qu’il n’est pas nécessaire de considérer la
consistance des chemins de longueur supérieure à 2 dans la représentation sous forme de graphe d’un
CSP où chaque sommet matérialise une variable.
Exemple 2.9. [Bes06] Soit un CSP P = (X, D, C) avec X = {X1 , X2 , X3 }, D = {D1 , D2 , D3 } où D1 =
D2 = D3 = {1, 2}, et C = {C1 , C2 } où C1 : X1 , X2 et C2 : X2 , X3 . P n’est pas chemin-consistant
puisque il n’y a aucune valeur de X2 qui satisfasse les contraintes de C pour les combinaisons de valeurs
(X1 , X3 ) ∈ {(1, 2), (2, 1)}. En revanche, le CSP P′ = (X, D, C ∪ {C3 : X1 = X3 }) est chemin-consistant.
De nombreuses évolutions de ces deux consistances ont ensuite vu le jour. On trouve aujourd’hui
un grand nombre de consistances discrètes différentes, chacune étant dotée d’une puissance de filtrage
relativement plus forte ou plus faible que les autres, et permettant de trouver le bonne équilibre selon les
cas entre puissance de filtrage et rapidité d’exécution7 .
f = (−1) s · m · 2E
où s est le bit de signe, m la mantisse et E l’exposant. On définit ainsi l’ensemble F des nombres flottants.
Remarquons l’inclusion F ⊂ R. Dans ces conditions, il n’y a pas d’autre choix que d’arrondir chaque réel
au flottant le plus proche, malgré les erreurs d’approximation que cela peut engendrer. L’exemple 2.10
montre les divergences qu’on peut observer en pratique entre le calcul sur les réels d’une part et le calcul
sur les flottants d’autre part, en fonction du nombre de bits utilisés pour représenter les nombres en
machine.
Exemple 2.10. [Rum88], [Che07] La fonction de R a pour expression :
x
f (x, y) = 333.75y6 + x2 (11x2 y2 − y6 − 121y4 − 2) + 5.5y8 +
2y
La valeur exacte de l’évaluation de f (77617, 33096) est −54767/66192, soit environ -0.8273960599. La
table ci-après contient les résultats obtenus pour f (77617, 33096) en fonction du nombre de bits utilisés
pour représenter les nombres, au moyen de l’outil MuPad 2.5.3 sur un processeur Intel Pentium M :
Précision Valeur
7 −2.47588 × 1027
16 −5.764075226 × 1017
24 −134217728.0
32 −0.84375
64 −0.8273960599 . . .
Comme on vient de le voir, un nombre réel n’est pas forcément représentable en machine, ce qui
introduit des erreurs d’évaluation dans les calculs. Cependant, il est possible de borner la valeur d’un
réel par deux flottants et d’effectuer les calculs sur de tels intervalles pour pallier le problème des erreurs
d’arrondis, grâce à l’arithmétique des intervalles.
domaine a été consacrée à l’étude des propriétés mathématiques de structures algébriques basées sur
la notion d’intervalles, à l’élaboration d’algorithmes de calcul d’approximations fiables et précises des
solutions faisables ou optimales de systèmes de contraintes linéaires et non-linéaires et à l’étude de la
convergence de ces algorithmes [BVH97]. Ainsi, depuis [Cle87], [OV90] et l’intégration efficace dans
les langages de programmation logique comme Prolog du calcul sur les intervalles mis au point quelques
décennies plus tôt, il est possible de raisonner sur des CSP dont les domaines des variables ne sont
pas discrets et énumérés, mais continus et seulement exprimés comme un encadrement par une borne
inférieure et une borne supérieure définies comme des nombres flottants.
Définition 2.7 (Intervalle à bornes flottantes). Étant donné l’ensemble F ∪ {−∞, +∞} des flottants au-
quel on a ajouté −∞ et +∞, et la relation ≤ sur cet ensemble, l’intervalle à bornes flottantes [a, b] est
l’ensemble des nombres réels compris entre a et b :
[a, b] = {x ∈ R : a ≤ x ≤ b} ∀a, b ∈ F ∪ {−∞, +∞}
On trouve dans [Lho93] les bases pour l’adaptation effective aux domaines continus des techniques
discrètes de filtrage par consistance, avec par exemple la définition d’une consistance d’arc aux bornes
des domaines. Puis, dans [BMVH94], l’arithmétique des intervalles —telle que définie par Moore dans
[Moo66]— est mise en lumière pour exhiber les principales caractéristiques utiles des intervalles à bornes
flottantes et de l’ensemble I de ces intervalles.
Définition 2.8 (Enveloppe). Soit r un sous-ensemble de R. L’enveloppe de r, notée (r), est le plus petit
intervalle à bornes flottantes contenant r :
∀ a ∈ R, (a) = ⌊a⌋F , ⌈a⌉F
∀ a ≤ b, a, b ∈ R, ([a, b]) = ⌊a⌋F , ⌈b⌉F
∀ a1 , . . . , an ∈ R, ({a1 , . . . , an }) = ([min(a1 , . . . , an ), max(a1 , . . . , an )])
où ⌈a⌉F (resp. ⌊a⌋F ) est le plus petit (resp. grand) élément de F ∪ {−∞, +∞} plus grand (resp. petit) ou
égal à a.
Étant donnée une opération arithmétique élémentaire ◦ ∈ {+, −, ×, /} et deux intervalles I1 , I2 , l’opé-
ration correspondante sur les intervalles est définie de la manière suivante :
I1 ◦ I2 = {x ◦ y : x ∈ I1 , y ∈ I2 }
Ce qui nous donne en pratique les règles de calcul suivantes, pour deux intervalles [a, b] et [c, d] :
En outre, l’arithmétique des intervalles nous permet de la même façon de définir une extension aux
intervalles des fonctions réelles :
Définition 2.9 (Extension aux intervalles). Étant donnée une fonction réelle f : Rn → R, la fonction sur
les intervalles F : In → I est une extension aux intervalles de f ssi ∀I ∈ In , { f (x) : x ∈ I} ⊆ F(I).
Notons qu’il existe de nombreuses extensions aux intervalles, basées sur des transformations symbo-
liques ou des relaxations numériques : les formes de B, les formes de H et l’extension de
T par exemple [Gou00], [BG06]. La Figure 2.1 donne un aperçu de la différence observable selon
l’extension utilisée pour l’évaluation d’une fonction.
20 É ’
Figure 2.1 – [Gou00] Évaluation de f (x) = (x6 /4 − 1, 9x5 + 0, 25x4 + 20x3 − 17x2 − 55, 6x + 48)/50 en
forme naturelle (gauche), en forme de B (droite), pour une largeur d’intervalle unitaire de 0,02.
Du point de vue de la résolution de contraintes, il est dès lors possible de raisonner sur les domaines
continus à la place de raisonner sur les domaines discrets, et d’utiliser l’arithmétique des intervalles
et notamment l’extension aux intervalles des fonctions réelles pour définir de nouvelles techniques de
filtrage par consistance, dédiées aux domaines continus cette fois.
Définition 2.10 (Consistance d’enveloppe). Etant donné un CSP P = (X, D, C), un domaine Di ∈ D est
enveloppe-consistant par rapport à une contrainte C ∈ C ssi Di = (Di ∩ {ai ∈ R : ∃a1 ∈ D1 . . . ∃ai−1 ∈
Di−1 ∃ai+1 ∈ Di+1 . . . ∃an ∈ Dn : C(a1 , a2 , . . . , an )}).
La consistance d’enveloppe, aussi appelée consistance 2B [Lho93], est une relaxation de la consis-
tance d’arc. Il s’agit d’une propriété locale, liée aux bornes du domaine d’une variable au niveau d’une
seule contrainte : une contrainte c est 2B consistante si pour n’importe quelle variable x de c, il y a
des valeurs dans les domaines de toutes les autres variables de c pour lesquelles c est satisfaite quand x
vaut la borne inférieure (resp. supérieure) de son domaine [CDR99]. Mais à cause des erreurs d’arrondis
inhérentes au calcul avec des nombres flottants, obtenir l’enveloppe d’un ensemble de réels est une tâche
difficile en elle-même, difficulté que l’algorithme de filtrage dédié HC3 parvient en partie à contourner
en rendant 2B consistant le système qu’on aura obtenu en décomposant les contraintes du problème en
contraintes aussi simples que possible10 : par exemple, la contrainte x + y · z = t pourra être décomposée
en deux contraintes y·z = α et x +α = t, introduisant une nouvelle variable α, et le filtrage sera réalisé sur
le système obtenu ainsi par décomposition [BGGP99]. Enfin, l’algorithme HC4, très proche de HC3, est
capable d’aboutir au même résultat mais sans avoir pour cela recours à l’introduction de telles variables
auxiliaires —l’exemple 2.11 ci-après montre comment cet algorithme HC4 utilise une contrainte pour
réduire les domaines des variables de cette dernière [Gou00]. Considérant une expression arithmétique
de la contrainte sous forme d’arbre où les feuilles sont des variables ou des constantes, les noeuds étant
des opérateurs arithmétiques et la racine un symbole de relation, le filtrage d’un pavé par HC4 est réalisé
en deux phases, l’évaluation ascendante et la propagation descendante :
Exemple 2.11. [Gou00] Filtrage par l’algorithme HC4 pour la contrainte 2x = z − y2 avec les domaines
D x = [0, 20], Dy = [−10, 10] et Dz = [0, 16] :
10
On utilisera le terme de contraintes primitives.
2.2. Techniques de résolution 21
Définition 2.11 (Consistance de pavé). Etant donné un CSP P = (X, D, C), un domaine Di ∈ D est
pavé-consistant par rapport à une contrainte C j ∈ C de la forme f j (x j1 , . . . , x jn ) = 0 ssi ∀i, Di = ({ai ∈
Di : 0 ∈ F j (D1 , D2 , . . . , Di−1 , ({ai }), Di+1 , . . . , Dm )}), où F j est une extension aux intervalles de f .
Notons que dans le cas où la variable considérée n’apparaît qu’une seule fois dans la contrainte, la
consistance d’enveloppe est équivalente et moins coûteuse à calculer que la consistance de pavé [CDR99].
En outre, la consistance de pavé est plus faible que la consistance d’arc, comme le montre l’exemple
suivant :
Exemple 2.12. [HMK97] Soit la contrainte x1 + x2 − x1 = 0 et les domaines D1 = D2 = [−1, 1]. D2
n’est pas arc-consistant puisqu’il n’existe aucune valeur n1 ∈ D1 telle que n1 + 1 − n1 = 0. Pourtant D2
est pavé-consistant puisque ([−1, 1] + [−1, −1+ ] − [−1, 1]) ∩ [0, 0] et ([−1, 1] + [1− , 1] − [−1, 1]) ∩ [0, 0]
ne sont pas vides.
De plus, la consistance de pavé est une relaxation de la consistance d’enveloppe dont le principe
est de remplacer le test de satisfaction de contraintes sur les réels par une procédure de réfutation sur
les intervalles, et pour atteindre la consistance de pavé, l’algorithme de filtrage BC3 cherche à trouver
des valeurs consistantes en bordure de Di [BMV94]. Comme le montre en pratique l’exemple 2.13,
l’implémentation standard de cet algorithme utilise une procédure de recherche dichotomique :
Exemple 2.13. [BG06] Soit y − x2 = 0 une contrainte avec D x = [0, 1] et Dy = [0, 4]. La borne
inférieure de Dy est pavé-consistante puisque 0 appartient à l’intervalle 0 − D2x = [−1, 0]. À l’inverse,
la borne supérieure est inconsistante. Le domaine de y peut alors être itérativement divisé en 2 pour y
trouver la valeur consistante la plus grande :
11
Cette définition ne fait référence qu’à des contraintes d’égalité, mais elle peut être ensuite facilement étendue aux inégalités
en considérant que f (x) ≤ 0 ⇔ f = z, z ∈ [−∞, 0] et f (x) ≥ 0 ⇔ f = z, z ∈ [0, +∞].
22 É ’
mixtes du deuxième type, qui utilisent des fonctions de transtypage réel vers entier, elles sont approxi-
mées par des contraintes numériques, e.g. on remplace i = ⌈r⌉ par r ≤ i et i ≤ r + 1, avant d’appliquer
les méthodes dédiées aux contraintes continues. Enfin, concernant les contraintes du troisième type, elles
sont traitées comme des contraintes continues dans lesquels les variables discrètes ont été au préalable
temporairement instanciées au point-médian de leur domaine courant. On reviendra sur ce principe en
fin de chapitre pour davantage d’explications —cf p.29.
2.2.1.4 Propagation
Quand un domaine a été réduit, la question se pose des répercussions sur les autres domaines avec
lesquels sont partagées des contraintes. Ceci est brièvement illustré par l’exemple introductif suivant :
Exemple 2.15. [Gra05] Soient les domaines D x = [0, 5] et Dy = [0, 4], et soient deux contraintes y = x+1
(c1 ) et y = x − 1 (c2 ). L’utilisation de c1 pour un filtrage par consistance d’enveloppe provoque les
réductions :
• Dy = [0, 4] ∩ [0, 5] + 1 = [1, 4]
• D x = [0, 5] ∩ [0, 4] − 1 = [0, 3]
Puis l’utilisation de c2 provoque les réductions :
• Dy = [1, 4] ∩ [0, 3] − 1 = [1, 2]
• D x = [0, 3] ∩ [1, 4] + 1 = [2, 3]
Enfin, l’utilisation à nouveau de c1 entraîne la réduction :
• Dy = [1, 2] ∩ [2, 3] + 1 = ∅
Le problème n’a aucune solution.
Ainsi, certaines réductions peuvent être très avantageusement enchaînées car si leur effet seul n’est
pas très important, l’effet conjugué de toutes peut avoir une capacité filtrante bien plus grande. Ce mé-
canisme d’itération du filtrage pose cependant un certain nombre de questions. Par exemple, l’itération
va-t-elle se terminer ? Les filtrages de domaines finissent-ils obligatoirement par ne plus rien ôter ? Com-
ment optimiser l’enchaînement des différents filtrages ? etc. Le cadre général des itérations chaotiques
va permettre de répondre à nombre de ces questions.
Cadre formel
A l’image de [GH88], [MR91] et [Fal94] notamment, nombre d’études théoriques ont déjà été réa-
lisées pour fournir à la propagation un cadre formel, avec théorèmes et propriétés prouvés, de façon à
pouvoir se doter d’un mécanisme réellement contrôlable et efficace. [Apt00] et [Bes06] sont eux aussi
le reflet d’une démarche de formalisation dont le but est de mieux connaître et maîtriser les tenants et
aboutissants de la propagation13 .
Définition 2.12 (Itération). Soit un ordre partiel (D, ⊑) avec un plus petit élément ⊥ et un ensemble fini
de fonction F = { f1 , . . . , fk } sur D. Une itération de F est une séquence infinie de valeurs d0 , d1 , . . .
définie inductivement par
d0 = ⊥
(
d j = fi j (d j−1 ), j , 0
où i j est un élément de [1..k]. Une séquence croissante d0 ⊑ d1 ⊑ d2 . . . finit par se stabiliser à d si
di+1 = di pour un certain j ≥ 0 avec i ≥ j.
13
La suite de cette section est tirée essentiellement de [Apt00].
24 É ’
Ce formalisme s’illustre ici de la façon suivante : l’ordre partiel s’applique aux CSP au regard d’une
consistance donnée, le plus petit élément ⊥ étant le CSP initial et les di étant obtenus par filtrages suc-
cessifs, et les fonctions sont les algorithmes de filtrage de domaine. Par la suite ne seront étudiées que
les itérations de fonctions ayant certaines propriétés :
Définition 2.13 (Monotonie). Soit un ordre partiel (D, ⊑) avec un plus petit élément ⊥ et une fonction f
sur D. f est monotone si x ⊑ y implique f (x) ⊑ f (y) ∀x ∀y.
Le rôle de la monotonie est éclairé par la simple observation suivante :
Lemme 2.1 (Stabilisation). Soit un ordre partiel (D, ⊑) avec un plus petit élément ⊥ et F un ensemble
fini de fonctions monotones sur D, et supposons qu’une itération de F finit par se stabiliser en un point
fixe d commun aux fonctions de F . Alors d est le plus petit point fixe commun aux fonctions de F .
Preuve. Considérons e un point fixe commun aux fonctions de F et prouvons que d ⊑ e. Soit d0 , d1 , . . .
l’itération en question. Pour un certain j ≥ 0, on a di = d pour i ≥ j. Il suffit de prouver par induction sur
i que di ⊑ d. Il est évident que c’est vrai pour i = 0 puisque d0 = ⊥. Supposons à présent que c’est vrai
pour un certain i ≥ 0. On a di+1 = f j (di ) pour un certain j ∈ [1..k]. Du fait de l’hypothèse d’induction et
de la monotonicité, il s’ensuit que f j (di ) ⊑ f j (e), et donc que di+1 ⊑ e puisque e est un point fixe de f j .
Fixons à présent un ordre partiel (D, ⊑) avec un plus petit élément ⊥, et un ensemble fini de fonction
F = { f1 , . . . , fk } sur D. On s’intéresse à calculer le plus petit point fixe commun aux fonctions de F .
C’est dans ce but qu’on étudie l’Algorithme 2.1, inspiré par un algorithme similaire tiré de [MR99].
Intuitivement, (i) signifie qu’il devra être impossible à la fonction de choix de retourner indéfiniment
une fonction qui ne réduit pas le domaine, et (ii) que l’invariant de boucle de l’algorithme est que toute
fonction h ∈ F − G est telle que h(d) = d. On garantit de cette façon une réduction optimale.
Théorème 2.1 (Propagation). L’algorithme de propagation est tel que :
(i) Toute exécution qui termine calcule dans d le plus petit point fixe commun des fonctions de F .
(ii) Toute exécution termine.
Preuve. (i) Montrons que le point fixe calculé est effectivement le plus petit des points fixes. Soit α un
point fixe commun. Prouvons par induction que α ∈ d à chaque pas de l’algorithme. C’est vrai pour le
2.2. Techniques de résolution 25
premier pas. Pour le k-ième pas, on a α ∈ dk par hypothèse et g(α) ∈ g(dk ) car g est monotone. Comme
α ∈ g(α) et dk+1 = g(dk ), il vient α ∈ dk+1 .
(ii) Soit l’ordre produit de (D, ⊐) et (N, <) défini sur les éléments de D × N par
où on utilise ⊐ l’ordre inverse défini par d1 ⊐ d2 ssi d2 ⊑ d1 ∨ d2 , d1 . Les couples (x, Card(G)) calculés
par l’itération sont strictement décroissants au sens de cet ordre produit, puisqu’à chaque fois soit le
domaine est réduit, soit le nombre d’éléments dans la liste G diminue. Or par hypothèse (D, ⊑) est fini
donc (D, ⊐) est bien fondé, ce qui implique la terminaison.
Ainsi, une forme possible de l’algorithme update, respectant les spécifications ci-dessus, peut être
celle de l’Algorithme 2.2 :
Ajoutons que la fonction update décrite ci-dessus peut être améliorée en tenant en compte de l’idem-
potence des fonctions de filtrage :
Définition 2.14 (Idempotence). Soit une fonction f sur un domaine D. f est idempotente si f ( f (x)) =
f (x) ∀x ∈ D.
Il ne sera donc pas judicieux d’ajouter à l’ensemble G des fonctions de filtrage idempotentes puisque,
quelle que soit la réduction qu’elles viennent de réaliser, les rappeler sur le domaine qu’elles ont elles-
mêmes réduit ne peut avoir aucun effet supplémentaire.
Toutes les combinaisons de valeurs sont générées et testées, intégralement, une par une.
Cette technique très naïve est facilement améliorable en ce que chaque combinaison testée est évi-
demment très similaire à la précédente. C’est de ce constat qu’est né le mécanisme de backtracking —
retour-arrière en français dans le texte— qui reprend la combinaison précédemment testée pour l’étendre
de manière incrémentale, tout en vérifiant à chaque étape la satisfaction des contraintes :
A chaque étage de l’arbre, on affecte une nouvelle valeur à une variable, on teste la satisfaction des
contraintes, puis on continue en profondeur avec la variable suivante, jusqu’à ce qu’il n’y ait plus de
variable à affecter et tant que les contraintes peuvent être satisfaites.
Le produit du déroulement de cet algorithme de recherche est appelé arbre de recherche. On peut en
outre se rendre compte du gain obtenu par filtrage des domaines en observant l’arbre de recherche cor-
respondant au même CSP que ci-dessus, mais lorsqu’un filtrage par consistance d’arc a été au préalable
appliqué aux domaines initiaux :
Exemple 2.18. Arbre de recherche pour le même CSP mais après filtrage par consistance d’arc :
2.2. Techniques de résolution 27
En effet, il n’y a plus que 2 valeurs à tester dans les domaines de X1 et X2 , une seule dans le domaine de
X4 .
Les deux paramètres qui régissent principalement le déroulement du backtracking concernent la mé-
thode de choix de variable et celle de choix de l’ordre d’instanciation de ces valeurs. Des stratégies de
choix de variables communément employées en matière de CSP discrets sont celle dite du fail-first, où
l’on choisit en priorité la variable de plus petit domaine, ou bien celle du most-constrained, où l’on pré-
fère choisir d’abord la variable la plus contrainte. L’idée qui gouverne ces choix est de chercher à éviter
d’explorer une longue branche de l’arbre se soldant par un échec. Pour cela, on privilégie les variables qui
sont les plus sujettes à provoquer des inconsistances pendant le filtrage. En ce qui concerne l’affectation
d’une valeur à la variable choisie, la règle générale consiste simplement à créer une branche pour chaque
valeur du domaine. Dans certains cas cependant, on peut définir un ordre sur les valeurs à affecter, selon
leur sémantique dans le problème considéré [Heu06]. On peut également limiter le nombre de branches
créées en un noeud, en n’énumérant qu’un nombre réduit de valeurs du domaines et en conservant les
autres dans une dernière branche.
On peut citer [DP88] comme exemple d’alternative au backtracking : une utilisation judicieuse d’une
consistance d’arc directionnelle alliée avec un bon choix de variable permet de résoudre certains CSP
sans retour-arrière. D’autres se sont également employés à circonscrire le besoin de backtrack aux plus
petites portions possibles du CSP, on peut citer notamment les méthodes de consistance adaptative utili-
sant des coupe-cycles [Dec92]. Enfin, de nombreuses améliorations de l’algorithme initial ont vu le jour
depuis son apparition dans les années 1970 et sont présentées dans l’Annexe A.
Branch-and-prune
L’idée de l’algorithme d’exploration dédié aux domaines continus est de procéder à une bissection
plutôt qu’à une énumération pour parcourir l’espace de recherche : à chaque noeud, on déclenche la
phase de propagation des réductions, puis lorsque les domaines ne peuvent plus être significativement
réduits, alors on choisira une variable dont le domaine va être coupé en deux parties. Chaque moitié
sera ensuite explorée dans une branche de l’arbre, de la même façon qu’avec un domaine discret chaque
branche correspond à une valeur du domaine initial —les autres domaines étant simplement recopiés à
l’identique. Cet algorithme dit de branch-and-prune que nous donnons ici —cf Algorithme 2.3— est
celui donné par V H dans [HMK97], évolution des travaux fondateurs de C vis-à-vis
de la résolution de problèmes à variables réelles [Cle87] mais aussi très nettement dans la continuité des
premières idées de M sur le sujet [Moo79].
28 É ’
À l’inverse du cas d’un CSP discret, une heuristique de choix de variables en usage en ce qui concerne
la résolution d’un CSP continu va donner la préférence à la variable dont le domaine est le plus large. Ceci
peut s’expliquer par le fait qu’un domaine large signifie qu’on a du mal à le réduire grâce aux contraintes,
et que le couper ainsi en deux permet de le réduire malgré l’inefficacité potentielle des contraintes. Notons
par ailleurs qu’une autre heuristique largement en usage consiste simplement à choisir les variables les
unes après les autres, en tourniquet —round-robin en anglais.
Exemple 2.19. [Cle87] Arbre de résolution pour le CSP de contrainte14 X · X = 2 avec X ∈ [1, 2] :
L’intervalle de départ est successivement réduit et coupé en deux jusqu’à ce qu’on obtienne un intervalle
de largeur minimale.
Ajoutons que la résolution d’un problème de satisfaction de contraintes doté d’une fonction-objectif
requiert de maintenir à jour des bornes fiables pour cette dernière [Han92]. De cette manière, on peut
14
Cet exemple est tiré des travaux fondateurs de C [Cle87], où le langage de modélisation n’englobe pas encore la
notion de puissance d’une variable.
2.2. Techniques de résolution 29
utiliser la fonction-objectif de la même manière que les autres contraintes pour filtrer les domaines. En
outre, de nouvelles bornes peuvent être obtenues de manière rigoureuse dans les régions satisfaisant les
contraintes et on utilise à cet effet des algorithmes de preuve d’existence [Kea96], l’idée étant de mettre
à jour l’encadrement de la valeur optimale connue à chaque fois qu’on découvre un point qui satisfait
toutes les contraintes du problème et dont la valeur de la fonction-objectif est meilleure que celle qu’on
avait pu obtenir jusque là. Enfin, l’ordre de traitement des pavés ainsi que la stratégie de branchement
sont deux paramètres importants pour l’efficacité du processus de résolution, puisqu’en effet le but est de
trouver le plus rapidement possible les pavés qui améliorent le plus l’encadrement de la fonction-objectif
et ainsi filtrer le plus efficacement possible l’espace de recherche [CR97]. Ce point sera abordé plus en
profondeur au chapitre suivant.
les autres variables discrètes et enfin les variables continues. En outre, l’utilisation d’un filtrage de type
forward checking permet d’améliorer les performances obtenues, spécialement quand il s’agit d’un pro-
blème dont les solutions sont réparties en un faible nombre de régions contiguës.
Comme le montre la figure ci-dessus, x1 ≤ 1 est une inégalité valide parce qu’elle n’ôte aucun point
entier de l’ensemble admissible. Cependant, il s’agit également d’une coupe pour la relaxation linéaire
du problème, qu’elle renforce.
III
M ODÉLISATION ET RÉSOLUTION
D ’ UN PROBLÈME DE CONCEPTION
EN ROBOTIQUE
Celui qui en toutes choses suit la Voie et règle ses principes sur la Voie,
Parce qu’il aspire à l’union suprême,
La Voie l’accueille avec joie.
Aussi sa conduite, ses projets, ses oeuvres ou ses abstentions,
Ont-ils d’heureux résultats.
3.1 Contexte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2 Modélisation du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2.1 Organe terminal sans jeu dans les articulations . . . . . . . . . . . . . . . . . 32
3.2.2 Modélisation des erreurs dues au jeu dans les articulations . . . . . . . . . . . 34
3.2.3 Erreur de pose maximum de l’organe terminal . . . . . . . . . . . . . . . . . . 35
3.3 Modification du cadre de résolution classique . . . . . . . . . . . . . . . . . . . . . . 35
3.4 Expériences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
Dans ce chapitre, nous modélisons et résolvons le problème de la prédiction des erreurs de position
et de rotation de l’organe terminal d’un système robotique mécanique, aussi appelées erreurs de pose, en
fonction du jeu dans ses articulations ou, plus précisément, du jeu introduit dans les articulations par des
imperfections de fabrication.
Après avoir montré comment le problème peut être exprimé sous la forme de deux problèmes conti-
nus d’optimisation sous contraintes, nous montrons comment le classique cadre continu de résolution
en programmation par contraintes peut être modifié pour nous permettre de gérer rigoureusement et
globalement ces problèmes d’optimisation difficiles. En particulier, nous présentons des expériences pré-
liminaires dans lesquelles notre optimiseur est très compétitif comparé aux méthodes les plus efficaces
présentées dans la littérature, tout en fournissant des résultats plus robustes.
31
32 Ḿ ́ ’ ̀
3.1 Contexte
La précision d’un système robotique mécanique est un aspect crucial pour l’exécution d’opérations
qui demandent de l’exactitude. Cependant, cette précision peut être affectée négativement par des erreurs
de positionnement et/ou d’orientation, aussi appelées erreurs de pose, de l’organe terminal du manipu-
lateur. Une principale source d’erreurs de pose réside dans le jeu qu’il y a dans les articulations, qui
introduit des degrés de liberté supplémentaires de déplacement entre les éléments pairés par les join-
tures du manipulateur. Il apparaît en outre que réduire ces déplacements augmente à la fois les coûts de
production et la difficulté d’assemblage du mécanisme. Actuellement, cette problématique est d’ailleurs
une tendance majeure de recherche en robotique [FS98, Inn02, MZL09]. Une façon de traiter le jeu
dans les articulations est de prévoir son impact sur les erreurs de pose. Dans ce but, les déplacements
impliqués ainsi qu’un ensemble de paramètres du système peuvent être modélisés par deux problèmes
d’optimisation sous contraintes, dont le maximum global va caractériser l’erreur de pose maximum. Il est
donc obligatoire d’obtenir un encadrement rigoureux de ce maximum global, ce qui élimine la possibi-
lité d’utiliser des optimiseurs locaux. Cette modélisation peut être vue comme un ensemble de variables
prenant leurs valeurs dans des domaines continus, accompagné d’un ensemble de contraintes et d’une
fonction-objectif. Il s’avère qu’un tel modèle est difficile à résoudre expérimentalement : les temps de
résolution peuvent augmenter de manière exponentielle avec la taille du problème, qui dans ce cas varie
avec le nombre d’articulations du manipulateur.
S’il existe des travaux de recherche traitant de la modélisation du jeu dans les articulations et la pré-
vision des erreurs dues à ce phénomène [FS98, Inn02], peu d’entre eux ont pour objet les techniques de
résolution [MZL09], et aucun ne fait appel à la programmation par contraintes. Dans ce chapitre, nous
investiguons donc l’usage de la programmation par contraintes pour résoudre de tels problèmes. En par-
ticulier, nous combinons l’algorithme classique de branch-and-bound avec l’arithmétique des intervalles
et les techniques de filtrage associées. L’algorithme de branch-and-bound nous permet de gérer la partie
optimisation du problème tandis que les temps de résolution peuvent être diminués grâce à ces tech-
niques de filtrage. En outre, la fiabilité du processus est garantie par l’utilisation de l’arithmétique des
intervalles, ce qui est intrinsèquement requis par l’application qui nous intéresse ici.
et n articulations, où la liaison 0 est la base fixe et la liaison n l’organe terminal. Ensuite, on définit des
repères F j d’origine O j et d’axes X j , Y j , Z j , le repère F j étant associé à la ( j − 1)ème liaison pour j variant
de 1 à n + 1. Le torseur suivant transforme F j en F j+1 :
Rj tj
" #
Sj = T , (3.1)
03 1
où R j est une matrice de rotation 3 × 3, t j ∈ R3 pointe de l’origine de F j vers F j+1 , et 03 est le vecteur
zéro à trois dimensions. De plus, S j peut s’exprimer de la façon suivante :
n
Y
P= S j, (3.3)
j=1
Cependant, nous devons inclure de petites erreurs dans l’équation (3.3) si nous voulons prendre en
compte le jeu dans les articulations.
Z3
Y3
a2
link 2
b2 Z2 end-effector
X3
Y2 F3
a1
X2
link 1 F2
b1 Z1
joint 2
X2
Y1
θ1
X1 joint 1
F1
Figure 3.1 – [BSG∗ 10] Manipulateur sériel composé de deux liaisons pivots.
34 Ḿ ́ ’ ̀
Du fait du jeu dans les articulations, le repère Fn+1 associé à l’organe terminal est transformé en Fn+1
′ .
D’après [MLS94], le déplacement amenant le repère F j en F j est donné par la matrice exponentielle de
′
δS j , eδS j . Il en résulte que le torseur représentant la pose de l’organe terminal décalé peut être calculé
tout au long de la chaîne cinématique par :
n
Y
P′ = eδS j S j , (3.10)
j=1
z zp
xp
y yp
Fp
Afin de mesurer l’erreur sur la pose de la plate-forme en mouvement, on calcule le torseur ∆P qui
transforme la pose nominale Fn+1 en la pose décalée Fn+1′ tout au long de la chaîne cinématique :
1
Y n
Y
−1 ′
∆P = P P = S−1
j eδS j S j . (3.11)
j=n j=1
D’après [CC09], il s’avère que ∆P peut aussi bien être représenté de manière vectorielle comme un
torseur de faible déplacement δp, de la façon suivante :
j
n Y
R j O3×3
X " #
−1
adj(Sk ) δs j , avec adj(S j ) ≡
δp =
T jR j R j
j=1 k=n
qui est la transformation adjointe du torseur S j et T j ≡ ∂(t j × x)/∂x la matrice produit vectoriel de t j .
1 L ← {B0 }
2 tant que L , ∅ et non critère_arrêt faire
3 (B, L) ← extract(L)
4 B ← contractC (B)
5 si splitable(B) alors
6 {B′ , B′′ } ← split(B)
7 L ← L ∪ {B′ , B′′ }
8 fin
9 fin
10 retourner L
taille des pavés dans L, l’algorithme s’arrêtant quand ils ont atteint une largeur inférieure à une précision
donnée. RealPaver [GB06] implémente un algorithme de ce type et nous nous en sommes servi de point
de départ.
Un algorithme de branch-and-prune peut être transformé en branch-and-bound capable de traiter des
problèmes de minimisation —les problèmes de maximisation étant traités d’une manière similaire. L’al-
gorithme 3.4 donne un exemple d’un tel algorithme de branch-and-bound. La fonction de coût est une
entrée supplémentaire, grâce à laquelle peut être maintenue une borne supérieure du minimum global
dans la variable m, initialisée à +∞ en début de la recherche (l.2). Cette borne supérieure est utilisée
périodiquement pour rejeter les parties de l’espace de recherche dont le coût lui est supérieur, et ce sim-
plement en ajoutant la contrainte f (x) ≤ m à l’ensemble des contraintes du problème (l.5). Finalement, la
borne supérieure est mise à jour à la ligne 6 (cf Algorithme 3.5), en cherchant des points faisables à l’inté-
1 L ← {B0 }
2 m ← +∞
3 tant que L , ∅ et non critère_arrêt faire
4 (B, L) ← extract(L)
5 B ← contractC∪{ f (X)≤m} (B)
6 (B, m) ← update(C, B, f, m)
7 si splitable(B) alors
8 {B′ , B′′ } ← split(B)
9 L ← L ∪ {B′ , B′′ }
10 fin
11 fin
12 retourner (m, L)
3.4. Expériences 37
rieur du pavé courant puis en évaluant la fonction de coût sur de tels points une fois trouvés. Dans notre
implémentation actuelle, des points sont générés aléatoirement dans le pavé courant et les contraintes
sont rigoureusement vérifiées en utilisant l’arithmétique des intervalles avant de mettre à jour la borne
supérieure. Remarquons que si des algorithmes de branch-and-bound plus élaborés font généralement
appel à de la recherche locale pour trouver de bons points faisables, les expériences que nous présentons
à la section suivante, et que nous avons réalisées en utilisant cette simple méthode de génération aléa-
toire de points, montrent déjà de bonnes performances. L’algorithme de branch-and-bound maintient un
encadrement du minimum global qui converge vers celui-ci à condition que des points faisables soient
trouvés pendant la recherche, ce qui est garanti dans le cas de contraintes d’inégalité qui ne sont pas
singulières.
Une implémentation efficace de l’algorithme de branch-and-bound requiert que la liste des pavés L
soit gérée avec attention : il faut la maintenir triée par rapport à la borne inférieure de l’objectif évalué
sur chaque pavé. Ainsi, à chaque fois qu’un pavé B est inséré dans L, l’évaluation aux intervalles f (B)
est calculée et la borne inférieure de cette évaluation est utilisée pour maintenir la liste triée. Ensuite, à la
ligne 4 on extrait le premier pavé de L de telle sorte que les régions les plus prometteuses de l’espace de
recherche soient explorées les premières, ayant pour effet des améliorations drastiques. Un autre avantage
de maintenir L triée est que le premier pavé de la liste contient la valeur la plus basse de borne inférieure
de l’objectif parmi tous les pavés. En conséquence, on peut utiliser cette valeur et mesurer sa distance
à la borne supérieure courante m pour arrêter l’algorithme quand l’encadrement du minimum global a
atteint la précision souhaitée.
L’algorithme de branch-and-bound que nous venons de décrire a été implementé dans RealPaver.
C’est cette implémentation que nous avons utilisée pour réaliser les expériences que nous présentons
dans la section suivante.
3.4 Expériences
Nous avons effectué un ensemble d’expériences pour analyser les performances de notre approche
pour résoudre le problème d’erreur de pose. Nous l’avons comparée avec GAMS/BARON [GAM] et
ECLi PSe [WNS97]. GAMS/BARON est un système de programmation mathématique largement utilisé
38 Ḿ ́ ’ ̀
de résolution de problèmes d’optimisation1 . Quant à ECLi PSe , c’est un des rares systèmes de program-
mation par contraintes capables de résoudre des problèmes continus d’optimisation.
Nous avons testé 8 modèles, 4 sur l’erreur maximum en translation (T) et 4 sur l’erreur maximum en
rotation (R). Pour chaque type de modèle, nous considérons à chaque fois des manipulateurs ayant de 2 à
5 joints. La Table 3.1 présente les résultats obtenus. Les colonnes 1 et 2 montrent le nombre de joints et le
type de problème. Les colonnes 3 à 5 donnent des informations pertinentes quant à l’instance considérée
(nombre de variables, de contraintes, d’opérations arithmétiques dans la fonction-objectif). Les colonnes
5 à 10 contiennent les temps de résolution observés en utilisant GAMS/BARON, RealPaver, et ECLi PSe .
Nos expériences ont été réalisées sur un Pentium D à 3 Ghz avec 2 Go de RAM. Le temps de résolution
présenté ici est à chaque fois le meilleur temps obtenu sur 5 lancements.
Problem Size
#joints Type GAMS RP ECLi PSe
#var #ctr #op
2 T 12 8 28 0.08 0.004 >60
2 R 6 4 18 0.124 0.004 >60
3 T 18 12 135 0.124 0.008 t.o.
3 R 9 6 90 0.952 0.004 t.o.
4 T 24 16 374 0.144 0.152 t.o.
4 R 12 8 205 2.584 0.02 t.o.
5 T 30 20 1073 0.708 >60 t.o.
5 R 15 10 480 9.241 0.26 t.o.
Les résultats montrent que notre approche est plus rapide dans presque tous les cas. Pour les ins-
tances de taille plus réduites comme 2R, 2T, 3R et 3T, RealPaver offre des performances remarquables
et peut être jusqu’à 100 fois plus rapide que GAMS/BARON. Un telle rapidité de convergence peut être
expliquée en partie par l’efficacité des techniques de filtrage, basées ici sur la consistance d’enveloppe
et bien optimisées ; mais surtout, par le fait qu’il n’y a pas besoin de calculs vraiment très précis pour
ce problème particulier. En fait, il n’est pas utile d’atteindre une précision de l’ordre de 10−8 , du fait
principalement de la limitation physique des outils de conception qui ne sont pas capable d’atteindre une
précision aussi grande.
Pour des instances plus grandes comme 4R et 5R, RealPaver reste nettement plus rapide. Cependant,
à mesure que la taille du problème augmente, et en particulier dès que le nombre de variables excède 20,
la convergence cesse d’être instantanée. Ceci est un phénomène commun qu’on peut expliquer à la fois
par le caractère exponentiel par rapport au nombre de variables de la complexité de ce type d’algorithme
de résolution, et aussi au fait que la fonction-objectif elle-même croît exponentiellement avec le nombre
de variables sur lesquelles elle porte. Néanmoins, nous estimons que les temps de résolution restent
raisonnables au vu de la complexité et de la taille des problèmes ainsi que des techniques utilisées2 .
Cependant, les temps de résolution ne sont pas le seul point à examiner dans cette situation. Au
contraire, il est également important de discuter la fiabilité des solutions calculées par BARON. En effet,
il a déjà été mentionné dans [LMR07] que celles calculées par BARON n’étaient pas fiables en général.
1
Dans notre version du logiciel, la résolution des problèmes linéaires (resp. non linéaires) est effectuée par CPLEX (resp.
MINOS).
2
Des expériences ultérieures ont en outre montré que notre approche permettait en l’état de résoudre l’instance 5T en un
temps de l’ordre de la seconde, à condition que les points tirés au sort par l’algorithme soit tels qu’ils permettent de mettre à
jour l’objectif de façon significative
3.5. Conclusion 39
Nous avons pu vérifier ici, en utilisant RealPaver, soit que l’encadrement optimal calculé par BARON
n’était pas faisable, soit qu’à l’inverse cet encadrement faisable n’était pas optimal.
Pour finir, remarquons le fait que bien que RealPaver et ECLi PSe utilisent des techniques de program-
mation par contraintes, le second est complètement dépassé même pour les plus petites de nos instances.
Ce n’est pas surprenant cependant, étant donné que le premier a été conçu pour résoudre des problèmes
continus non linéaires alors qu’ECLi PSe est un outil plus générique et extensible.
3.5 Conclusion
Dans ce premier chapitre, nous avons montré comment utiliser la programmation par contraintes pour
résoudre le problème du calcul de l’erreur de pose maximum d’un système robotique mécanique due au
jeu dans les articulations. Nous avons commencé par modéliser le problème comme un CSP continu doté
d’une fonction-objectif, avant de montrer comment transformer le schéma branch-and-prune de résolu-
tion qu’emploie classiquement la programmation par contraintes, en un algorithme branch-and-bound
dédié à la résolution rigoureuse de problèmes d’optimisation sous contraintes. Enfin, nous avons présenté
des résultats expérimentaux, obtenus à partir de plusieurs instances du problème initial, et montré que
notre approche était très compétitive par rapport à des outils à la pointe des techniques d’aujourd’hui.
Il est important de noter que cette approche n’est pas spécifique à la robotique mais est applicable à
n’importe quel problème où sont requis des calculs rigoureux, gage de fiabilité des résultats numériques
obtenus. De plus, des travaux futurs devraient permettre de mettre en oeuvre d’autres critères de filtrage
utilisant la fonction-objectif et d’ainsi accélérer encore la convergence de l’algorithme d’optimisation,
permettant de traiter efficacement des problèmes de plus grande taille. Par ailleurs, nous avons déjà
commencé à explorer des pistes pour l’amélioration des performances de notre optimiseur global sous
contraintes d’inégalités, et notamment l’utilisation de la contrainte additionnelle f ′ (B) = 0, qui permet
sur des exemples triviaux d’effectuer des réductions de domaine drastiques, mais qui pour le problème
traité ici ne s’est pas montrée d’une utilité particulière. D’autres expériences devraient nous permettre de
mieux comprendre cet état de fait.
Notons également qu’au niveau de l’exploitation des résultats du point de vue de la robotique et non
de l’informatique cette fois, l’enjeu sera pour nous de faire l’application de la méthode que nous avons
40 Ḿ ́ ’ ̀
décrite ici afin de calculer l’erreur de pose maximum pour un échantillon de points de l’espace de travail
du robot. Cela nous permettra de produire la figure dite des isocontours de l’erreur maximum en fonction
de la position du robot dans l’espace de travail, d’une nature similaire à celle présentée en Figure 3.3,
et qui permettra de mieux connaître l’impact du jeu dans les articulations sur l’erreur de pose du robot.
Il pourra être intéressant au passage de comparer les performances de nos algorithmes avec celles des
algorithmes génétiques habituellement utilisés à cet effet. Enfin, d’autres travaux sont déjà en cours pour
également prendre en compte les imperfections d’usinage des pièces du robot.
IV
C OLLABORATION DE DEUX
SOLVEURS DISCRET ET CONTINU
Dans ce chapitre, nous mettons en oeuvre une collaboration de deux solveurs de contraintes, l’un
dédié aux contraintes discrètes, l’autre dédié aux contraintes continues. Notre idée directrice est de partir
des outils existants pour voir ce qui manque vraiment à chacun pour traiter efficacement les problèmes
mixtes, que ce soit sur le plan de la modélisation ou bien des techniques de résolution, et d’interfacer les
outils entre eux pour que l’un comble les lacunes de l’autre.
Nous commençons par une présentation détaillée des capacités de modélisation des principaux outils
existants, dans le but de choisir lesquels faire collaborer. Ensuite, nous décrivons la manière dont nous
avons pu mettre en oeuvre la collaboration entre ces outils. Finalement, nous présentons les résultats
expérimentaux que nous avons obtenus sur un ensemble de problèmes de test.
41
42 C
• des outils boîte-noire qui prennent en entrée un fichier contenant le modèle écrit dans un langage
dédié uniquement à la modélisation ;
• des interpréteurs CLP, qui travaillent avec des modèles définis sous la forme de programmes écrits
en langage déclaratif de type Prolog étendu ;
• des bibliothèques d’objets et de fonctions dédiés aux CSP écrits dans un langage à objets type C++
ou Java, qui permettent d’exprimer les modèles sous la forme de programmes compilables dans ce
langage.
Nous ne cherchons pas à être exhaustifs face au grand nombre d’outils non commerciaux disponibles
et ne référençons plutôt que les outils les plus connus, dans chacun des différents panels correspondant
à chaque type d’outil existant. Par contre, nous nous efforcerons de l’être concernant la mise-en-avant
des propriétés de chacun, au regard de notre problématique de résoudre des problèmes mixtes, pour
proposer un aperçu le plus complet possible des différentes possibilités offertes par chacun en termes de
modélisation. On trouvera dans la Table 4.1 un condensé de l’ensemble de nos observations à ce sujet.
Méthode utilisée
Pour mieux cerner les capacités qu’offre chacun des outils, nous proposons d’étudier la possibilité
de modéliser trois CSP de trois types différents : discret, continu, mixte. Le premier est un CSP pure-
ment discret bien connu dans la communauté de la programmation par contraintes, il s’agit de résoudre
l’équation SEND+MORE=MONEY où chaque lettre correspond à un chiffre distinct :
Exemple 4.1. Le problème SEND+MORE=MONEY est modélisé à l’aide du CSP discret suivant :
X = S , E, N, D, M, O, R, Y
D = DS = DE = DN = DD = D M = DO = DR = DY = {0, . . . , 9}
S , 0,
M 0,
C =
,
1000S + 100E + 10N + D + 1000M + 100O + 10R + E
= 10000M + 1000O + 100N + 10E + Y
Le deuxième modèle que nous avons cherché à résoudre est le CSP continu qui correspond à la re-
cherche de la configuration 3D d’une molécule de cyclohexane, décrite par un système de trois équations
non linéaires :
Exemple 4.2. Le problème cyclohexane est modélisé à l’aide du CSP continu suivant :
4.1. Présentation des outils existants 43
X = x, y, z
D = D x = Dy = Dz = [−1 × 108 , 1 × 108 ]
2
y (1 + z2 ) + z(z − 24y) = −13,
x2 (1 + y2 ) + y(y − 24x) = −13,
C =
z2 (1 + x2 ) + x(x − 24z) = −13
Quant au troisième, il s’agira d’un CSP mixte, dont l’objet est le déplacement dans son espace de
travail d’un bras de robot en fonction de ses caractéristiques :
Exemple 4.3. Le problème robot est modélisé à l’aide du CSP mixte suivant :
i, j, a, b, α, β, x, y, d
X =
D =
Di = D j = {2, . . . , 5}, Da = Db = [2, 8], D x = [0, 10], Dy = [0, 6], Dα = Dβ = [0, π],
Dd = [1.99, 2.01]
a · sin(α) = b · sin(β − α) + y,
a · cos(α) = x − (b · cos(β − α)),
a cos(α) 10,
· ≤
a sin(α) 6,
C = · ≤
2 + (y − 4)2 ≤ 4,
(x 8)
−
a (i 1)d,
= −
b = ( j − 1)d
Minion
Minion est un solveur de contraintes de type boîte-noire dont le langage de modélisation utilisé par
défaut est basé sur la notion de représentation matricielle d’un CSP pour pouvoir définir très simplement
des contraintes sur les lignes, colonnes ou plans [GJM06].
Plus récemment, il a été interfacé avec le langage haut niveau de modélisation de problèmes combina-
toires Essence [FHJ∗ 05] —un sous-ensemble de celui-ci nommé Essence’ pour être exact— au moyen du
traducteur de modèles Taylor [GMR07]. C’est sous cette dernière forme que nous donnons ici le modèle
du problème discret S END + MORE = MONEY (cf Prog. 4.1).
Comme on peut s’en rendre compte à la lecture du modèle, le langage de modélisation Essence’
utilisé pour exprimer des modèles à résoudre par Minion est intuitif dans son utilisation et n’est pas
difficile à comprendre. On s’aperçoit ainsi que Minion permet la déclaration de contraintes globales et de
contraintes arithmétiques sur des variables entières.
Cependant, le solveur ne permet pas la déclaration de variables de type réel. Dans l’optique de ré-
soudre des problèmes mixtes, il faut donc y incorporer la notion de domaine de type intervalle à bornes
flottantes, ainsi que tout ce qui permet une gestion efficace de ce type de domaine : arithmétique des
44 C
intervalles, possibilité d’utiliser des expressions arithmétiques sur des variables réelles pour filtrer les
domaines, bissection de domaines et heuristique de branchement associée.
P. 4.1 : Modèle Minion pour le problème discret S END + MORE = MONEY.
RealPaver
RealPaver [GB06] est lui aussi un outil de résolution qui traite des modèles exprimés dans un langage
de haut niveau, dédié à la modélisation de problèmes de satisfaction de contraintes.
Comme un peut s’en apercevoir avec le Programme 4.2, et comme le signifie en outre le nom de
RealPaver, c’est un solveur de contraintes capable de paver des domaines réels étant données des équa-
tions/inéquations restreignant les valeurs admissibles pour les variables : utilisant la notion d’intervalle
à bornes flottantes, il lui est possible d’effectuer des calculs corrects sur les réels.
Quant au Programme 4.3, il représente la formulation du problème discret S END + MORE =
MONEY. Comme on peut le voir dans cet exemple, RealPaver ne dispose pas de contraintes globales
ni d’un opérateur de différence. S’il est tout de même possible d’exprimer ces contraintes en utilisant
l’opérateur de valeur absolue et en contraignant son résultat pour disposer d’une contrainte de différence,
et d’utiliser une clique de telles contraintes pour modéliser un alldifferent1 , on s’aperçoit ici que
c’est le manque de contraintes globales, discrètes, qui fait le point faible de RealPaver dans l’optique de
résoudre efficacement des problèmes mixtes.
Finalement, le Programme 4.4 montre la modélisation du problème mixte du positionnement d’un
bras de robot. Comme dans l’exemple précédent, on constate que RealPaver est capable de gérer l’inté-
gralité de tout ou partie de ses variables et permet de déclarer des variables discrètes dont le domaine est
un intervalle d’entiers.
Ainsi, RealPaver est déjà capable de traiter des problèmes continus, discrets ou mixtes. Cependant,
l’absence de contraintes globales est sa principale faiblesse au vu de notre problématique, que ce soit en
terme d’aisance de modélisation ou d’efficacité de résolution.
1
Ce point sera d’ailleurs abordé plus en profondeur au chapitre suivant.
4.1. Présentation des outils existants 45
• l. 1 : entête du modèle
• l. 2-7 : déclaration d’une constante p et de variables réelles pour lesquelles on cherche à
atteindre la précision p
• l. 8-11 : ajout des contraintes de configuration de la molécule
• l. 12 : pas de variables auxiliaires, on lance donc la résolution sur toutes les variables puisque
solve ∗ est identique à solve x, y, z
P. 4.3 : Modèle RealPaver pour le problème discret S END + MORE = MONEY.
• l. 1 : entête du modèle
• l. 2-6 : déclaration de variables entières
• l. 9 : ajout des contraintes empêchant S et M de valoir 0
• l. 11-20 : ajout d’une clique de contraintes binaires pour que toute paire de variables prenne
un couple de valeurs différentes
• l. 22-24 : ajout de la contrainte arithmétique spécifiant la relation de somme entre les digits
• l. 26 : pas de variables auxiliaires, on lance la résolution sur toutes les variables
46 C
• l. 1 : entête du modèle ;
• l. 3-5 : déclarations de constantes ;
• l. 8-12 : déclarations des variables du modèle, certaines réelles et d’autres entières, en utilisant
les constantes définies en tête de modèle ou prédéfinies comme Pi ;
• l. 15-21 : ajout des contraintes du modèle ;
• l. 23 : pas de variables auxiliaires, on lance la résolution sur toutes les variables.
ECLi PSe
ECLi PSe est un langage de programmation logique à contraintes dont l’origine remonte au début des
années 1990 [AW07]. Aussi, exprimer un modèle pour ECLi PSe consiste à écrire un prédicat en langage
déclaratif, semblable à n’importe quel autre prédicat ou presque.
Le Programme 4.5 correspond au modèle discret S END + MORE = MONEY. Le langage utilisé
ici est très proche d’un langage de modélisation dédié même si cela reste la syntaxe d’un langage de
la famille Prolog qui a été élargie, pour accueillir l’expression de contraintes arithmétiques notamment.
ECLi PSe dispose en outre d’un catalogue de contraintes globales comme cet alldifferent qu’on utilise
ici pour contraindre les différents digits à prendre des valeurs distinctes.
Mais ECLi PSe offre également la possibilité de déclarer des variables réelles et des contraintes conti-
nues, comme le montre le Programme 4.6 : le symbole $ permet de faire la distinction entre domaine
discret et domaine continu, entre contrainte discrète et contrainte continue.
Ainsi, il est tout à fait possible d’exprimer et de résoudre grâce à ECLi PSe des problèmes mixtes.
Le Programme 4.7 montre d’ailleurs l’exemple du problème mixte du bras de robot. Même si dans ce
cas précis n’interviennent pas simultanément contraintes arithmétiques continues et contraintes globales
4.1. Présentation des outils existants 47
P. 4.5 : Modèle ECLi PSe pour le problème discret S END + MORE = MONEY.
1 sendmore(Digits) :− 8 1000*S+100*E+10*N+D
2 Digits = [S,E,N,D,M,O,R,Y], 9 + 1000*M+100*O+10*R+E
3 Digits : : [0. .9], 10 #= 10000*M+1000*O+100*N+10*E+Y,
4 11
5 S #\= 0, 12 alldifferent(Digits),
6 M #\= 0, 13
7 14 labeling(Digits).
• l. 1 : entête du modèle
• l. 2-3 : déclarations des variables, entières, et de leurs domaines initiaux
• l. 5-6 : ajout des contraintes empêchant S et M de valoir 0
• l. 8-10 : ajout de la contrainte arithmétique spécifiant la relation de somme entre les digits
• l. 12 : ajout d’une contrainte globale alldifferent sur les variables du problème
• l. 14 : lancement de l’exploration de l’espace de recherche en énumérant les domaines des
variables
discrètes, il est tout à fait possible de modéliser de tels problèmes et la force d’ECLi PSe vis-à-vis de notre
problématique est bel-et-bien l’étendue de ses capacités de modélisation.
Cependant, si le moteur de résolution d’ECLi PSe est parfaitement opérationnel en ce qui concerne le
traitement des problèmes discrets, il n’en est pas de même quant à la résolution de problèmes continus ou
mixtes. Comme nos expériences le mettront en lumière plus loin dans ce chapitre, l’efficacité d’ECLi PSe
est très rapidement mise à mal par la présence de variables continues dans le modèle2 .
1 cyclohexane(Vars) :− 6 13+Y^2*(1+Z^2)+Z*(Z−24*Y)$=0,
2 7 13+Z^2*(1+X^2)+X*(X−24*Z)$=0,
3 Vars = [X,Y,Z], 8 13+X^2*(1+Y^2)+Y*(Y−24*X)$=0,
4 Vars $ : : [−1e8. .1e8], 9
5 10 locate([X,Y,Z],[X,Y,Z],1e−6,lin).
• l. 1 : entête du modèle
• l. 2-3 : déclarations des variables, réelles (opérateur $ ::) et de leurs domaines initiaux
• l. 5-7 : ajout des contraintes du modèle
• l. 10 : lancement de l’exploration de l’espace de recherche en bissectant les domaines jusqu’à
une largeur finale de 1 × 10−6
2
Il s’agit en effet d’un programme que son ancienneté a rendu très complet, mais qui, bien que fournissant de très bons
résultats en terme de correction de la solution, est hélas loin d’être optimal quant à la rapidité d’obtention de cette solution, en
particulier en ce qui concerne les problèmes continus ou mixtes.
48 C
P. 4.7 : Modèle ECLi PSe pour le problème mixte du bras de robot.
1 arm(Vars) :− 13 A * sin(Alpha)
2 14 $= B * sin(Beta − Alpha) + Y,
3 Tx is 10.0, Ty is 6.0, 15 A * cos(Alpha)
4 X0 is 8, Y0 is 4, 16 $= X − (B * cos(Beta − Alpha)),
5 D is 1.99 2.01, 17 A * cos(Alpha) $=< Tx,
6 18 A * sin(Alpha) $=< Ty,
7 Vars = [A,B,Alpha,Beta,X,Y,I,J], 19 (X − X0)^2 + (Y − Y0)^2 $=< 4,
8 A $ : : [2.0. .8.0], B $ : : [2.0. .8.0], 20 A $= (I − 1) * D,
9 Alpha $ : : [0. . pi], Beta $ : : [0. . pi], 21 B $= (J − 1) * D,
10 X $ : : [0. .Tx], Y $ : : [0. .Ty], 22
11 I : : [2. .5], J : : [2. .5] ; 23 locate([A,B,Alpha,Beta,X,Y,I,J],
12 24 [A,B,Alpha,Beta,X,Y,I,J],1e−6,lin).
• l. 1 : entête du modèle
• l. 3-5 : déclarations de constantes
• l. 7-11 : déclarations des variables, réelles ou entières, et de leurs domaines initiaux
• l. 13-21 : ajout des contraintes du modèle
• l. 23-24 : lancement de l’exploration de l’espace de recherche en bissectant les domaines
jusqu’à une largeur finale de 1 × 10−6
CHR
L’autre outil que nous présentons dans cette catégorie des outils de résolution traitant des modèles
exprimés en langage déclaratif est Constraint Handling Rules (CHR). Il existe des implémentations de ce
langage dans différents langages hôtes comme Java, Prolog et même Haskell3 . Difficile à classer, CHR
est un langage à écrire des solveurs plutôt qu’un solveur à proprement parler [Rai08].
Le Programme 4.8 correspondant à un modèle CHR pour résoudre le problème discret S END +
MORE = MONEY. La syntaxe déclarative rappelle fortement Prolog. En fait, le langage CHR ne sert
pas tant à exprimer le modèle qu’à écrire le solveur pour le résoudre, un solveur par type de modèle pour
ainsi dire.
Si CHR est capable en théorie de traiter les contraintes globales et les domaines autres que discret —
cette flexibilité est même une des raisons d’être de CHR— il faut d’abord trouver un solveur adéquat et à
défaut l’écrire. On trouve par exemple dans [Rai08] une méthode de génération automatique de solveurs
CHR pour les contraintes globales. Ainsi, au regard de notre problématique de la résolution efficace des
problèmes mixtes, l’essentiel du travail à réaliser ici consisterait à faire un état de l’art en profondeur du
millier d’articles traitant de CHR référencés à ce jour, et à fusionner tout ce qui a été déjà fait en CHR
avant de compléter les manques grâce à d’autres outils le cas échéant4 .
3
Pour les besoins de ce comparatif, nous nous sommes basés sur l’implémentation fournie avec YAP Prolog [SCDRA].
4
Nos travaux sont antérieurs à la sortie du livre [Frü09] qui constitue à n’en pas douter un parfait point de départ dans cette
direction.
4.1. Présentation des outils existants 49
P. 4.8 : Modèle CHR pour le problème discret S END + MORE = MONEY.
1 sum(1000,S,100,E,1,H1), 14 neq(S,0),neq(M,0),
2 sum(1,H1,10,N,1,H2), 15
3 sum(1,H2,1,D,1,H3), 16 neq(S,E), neq(S,N), neq(S,D), neq(S,M),
4 sum(1,H3,1000,M,1,H4), 17 neq(S,O), neq(S,R), neq(S,Y), neq(E,N),
5 sum(1,H4,100,O,1,H5), 18 neq(E,D), neq(E,M), neq(E,O), neq(E,R),
6 sum(1,H5,10,R,1,H6),
19 neq(E,Y), neq(N,D), neq(N,M), neq(N,O),
7 sum(1,H6,1,E,1,Z),
20 neq(N,R), neq(N,Y), neq(D,M), neq(D,O),
8 sum(10000,M,1000,O,1,G1),
21 neq(D,R), neq(D,Y), neq(M,O), neq(M,R),
9 sum(1,G1,100,N,1,G2),
10 sum(1,G2,10,E,1,G3), 22 neq(M,Y), neq(O,R), neq(O,Y), neq(R,Y),
11 sum(1,G3,1,Y,1,Z). 23
12 24 fd(S,0,9),fd(E,0,9),fd(N,0,9),fd(D,0,9),
13 25 fd(M,0,9),fd(O,0,9),fd(R,0,9),fd(Y,0,9).
• l. 1-11 : ajout de la contrainte arithmétique spécifiant la relation de somme entre les digits
• l. 14 : ajout des contraintes empêchant S et M de valoir 0
• l. 16-22 : ajout d’une clique de contraintes binaires pour que toute paire de variables prenne un
couple de valeurs différentes
• l. 24-25 : lancement de l’exploration de l’espace de recherche en donnant les domaines initiaux
des variables
Gecode
Commençons avec Gecode, une bibliothèque de résolution de problèmes sous contraintes écrite en
C++. Tout comme Minion, son efficacité à résoudre des problèmes, même imposants en taille, n’est plus
à démontrer [Tac09].
Le Programme 4.9 donne un bon aperçu de ce à quoi ressemble un modèle au format tel que le
traite Gecode. Il s’agit ainsi d’une classe C++ et d’une fonction main. Dans la première, les variables et
les contraintes sont déclarées et instanciées, et la résolution paramétrée. Dans la deuxième, on lance la
méthode de résolution associée à la classe qu’on vient de définir.
Pour autant efficace qu’il soit toutefois, Gecode ne propose pas le type de variable réel. Il n’est donc
pas possible de déclarer de domaines réels ni de modéliser des contraintes continues ou mixtes. Dans
l’optique de résoudre des problèmes mixtes, il faudrait donc y inclure la notion d’intervalle à bornes
flottantes, de même que tout ce qui permet de gérer efficament ce type de domaine, comme l’arithmétique
des intervalles, la possibilité de filtrer les domaines en utilisant des expressions arithmétiques pour filtrer
les domaines, le mécanisme de bissection de domaines et l’heuristique de branchement associée.
50 C
P. 4.9 : Modèle Gecode pour le problème discret S END + MORE = MONEY.
• l. 1-8 : déclaration d’une classe C++ dédiée au problème, et instanciation des paramètres hérités
comme le nombre de variables nl
• l. 10-12 : déclaration des variables
• l. 14-15 : ajout des contraintes empêchant S et M de valoir 0
• l. 17-20 : ajout de la contrainte arithmétique spécifiant la relation de somme entre les digits
• l. 22 : ajout d’une contrainte de type alldifferent portant sur l’ensemble des digits
• l. 24 : définition de la stratégie d’exploration
• l. 29-38 : initialisation et lancement de la résolution
Choco
Choco est une bibliothèque écrite en Java qui permet la modélisation et la résolution de problèmes
sous contraintes. L’implémentation a été faite avec deux objectifs : d’une part offrir une architecture
logicielle modulaire qui accepte aisément des extensions, et d’autre part, mettre en oeuvre une gestion
optimisée des événements de propagation [Lab00]5 .
Le Programme 4.10 correspond à la traduction du modèle S END+ MORE = MONEY. Remarquons
que la déclaration d’expressions arithmétiques pour les contraintes se fait de manière préfixée, à l’aide
de méthodes de la classe Problem : le langage Java ne permet pas de surcharger ses opérateurs et il est
donc impossible d’adopter une syntaxe plus naturelle comme le permet par exemple Gecode, écrit lui en
C++.
En outre, C permet la modélisation de problèmes continus en proposant d’instancier des va-
riables de type réel comme on peut le voir dans le Programme 4.11. On s’aperçoit même finalement
5
Les présents travaux ont été réalisés avant la récente refonte de la bibliothèque présentée dans [RJL∗ 08] et se basent sur la
version 1.2.06 qui lui est antérieure.
4.1. Présentation des outils existants 51
que Choco est déjà, dans l’état actuel, tout à fait capable de traiter des problèmes mixtes comme notre
problème du bras de robot (cf Prog. 4.12).
Cependant, ses limites se trouvent principalement dans les opérations arithmétiques proposées (peu
de non-linéaire, erreurs de calculs dans les fonctions trigonométriques) et les types de contraintes numé-
riques disponibles (égalité seulement). En outre, le langage de modélisation proposé est difficile à utili-
ser et les modèles obtenus sont peu lisibles, comme on peut déjà s’en apercevoir à travers les quelques
exemples présentés ici.
P. 4.10 : Modèle Choco pour le problème discret S END + MORE = MONEY.
• l. 1-7 : déclaration d’une classe Java et de sa fonction main dédiée à la résolution du problème
• l. 9-16 : déclaration des variables correspondant aux différents digits, de type entier, ayant leur
domaine énuméré, de domaine initial {0, . . . , 9}
• l. 18-23 : déclaration des variables correspondant à chaque ligne de la somme, de type entier,
ayant leur domaine borné, de domaine initial [0, 106 ]
• l. 26-27 : ajout d’une contrainte alldifferent portant sur l’ensemble des digits
• l. 29-30 : ajout des contraintes S , 0 et M , 0
• l. 32-42 : ajout des contraintes qui relient chaque variable correspondant à un des trois termes
de la somme aux digits qui le constitue, e.g. S END = S × 1000 + E × 100 + N × 10 + D
• l. 44-46 : ajout de la contrainte S END + MORE = MONEY
• l. 38 : lancement de la résolution
52 C
• l. 1-7 : déclaration d’une classe Java et de sa fonction main dédiée à la résolution du problème
• l. 8 : paramétrage de la précision à atteindre sur le domaine des variables réelles
• l. 10-12 : déclaration des variables du modèle, de type réel, de domaine initial [−1 × 108 , 1 × 108 ]
• l. 14-24 : constructions des expressions arithmétiques qui forment le corps des contraintes, e.g.
y2 (1 + z2 ) + z(z − 24y)
• l. 26-39 : ajout des contraintes du modèle en utilisant les expressions instanciées ci-dessus, en
spécifiant grâce à la fonction addBoxedVar le fait que les variables devront être réduite par
consistance de pavé
• l.41-46 : configuration du solveur associé au modèle sur la façon dont prendre en compte les
variables réelles dans la stratégie d’exploration
• l. 47 : lancement de la résolution
4.1. Présentation des outils existants 53
• l. 1-5 : déclaration d’une classe Java et de sa fonction main dédiée à la résolution du problème
• l. 7-14 : déclaration des variables continues du modèle, de type réel, avec leurs domaines
initiaux respectifs
• l. 16-17 : déclaration des variables discrètes du modèle, ayant leurs domaine bornés, de type
entier, avec leurs domaines initiaux respectifs
• l. 18-19 : déclaration de variables continues correspondant aux variables discrètes du modèle
• l. 20-21 : ajout de contraintes reliant chaque variable discrète à sa variable continue équivalente
• l. 23-33 : constructions des expressions arithmétiques qui forment le corps des contraintes, e.g.
(x − 8)2 + (y − 4)2 ≤ 4
• l. 35-37 : ajout des contraintes du modèle en utilisant les expressions instanciées ci-dessus
• l. 39-46 : ajout des autres contraintes du modèle en utilisant les variables continues homologues
des variables discrètes
• l. 48-54 : configuration du solveur associé au modèle sur la façon dont prendre compte les
variables réelles dans la stratégie d’exploration
• l. 56 : lancement de la résolution
54 C
4.1.4 Synthèse
La Table 4.1 récapitule l’ensemble de ce qui a été montré des capacités de modélisation6 de chacun
des solveurs. Quant à la Figure 4.1, elle représente sous une forme plus synthétique ce qui manque à
chaque outil pour résoudre efficacement les problèmes mixtes.
Au regard de nos observations, nous avons choisi en conséquence de faire collaborer Choco et Real-
Paver. En effet, il ne manque quasiment rien au premier pour être complètement opérationnel sur les
problème mixtes et c’est bel-et-bien le deuxième, RealPaver, qui est capable de lui apporter le peu qui
lui manque le plus efficacement7 parmi les outils à notre disposition. Ainsi, nous allons confier à Choco
le rôle principal puisqu’il est en effet déjà très complet, et nous baser sur RealPaver pour combler sa
principale faiblesse : la gestion de la partie continue ou mixte du modèle.
Figure 4.1 – Schéma de synthèse, la taille de la police utilisée est proportionnelle à la quantité estimée de
travail à fournir pour compléter l’outil en question : Choco, CHR et ECLi PSe sont déjà assez bien dotés,
RealPaver, Minion et Gecode sont quant à eux plus spécifiques à un type de problème.
RealPaver
RealPaver est un outil logiciel pour la modélisation et la résolution de systèmes de contraintes non
linéaires. La version 1.0 que nous avons utilisée est utilisable sous la forme d’une bibliothèque C/C++
avec une API qui lui permet de s’intégrer facilement dans d’autres applications.
L’implémentation du parseur, des expressions arithmétiques ainsi que leur évaluation en utilisant
l’arithmétique des intervalles et la gestion en mémoire des boîtes sont réalisées en C. Ainsi, que ce soit
par exemple le type des variables ou des noeuds des arbres représentant les opérations arithmétiques,
cela n’est pas matérialisé par des classes particulières d’objets mais par des valeurs spécifiques pour un
champ donné d’une structure de donnée générique. À l’inverse, les opérateurs de réduction et tout ce qui
concerne la propagation mais aussi les stratégies d’exploration, ainsi que les autres structures de données
haut-niveau du solveur, sont réalisés en C++ et utilisent le concept d’héritage pour se spécialiser le cas
échéant. Les différentes possibilités de configuration du solveur sont ainsi exprimées sous la forme d’une
hiérarchie de classes au sein de laquelle on peut choisir celle qui correspond au paramétrage désiré.
Choco
Choco est une bibliothèque Java initialement dédiée à la résolution de problèmes discrets. Puis elle
a été étendue aux problèmes continus ou mixtes par l’adjonction du type de variable réel ainsi que du
type de domaine intervalle à bornes flottantes et des règles d’évaluation des expressions arithmétiques se
basant sur l’arithmétique des intervalles.
56 C
Contrairement à RealPaver, Choco est foncièrement et intégralement orienté objet. Ainsi, une hiérar-
chie de contraintes abstraites les répartit selon leur genre —globale, arithmétique, booléenne, etc.— et
le type des variables sur lesquelles elles portent —entière, réelle, booléenne, ensembliste, etc. Idem pour
tout ce qui concerne les stratégies d’exploration et toutes les structures de données haut-niveau du solveur
qui sont représentées par une hiérarchie ad hoc. En outre, un mécanisme basé sur la notion d’événement
déclenche automatiquement la méthode adéquate comme par exemple lors d’un changement de borne
inférieure qui entraîne la réévaluation de toutes les contraintes utilisant la variable dont le domaine a été
modifié.
des réductions de domaines pendant l’exploration de l’espace de recherche par Choco, après par exemple
qu’une ou des variables sur lesquelles ce propagateur porte aient vu leurs domaines respectifs réduits
par l’application du filtrage d’une autre contrainte Choco, la méthode de calcul de point du fixe du
propagateur est appelée avec en paramètre la liste des domaines —tel qu’actuellement réduits— des
variables sur lesquelles il porte. Ces domaines sont recopiés dans RealPaver avant d’être réduits grâce
aux algorithmes de filtrage de ce dernier. Enfin, quand le point-fixe a été atteint, les nouveaux domaines
sont retransmis à Choco afin de poursuivre la résolution au point où elle en était rendue, à savoir ou
bien mener à son terme la propagation des réductions de domaine, ou bien choisir une variable pour
énumérer/bissecter son domaine et ainsi poursuivre l’exploration de l’espace de recherche.
P. 4.13 : Déclaration en Java des fonctions C++ appelables via JNI.
l’essentiel du travail à fournir —du point de vue de Choco— pour s’insérer correctement dans le proces-
sus de résolution. Elle est dotée d’une méthode propagate dans laquelle on va concrètement faire appel
à RealPaver via des méthodes C++ dédiées que nous avons implémentées.
Mis à part cela, il suffit d’étendre le langage de modélisation de façon à pouvoir faire appel à une mé-
thode d’ajout de contrainte au propagateur et faire le lien, à l’aller comme au retour, entre les domaines
traités par Choco et ceux retournés par RealPaver.
P. 4.14 : Extrait d’une fonction C++ faisant appel à des éléments écrits en Java.
• l. 1-2 : entête de la fonction, généré automatiquement par l’exécutable jni à partir d’un
fichier java faisant usage de la bibliothèque (cf Prog. 4.13)
• l. 5 : récupération du premier élément du tableau vars
• l. 6 : instanciation d’un objet représentant la classe de cet objet
• l. 7 : instanciation d’un objet représentant la méthode toString de cette classe
• l. 9 : boucle de traitement de chacun des éléments du tableau vars
• l. 11 : récupération du i-ème élément du tableau vars
• l. 12 : appel à la méthode toString de cet élément et stockage dans une chaîne de la valeur
retournée
4.3 Expériences
Dans cette section nous décrivons les expériences que nous avons réalisées pour évaluer le compor-
tement et les performances de la collaboration que nous avons mise en oeuvre.
MINLPlib
MINLPlib est une collection de problèmes d’optimisation mixtes, c’est-à-dire qui font intervenir à la
fois des variables discrètes —parmi lesquelles certaines peuvent être booléennes et d’autres entières— et
des variables continues [BDM03]. Une petite moitié des problèmes que nous utilisons dans cette section
provient de cette bibliothèque :
• synthes3 : un problème simplifié de synthèse de procédé qui vise à déterminer la structure opti-
male et les conditions d’opérations d’un procédé qui doit satisfaire des spécifications de conception
données [DG86] ;
• st_miqp5 : un problème de programmation mixte quadratique indéfinie MIQP [TS02] ;
• ex1263 : un problème d’optimisation de découpe de papier cherchant à optimiser le profit en
minimisant les pertes de matière première [HWPS98].
Placement d’antennes
Un autre problème consiste à placer des antennes sur une surface rectangulaire de façon à maximi-
ser la distance entre les antennes, étant données des contraintes de distance minimimum et maximum
possibles entre deux antennes et matérialisant des limitations techniques8 .
En outre, nous avons considéré une contrainte supplémentaire qui consiste à devoir placer certaines
antennes à des positions discrètes et les autres à des positions quelconques. Nous testons deux instances
de ce problème pour sa version à 4 antennes, l’une où une seule des antennes doit être placée à des
coordonnées entières, l’autre où trois des antennes le doivent.
Fractions
Nous avons également testé le comportement de la collaboration sur des instances d’un problème
faisant intervenir non seulement des contraintes arithmétiques mais aussi des contraintes globales dis-
crètes. Il s’agit d’une variante du problème classique de Y qui consiste à résoudre l’équation
suivante :
A D G
+ + =1
BC EF HI
où chaque lettre est associée à un digit distinct et différent de zéro. Nous avons également étendu ce
problème en une deuxième instance, plus difficile à résoudre, en ajoutant des variables J, K, L et un
J
quatrième terme KL à l’équation. En outre, afin de garder le problème consistant, nous avons étendu
les domaines initiaux jusqu’à 12 au lieu de 9 i.e. les digits s’expriment en base 13 au lieu de la base 10.
Notons enfin que pour chacune des instances nous avons également posé des contraintes de façon à briser
A D A
les symétries, e.g. BC ≤ EF et 4 BC ≥ 1.
4.3.2 Résultats
Les résultats expérimentaux sont présentés dans la Table 4.2. Tous les calculs ont été réalisés sur
un Intel Xeon 3 Ghz avec 16 Go de RAM. Nous avons comparé la collaboration des solveurs Choco
1.2.06 et RealPaver 1.0 avec les outils RealPaver 1.0 seul, Choco 1.2.06 seul et ECLi PSe 5.10. Chaque
colonne correspond à un outil différent RP, COOP, Choco et ECLi PSe . Une colonne #var renseigne pour
chaque instance le rapport du nombre de variables discrètes sur le nombre de variables continues. Les
8
Les contraintes d’allocation de fréquences ne sont pas prises en compte ici.
4.3. Expériences 61
temps de résolution en secondes sont donnés sur la première ligne et le nombre de branchements sur la
deuxième, exprimés en milliers pour les variables discrètes et pour les variables continues. Un nombre
de branchements valant ∅ signifie que ce nombre est sans objet pour le problème et le type de variable
donnés. Un point d’interrogation marque un temps de résolution dépassant les 10000 secondes. Une barre
horizontale représente une impossibilité de modéliser le problème avec l’outil donné.
Constatons pour commencer que la collaboration des solveurs est la plus efficace pour traiter les
instances du problème fractions, à savoir les problèmes faisant intervenir à la fois des contraintes
arithmétiques nécessitant d’être évaluées grâce à l’arithmétique des intervalles, déléguées à RealPaver,
et des contraintes globales traitées par les algorithmes dédiés et state-of-the-art implémentés dans Choco,
pour filtrer ici le alldifferent. Ce type de problème mixte est donc logiquement celui qui est à même
de profiter le plus de ce type de collaboration de solveurs et matérialise le mieux le bien-fondé de l’idée
directrice de ces travaux consistant à utiliser les forces de chacun des outils pour combler les faiblesses
de l’autre. En pratique, on peut observer ici que la collaboration est quinze fois plus rapide que RealPaver
seul, ce qui s’explique aisément par le fait que ce dernier utilise une clique de diségalités binaires pour
représenter le alldifferent et perd ainsi le double avantage d’une information globale à toutes les
variables traitée qui plus est de façon optimale par les algorithmes de filtrage de Choco. En outre, la
collaboration est huit fois plus rapide que Choco seul, qui pèche de par sa capacité nettement moins
bonne à utiliser les contraintes arithmétiques pour filtrer les domaines et doit donc avoir recours à un
nombre de branchements deux fois plus élevé. Enfin, remarquons qu’elle est quatre fois plus rapide
qu’ECLi PSe et que cela semble dû à l’efficacité de RealPaver pour le calcul sur les intervalles, puisqu’on
observe le même nombre de branchements mais pour un temps de résolution bien supérieur.
En outre, la coopération que nous avons mise en oeuvre est largement plus rapide que Choco seul
dans la très grande majorité des cas. En effet, RealPaver est nettement plus efficace que Choco pour
traiter les contraintes arithmétiques, ce qui peut très certainement s’expliquer par le fait que RealPaver
antennes_41 2/6 ? ? ? ?
260.9 381.3 1862.4 839.41
antennes_43 6/2 274.4 275.1 174.6 146.0
102.8 108.3 66.1 186.7
2.0 0.91 2.5 1.42
fractions_10 9/0 6.7 1.3 2.0 0.86
∅ ∅ ∅ ∅
3351.6 220.9 1707.0 848.2
fractions_13 12/0 667.2 484.7 1033.7 485.0
∅ ∅ ∅ ∅
Table 4.2 – Temps de résolution en secondes (première ligne) et nombre de branchement sur entier/réel
en milliers (deuxième ligne).
62 C
a été conçu dans ce but, autour d’un cœur optimisé en C et dédié au calcul sur les intervalles à bornes
flottantes, alors que Choco, codé en Java, a simplement été étendu après sa conception. L’exception
à cette tendance générale, qu’on observe pour ex1263, tient au nombre élevé de variables utilisées et
montre une faiblesse de notre modèle de collaboration : le temps passé à recopier les domaines d’un
solveur à l’autre solveur. Pour cet exemple, on a en effet relevé un temps de transfert s’élevant à plus de
50 secondes. De plus, on voit bien dans cet exemple notamment que Choco, moins efficace pour réduire
les domaines, doit effectuer plus de branchements.
À l’opposé, RealPaver est plus rapide que la collaboration dans la majorité des cas —toutes les
instances des problèmes qui ne font pas intervenir de contraintes globales. Le gain procuré par la stratégie
d’exploration de Choco, énumérant intégralement les domaines discrets, n’est pas compensé par les
temps de transfert des domaines. Mais ces temps de transfert, s’il est indéniable qu’ils participent à la
différence observée, n’expliquent pas toute la différence et on doit aussi considérer le fait que le point
fixe est calculé dans la coopération sur toutes les contraintes à chaque fois, à chaque modification d’une
variable, alors que dans RealPaver ce rappel des contraintes est plus souple et seules les contraintes
impliquant les variables modifiées sont propagées. Ainsi la différence est expliquée à la fois par les
temps de transfert des domaines mais aussi par le caractère plus lourd, moins fin, du calcul du point fixe
par Choco à l’aide de la contrainte utilisant RealPaver.
Notons pour finir que ECLi PSe apparaît ici bien plus performant sur les problèmes avec un rapport
élevé du nombre de variables discrètes sur le nombre de variables continues, comme le montre ex1263
par opposition à st_miqp5 par exemple : dans une situation du genre du premier cas, le solveur est très
compétitif pour résoudre un problème où les trois quarts des variables sont discrètes —binaires qui plus
est— alors que dans le second cas, il est largement dépassé par au moins l’un des autres outils, sinon une
majorité d’entre eux.
4.4 Conclusion
Dans ce chapitre nous avons présenté nos travaux concernant la réalisation d’une collaboration entre
deux solveurs, l’un dédié aux domaines discrets, l’autre dédié aux domaines continus. Après avoir étudié
en détails les capacités de modélisation en termes de problèmes mixtes des principaux outils non com-
merciaux disponibles actuellement, nous avons entrepris d’utiliser les capacités de propagation sur les
contraintes arithmétiques d’un solveur continu pour améliorer les capacités dans ce domaine d’un solveur
discret. Nous avons alors pu mesurer l’efficacité de la collaboration que nous avons mise en oeuvre en
la comparant à différents outils vis-à-vis de la résolution d’une sélection de problèmes mixtes. À travers
ces travaux, nos principales contributions sont les suivantes :
• Nous avons mis au point une plateforme logicielle doté d’un langage de modélisation facilitant
l’expression des problèmes mixtes en fusionnant les hautes capacités d’expression de contraintes
arithmétiques du langage de modélisation dédié d’un solveur continu, avec les nombreuses primi-
tives de modélisation de contraintes globales d’un solveur discret, améliorant en cela les capacités
de chacun des deux outils pris individuellement.
• Nous avons également mis en évidence la nette supériorité en termes d’efficacité de résolution
d’une collaboration d’outils pour résoudre les problèmes mixtes qui font intervenir à la fois des
contraintes arithmétiques, nécessitant une évaluation par l’arithmétique des intervalles, et des con-
traintes globales, requérant des algorithmes de filtrage dédiés et optimisés.
4.4. Conclusion 63
En outre, des travaux futurs pourraient permettre de profiter encore davantage des forces d’une telle
collaboration en tentant, sinon de s’en affranchir, au moins d’en améliorer les points faibles majeurs que
sont actuellement les temps de recopie de domaines d’un solveur à l’autre et le manque de souplesse
du propagateur utilisé pour calculer un point-fixe des domaines à partir des contraintes arithmétiques du
modèle. En ce qui concerne le premier point, on peut en effet concevoir une symbiose plus profonde
entre les deux outils, qui pourrait amener par exemple à proposer une gestion plus intelligente de la
recopie des domaines, au moins dans le cas particulier où il n’y a pu se produire aucune réduction et où
la recopie en retour est inutile. Quant au deuxième point, on peut également envisager une collaboration
plus en profondeur qui permettrait une sélection plus fine des opérateurs de réduction à rappeler en cas
de modification de domaines, par exemple en offrant la possibilité de n’effectuer le calcul de point-fixe
que sur un sous-ensemble des contraintes arithmétiques déléguées au solveur continu.
V
U TILISATION DE LA CONTRAINTE
DISCRÈTE A LL D IFFERENT DANS
UN SOLVEUR CONTINU
5.1 Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.2 Modéliser des MINLP avec la contrainte AllDifferent . . . . . . . . . . . . . . . . . . 66
5.2.1 Exemples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.2.2 Synthèse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.3 Traiter le AllDifferent dans un solveur continu . . . . . . . . . . . . . . . . . . . . . . 71
5.3.1 Reformulations de la contrainte . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.3.2 Implémentation de filtrages dédiés . . . . . . . . . . . . . . . . . . . . . . . . 72
5.4 Expériences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.4.1 Description des problèmes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.4.2 Résultats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
Ce chapitre présente une étude expérimentale investigant le champ des techniques de programmation
par contraintes dédiées à la résolution de problèmes mixtes, un champ où peu de travaux ont déjà été
réalisés [Gel98], [GF03]. Dans ce chapitre, nous nous intéressons à la recherche d’une manière simple
de permettre aux solveurs continus de gérer efficacement les contraintes globales et de tirer avantage de
cette force pour résoudre des problèmes mixtes. En particulier, nous étudions la possibilité d’utiliser la
contrainte alldifferent pour modéliser des problèmes mixtes du monde réel et les résoudre dans un
solveur continu. Cette contrainte est en effet l’une des plus utilisées et des mieux étudiées.
Le chapitre est organisé comme suit. D’abord, nous présentons ce qui nous motive à aller dans cette
direction. Ensuite, nous montrons comment utiliser la contrainte alldifferent pour la modélisation
de quelques problèmes mixtes du monde réel. Nous montrons alors les solutions possibles pour la for-
mulation de la contrainte et son filtrage dans un solveur continu. Enfin, nous présentons les résultats des
expériences que nous avons réalisées pour comparer leur comportement, avant de conclure le chapitre.
65
66 U AD
5.1 Motivations
Les problèmes d’optimisation non linéaire mixtes-entiers (MINLP) sont des problèmes d’optimisa-
tion sous contraintes qui font intervenir à la fois des variables discrètes et des variables continues, dotés
d’une fonction-objectif éventuellement non convexe et de contraintes qui peuvent être non linéaires. Des
problèmes de ce genre apparaissent dans de nombreux domaines d’application comme par exemple la
bioinformatique, l’optimisation de portefeuilles d’actions et le design en ingénierie. La plupart de ces
méthodes de résolution sont basées sur le cadre classique du branch-and-bound et ont recours à des
techniques de différentes sortes, comme les relaxations linéaires, les méthodes de coupes ou encore la
reformulation symbolique.
Récemment, une tendance majeure a émergé dans le champ de la recherche opérationnelle, la conver-
gence avec la programmation par contraintes. Elle est dédiée à la recherche de moyens pour la première
de tirer avantage de méthodes issues de la seconde [Hoo07]. En fait, il se trouve qu’intégrer la program-
mation par contraintes et la recherche opérationnelle est une tendance très naturelle puisqu’elles ont des
buts très similaires. L’utilisation de l’arithmétique des intervalles en programmation par contraintes est
par exemple un moyen puissant de résoudre les MINLP de façon rigoureuse et complète, i.e. de les ré-
soudre avec les deux intéressantes caractéristiques de calculer des solutions garanties et de ne perdre
aucune solution. Cela devient très intéressant quand on considère le fait qu’un outil de dernière généra-
tion bien connu comme BARON n’est pas digne de confiance en général [BSG∗ 10], [LMR07].
Mais ce qui nous intéresse particulièrement ici est l’utilisation de contraintes globales en recherche
opérationnelle. Les contraintes globales sont de puissantes primitives de modélisation qui sont très lar-
gement utilisées dans le domaine de la programmation par contraintes. Les utilisateurs disposent d’un
catalogue de plusieurs centaines de contraintes préconstruites parmi laquelle ils n’ont qu’à choisir celle
qui convient le mieux au problème qu’ils veulent modéliser [BCDP07]. En outre, chaque contrainte glo-
bale a son propre algorithme de filtrage dédié de façon à ce que le processus de résolution s’adapte le
mieux au problème modélisé et rejette facilement de larges parties de l’espace de recherche. L’utilisation
de contraintes globales en recherche opérationnelle a ainsi deux grands avantages [Hoo07] :
• rendre les modèles plus faciles à construire et à mettre à jour ;
• révéler la structure du problème au solveur pour permettre de le résoudre plus efficacement.
5.2.1 Exemples
Ordonnancement robuste à une machine
Pour commencer, nous présentons la formulation en programmation mixte-entière de M
pour le problème d’ordonnancement robuste à une machine avec des données intervalles qu’on peut
5.2. Modéliser des MINLP avec la contrainte AllDifferent 67
trouver dans [Mon07], où l’on considère pour chaque tâche un intervalle de temps de traitements pos-
sibles. L’optimisation est réalisée selon un critère de robustesse, i.e. on cherche l’ordonnancement qui se
comporte le mieux dans le pire des cas.
Chacune de N tâches données doit se voir affecter un ordre de traitement distinct, ce qui en pro-
grammation non linéaire mixte-entière est normalement modélisé en utilisant des variables binaires, e.g.
xik est vraie quand la tâche i est traitée en kème position, et sous la contrainte qu’à chaque tâche (resp.
position) soit affectée exactement une position (resp. tâche) :
N
X
xik = 1,
i=1
XN
xik = 1.
k=1
alldifferent([t1 , . . . , tN ]).
Le modèle que nous fournissons ici pour le problème d’ordonnancement robuste résulte d’une tech-
nique de reformulation bien connue qui exploite la propriété d’unimodularité du sous-problème de maxi-
misation [Law76] :
N
X N
X
min ηi + τk
i=1 k=1
s.c.
(1) alldifferent(t1 , . . . , tN ),
X N
(2) ti = xik · k,
k=1
N
X
(3) xik = 1,
k=1
k
X N
X
(4) ηi + τk ≥ pi (k − j)xi j + pi (k − j)xi j ,
j=1 j=k
xik ∈ {0, 1},
ti ∈ {1, . . . , N},
ηi ∈ R,
τk ∈ R.
La contrainte (1) exprime que chaque tâche doit avoir un ordre de traitement différent ; (2) transmet
cette valeur aux variables binaires, ce que requiert la contrainte (4). La contrainte (3) empêche cette
transmission d’impliquer deux variables binaires xik simultanément vraies pour un même i. Enfin, (4)
fait simplement le lien entre les variables dont on maximise la somme et l’expression arithmétique du
coût de regret pour le pire des scénarios associé à cet ordonnancement. De plus amples détails sont
disponibles dans [Mon07].
68 U AD
Ici aussi, étant donné un ensemble de N variables entières ni prenant leur valeur dans {1, . . . , N}, les
deux contraintes ci-dessus peuvent être reformulées en utilisant la contrainte globale alldifferent :
alldifferent([n1 , . . . , nN ]).
s.c.
(1) alldifferent(n
X 1 , . . . , nN ),
(2) ni = Zi j · j,
X j
(3) Zi j = 1,
j
(4) di T c ≤ G
X Gi αβX
i γi T Pi +X
i
i
(1 − e−βi T PXi ),
X
(5) Tc = T Pi + τi j Zi j + sli j ,
i i j i j
(6) −(1 − Zi j )T c ≤ T S j − (T S i + T Pi + τi j + sli j ),
(7) (1 − Zi j )T c ≥ T S j − (T S i + T Pi + τi j + sli j ),
T ci , T Pi , T S i , di , τi j , sli j ∈ R.
ni ∈ {1, . . . , N},
Zi j ∈ {0, 1},
La fonction objectif à optimiser est composée de la somme des coûts par cycle du stockage, des
matières premières et de l’entretien. Comme dans l’exemple précédent, la contrainte (1) exprime que
chaque tâche doit avoir un ordre de traitement distinct, (2) transmet cette valeur aux variables binaires
comme le nécessite la contrainte (5). La contrainte (3) empêche cette transmission de faire intervenir
5.2. Modéliser des MINLP avec la contrainte AllDifferent 69
deux variables Zi j simultanément vraies pour un i donné. La contrainte (4) exprime que pour chaque
produit, on impose que la quantité produite soit supérieure à la demande pendant le cycle, tandis que
la contrainte (5) impose que la durée du cycle soit supérieure à la somme des temps de traitement et
de transition entre produits. Pour finir, les contraintes (6) et (7) déclarent que chaque produit ne peut
commencer à être produit qu’après le précédent traitement suivi d’une réinitialisation de la chaîne de
production pour le bon produit. De plus amples détails sont disponibles dans [CLPP05].
4 → 7 → 3 → 1
2 → 8 → 10 → 12
5 → 11 → 6 → 9
Ci-dessus on peut voir une trajectoire possible, d’une durée de 4 cycles, pour un réacteur de 12
noeuds et dans lequel 3 barres de combustible —soit un quart du nombre total de barres— sont retirées
du cœur à chaque fin de cycle. Chaque nombre y réfère à un noeud différent du réacteur. Selon cette
matrice, les barres placées aux noeuds 1, 12 et 9 sont retirées à chaque rechargement pendant que des
barres fraîches sont insérées dans les noeuds 4, 2 et 5. La barre au noeud 7 est déplacée au noeud 3, la
barre au noeud 10 est déplacée au noeud 12, etc.
Comme précédemment, l’approche MINLP classique pour modéliser l’affectation d’un noeud à une
cellule distincte de la matrice et réciproquement est d’utiliser des variables binaires : xi,l,m est valuée à
vrai lorsque que la barre au noeud i est affectée à la colonne l et à la ligne m dans la matrice de trajectoire.
On doit ensuite contraindre ces variables de manière à ce qu’une exactement soit mise à vrai pour un i
donné puisque chaque noeud doit apparaître dans la matrice une fois et une seule. Réciproquement, une
variable et une seule doit être mise à vrai pour une position donné (l, m) de la matrice puisque chaque
cellule doit se voir affecter exactement un noeud. Ainsi on a :
N
X
xi,l,m = 1,
i=1
X L X M
xi,l,m = 1.
l=1 m=1
Il s’avère que le problème d’affectation des noeuds du réacteur aux cellules de la matrice de trajec-
toire peut s’exprimer au moyen d’une contrainte alldifferent. En effet, considérons la numérotation
suivante des cellules :
70 U AD
1 2 3 4
5 6 7 8
9 10 11 12
Ainsi, à chaque noeud i doit être affecté un numéro distinct ni , compris dans {1, . . . , 12}, ce qui une
fois encore correspond parfaitement avec la signification de la contrainte alldifferent. La contrainte
d’affectation de la matrice de trajectoire est donc simplement modélisée comme suit :
alldifferent([n1 , . . . , nN ]).
Nous donnons à présent un aperçu du modèle complet du problème de motif optimal de rechargement
du réacteur nucléaire où chaque cycle est divisé en un nombre discret d’étapes :
ef f
max kT
x,k∞ ,R,ke f f
s.c.
(1) alldifferent(n1 , . . . , nN ),
X
(2) ni = xi,l,m · (l + (m − 1) · L),
XX l,m
(3) xi,l,m = 1,
l m
f lim
(4) ∞R
ki,t i,t ≤ X N X X
(5) ∞
ki,1 = xi,l,m x j,l−1,m k∞
j,T
l>1Xm j
+( xi,1,m )k f resh ,
m
(6) X ki,t+1 = ki,t − (αPc ∆t )ki,t Ri,t ,
∞ ∞ ∞
∞
(7) ki,t Ri,t = 1
i
ef f
X
(8) kT Ri, j = Gi, j k∞
j,T R j,t ,
j
ef f
kt ≥ 0,
∞
ki,t ≥ 0,
Ri,t ≥ 0,
ni ∈ {1, . . . , N},
xi,l,m ∈ {0, 1},
Les contraintes (1) à (3) modélisent le problème d’affectation de la matrice de trajectoire : (1) ex-
prime que chaque noeud doit y avoir une position distincte tandis que (2) transfère cette position aux
variables binaires pour que la contrainte (5) puisse accéder à cette information et sache quels noeuds se
sont vus affecter une barre fraîche ; (3) empêche ce transfert d’impliquer deux variables xi,l,m mises à vrai
pour un même i donné. Les contraintes (4) à (7) décrivent le comportement des barres de combustible
en termes de lois physiques : (4) requiert que la puissance maximum par noeud ne soit pas trop grande,
pour des raisons évidentes de sécurité ; (5) définit la valeur initiale du coefficient de multiplication de
puissance de chaque noeud selon qu’il contient ou non une barre fraîche ; (6) décrit la décroissance de
la réactivité locale ; (7) sert à normaliser la puissance en chaque noeud. Enfin, la contrainte (8) exprime
la réactivité moyenne du cœur en fonction de la réactivité locale à chaque noeud et de la probabilité Gi j
qu’un neutron émis au noeud i soit absorbé au noeud j.
5.3. Traiter le AllDifferent dans un solveur continu 71
5.2.2 Synthèse
Ces trois exemples nous ont permis de bien comprendre la façon dont on peut utiliser la contrainte
discrète alldifferent pour modéliser des problèmes MINLP issus du monde réel. L’idée est d’être ca-
pable de répérer si une partie du modèle ne correspond pas à un problème d’affectation et/ou d’identifier,
dans les modèles de problèmes de type ordonnancement/affectation sous contraintes, s’il est possible de
trouver une partie avec la forme générique suivante :
Ni
1
Ni p
X X
xi1 ,...,i p , j1 ,..., jq = 1,
...
i1 =1 i p =1
N j1 N jq
X X
xi1 ,...,i p , j1 ,..., jq = 1,
...
j1 =1 jq =1
Une variable yi j valant vrai signifie que la ième variable prend pour valeur la jème valeur de son
domaine. De nombreux problèmes d’affectation utilisent ce genre de modélisation. Quand il y a autant
de valeurs dans le domaine que de variables à contraindre par le alldifferent, on remplace par un
égal l’inférieur-ou-égal dans la seconde équation :
Xn
yi j = 1, i = 1, . . . , n
j=1
Xn
yi j = 1, j = 1, . . . , n
i=1
Pour finir, remarquons qu’il est possible au besoin de faire le lien pour un i donné entre les variables
booléennes yi j et la valeur entière affectée à la variable xi . Considérant par exemple xi ∈ {v1 , . . . , vn }, on
a :
m
X
xi = v j yi j , i = 1, . . . , n
j=1
Heureusement, on peut savoir quelles sont les contraintes de la relaxation qui sont le plus susceptibles
de pouvoir filtrer les domaines étant donné un ensemble de domaines à réduire :
k k
X k(k + 1) X k(k + 1)
x k N tel que xi <
i ≥ ∀ <
2 2
i=1
i=1
N N
k(k + 1) N(N + 1) k(k − 1)
X X
x k 1 tel que xi >
≤ ∀ > −
i
2 2 2
i=k i=k
à condition de classer les variables en ordre croissant (resp. décroissant) de la borne inférieure (resp.
supérieure) de leur domaine i.e. xi ≤ x j ∀i ≤ j (resp. xi ≥ x j ∀i ≥ j ).
La formule ci-dessus s’implémente facilement en un algorithme de filtrage de bornes, en se servant
de la première ligne pour filtrer les bornes inférieures et de la seconde pour filtrer les bornes supérieures
1 B′ ← sort(B)
2 sum ← 0
3 pour k ← 1 à |B′ | − 1 faire
4 si |B′k | = |B′k+1 | = 1 et B′k = B′k+1 alors
5 retourner f alse
6 fin
7 sum ← sum + B′k
k(k+1)
8 si sum < alors
2
k
X k(k + 1)
B′ ← reduce B′ , xi ≥
9 i=1
2
10 si empty(B′ ) alors
11 retourner f alse
12 fin
13 fin
14 fin
15 retourner true
des domaines des variables. Nous ne donnons ici que l’algorithme pour filtrer les bornes inférieures —cf
Algo. 5.1— puisque l’autre est trivial à obtenir par symétrie. Réduire un pavé consiste alors à utiliser en
même temps l’égalité (1) de la formule de la relaxation, et les inégalités de la forme ci-dessus, que nous
qualifierons d’ad hoc, comme le montre l’Algorithme 5.2.
Pour finir, signalons que nous voulons également tester expérimentalement une deuxième façon d’uti-
liser cette relaxation dont nous venons de décrire l’usage communément pratiqué. Dans la section sui-
vante, nous comparerons deux versions du filtrage :
• en utilisant l’égalité (1) et les inégalités ad hoc ;
• en utilisant l’égalité (1) seule.
# D x1 ∪ · · · ∪ D xn ≥ n
où D xi est le domaine de xi pour tout i. Ceci est en fait une condition nécessaire bien connue, basée sur
la notion d’intervalle de H —notion qu’on détaillera davantage dans la sous-section suivante— en
particulier l’intervalle de H associé à tout l’ensemble de variables [Hal35].
L’algorithme de filtrage dédié que nous avons implémenté est présenté comme l’Algorithme 5.3. En
premier lieu, si l’opérateur de cardinalité de domaine calcule un entier inférieur à n le nombre de va-
riables, alors la contrainte ne tient pas et retourne une inconsistance. Deuxièmement, quand une variable
est instanciée à une valeur, cette valeur est retirée des domaines de toutes les variables pour lesquelles il
s’agit d’une valeur en borne de domaine. Cette technique correspond d’ailleurs au traitement classique
d’une contrainte de diségalité : étant donnés x , y et x = a alors a est retiré du domaine de y.
5.3. Traiter le AllDifferent dans un solveur continu 75
Nous pouvons également implémenter un algorithme dédié de filtrage par consistance de bornes
comme celui donné dans [LOQTvB03], ou bien encore un algorithme classique de matching comme
celui décrit dans [Rég94]. Par ailleurs, dans le cadre d’une étude expérimentale utilisant un solveur
travaillant sur des intervalles, nous nous focalisons sur le cas des domaines intervalles d’entiers sans
trous. Puisque nous savons déjà que l’algorithme de matching de Ŕ n’est pas le plus efficace sur ce
type de domaines [vH01], nous nous intéresserons plutôt à l’algorithme de consistance de bornes donné
dans [LOQTvB03], i.e. la version la plus récente et la plus efficace des évolutions de la version initiale
de l’algorithme bien connu de P [Pug98].
L’idée clé sur laquelle se base cet algorithme est la recherche d’intervalles de H, c’est-à-dire
des intervalles dont la taille est égale au nombre de variables dont ils peuvent contenir le domaine. Par
exemple, étant données quatre variables contraintes par une contrainte alldifferent, e.g. x1 ∈ {1, 2},
x2 ∈ {2, 3}, x3 ∈ {1, 2, 3} et x4 ∈ {1, . . . , 5}, l’intervalle d’entiers IH = {1, 2, 3} est un intervalle de H
76 U AD
1 pour i ← 1 à nb + 1 faire
2 t[i] ← h[i] ← i − 1
3 d[i] ← bounds[i] − bounds[i − 1]
4 fin
5 pour i ← 0 à niv − 1 faire
6 x ← maxsorted[i].minrank
7 y ← maxsorted[i].maxrank
8 z ← pathmax(t, x + 1)
9 j ← t[z]
10 d[z] ← d[z] − 1
11 si d[z] = 0 alors
12 t[z] ← z + 1
13 z ← pathmax(t, t[z])
14 t[z] ← j
15 fin
16 pathset(t, x + 1, z, z)
17 si d[z] < bounds[z] − bounds[y] alors
18 retourner f alse
19 fin
20 si h[x] > x alors
21 w ← pathmax(h, h[x])
22 Bk ← maxsorted[i].min ← bounds[w]
23 pathset(h, x, w, w)
24 fin
25 si d[z] = bounds[z] − bounds[y] alors
26 pathset(h, h[y], j − 1, y)
27 h[y] ← j − 1
28 fin
29 fin
30 retourner true
par rapport à l’ensemble de variables {x1 , x2 , x3 } et il est impossible pour la variable x4 de se voir affecter
une des valeurs contenue dans IH : il n’y a en effet pas d’autre choix que celui d’affecter ces valeurs à
x1 , x2 et x3 —peu importe quelle variable prend quelle valeur— et ces valeurs peuvent alors être retirées
du domaine initial de x4 qui est ainsi réduit à {4, 5}.
L’algorithme complet pour le filtrage des bornes inférieures est donné ici bien qu’il soit un peu
difficile à comprendre de prime abord. Il s’agit de l’Algorithme 5.4 et il se dérive symétriquement en
un algorithme de filtrage des bornes supérieures que par contre nous ne donnerons pas ici. Ces deux
algorithmes font usage d’un certain nombre de structures de données dédiées, triées, qui ont déjà été
optimisées plusieurs fois depuis la première parution de l’algorithme dans [Pug98]. De plus, nous n’ex-
pliquerons pas son fonctionnement en détails ici car cela prendrait une place importante et là n’est pas
le but de cette section. Cela nous paraît par contre important de montrer cet algorithme ici, de façon à
ce que chacun puisse le comparer avec les précédents algorithmes présentés dans cette section et puisse
réaliser qu’il a un degré beaucoup plus élevé de complexité dans son principe, en particulier quand on
le compare à l’Algorithme 5.3. En conséquence, des lecteurs intéressés qui auraient besoin de plus d’ex-
plications ne devront pas hésiter à se référer à l’abondante littérature couvrant ce sujet, à commencer
par [Pug98], [LOQTvB03] et [vH01].
5.4 Expériences
Notre but est d’étudier l’impact sur la résolution du choix d’une modélisation et d’une méthode
de filtrage pour le alldifferent. Comme on l’a déjà mentionné à plusieurs reprises, peu ou pas de
travaux ont déjà été réalisés dans ce domaine de la programmation par contraintes pour les problèmes
mixtes, faisant intervenir à la fois des contraintes globales, discrètes, et des contraintes non linéaires,
continues.
La Table 5.1 résume les caractéristiques principales des instances que nous avons résolues de chacun
de ces types de problèmes.
Table 5.1 – Pour chaque problème, #var donne le nombre variables discrètes/continues dans le modèle
utilisant la contrainte alldifferent , et #sol représente le nombre de solutions au problème.
puisqu’en effet ils font intervenir à la fois des contraintes globales et des contraintes arithmétiques né-
cessitant une évaluation aux intervalles. Cette sous-section présente brièvement les problèmes que nous
avons utilisés à cet effet.
Sevvoth
Le problème Sevvoth [Kei] est un problème du type SEND+MORE=MONEY, dans lequel les mots à
additionner sont ceux du poème ci-après, le terme résultant étant le mot SEVVOTH 2 lui-même. Comme
dans tous les problèmes de ce genre, un chiffre est associé à chaque lettre et tous doivent être distincts les
uns des autres —ce qu’on modélise bien sûr au moyen d’une contrainte alldifferent. On ajoute alors
la contrainte arithmétique qui représente la somme globale, chaque mot MOT générant un terme comme
M ∗ 100 + O ∗ 10 + T . Le poème utilisé est le suivant :
2
Il s’agit du nom d’une île mythique quelque part en Mer du Nord.
5.4. Expériences 79
Fractions
Le problème fractions est une extension du problème classique de Y qui vise à résoudre
l’équation suivante :
A D G
+ + =1
BC EF HI
où chaque lettre est associée à un digit distinct et différent de zéro. Pour compliquer le problème, nous
J
avons comme au chapitre précédent ajouté trois variables J, K, L et un quatrième terme KL à l’équation.
Pour garder le problème consistant, nous avons également dû étendre les domaines initiaux de façon à ce
que chaque digit s’exprime en base 13. Nous avons également imposé les contraintes suivantes, de façon
à casser les symétries et limiter le nombre de solutions :
A
≤ D ≤ G ≤ J ,
BC A EF HI KL
4∗ ≥ 1,
4 ∗ BC
J
≤ 1.
KL
Carré magique
Le problème du carré magique est un puzzle mathématique bien connu où l’on doit remplir une
grille de dimensions n × n avec les n2 premiers nombres entiers de façon à ce que la somme des cases de
chaque ligne, colonne ou diagonale soit égale à n2 . Brièvement, les variables du problèmes sont bien sûr
les n2 cases de la grille, de domaine initial {1, . . . , n2 }, soumises à une contrainte alldifferent pour
qu’elles prennent toutes des valeurs deux à deux distinctes ; les contraintes restantes sont des contraintes
arithmétiques matérialisant mathématiquement le fait que la somme des cases de chaque ligne, colonne
ou dialogue doit être égale à n2 .
N reines
Pour finir, le problème des N reines est le casse-tête bien connu qui vise comme sont nom l’indique
à placer n reines sur un échiquier de dimensions n × n de façon à satisfaire la contrainte que deux reines
ne peuvent pas être prises l’une par l’autre, i.e. toute paire de reines ne peut pas être sur une seule ligne,
une seule colonne ou une seule diagonale. Une modélisation possible consiste à affecter à chaque reine
une variable Ri représentant sa ligne, la reine Ri se trouvant dans la ième colonne. On ajoute alors trois
contraintes alldifferent :
• une sur Ri , ∀i : deux reines ne peuvent pas être sur la même ligne ;
• une sur Ri + i, ∀i : deux reines ne peuvent pas être sur la même diagonale NO-SE ;
• une sur Ri − i, ∀i : deux reines ne peuvent pas être sur la même diagonale NE-SO.
Le premier de ces problèmes, int-frac, est un problème purement discret doté d’une structure inspirée
du problème fractions :
alldifferent(x1 , . . . , xn ),
x1 xn−1 1 n−1
+ ... + = + ... + ,
x2 xn 2 n
xi + 1 ≤ xi+2 , 1 ≤ i ≤ n, i odd,
xi ∈ {1, . . . , n}, 1 ≤ i ≤ n.
Le deuxième problème, mixed-sin, est un problème mixte qui fait appel à la fois à des contraintes
linéaires et non linéaires :
alldifferent(p
X 1 , . . . , pn ),
xi pi + x j = i, 1 ≤ i ≤ n,
j,i
sin(πxi ) = 0, 1 ≤ i ≤ n,
pi ∈ {1, . . . , n}, 1 ≤ i ≤ n,
xi ∈ [−100, 100], 1 ≤ i ≤ n.
Quant au troisième problème, hard-mixed-sin, il s’agit d’un problème mixte inspiré du précédent
mais dans lequel on a durci les contraintes non linéaires :
alldifferent(pX
1 , . . . , pn ),
sin(xi pi ) + x j = i, 1 ≤ i ≤ n,
j,i
sin(πxi ) = 0, 1 ≤ i ≤ n,
pi ∈ {1, . . . , n}, 1 ≤ i ≤ n,
xi ∈ [−10, 10], 1 ≤ i ≤ n.
5.4.2 Résultats
Les résultats expérimentaux sont présentés dans les Tables 5.2, 5.3 et 5.4. Toutes les expériences
ont été conduites sur un Intel Xeon à 3 GHz avec 16 Go de RAM, en utilisant le solveur RealPaver
1.0 [GB06] étendu pour l’occasion de façon à pouvoir traiter des contraintes alldifferent de chacune
des façons que nous avons décrites en détails à la Section 5.3. Les six premières colonnes correspondent
aux différentes méthodes ou combinaisons de méthodes que nous avons voulu comparer : MILP pour la
5.4. Expériences 81
formulation MILP de la contrainte ; , pour la formulation en utilisant seulement les contraintes de disé-
galités ; pour ces deux formulations, une colonne étiquetée +≥≤ (resp. +SUM) indique l’usage supplémen-
taire d’une relaxation convexe d’enveloppe utilisant les inégalités ad hoc (resp. utilisant l’égalité-somme
seule). Les deux dernières colonnes CARD et ALLDIF correspondent naturellement à la formulation de la
contrainte utilisant l’opérateur de cardinalité et à la formulation utilisant l’algorithme dédié de filtrage par
consistance de borne. Les temps de résolution en secondes sont donnés sur la première ligne, le nombre
de branchements sur la seconde et exprimé en milliers. Un point d’interrogation indique une méthode qui
n’a pas pu résoudre le problème donné en un temps raisonnable, c’est-à-dire inférieur à 10000 secondes.
Un nombre de branchements valant ∅ signifie qu’il est impossible de brancher sur ce type de variables
pour le problème considéré.
Pour commencer, comparant les relaxations SUM et ≥≤ on peut s’apercevoir que la première est plus
efficace que la seconde dans presque tous les cas et quel que soit le type de problème. À chaque fois, on
constate que le nombre de branchements est à peu près le même, montrant que la contrainte de somme
sur toutes les variables réalise quasiment la totalité du travail de filtrage. Ce constat est d’ailleurs assez
bien reflété par le graphique présenté en Figure 5.1, qui montre bien l’impact qu’a cette contrainte sur le
temps de résolution quelles que soient les contraintes qu’on lui rajoute dans la relaxation. Ainsi, l’effort
additionnel fourni par les inégalités ad hoc ne génère pas suffisamment de rejets de valeurs pour valoir
la peine d’être réalisé dans le cas général, et peut même aller jusqu’à faire accroître significativement les
temps de résolution —jusqu’à 200% par rapport au temps standard MILP ou ,. Un contre-exemple existe
cependant, quand on résout les instances du problème purement discret int-frac en utilisant la formulation
MILP : dans ce cas, utiliser les inégalités ad hoc aide à réduire le temps de résolution d’un dixième et le
nombre de branchement d’un tiers. Une explication pour ce phénomène peut être qu’il existe des valeurs
du domaine qui sont retirées du domaine par les contraintes , mais pas par les contraintes MILP, et que
MILP MILP , ,
Pb MILP , CARD ALLDIF
+≥≤ +SUM +≥≤ +SUM
16.6 10.7 9.7 14.3 10.0 8.5 5.4 5.1
sevvoth 4.7 2.1 2.1 5.5 3.0 3.0 2.1 1.7
∅ ∅ ∅ ∅ ∅ ∅ ∅ ∅
1308.9 1074.3 917.8 366.0 193.0 125.8 20.3 15.7
fractions 590.0 291.1 292.6 538.1 163.0 163.7 87.4 68.6
∅ ∅ ∅ ∅ ∅ ∅ ∅ ∅
? ? ? 2134.4 5989.7 2227.4 71.9 49.4
m-square 1859.8 1856.0 1859.8 155.5 98.0
∅ ∅ ∅ ∅ ∅
? 1536.2 ? 94.3 116.1 59.1 10.1 5.8
queens-10 163.0 49.8 27.7 27.9 18.2 8.0
∅ ∅ ∅ ∅ ∅ ∅
? ? ? 675.9 735.7 362.5 52.9 24.9
queens-11 310.1 142.2 143.1 91.4 36.9
∅ ∅ ∅ ∅ ∅
? ? ? 5200.1 4805.6 2289.1 260.3 133.5
queens-12 2039.2 782.6 785.3 377.9 187.7
∅ ∅ ∅ ∅ ∅
Table 5.2 – Temps de résolution en secondes (première ligne) et nombre de branchements sur variable
entière/réelle en milliers (seconde ligne).
82 U AD
80
70
60
50
40
30
20
10
0
0 100 200 300 400 500
Figure 5.1 – Temps de résolution en fonction du nombre de contraintes d’une relaxation dont les
contraintes ont été tirées aléatoirement, pour le problème Sevvoth. Le groupe des points correspondant
aux meilleurs temps de résolution quelque soit la taille correspond uniquement à des relaxations utilisant
la contrainte de somme sur toutes les variables.
ces valeurs peuvent aussi être retirées par les inégalités mais pas par la contrainte d’égalité sur toutes
les variables. En effet, pour le problème considéré on n’observe pas ce phénomène quand on utilise
la formulation , pour laquelle les valeurs filtrées par les inégalités semblent être déjà filtrées par les
diségalités, rendant les inéquations inutiles.
Concernant à présent les modélisations spécialisées CARD et ALLDIF, dans le cas particulier des pro-
blèmes purement discrets on s’aperçoit sans surprise qu’il y a une amélioration des temps de résolution
d’un facteur allant jusqu’à 10 (resp. 100) provenant de l’utilisation de ces formulations et des algorithmes
de filtrage dédiés en comparaison avec la formulation , (resp. MILP) —cf Table 5.2. De plus, dans ce
cas on s’aperçoit qu’il y a seulement un facteur inférieur à 2 entre CARD et ALLDIF en termes de temps
de résolution, tandis que dans le cas plus général des problèmes mixtes —cf Tables 5.3 et 5.4— et plus
particulièrement en ce qui concerne les problèmes issus du monde réel, il n’y a aucune différence signi-
ficative entre CARD et ALLDIF. Une première raison à cela peut être trouvée dans le fait qu’une grosse
partie du travail de filtrage a été réalisée par le processus d’exploration, au moyen des branchements,
et non par les algorithmes de filtrage eux-mêmes, signifiant finalement que plus le filtrage est simple et
plus la résolution est rapide : on s’aperçoit en effet qu’il y a une différence à peine significative vis-à-vis
de la formulation , pour les problèmes issus du monde réel. Mais il apparaît également dans ce dernier
cas que 90 % du travail de résolution est passé à traiter la partie continue/mixte des modèles. En effet,
le rapport entre les nombres de branchements sur variables entières et sur variables réelles nous prouve
qu’un dixième seulement du temps de résolution est consacré à l’instanciation des variables entières,
d’autant plus que la stratégie de branchement que nous avons utilisée consiste à toujours choisir d’abord
les variables entières. Ceci explique également pourquoi on n’observe pas ici une différence plus impor-
tant avec la formulation MILP, nettement moins performante pourtant sur les problèmes discrets ou les
problèmes mixtes artificiels, les seconds ayant une partie continue d’un volume bien inférieur à celle des
problèmes issus du monde réel.
Ainsi, nos expériences montrent que la formulation CARD est une bonne approche pour étendre un
solveur fonctionnant sur les intervalles de façon à gérer plus efficacement les problèmes mixtes faisant
usage de la contrainte alldifferent :
• le concept de cardinalité d’ensemble est une notion mathématique très utilisée, e.g. il y a déjà un
mot-clé card dans le langage de modélisation AMPL ;
5.4. Expériences 83
• l’algorithme de filtrage que nous avons implémenté est vraiment très simple à mettre en oeuvre ;
• cet algorithme peut drastiquement réduire les temps de résolution face à une modélisation classique
MILP ou même ,, même dans le cas particulier de problèmes purement discrets ;
• il est très compétitif par rapport à l’algorithme state-of-the-art de filtrage de la contrainte alldif-
ferent, spécialement quand il s’agit de problèmes mixtes issus du monde réels et dotés d’une
partie continue/mixte très importante.
D’autre part, nos expériences montrent clairement une nouvelle utilisation possible pour la relaxation
convexe d’enveloppe, qui non seulement se comporte aussi bien que l’utilisation standard dans le cas gé-
néral des problèmes mixtes non linéaires, mais qui en plus est très significativement meilleure en termes
de performances dans certains cas. Trouver quelles inégalités sont les plus à propos parmi toutes celles
que contient l’ensemble de contraintes de la relaxation, de taille exponentielle en nombre de variables,
ne semble pas valoir l’effort au vu du gain de performance. Utiliser la contrainte de somme sur toutes les
variables est suffisant pour atteindre bien plus rapidement une puissance de filtrage au moins équivalente
dans presque la totalité des cas.
MILP MILP , ,
Pb MILP , CARD ALLDIF
+≥≤ +SUM +≥≤ +SUM
817.8 222.2 252.6 33.2 16.8 9.4 1.7 1.0
int-frac-12 394.6 80.9 115.0 51.3 14.7 14.7 8.7 6.0
∅ ∅ ∅ ∅ ∅ ∅ ∅ ∅
? 6989.5 7878.3 996.0 416.0 225.3 29.9 19.1
int-frac-14 1824.4 2771.4 1160.2 263.3 263.5 149.1 102.0
∅ ∅ ∅ ∅ ∅ ∅ ∅
? ? ? ? ? 6189.1 725.7 441.8
int-frac-16 5607.6 3024.5 2082.7
∅ ∅ ∅
158.7 167.8 162.0 9.3 8.9 8.1 6.1 5.8
mixed-sin-5 0.1 0.1 0.1 7.9 6.9 6.9 6.3 5.3
271.2 271.2 271.2 11.3 9.3 9.3 8.4 8.4
6209.6 6612.4 6320.8 197.0 148.1 138.9 99.5 96.0
mixed-sin-6 0.7 0.7 0.7 123.0 83.7 83.7 72.6 59.1
8143.2 8143.5 8143.5 184.4 118.7 118.7 105.0 105.1
? ? ? ? 8866.2 8085.6 5530.5 5290.3
mixed-sin-7 3617.3 3621.0 3108.1 2558.0
6638.9 6639.3 5806.6 5792.3
hard- 4.5 4.5 4.3 0.9 0.8 0.9 0.8 0.8
0.0+ 0.0+ 0.0+ 0.3 0.2 0.2 0.2 0.1
mixed-sin-4 9.7 9.7 9.7 1.6 1.4 1.4 1.4 1.4
hard- 542.3 649.3 554.2 90.2 91.8 79.6 69.3 67.9
0.1 0.1 0.1 25.1 12.1 12.1 10.3 8.3
mixed-sin-5 956.5 956.5 956.5 126.3 108.8 108.8 106.8 106.7
Table 5.3 – Temps de résolution en secondes (première ligne) et nombre de branchements sur variable
entière/réelle en milliers (seconde ligne).
84 U AD
MILP MILP , ,
Pb MILP , CARD ALLDIF
+≥≤ +SUM +≥≤ +SUM
74.3 74.0 74.3 75.6 75.3 73.5 73.2 73.2
nuclear-A 0.0+ 0.0+ 0.0+ 0.0+ 0.0+ 0.0+ 0.0+ 0.0+
6.4 6.4 6.4 6.4 6.4 6.4 6.4 6.4
624.6 628.3 633.1 626.7 631.8 631.4 626.6 624.9
nuclear-B 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
12.6 12.6 12.6 12.6 12.6 12.6 12.6 12.6
850.9 1230.6 586.8 9.9 9.8 10.1 9.3 9.2
schedu-9 9.1 0.0+ 8.9 0.0+ 0.0+ 0.0+ 0.0+ 0.0+
160.5 261.7 108.8 2.0 2.0 2.0 2.0 2.0
1444.9 2545.5 1136.6 43.1 42.9 42.9 40.3 40.1
schedu-10 63.6 0.0+ 63.0 0.0+ 0.0+ 0.0+ 0.0+ 0.0+
165.4 377.7 120.7 6.6 6.6 6.6 6.6 6.6
? ? ? 231.7 227.2 229.6 215.9 213.9
schedu-11 0.0+ 0.0+ 0.0+ 0.0+ 0.0+
25.5 25.5 25.5 25.5 25.5
210.9 212.1 210.1 9.3 9.6 9.5 9.0 9.0
ELSP-8 0.0+ 0.0+ 0.0+ 0.0+ 0.0+ 0.0+ 0.0+ 0.0+
45.5 45.5 45.5 2.0 2.0 2.0 2.0 2.0
Table 5.4 – Temps de résolution en secondes (première ligne) et nombre de branchements sur variable
entière/réelle en milliers (seconde ligne).
5.5 Conclusion
Dans ce chapitre, nous avons présenté une étude expérimentale à propos de la résolution dans un sol-
veur continu de problèmes mixtes utilisant la contrainte discrète alldifferent . Nous avons d’abord
montré comment faire usage de cette contrainte pour modéliser une certaine catégorie de problème
MINLP issus du monde réel et qu’on peut trouver dans la littérature. Nous avons aussi passé en re-
vue les principales façons dont on peut traiter une contrainte alldifferent dans un solveur continu,
basé sur la notion d’intervalle, en reformulant la contrainte ou bien en implémentant des algorithmes de
filtrage dédiés. Finalement, nous avons réalisé des expériences sur différents types de problèmes de façon
à étudier l’impact du choix de la méthode de traitement sur la résolution proprement dite. À travers ces
travaux, nos principales contributions tiennent en deux points :
• Nous avons montré qu’un algorithme de filtrage dédié mais très simple à implémenter, comme
celui que nous avons réalisé pour la contrainte CARD, était très compétitif par rapport au très opti-
misé algorithme state-of-the-art de filtrage en ce qui concerne la résolution de problèmes mixtes
très contraints issus du monde réel, tout en restant relativement compétitif dans le cas particulier
de problèmes purement discret, ce qui n’est pas du tout le cas des formulations , et surtout MILP,
cette dernière étant pourtant la formulation d’usage ;
5.5. Conclusion 85
• Nous avons découvert que la relaxation convexe d’enveloppe telle qu’en usage actuellement pou-
vait être dépassée en termes de performances dans notre solveur à intervalles par un ensemble de
contraintes beaucoup plus facile à construire et à utiliser, et qui pouvait permettre une résolution
jusqu’à deux fois plus rapide pour résoudre certaines instances des problèmes mixtes issus du
monde réel.
Ce chapitre est dédié à la spécialisation aux contraintes mixtes des techniques de filtrage de domaines
basées sur l’arithmétique des intervalles. Après la présentation des travaux de A et Z à propos
de l’utilisation d’une arithmétique des intervalles d’entiers pour la résolution de problèmes purement
discrets, nous montrons comment le cadre classique de résolution dédié aux problèmes réels peut être
modifié pour mieux prendre en compte la notion d’intégralité des variables discrètes. Nous donnons en-
suite quelques précisions à propos de l’implémentation de ces modifications dans un solveur à intervalles,
puis nous réalisons des expériences pour en étudier le comportement en pratique.
87
88 Ś ’ ́ ’́
Les travaux présentés dans ce chapitre prolongent ceux de A et Z dans la direction d’une
meilleure connaissance des avantages en pratique de l’arithmétique des intervalles d’entiers par rapport à
l’arithmétique des intervalles de réels. En particulier, les auteurs n’ont pas été intéressés par une compa-
raison du comportement ni en théorie ni en pratique des méthodes de résolution basées sur l’arithmétique
des intervalles réels et celles basées sur l’arithmétique des intervalles d’entiers. D’autre part, nous élargis-
sons le spectre des opérations étendues aux intervalles d’entiers et enfin, notre approche n’est pas limitée
aux problèmes purement discrets : nous nous intéressons aux problèmes mixtes.
6.2.1.1 Logarithme
Soit x une variable entière. Supposons que dx = [a, b] contient strictement 0, i.e. a < 0 < b. Comme
le montre la Figure 6.2, il n’existe aucun nombre réel inférieur à 0 qui puisse être obtenu par log(x)
quand x est un entier. En effet, les seuls nombres dont le logarithme est un nombre réel inférieur à 0 sont
des nombres réels compris entre 0 et 1. Ainsi, effectuer l’opération de logarithme sur des entiers n’offre
aucune possibilité d’obtenir un nombre inférieur à 0. Ce qui nous donne l’encadrement suivant :
Comparons cette nouvelle définition pour l’opération de logarithme sur les entiers avec la définition
originale de cette opération telle qu’utilisée sur les réels. Quand la contrainte d’intégralité du domaine
log(x)
y ∈ [0, log(b)]
est relâchée et que le calcul est effectué sur un intervalle de réels, c’est la définition suivante qui nous
donne le résultat de l’encadrement :
Ainsi, puisque la fonction log(x), qui tend vers −∞ en 0, ne peut atteindre aucune valeur inférieure
à 0 quand x est un entier, on peut redéfinir l’opération correspondante sur les intervalles quand on traite
une variable dont le domaine est celui des entiers. Utiliser la définition (6.1) nous permettra dans ce cas
d’effectuer des évaluations plus précises des expressions arithmétiques, et cela sans perdre de solution
potentielle, puisque la contrainte y = log(x) avec x ∈ N n’admet aucune solution quand dy ∩ [0, +∞] = ∅.
6.2.1.2 Inverse
À présent, soit x une variable entière et supposons à nouveau que dx = [a, b] contienne strictement 0,
i.e. a < 0 < b. D’une façon similaire à ce que nous avons observé pour le logarithme, on s’aperçoit qu’il
n’existe aucun nombre réel inférieur à -1 ou supérieur à 1 qui puisse être obtenu par 1/x quand x est un
entier. En effet, les seuls nombres dont l’inverse est un nombre réel inférieur à -1 ou supérieur à 1 sont
des nombres réels compris entre -1 et 1 (voir Figure 6.3). Aussi, effectuer l’opération d’inverse sur des
entiers n’a aucune chance de retourner un nombre inférieur à -1 ou supérieur à 1. En conséquence, on a
l’encadrement suivant :
1/x ∈ [−1, 1/a] ∪ [1/b, 1] (6.3)
Comme précédemment pour le logarithme, comparons cette nouvelle définition de l’opération d’in-
verse sur les entiers avec la définition originale de cette opération qu’on utilise sur les réels. Lorsque
qu’on relâche la contrainte d’intégralité et qu’on effectue le calcul sur un intervalle de réels, le résultat
de l’encadrement nous est donné par la définition suivante :
y ∈ [ 1b , 1]
1
x x
1
y ∈ [−1, a]
Il est clair que la définition (6.3) permet, en resserrant les bornes des valeurs possibles, d’obtenir
un encadrement plus précis que la définition (6.4). En fait, il n’y a aucun nombre réel dans [−∞, −1) ∪
(1, +∞] qui puisse être le résultat de 1/x quand x est un entier car la fonction 1/x, qui tend vers −∞
en 0− et vers +∞ en 0+ , ne peut atteindre aucune valeur inférieure à -1 ou supérieure à 1 quand x
est un entier. Aussi, on peut rédéfinir l’opération correspondante sur les intervalles quand une traite
une variable dont le domaine est celui des entiers. L’utilisation de la définition (6.3) va donc permettre
d’effectuer le cas échéant des évaluations plus fines des expressions arithmétiques. Aucune solution
potentielle ne sera perdue puisque la contrainte y = 1/x avec x ∈ N n’admet aucune solution quand
dy ∩ ([−1, 1/a] ∪ [1/b, 1]) = ∅.
6.2.2 Application
Notre but ici n’est pas de présenter une arithmétique des intervalles complètement étendue aux pro-
blèmes mixtes. Cependant, il nous faut remarquer que la division est particulièrement utile dans le
contexte de la satisfaction de contraintes arithmétiques. En outre, elle est fortement mise à contribu-
tion au cours de la résolution de problèmes polynomiaux. Comme nous allons pouvoir le constater d’ici
la fin de cette section, l’utilisation d’une définition plus précise du calcul sur les intervalles va permettre
d’améliorer le filtrage d’une manière très significative.
Soit x ∈ [−10, 10] une variable entière et soit y ∈ [−10, 10] une variable réelle. Considérons le
système suivant :
(
x·y = 1
x+y = 2
Mais cette réduction du domaine de la variable x est trop faible pour avoir une incidence en retour
sur le domaine de y en utilisant la définition (6.4). Les domaines obtenus en fin de propagation restent
donc dx = dy = [−8, 10].
92 Ś ’ ́ ’́
dy := dy ∩ (1/dx ) := [ 13 , 1]
Enfin, en propageant à nouveau dans la deuxième contrainte1 puis dans la première, on obtient fina-
lement les réductions suivantes :
dx := dx ∩ (2 − dy ) := 1
dy := dy ∩ (1/dx ) := 1
De par l’utilisation de la nouvelle définition, on constate que les domaines obtenus en fin de pro-
pagation sont cette fois l’unique solution du système. Ainsi, la prise en compte du caractère entier des
variables pendant l’évaluation aux intervalles nous permet d’obtenir un filtrage bien plus fort.
au moins une variable entière x, de domaine [a, b] enveloppe ou pavé-consistant, il n’est pas possible
de trouver une valeur v ∈ [a, ⌈a⌉] (resp. [⌊b⌋, b]) telle que C(. . . , v, . . .) est satisfaite. Ainsi, une façon
de faire couramment mise en oeuvre est d’implémenter une procédure de filtrage dédiée à la contrainte
d’intégralité, qui a pour but de tronquer un domaine à ses bornes entières de la manière suivante :
dx := (int(dx ))
où int est l’opérateur qui tronque à ses bornes entières un intervalle à bornes flottantes, i.e.
(
[⌈a⌉, ⌊b⌋] si a ≤ b + 1
int([a, b]) :=
∅ sinon
dx := (dx ∩ F(dy ))
(6.5)
dx := (int(dx ))
où le traitement de la contrainte d’intégralité est réalisé après le filtrage du domaine. Mais nous allons voir
à présent qu’il est en fait possible d’effectuer la troncature des domaines de variables entières directement
à l’intérieur des opérateurs de filtrage associées aux contraintes arithmétiques elles-mêmes, obtenant de
cette manière une amélioration radicale des opérations de filtrage.
où la troncature aux bornes entières est réalisée avant le calcul de l’enveloppe de dx ∩F(dy ). L’exemple 6.5
va d’ailleurs nous donner un aperçu de combien il peut être intéressant de prendre ainsi en compte la
contrainte d’intégralité le plus tôt possible au cours des calculs sur les intervalles.
Exemple 6.5. Soit la contrainte y = x2 où x est une variable entière avec D x = [−2, 2] et y une variable
réelle avec Dy = [2, 3]. Le domaine de x est calculé au moyen de l’opération inverse du carré, appliquée
94 Ś ’ ́ ’́
√ √ √ √
à y et à son domaine. Il s’ensuit que x doit appartenir à l’intervalle [− 3, − 2] ∪ [ 2, 3]. Appliquant
la méthode (6.5), on obtient :
√ √ √ √
dx := (int(([− 3, − 2] ∪ [ 2 3]))
dx := (int([−1.732 . . . , 1.732 . . .])
dx := ([−1, 1]))
dx := [−1, 1]
L’exemple précédent met en évidence le gain potentiel qu’on peut attendre de cette nouvelle manière
de prendre en compte les contraintes d’intégralité des variables. Posons la proposition suivante :
Proposition 6.1. Soient d1x et d2x les domaines obtenus respectivement en utilisant les méthodes de cal-
cul (6.5) et (6.6). Alors on a d2x ⊆ d1x .
d’après est une propriété élémentaire de l’opération d’enveloppe. Il vient ensuite naturellement que
6.4 Implémentation
Dans cette section, nous donnons des détails à propos de l’implémentation que nous avons pu réaliser
des améliorations du filtrage dédiées aux contraintes arithmétiques mixtes, décrites dans les sections
précédentes.
6.4. Implémentation 95
Plate-forme utilisée
Nous avons de nouveau utilisé le solveur RealPaver [GB06] pour effectuer l’implémentation des opti-
misations du filtrage. Ainsi qu’on a pu s’en apercevoir précédemment, il s’agit d’un solveur de contraintes
continues, linéaires ou non linéaires, doté d’une API C++ construite autour d’un cœur en C. Il peut ainsi
être utilisé comme une application indépendante ou bien encapsulé dans une autre application. Nous
rappelons brièvement quelques unes de ses principales caractéristiques :
• les modèles peuvent être exprimés en utilisant un langage de modélisation dédié (voir Fig. 6.1) ;
• les domaines en entrée peuvent être des intervalles ou unions d’intervalles ;
• les variables déclarées peuvent être réelles ou entières ;
• les contraintes arithmétiques linéaires et non linéaires sont autorisées ;
• les autres contraintes à disposition incluent les contraintes en extension, les contraintes condition-
nelles et les contraintes par morceaux.
Comme on l’a déjà partiellement vu au Chapitre 3, la représentation des domaines, les structures
arborescentes utilisées pour représenter les expressions arithmétiques et les évaluer, ainsi que le parsing
des modèles, sont écrits en C ; les opérateurs de réduction de domaines, la propagation de contraintes,
les heuristiques de sélection de variable et de découpage de domaine sont implémentés en C++.
Les opérateurs de réduction de domaine implémentés incluent le filtrage par consistance d’enveloppe,
de pavé et 3B. Il y aussi de multiples heuristiques de choix de variables : priorité aux variables les plus
contraintes, priorité aux domaines entiers les plus fins (resp. réels les plus larges), priorité aux variables
de décision, tourniquet, etc. Enfin, les stratégies de découpage de domaine implémentés incluent la bis-
section et l’énumération.
Le plus gros du travail que nous avons réalisé ici concerne la partie bibliothèque de calcul sur les
intervalles du logiciel. C’est là que nous avons modifié les méthodes existantes et rajouté de nouvelles
P. 6.1 : Modèle RealPaver pour un problème d’intersection entre un cercle et une parabole.
1 problem circle
2
3 num p := 1.0e−8,
4 D := 2.75,
5 x0 := 2,
6 y0 := 1,
7 R := 1.25 ;
8
9
10 var x : real / p ˜ [−oo,+oo],
11 y : real / p ˜ [−100,100] ;
12
13 ctr (x−x0)^2 + (y−y0)^2 = R^2,
14 x^2 + y^2 = D^2 ;
15
16 solve * ;
17
18 end
96 Ś ’ ́ ’́
méthodes pour prendre en compte les nouvelles opérations arithmétiques et la gestion au plus tôt des
contraintes d’intégralité.
L’ajout des nouvelles opérations dans la bibliothèque est quasi immédiat étant données les règles de
calcul que nous avons données en Section 6.2. Nous avons ainsi implémenté deux nouvelles méthodes
effectuant respectivement l’opération de logarithme et la division d’un intervalle d’entiers. En ce qui
concerne ce dernier, remarquons qu’il n’y a aucun changement à faire quand il s’agit de traiter des
domaines qui ne contiennent pas 0.
Étant donnés une expression arithmétique et un pavé correspondant aux domaines courants des va-
riables, l’évaluation aux intervalles consiste à appliquer les règles de l’arithmétique des intervalles afin
d’obtenir un intervalle correspondant à l’encadrement du résultat. Concernant l’implémentation de nos
modifications, nous devons ici définir les bonnes méthodes à appeler e.g. c’est là que nous spécifions
le fait que pour une variable entière dont on veut effectuer le logarithme, on fait appel à la nouvelle
méthode.
Enfin, rappelons que la propagation descendante consiste à inférer de nouveaux domaines pour les
variables, étant donnés leur expression arithmétique et un pavé correspondant aux domaines courants des
autres variables. Aussi, nous avons dû modifier les appels aux méthodes de projection dans le cas des
noeuds ayant pour fils une variable entière —et récursivement les noeuds ayant pour fils un noeud ayant
pour fils une variable entière— de façon à faire appel aux bonnes opérations. Notons qu’il faut également
modifier les appels aux méthodes de projection concernant l’opérateur de multiplication, puisque x· y = z
implique x = z/y et y = x/z comme nous l’avons mis en évidence à la section précédente.
Figure 6.4 – Propagation descendante utilisant la méthode normale (gauche), utilisant la méthode de
troncature au plus tôt (droite).
De cette façon, on peut effectivement calculer l’enveloppe du résultat de la division, par exemple, après
seulement avoir replongé dans le domaine des entiers le résultat obtenu sur les réels. C’est cet algorithme
que nous avons implémenté dans RealPaver.
L’étape suivante consiste à remplacer chaque appel à la méthode d’enveloppe standard, après une
opération sur des variables entières capable de générer des trous dans les domaines, comme la division
ou l’élévation à une puissance paire, par un appel à cette nouvelle méthode. Dans chaque opérateur de
projection vers une variable entière capable de créer un trou dans le domaine et de retourner une union
d’intervalles —e.g. élévation/racine avec puissance paire ou multiplication/division— il faut remplacer
l’appel à la méthode standard par un appel à la nouvelle méthode. Nous avons ainsi défini de nouvelles
fonctions de projection, qui sont appelées à la place des anciennes, dans le cas d’une projection vers une
variable entière (voir Figure 6.4). De cette façon, on s’assure de ne plus surestimer les domaines qui ne
contiennent aucun entier mais dont l’enveloppe en contient.
6.5 Expériences
Les expériences présentées ici ont été réalisées sur un Intel Xeon à 3 GHz avec 16 Go de RAM. Nous
avons utilisé la version 1.0 de RealPaver ayant été modifiée de la façon que nous avons décrite dans la
section précédente. Dans le but de mesurer les apports de nos optimisations du filtrage, nous avons utilisé
deux problèmes mixtes dimensionnables que par simplicité nous nommerons problème A et problème B.
Problème A
Soit x un vecteur de variables réelles, de dimension n, et y un vecteur de variables entières, de dimen-
sion n également. Supposons que chaque variable ait pour domaine initial [−108 , 108 ] et considérons le
système d’équations suivant :
xi yi = 1, i = 1, . . . , n
n
X
x2i = 2n.
i=1
6.5. Expériences 99
Les résultats expérimentaux pour ce problème sont présentés dans la Table 6.1 et les Figures 6.5
et 6.6. Nous avons mesuré le temps de résolution, en secondes, et le nombre de branchements, pour une
dimension n du problème variant de 3 à 20. La première colonne correspond à la dimension de l’instance
testée. La deuxième colonne montre les résultats obtenus avec l’implémentation standard. La troisième
colonne présente quant à elle les résultats obtenus avec la nouvelle implémentation.
180 1.2e+06
std std
new new
160
1e+06
140
120 800000
100
600000
80
60 400000
40
200000
20
0 0
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16 18 20
Figure 6.5 – Pour le Problème A, temps standard face à nouveau temps (gauche), nombre de découpages
standard face à nouveau nombre de découpages (droite).
Pour ce problème, la nouvelle méthode se comporte mieux que la méthode standard quelle que soit la
dimension du problème. De plus, le gain en temps et en nombre de branchements est exponentiel comme
le montre la Figure 6.6. En effet, la méthode standard ne permet de réaliser que de faibles réductions
de domaines en utilisant la première contrainte, si bien que pour prouver l’inconsistance du problème,
la propagation de contraintes doit être combinée avec une technique de découpage de domaine. Ce-
pendant, cela requiert un nombre exponentiel de découpages, exactement 2n−1 . Mais avec la méthode
1e+07
time
#split
1e+06
100000
10000
1000
100
10
1
2 4 6 8 10 12 14 16 18 20
Figure 6.6 – Gains en temps et en nombre de découpages, en fonction du nombre de variables de l’ins-
tance résolue du Problème A (échelle logarithmique).
Problème B
Posons à présent le problème B. Soit x un vecteur de variables réelles, de dimension n, et y un vecteur
de variables entières, de dimension n également. Supposons que chaque variable a pour domaine inital
[−108 , 108 ] et considérons le système d’équations suivant :
(xi − 1)yi = 1, i = 1, . . . , n
n
Y
log( yi ) = 1.
i=1
Les résultats expérimentaux pour ce problème sont présentés dans la Table 6.2 et les Figures 6.7
et 6.8. À nouveau, nous avons mesuré le temps de résolution, en secondes, et le nombre de découpages,
pour une dimension n du problème variant de 3 à 20. La première colonne correspond à la dimension de
l’instance testée. La deuxième colonne montre les résultats obtenus avec l’implémentation standard. La
troisième colonne présente quant à elle les résultats obtenus avec la nouvelle implémentation.
Comme dans le problème précédent, la nouvelle méthode se comporte mieux que la méthode standard
quelle que soit la taille du problème. En outre, le gain est cette fois linéaire par rapport à la dimension du
problème, en temps comme en nombre de branchements —cf Figure 6.8— i.e. la résolution du problème
de dimension n avec la nouvelle méthode est n fois plus rapide et le nombre de branchements est n fois
inférieur. Ceci s’explique par le fait que de nombreuses branches de l’arbre de recherche sont coupées
après que la propagation de contraintes ait pu réussir à prouver l’inconsistance du pavé courant grâce aux
nouvelles méthodes de filtrage que nous avons implémentées spécialement pour les domaines entiers. En
6.5. Expériences 101
conséquence, ces branches inconsistantes ne sont pas explorées inutilement et la résolution s’en trouve
accélérée d’autant.
Ainsi, au vu des résultats mesurés sur ce second problème, les techniques que nous présentons dans
ce chapitre apparaissent de nouveau très prometteuses. Non seulement elles étaient déjà intéressantes
d’un point de vue théorique, montrant comment l’arithmétique des intervalles peut être améliorée afin
de gérer les domaines entiers d’une meilleure façon, avec plus de précision, mais il s’avère également
qu’elles sont très efficaces d’un point de vue pratique pour les problèmes considérés, quel que soit leur
dimension.
2000 1e+07
std std
1800 new 9e+06 new
1600 8e+06
1400 7e+06
1200 6e+06
1000 5e+06
800 4e+06
600 3e+06
400 2e+06
200 1e+06
0 0
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16 18 20
Figure 6.7 – Pour le Problème B, temps standard face à nouveau temps (gauche), nombre de découpages
standard face à nouveau nombre de découpages (droite).
20
time
18 #split
16
14
12
10
0
2 4 6 8 10 12 14 16 18 20
Figure 6.8 – Gains en temps et en nombre de découpages, en fonction du nombre de variables de l’ins-
tance résolue du Problème B.
6.6 Conclusion
Dans ce chapitre, nous avons présenté des améliorations pour les techniques de filtrage de domaine
basées sur l’arithmétique des intervalles. Ces améliorations dédiées à une meilleure utilisation des con-
traintes mixtes, i.e. des contraintes qui font intervenir à la fois des variables discrètes et des variables
continues. Prolongeant les travaux antérieurs de A et Z à propos d’une arithmétique des inter-
valles d’entiers, nous avons introduit de nouvelles règles de calcul sur les intervalles de certaines opéra-
tions arithmétiques. Celles-ci nous permettent de faire des encadrements plus fins du résultat de l’éva-
luation aux intervalles d’une expression arithmétique, en prenant en compte le caractère éventuellement
entier des variables sur lesquelles portent ces opérations. Nous avons également mis en lumière une nou-
velle façon de prendre en compte les contraintes d’intégralité au cours de la propagation de contraintes,
plus tôt dans le processus de propagation. Nous avons ensuite montré comment nous avons pu implémen-
ter ces améliorations dans un solveur à intervalles. Enfin, nous avons présenté des expériences dans le
but d’observer le comportement de nos améliorations vis-à-vis de la résolution de problèmes mixtes de
dimension variable. À travers ces travaux, nos contributions sont ainsi les suivantes :
• Nous donnons de nouvelles définitions pour des opérations typiques de l’arithmétique des inter-
valles, dans le cas particulier où leurs opérandes sont des variables entières, et montrons que cela
permet, à travers un augmentation de la force du filtrage obtenu, d’accélérer de manière très signi-
ficative la résolution des problèmes ;
• Nous proposons une nouvelle façon de gérer les contraintes d’intégralité des variables pendant le
processus de résolution : prendre en compte ces contraintes le plus tôt possible au cours de la pro-
pagation, c’est-à-dire non pas comme un post-traitement au filtrage des domaines mais au moment
même de celui-ci, permet d’accélérer de manière drastique la résolution des problèmes.
De futurs travaux pourraient être consacrés à mieux comprendre l’impact de chacune des amélio-
rations du filtrage sur la résolution des problèmes, en réalisant des observations plus précises comme
compter le nombre d’opérations arithmétiques effectuées, ou activer/désactiver chacune des optimisa-
6.6. Conclusion 103
tions, avant de comparer les résultats obtenus. En outre, il pourrait être intéressant de réaliser une étude
plus exhaustive de l’arithmétique des intervalles d’entiers, dans le but de trouver éventuellement d’autres
opérations arithmétiques qui pourraient bénéficier d’une telle spécialisation aux intervalles d’entiers, ou
bien le cas échéant de prouver que cela n’est pas possible. Pour finir, il faudrait étudier expérimentalement
l’impact des améliorations que nous présentons ici dans le cadre d’applications de la programmation par
contraintes issues du monde réel.
.
VII
C ONCLUSION ET PERSPECTIVES
105
106 C
par rapport à des outils à la pointe des techniques d’aujourd’hui. D’autre part, nous avons pu constater
par nous-mêmes une divergence assez significative dans les résultats obtenus par les différentes méthodes
que nous comparions. Nous avons suggéré que cette divergence pouvait être mise sur le compte de la
non-fiabilité des outils d’optimisation non basés sur l’arithmétique des intervalles, un état de fait déjà
mentionné dans la littérature [LMR07].
Ainsi, nous avons eu l’opportunité de faire l’application des techniques de programmation par con-
traintes pour la modélisation et la résolution d’un problème de conception en robotique pour lequel la
PPC n’avait jamais été mise à contribution jusqu’ici. L’intégration de méthodes orientées intervalle pour
l’optimisation permet d’acquérir une robustesse des solutions que les méthodes numériques classiques ne
permettent pas d’atteindre à cause des erreurs d’arrondis. Nos expériences ont par ailleurs montré qu’en
ce qui concerne le problème considéré, ce gain en robustesse s’accompagnait d’un gain en efficacité de
résolution, en temps de calcul.
offrant la possibilité de s’abstraire d’un mécanisme de modélisation rigide et difficile à utiliser. En outre,
cette collaboration nous permet de résoudre plus efficacement certains types de problèmes mixtes. À tra-
vers ces travaux, nos principales contributions sont ainsi les suivantes
• Nous avons mis au point une plateforme logicielle doté d’un langage de modélisation facilitant
l’expression des problèmes mixtes en fusionnant les hautes capacités d’expression de contraintes
arithmétiques du langage de modélisation dédié d’un solveur continu, avec les nombreuses primi-
tives de modélisation de contraintes globales d’un solveur discret, améliorant en cela les capacités
de chacun des deux outils pris individuellement.
• Nous avons également mis en évidence la nette supériorité en termes d’efficacité de résolution
d’une collaboration d’outils pour résoudre les problèmes mixtes qui font intervenir à la fois des
contraintes arithmétiques, nécessitant une évaluation par l’arithmétique des intervalles, et des con-
traintes globales, requérant des algorithmes de filtrage dédiés et optimisés.
férentes modélisations et algorithmes de filtrage associés. Il ressort de nos expériences que même s’il
reste très avantageux en général d’opter pour une modélisation globale de la contrainte plutôt que son
éclatement en une clique de contraintes binaires, il n’est pas intéressant d’opter pour un algorithme de
filtrage trop coûteux, comme c’est le cas de l’algorithme de filtrage de type consistance de bornes dans
sa version state-of-the-art, en ce qui concerne la résolution des problèmes difficiles issus du monde réel.
À travers ces travaux, nos principales contributions tiennent en deux points :
• Nous avons montré qu’un algorithme de filtrage global, dédié mais très simple à implémenter
comme celui que nous avons réalisé pour la contrainte CARD, était très compétitif par rapport au
très optimisé algorithme de filtrage en ce qui concerne la résolution de problèmes mixtes très
contraints issus du monde réel, tout en restant relativement compétitif dans le cas particulier de
problèmes purement discret, ce qui n’est pas du tout le cas des formulations , et surtout MILP,
cette dernière étant pourtant la formulation d’usage ;
• Nous avons découvert que la relaxation convexe d’enveloppe telle qu’en usage actuellement pou-
vait être dépassée en termes de performances dans notre solveur à intervalles par un ensemble de
contraintes beaucoup plus facile à construire et à utiliser, et qui pouvait permettre une résolution
jusqu’à deux fois plus rapide pour résoudre certaines instances des problèmes mixtes issus du
monde réel.
à-vis de la résolution de problèmes mixtes de dimension variable. La mise en évidence expérimentale des
effets sur la résolution de problèmes mixtes des nouvelles techniques mises en oeuvre nous a ainsi permis
de constater l’impact drastique que peuvent avoir en pratique ces optimisations naturelles du calcul sur
les temps effectifs de résolution des problèmes, au delà de l’efficacité théorique que nous avions déjà pu
remarquer. Ainsi, nos contributions à travers ces travaux sont les suivantes :
• Nous donnons de nouvelles définitions pour des opérations typiques de l’arithmétique des inter-
valles, dans le cas particulier où leurs opérandes sont des variables entières, et montrons que cela
permet, à travers un augmentation de la force du filtrage obtenu, d’accélérer de manière très signi-
ficative la résolution des problèmes ;
• Nous proposons une nouvelle façon de gérer les contraintes d’intégralité des variables pendant le
processus de résolution : prendre en compte ces contraintes le plus tôt possible au cours de la pro-
pagation, c’est-à-dire non pas comme un post-traitement au filtrage des domaines mais au moment
même de celui-ci, permet d’accélérer de manière drastique la résolution des problèmes.
Des travaux futurs pourraient permettre d’améliorer les points faibles majeurs d’une collaboration de
solveurs telle que nous l’avons mise en oeuvre. En effet, on peut imaginer pouvoir profiter davantage des
forces d’une telle collaboration en s’affranchissant de ses principaux inconvénients que sont pour nous
les temps de recopie de domaines d’un solveur à l’autre, ainsi que le manque de souplesse du propaga-
teur que nous avons utilisé pour calculer le point-fixe des domaines à partir des contraintes arithmétiques
continues et mixtes du modèle :
• En ce qui concerne le premier point, on peut en effet concevoir une symbiose plus profonde entre
les deux outils, qui pourrait amener par exemple à proposer une gestion plus intelligente de la
recopie des domaines, au moins dans le cas particulier où il n’y a pu se produire aucune réduction
et où la recopie en retour est inutile.
110 C
• Quant au deuxième point, on peut également envisager une collaboration plus en profondeur qui
permettrait une sélection plus fine des opérateurs de réduction à rappeler en cas de modification
de domaines, par exemple en offrant la possibilité de n’effectuer le calcul de point-fixe que sur un
sous-ensemble des contraintes arithmétiques déléguées au solveur continu.
• Il est également intéressant d’envisager une collaboration de solveurs centrée cette fois sur le sol-
veur continu. En effet, selon le type de problèmes il peut y avoir ou bien une majorité de contraintes
arithmétiques, ou bien une majorité de contraintes globales discrètes. Dans le premier cas, il est
sûrement avantageux de confier le rôle principal au solveur continu et de ne déléguer au solveur
discret que le filtrage des contraintes continues.
• Il semble important de trouver davantage de problèmes mixtes, issus du monde réel ou bien généré
d’une façon ou d’une autre, mais ayant un rapport plus équilibré entre le nombre de contraintes
discrètes et le nombre de contraintes continues, de façon à mieux connaître les comportements des
différentes modélisations selon les caractéristiques des problèmes résolus.
• À l’heure actuelle, aucun solveur de la platforme de résolution AMPL n’est en mesure de traiter le
mot-clé card, pourtant déjà une primitive de son langage de modélisation. Au vu de nos résultats,
il semblerait très utile d’y implémenter une méthode de filtrage associée qui pourrait permettre de
modéliser simplement et efficacement la contrainte alldifferent.
• Au delà des travaux que nous présentons ici sur la contrainte alldifferent, on peut tout à fait
imaginer réaliser un travail similaire à propos d’autres contraintes globales. Ainsi, on pourrait
même à terme se doter d’un ensemble minimal de primitives de type card, permettant d’expri-
mer un ensemble maximal de contraintes globales et implémentable facilement dans un solveur de
contraintes continues.
Nous voyons ici deux pistes principales pour de futures recherches sur le thème d’une meilleure maî-
trise d’une arithmétique des intervalles mixte :
• De futurs travaux pourraient être consacrés réaliser une étude plus exhaustive de l’arithmétique des
intervalles, dans le but de trouver éventuellement d’autres opérations arithmétiques qui pourraient
7.2. Perspectives de recherche 111
bénéficier d’une telle spécialisation aux intervalles d’entiers, ou bien le cas échéant de prouver que
cela n’est pas possible.
• En outre, il pourrait être intéressant de mieux comprendre l’impact sur la résolution des problèmes
de chacune des améliorations du filtrage que nous avons proposées ici, en réalisant des obser-
vations plus précises comme compter le nombre d’opérations arithmétiques effectuées, ou acti-
ver/désactiver chacune des optimisations, avant de comparer les résultats obtenus.
Mais au-delà des problèmes mixtes, il nous semble intéressant, pour conclure cette thèse, d’évoquer
des pistes dans le cadre plus général de la résolution de tout problème à variables continues, qu’il im-
plique également ou non des variable discrètes. En effet, nous avons eu l’opportunité au cours de cette
thèse de faire l’application des techniques de programmation par contraintes pour la modélisation et
la résolution d’un problème de conception en robotique, pour lequel la programmation par contraintes
n’avait jamais été mise à contribution jusqu’ici. Pour ce faire, nous avons eu à modifier le cadre classique
de résolution d’un solveur de contraintes continues, destiné à la satisfaction de contraintes uniquement,
pour y intégrer la notion de recherche rigoureuse d’optimum. En conclusion de ces travaux, nous sont
apparues des directions de recherche intéressantes que nous partageons ici :
• Nous avons non seulement pu constater que pour le problème considéré, les méthodes de pro-
grammation par contraintes sont très compétitives par rapport aux méthodes de programmation
mathématique, mais surtout que les résultats obtenus par l’une et l’autre méthode divergent sensi-
blement. Ainsi, il semble primordial d’élucider avec certitude la cause de cette divergence afin de
chercher à corriger ce qui en est la source.
• De plus, des travaux futurs devraient permettre de mettre en oeuvre d’autres critères de filtrage
utilisant la fonction-objectif, par exemple en se basant sur de la recherche locale, et d’ainsi ac-
célérer encore la convergence de notre algorithme d’optimisation de façon à permettre de traiter
efficacement des problèmes de taille plus importante.
• Notons également qu’au niveau de l’exploitation des résultats d’un point de vue purement robo-
tique, nous avons encore à faire l’application de la méthode que nous avons décrite ici afin de
tracer les isocontours de l’erreur maximum dans l’espace de travail, à partir du calcul de l’erreur
de pose maximum pour un échantillon de points de l’espace de travail du robot et de façon à mieux
connaître l’impact du jeu dans les articulations sur l’erreur de pose.
• En outre, d’autres travaux sont déjà en cours pour également prendre en compte les imperfections
d’usinage des pièces du robot. L’enjeu est ici de chercher à compenser le jeu dans les articulations
par ces imperfections et ainsi réduire drastiquement les coûts d’assemblage.
reuses dans des outils de modélisation mathématique comme MatLab. Des travaux sont déjà en
cours au sein de l’équipe Méthodes Ensemblistes pour l’Optimisation, du Laboratoire d’Informa-
tique de Nantes Atlantique, au sein de laquelle a été réalisée cette thèse, à propos de l’inclusion de
telles méthodes dans ces outils.
.
Quand les lettrés supérieurs ont entendu parler de la Voie,
Ils la pratiquent avec zèle.
Quand les lettrés du second ordre ont entendu parler de la Voie,
Tantôt ils la conservent, tantôt ils la perdent.
Quand les lettrés inférieurs ont entendu parler de la Voie,
Ils la tournent en dérision.
S’ils ne la tournaient pas en dérision,
Elle ne mériterait pas le nom de Voie !
Généralisations
Dans [Fre85] a été élaborée une généralisation des consistances d’arc et de chemin, la notion de
(i, j)-consistance. Grâce à ce nouveau formalisme, on peut voir la consistance d’arc comme étant la (1,1)-
consistance, où toute valeur doit, pour être conservée, pouvoir être complétée par la valeur d’une variable
adjacente dans le graphe. De même, la consistance de chemin est la (2,1)-consistance, où tout couple de
valeurs doit pouvoir être complété par la valeur d’une variable adjacente.
Définition A.1 ((i,j)-consistance). Un CSP discret P = (X, D, C) est (i,j)-consistant ssi toute instancia-
tion localement consistante de i variables peut être étendue à j variables.
F a auparavant élaboré dans [Fre78] la notion de k-consistance, une autre généralisation des
consistances d’arc et de chemin. L’idée est de considérer les inconsistances par ensemble de k variables,
et de pouvoir concevoir des consistances plus fortes que la consistance de chemin. De cette façon, la
consistance d’arc est la 2-consistance, et la consistance de chemin la 3-consistance. Mais s’il est dans
l’idée possible d’aller au delà de 3, on conçoit aisément que la complexité algorithmique d’un tel filtrage
puisse le rende inutilisable en pratique.
Définition A.2 (k-consistance). Un CSP discret P = (X, D, C) est k-consistant ssi toute instanciation
consistante de k − 1 variables peut être étendue à k variables.
Définition A.3 (k-consistance forte). Un CSP discret P = (X, D, C) est fortement k-consistant ssi il est
i-consistant pour tout i ≤ k.
Si on peut dès lors parler de PC forte, quand un CSP est à la fois 2-consistant et 3-consistant, il reste
que la PC est algorithmiquement prohibitive à assurer. Pour pallier cet inconvénient tout en cherchant à
filtrer plus de valeur que l’AC, diverses consistances ont progressivement été proposées, des consistances
intermédiaires entre AC et PC forte. Notons que ces nouvelles consistances sont dites de domaine, c’est-
à-dire qu’il ne s’agira pas comme pour la PC de retirer simultanément d’un ensemble de variables des
n-uplets inconsistants, mais uniquement de retirer du domaine d’une variable les valeurs ne respectant
pas une certaine propriété.
119
120 ANNEXE A
Consistances restreintes
Dans [Ber95], on trouve ainsi une première restriction de la PC, de façon à accélérer le filtrage tout
en conservant un important pouvoir filtrant :
Définition A.4 (Consistance de chemin restreinte). Etant donné un CSP P = (X, D, C), le domaine D
d’une variable X est chemin-consistant restreint ssi il est arc-consistant et si, pour toute valeur u de D
n’ayant qu’un seul support v dans le domaine d’une variable Y pour une contrainte C XY donnée, chaque
variable Z telle qu’il existe à la fois une contrainte C XZ portant sur X et Z et une contrainte CYZ portant
sur Y et Z contient dans son domaine une valeur w telle que (u, w) et (v, w) satisfont respectivement C XZ
et CYZ – i.e. la paire (u, v) est chemin-consistante.
Le terme anglo-saxon original est restricted path consistency (RPC). Cette consistance est destinée à
être un intermédiaire entre AC et PC. En effet, la RPC est moins forte que la PC puisqu’elle ne s’intéresse
qu’aux valeurs n’ayant qu’un seul support pour une contrainte, et non à toutes les valeurs comme c’est
le cas pour la PC. Elle reste cependant plus forte que l’AC puisqu’elle considère des triplets de variables
et non juste des couples.
On trouve aujourd’hui un algorithme de filtrage par RPC, dans le pire des cas de complexité tempo-
relle en O(en + ed2 + cd2 ) et spatiale en O(ed + cd), avec d la taille du plus grand domaine, n le nombre de
variables, e le nombre de contraintes et c le nombre de cliques1 de trois sommets dans le graphe [DB97a].
Exemple A.6. [Ber95] Filtrage par RPC :
Définition A.5 (Consistance de chemin k-restreinte). Etant donné un CSP P = (X, D, C), le domaine D
d’une variable X est chemin-consistant k-restreint ssi il est arc-consistant et si, pour toute valeur u de D
ayant au plus k supports vi dans le domaine d’une variable Yi pour une contrainte C XYi donnée, chaque
variable Z telle qu’il existe à la fois une contrainte C XZ portant sur X et Z et une contrainte CYi Z portant
sur Yi et Z contient dans son domaine une valeur w telle que (u, w) et (vi , w) satisfont respectivement C XZ
et CYi Z i.e. chaque paire (u, vi ) est chemin-consistante.
Il s’agit d’une extension de la RPC aux valeur ayant non pas 1 mais au plus k supports. Elle est donc
logiquement plus forte que la RPC puisqu’elle examine plus de valeurs. Poursuivant les recherches dans
1
Une clique est un graphe dans lequel tout sommet est relié par un arc à chacun des autres sommets.
ANNEXE A 121
la direction du meilleur compromis entre filtrage et complexité, Debruyne propose dans [DB97a] une
autre évolution de la RPC :
Définition A.6 (Consistance de chemin max-restreinte). Etant donné un CSP discret P = (X, D, C), le
domaine D d’une variable X est chemin-consistant max-restreint ssi toute valeur de D possède un support
chemin-consistant dans chaque domaine de variable Y pour laquelle il existe une contrainte portant à la
fois sur X et sur Y.
La max-RPC est plus forte que la RPC, car ici on s’intéresse à toutes les valeurs —non pas juste à
celles ayant un unique support— et plus forte également que la k-RPC, puisqu’on ne se restreint pas à
examiner les valeurs qui satisfont un nombre de support d’au plus k, toutes les valeurs sont considérées.
Pourtant, la max-RPC peut être plus facile à obtenir que n’importe quelle k-RPC : en effet, il n’y a pas
à chercher dans le CSP quelles valeurs satisfont la condition d’avoir au plus k supports, ce qui pouvait
s’avérer très coûteux quand k est grand.
Dans sa version optimale, l’algorithme a dans le pire des cas une complexité temporelle en O(en +
ed2 + cd3 ) et spatiale en O(ed + cd), avec d la taille du plus grand domaine, n le nombre de variables, e
le nombre de contraintes et c le nombre de 3-cliques dans le graphe [DB97a].
Exemple A.7. [Deb01] Filtrage par max-RPC :
Consistances inverses
Toujours dans cette optique de trouver un compromis optimal entre complexité et efficacité de filtrage,
dans [FE96] est élaborée la notion de consistance de chemin inverse :
Définition A.7 (Consistance de chemin inverse). Etant donné un CSP P = (X, D, C), le domaine D d’une
variable X est chemin-consistant inverse ssi toutes les valeurs des domaines de toutes deux variables Y
et Z sont consistantes pour toute valeur de D par rapport à toute contrainte portant sur un sous-ensemble
de {X, Y, Z}.
Selon le formalisme de F, il s’agit de la (1,2)-consistance, toute valeur peut être étendue à tout
couple de valeurs supplémentaire. On parle de PC inverse puisque la PC est la (2,1)-consistance, où tout
couple de valeur peut être étendu à toute troisième valeur.
Dans sa version optimale, l’algorithme a dans le pire des cas une complexité temporelle en O(en +
ed2 + cd3 ) et spatiale en O(ed + cd), avec d la taille du plus grand domaine, n le nombre de variables, e
le nombre de contraintes et c le nombre de 3-cliques dans le graphe [Deb00].
122 ANNEXE A
Définition A.8 (Consistance inverse de voisinage). Etant donné un CSP P = (X, D, C), le domaine D
d’une variable X est inverse voisinage-consistant ssi pour toute valeur v de D, il existe des valeurs wi
dans les domaines de chacune des variables Yi pour lesquelles existe une contrainte C j portant sur X et
sur une des Yi au moins, wi telles que toutes les C j sont satisfaites.
Le terme original est neighborhood inverse consistency. Les deux idées de bases sont la consistance
inverse, développée par F dès 1985 sans être implémentée, et la consistance de voisinage, qui
insiste sur le sous-problème que forme une valeur et les contraintes qu’elles imposent autour d’elle dans
le graphe [FE96].
L’algorithme de filtrage par NIC a aujourd’hui dans le pire des cas une complexité temporelle en
O(g2 (n + ed)dg+1 ) et spatiale en O(n), avec d la taille du plus grand domaine, n le nombre de variables, e
le nombre de contraintes et g le degré maximum d’un noeud du graphe [FE96].
Exemple A.9. [FE96] Filtrage par NIC :
Définition A.9 (Consistance de singleton). Etant donné un CSP P = (X, D, C) et une consistance discrète
CONS, le domaine D d’une variable X est singleton CONS-consistant ssi D n’est pas vide et pour toute
valeur v de D, le CSP dans lequel on a remplacé X par v n’est pas CONS-inconsistant.
ANNEXE A 123
Ici, on est sûr qu’il existe dans chaque domaine au moins une valeur qui participe avec v à une solu-
tion. En outre, en fonction de la consistance sur laquelle se base la consistance de singleton, on obtient
une consistance différente : consistance d’arc singleton (SAC) pour la consistance d’arc, consistance de
chemin restreinte singleton (SRPC) pour la RPC, etc.
L’algorithme de filtrage par SAC a dans le pire des cas une complexité temporelle en O(en2 d4 ) et
spatiale en O(ed), tandis que l’algorithme de filtrage par SRPC a lui dans le pire des cas une complexité
temporelle en O(en + n2 (e + c)d4 ) et spatiale en O(ed + cd) –avec d la taille du plus grand domaine, n le
nombre de variables, e le nombre de contraintes et c le nombre de 3-cliques du graphe [DB97b].
Exemple A.10. [Deb05] Filtrage par SRPC :
Synthèse
Pour conclure sur ce sujet, remarquons que ces consistances peuvent ainsi être mises en relation les
unes par rapport aux autres en fonction de la quantité de valeurs qu’elles filtrent. Par exemple, RPC est
plus forte qu’AC mais plus faible que PIC, ou encore SAC est plus forte que max-RPC. Le résultat de
cette comparaison exhaustive est la figure A.1 page 123, extraite de [DB97b] et [Deb01], où est représenté
cet ordre sur les consistances discrètes. On notera que certaines sont incomparables : cet ordre n’est pas
total.
Consistances fortes
Une limitation de la consistance d’enveloppe, que partage également la consistance de pavé, est la
localité de son application i.e. elle prend en compte les contraintes séparément. Or dès les premières
tentatives d’adapter aux domaines continus les techniques de filtrage par consistance, L définit
dans [Lho93] une consistance plus forte, capable de filtrer les domaines en utilisant toutes les contraintes
à la fois :
Définition A.10 (Consistance 3B). Etant donné un CSP P = (X, D, C), un domaine Di = [a, b] de D est
3B-consistant ssi les CSP P′ et P′′ dans lesquels on a remplacé Di respectivement par [ai , a+i ] et [b−i , bi ]
ne sont pas inconsistants après filtrage par consistance d’enveloppe.
Il s’agit effectivement d’une consistance de plus grande force puisqu’elle garantit qu’un domaine ne
puisse plus être réduit par aucune contrainte. De plus, il s’agit en fait d’une famille de consistances plutôt
que d’une consistance, un peu à l’image de la consistance discrète de singleton décrite précédemment
dans cette section.
Exemple A.11. [BG06] Soient les contraintes x + y = 0 et x − y = 0 avec D x = Dy = [−1, 1]. Ces
domaines sont irréductibles au sens de la consistance d’enveloppe. Cependant, la consistance 3B permet
de filtrer des valeurs. En effet, pour x = −1 par exemple, les contraintes impliquent à la fois y = 1 et
y = −1, ce qui provoque le rejet de cette valeur de x.
Notons en outre que la définition de la consistance 3B peut être réécrite pour n’importe quelle consis-
tance, en l’occurrence la consistance de pavé, comme dans [VHMD97] :
Définition A.11 (Bound-consistance). Etant donné un CSP P = (X, D, C), un domaine Di = [a, b] de D
est Bound-consistant ssi les CSP P′ et P′′ dans lesquels on a remplacé Di respectivement par [a, a+ ] et
[b− , b] ne sont pas inconsistants après filtrage par consistance de pavé.
Plus récemment a été introduite dans [CTN04] une nouvelle consistance ayant pour objet les do-
maines continus. Elle se base sur la notion de sous-ensemble arc-consistant d’un CSP, et est plus forte
que l’AC.
Définition A.12 (Consistance d’ensemble de pavés). Un CSP P = (X, D, C) est ensemble de pavés-
consistant ssi D est l’ensemble {B′ }, B′ ∈ D, des pavés maximaux tels que (X, {B′ }, C) est arc-consistant.
Cette nouvelle consistance est plus forte que l’AC, mais reste plus faible qu’une résolution globale
du CSP i.e. elle laisse dans les domaines certaines valeurs qui ne sont pas des solutions du CSP. Ces deux
points sont illustrés respectivement par les deux exemples A.12 et A.13 :
Exemple A.12. [CTN04] Soient les domaines D x = Dy = Dz = [−2, 2] et Dw = [1, 4]. Soient les
contraintes x2 = w, x = y, y = z et x = −z. La consistance d’arc est atteinte pour les domaines
ANNEXE A 125
D x = Dy = Dz = [−2, −1] ∪ [1, 2] et Dw = [1, 4], alors que le consistance d’ensemble de pavé élimine
le pavé tout entier puisque, dans le domaine de x, ni [−2, −1] ni [1, 2] n’appartiennent à un sous-pavé
arc-consistant.
Exemple A.13. [CTN04] Soient les domaines D x = Dy = Dz = [0, 2]. Soient les contraintes x = y,
x + z = 2 et y = z. Comme le pavé initial est arc-consistant, il est ensemble de pavés-consistant. Or la
vraie solution est (1, 1, 1).
Algorithmes prospectifs
Les algorithmes prospectifs, tout d’abord, comme leur nom l’indique ont la tête tournée vers le
futur. Il s’agit d’augmenter le niveau de filtrage à chaque noeud de façon à limiter les explorations
infructueuses.
Un premier algorithme de ce type est le forward checking : après chaque instanciation, d’une variable
x, on élimine du domaine de chaque variable y non encore instanciée les valeurs inconsistantes avec les
contraintes contenant x et y [HE80].
Exemple A.14. [Bar] Forward checking pour le problème des 4 reines3 :
Ensuite est apparu dans [SF94] MAC pour maintaining arc consistency. Le CSP est initialement
rendu arc-consistant et, à chaque nouvelle variable instanciée, on en propage les effets jusqu’à ce que le
2
Ce terme désigne un comportement perfectible du backtracking qui consiste à traverser inutilement des portions de l’arbre
: quand une valeur affectée à X1 est incompatible avec toute valeur de X3 , chercher entre temps quelles valeurs de X2 sont
compatibles avec celle de X1 est une perte de temps.
3
Le problème des n reines est un classique de la programmation par contraintes : il s’agit de placer n reines sur un échiquier
de telle sorte qu’aucune ne soit en position d’être prise par une autre.
126 ANNEXE A
CSP soit de nouveau arc-consistant : après chaque instanciation d’une variable, on vérifie, pour toutes les
valeurs de chaque domaine de variable non encore instanciée, qu’il existe un support pour les contraintes
qui la contiennent.
Algorithmes rétrospectifs
La deuxième catégorie d’algorithmes de recherche est celle des algorithmes rétrospectifs qui, eux,
ont la tête tournée vers le passé et regardent dans les anciens échecs d’instanciation.
Un algorithme de ce type est par exemple le conflict-directed backjumping (CBJ) apparu dans [Pro93]
: chaque échec sur le choix d’une valeur est expliqué par les précédents choix qui entrent en conflit, et en
cas d’impossibilité d’affecter une variable –quand aucune valeur potentielle n’est consistante– on revient
en arrière jusqu’à la variable la plus récemment instanciée dont la valeur provoque une incompatibilité.
On mémorise pour chaque reine la (les) ligne(s) de la (des) reine(s) qui bloque(nt) une (des) case(s).
Ainsi, quand à la ligne 6 il n’y a plus d’alternative à tester, on revient en arrière jusqu’au choix le plus
récent, c’est-à-dire 4.
Un autre exemple d’algorithme rétrospectif est le dynamic backtracking (DBT) [Gin93]. Comme pour
le CBJ, l’idée est de mémoriser et d’utiliser les causes des échecs d’instanciation de valeurs. Mais cette
fois, on ne revient que sur la valeur incompatible la plus récente et on ne remet pas en cause les autres
variables déjà instanciées. Ainsi, il s’agit plus d’une réparation que d’un retour arrière.
Affectant 1 à A, 2 à B, 3 à C puis 3 à D, il n’y a plus aucune valeur possible pour E. D étant la variable
la plus récemment affectée, on revient sur sa valeur 3, rendue inconsistante du fait des valeurs de A
et B. A nouveau il n’y a plus aucune valeur possible pour D –puisque 3 vient de devenir impossible à
affecter– on revient donc sur B la variable explicative la plus récemment affectée, pour laquelle la valeur
2 devient inconsistante. Au prochain pas, on devra lui affecter une autre valeur sans revenir entre temps
sur la valeur de C.
Algorithmes rétro-prospectifs
Enfin, la troisième catégorie d’algorithme de recherche est celle des algorithmes rétro-prospectifs,
qui conjuguent les propriétés des deux précédentes catégories. L’algorithme MAC-CBJ, par exemple,
intègre au CBJ le maintien de la consistance d’arc, mais en pratique il est rarement plus efficace que
MAC. Cela n’est pas le cas de l’algorithme MAC-DBT qui adjoint le maintien de la consistance d’arc au
retour-arrière dynamique [JDB00] : cet algorithme se base sur la notion d’explication —au sens de cause
d’une inconsistance— qui y semble très prometteuse en ce qui concerne l’amélioration de l’efficacité des
stratégies de recherche.
B IBLIOGRAPHIE
[Apt00] A K. R. : The role of commutativity in constraint propagation algorithms. ACM Tran-
sactions on Programming Languages and Systems 22 (2000), 1002–1036.
[Apt03] A K. : Principles of Constraint Programming. Cambridge University Press, New York,
NY, USA, 2003.
[AW07] A K. R., W M. : Constraint Logic Programming using Eclipse. Cambridge
University Press, New York, NY, USA, 2007.
[AZ03] A K. R., Z P. : A comparative study of arithmetic constraints on integer inter-
vals. In Apt, Fages, Rossi, Szeredi, Váncza (Eds.) Recent Advances in Constraints, LNAI
3010 (2003), Springer-Verlag, pp. 1–24.
[AZ07] A K. R., Z P. : An analysis of arithmetic constraints on integer intervals.
Constraints 12, 4 (2007), 429–468.
[Bar] B R. : Online guide to constraint programming (1998). http ://kti.mff.cuni.cz/
∼bartak/constraints/. Visité le 15 Mai 2010.
[Bar99] B R. : Constraint programming : In pursuit of the holy grail. In Proceedings of
WDS99 (invited lecture) (1999), pp. 555–564.
[BCDP07] B N., C M., D S., P T. : Global constraint catalogue : Past,
present and future. Constraints 12, 1 (2007), 21–62.
[BDM03] B M. R., D A. S., M A. : Minlplib - a collection of test models for
mixed-integer nonlinear programming. INFORMS Journal on Computing 15, 1 (2003),
114–119.
[Ber95] B P. : Improving domain filtering using restricted path consistency. In Procee-
dings of IEEE CAIA-95 (1995).
[Bes94] B̀ C. : Arc-consistency and arc-consistency again. Artificial Intelligence 65 (1994),
179–190.
[Bes06] B̀ C. : Constraint propagation. In Handbook of Constraint Programming (2006),
Elsevier.
[BFR99] B C., F E. C., R J. : Using constraint metaknowledge to reduce arc
consistency computation. Artificial Intelligence 107 (1999), 125–148.
[BG06] B F., G L. : Continuous and interval constraints. In Handbook of
Constraint Programming (2006), Elsevier, pp. 592–624.
[BGGP99] B F., G F., G L., P J.-F. : Revising hull and box consis-
tency. In International Conference on Logic Programming (ICLP’99) (1999).
[BLL∗ 09] B P., L J., L L., M F., W A. : Branching and bounds tighte-
ningtechniques for non-convex minlp. Optimization Methods Software 24, 4-5 (2009),
597–634.
129
130 BIBLIOGRAPHIE
[BMV94] B F., MA D., V H P. : Clp(intervals) revisited. In ILPS’94
(International Symposium on Logic Programming) (Ithaca, NY, USA, 1994), The MIT
Press.
[BMVH94] B F., MA D., V-H P. : Clp(intervals) revisited. In Procee-
dings of the International Symposium on Logic Programming (1994), pp. 124–138.
[BSG∗ 10] B N., S R., G A., C S., C P. : Finding the maximal pose error
in robotic mechanical systems using constraint programming. In Proceedings of IEAAIE
2010 (2010).
[BVH97] B F., V H P. : Introduction to the special issue on interval
constraints. Constraints 2 (1997), 107–112.
[CC09] C P., C S. : The Kinematic Sensitivity of Robotic Manipulators to Manufactu-
ring Errors. Tech. rep., IRCCyN, Internal Report No RI2009_4, 2009.
[CDR99] C H., D F., R M. : Comparing Partial Consistencies. Reliable Com-
puting 5, 3 (1999), 213–228.
[Che07] C R. : Résolution par satisfaction de contraintes appliquée à l’aide à la dé-
cision en conception architecturale. PhD thesis, École Nationale Supérieure d’Arts et
Métiers de Bordeaux, 2007.
[Cle87] C J. G. : Logical arithmetic. Future Computing Systems 2, 2 (1987), 125–147.
[CLPP05] C-L J., P J. M., P L. G. : Mixed integer optimization for cyclic
scheduling of multiproduct plants under exponential performance decay. Chemical engi-
neering research & design 83, 10 (2005), 1208–1217.
[CR97] C T., R D. : Subdivision direction selection in interval methods for global opti-
mization. SIAM Journal on Numerical Analysis 34 (1997), 922–938.
[CTN04] C G., T G., N B. : New light on arc consistency over continuous
domains. In International workshop on Constraint Propagation and Implementation. At
CP conference (2004).
[CTN05] C G., T G., N B. : Box-set consistency for interval-based
constraint problems. In SAC (2005), pp. 1439–1443.
[DB97a] D R., B̀ C. : From restricted path consistency to max-restricted path
consistency. In Proceedings of Third International Conference on Principles and Prac-
tice of Constraint Programming (1997), pp. 312–326.
[DB97b] D R., B̀ C. : Some practicable filtering techniques for the constraint satis-
faction problem. In Proceedings of IJCAI-97 (1997), pp. 412–417.
[Deb00] D R. : A property of path inverse consistency leading to an optimal pic algorithm.
In Proceedings ECAI’00 (2000), pp. 88–92.
[Deb01] D R. : Domain filtering consistencies. Journal of Artificial Intelligence Research
14 (2001), 205–230.
[Deb05] D R. : Cours de Master sur les réseaux de contraintes à domaines discrets, Faculté
des sciences et techniques, Nantes, 2005.
[Dec92] D R. : Constraint networks. In Encyclopedia of Artificial Intelligence (1992),
Wiley and Sons.
BIBLIOGRAPHIE 131
[GH88] G̈ H.-W., H J. : Some fundamental properties of local constraint propaga-
tion. Artificial Intelligence 36 (1988), 237–247.
[Gin93] G M. L. : Dynamic backtracking. Journal of Artificial Intelligence Research 1, 1
(1993), 25–46.
[GJM06] G I. P., J C., M I. : Minion : A fast scalable constraint solver. In Procee-
dings of ECAI 2006 (2006), pp. 98–102.
[GMR07] G I. P., M I., R A. : Tailoring solver-independent constraint models : A case
study with ESSENCE’ and Minion. In Proceedings of the 7th International Symposium
on Abstraction, Reformulation and Approximation (2007), pp. 18–21.
[Gou00] G F. : Langages et Environnements en Programmation par Contraintes d’Inter-
valles. PhD thesis, Ecole Doctorale : Sciences pour l’Ingénieur, Nantes, 2000.
[Gra05] G L. : Cours de Master sur les contraints numériques, Faculté des sciences et
techniques, Nantes, 2005.
[GW99] G I., W T. : CSPLib : a benchmark library for constraints. Tech. rep., Techni-
cal report APES-09-1999, 1999. Available from http ://csplib.cs.strath.ac.uk/. A shorter
version appears in the Proceedings of the 5th International Conference on Principles and
Practices of Constraint Programming (CP-99).
[Hal35] H P. : On representatives of subsets. Journal of the London Mathematical Society 10
(1935), 26–30.
[Han92] H E. : Global Optimization Using Interval Analysis. Marcel Dekker, New York,
1992.
[HDmT92] H P. V., D Y., T C. : A generic arc-consistency algorithm and its
specializations. Artificial Intelligence 57 (1992), 291–321.
[HE80] H R., E G. : Increasing tree search efficiency for constraint satisfaction
problems. Artificial Intelligence, 14 (1980), 263–313.
[Heu06] H M. : Modélisation et résolution d’une application d’aide au déploiement d’an-
tennes radio en programmation par contraintes sur le discret et le continu. PhD thesis,
UFR Sciences et Techniques, Université de Nantes, 2006.
[HMK97] H P. V., MA D., K D. : Solving polynomial systems using a
branch and prune approch. SIAM Journal on Numerical Analysis 34, 2 (1997), 797–827.
[Hoo02] H J. N. : Logic, optimization and constraint programming. INFORMS Journal on
Computing 14 (2002), 295–321.
[Hoo06] H J. : Operations research methods in constraint programming. In Handbook of
Constraint Programming (2006), Elsevier, pp. 521–564.
[Hoo07] H J. N. : Integrated Methods for Optimization, vol. 100 of International Series in
Operations Research & Management Science. Springer, 2007.
[HS05] H H. H., S̈ T. : Stochastic Local Search : Foundations and Applications. Mor-
gan Kaufmann Publishers, 2005.
[HSG01] H P., S D., G E. : A framework for cooperating solvers - a pro-
totypic implementation. In Proceedings of the Workshop on Cooperative Solvers in
Constraint Programming - CoSolv. At the Seventh International Conference on Principles
and Practice of Constraint Programming (2001).
BIBLIOGRAPHIE 133
[HT06] H H. H., T E. : Local search methods. In Handbook of Constraint Programming
(2006), Elsevier, pp. 135–167.
[HWPS98] H I., W T., P̈ R., S H. : Different transformations for
solving non-convex trim-loss problems by MINLP. European Journal of Operational
Research 105, 3 (1998), 594–603.
[IEE85] IEEE T P754 : ANSI/IEEE 754-1985, Standard for Binary Floating-Point Arithmetic.
IEEE, New York, Aug. 12 1985.
[Inn02] I C. : Kinematic Clearance Sensitivity Analysis of Spatial Structures With Re-
volute Joints. ASME Journal of Mechanical Design 124 (2002), 52–57.
[JDB00] J N., D R., B P. : Maintaining arc-consistency within dynamic
backtracking. In Principles and Practice of Constraint Programming (2000), pp. 249–
261.
[JL87] J J., L J.-L. : Constraint logic programming. In POPL (1987), pp. 111–119.
[Kea96] K R. : Rigorous Global Search : Continuous Problems. Kluwer Academic Publi-
shers, Dordrecht, Netherlands, 1996.
[Kei] K M. : The alphametics page. http ://www.mathematik.uni-bielefeld.de/∼sillke/
puzzles/alphametic/alphametic-mike-keith.html. Visité le 20 avril 2010.
[KK94] K B. K., K S. N. : An augmented lagrange multiplier based method for
mixed integer discrete continuous optimization and its applications to mechanical design.
ASME Journal of Mechanical Design (1994), 116–318.
[Lab00] L F. : Choco : Implementing a cp kernel. In CP00 Post Conference Workshop on
Techniques for Implementing Constraint programming Systems (TRICS) (2000).
[Law76] L E. : Combinatorial Optimization : Networks and Matroids. Holt, Rinehart and
Winston, New-York, 1976.
[Lho93] L O. : Consistency techniques for numeric csp. In IJCAI 1993 (1993).
[LMR07] L Y., M C., R M. : An Efficient and Safe Framework for Solving Opti-
mization Problems. Journal of Computational and Applied Mathematics 199, 2 (2007),
372–377. Special Issue on Scientific Computing, Computer Arithmetic, and Validated
Numerics (SCAN 2004).
[LOQTvB03] Ĺ-O A., Q C.-G., T J., B P. : A fast and simple algorithm
for bounds consistency of the alldifferent constraint. In Proceedings of the Eighteenth
International Joint Conference on Artificial Intelligence (2003), pp. 245–250.
[Mac77] M A. K. : Consistency in networks of relation. Artificial Intelligence 8 (1977),
99–118.
[MF85] M A. K., F E. C. : The complexity of some polynomial network consis-
tency algorithms for constraint satisfaction problems. Artificial Intelligence 25 (1985),
65–74.
[MH86] M R., H T. C. : Arc and path consistency revisited. Artificial Intelligence 28
(1986), 225–233.
[MLS94] M R., L Z., S S. : A Mathematical Introduction to Robotic Manipulation.
CRC, Boca Raton, Fl., 1994.
134 BIBLIOGRAPHIE
5.1 Pour chaque problème, #var donne le nombre variables discrètes/continues dans le mo-
dèle utilisant la contrainte alldifferent , et #sol représente le nombre de solutions
au problème. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.2 Temps de résolution en secondes (première ligne) et nombre de branchements sur va-
riable entière/réelle en milliers (seconde ligne). . . . . . . . . . . . . . . . . . . . . . . 81
5.3 Temps de résolution en secondes (première ligne) et nombre de branchements sur va-
riable entière/réelle en milliers (seconde ligne). . . . . . . . . . . . . . . . . . . . . . . 83
5.4 Temps de résolution en secondes (première ligne) et nombre de branchements sur va-
riable entière/réelle en milliers (seconde ligne). . . . . . . . . . . . . . . . . . . . . . . 84
137
TABLE DES FIGURES
2.1 [Gou00] Évaluation de f (x) = (x6 /4 − 1, 9x5 + 0, 25x4 + 20x3 − 17x2 − 55, 6x + 48)/50 en
forme naturelle (gauche), en forme de B (droite), pour une largeur d’intervalle
unitaire de 0,02. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.1 Schéma de synthèse, la taille de la police utilisée est proportionnelle à la quantité estimée
de travail à fournir pour compléter l’outil en question : Choco, CHR et ECLi PSe sont
déjà assez bien dotés, RealPaver, Minion et Gecode sont quant à eux plus spécifiques à
un type de problème. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.2 Principe de la collaboration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.1 Temps de résolution en fonction du nombre de contraintes d’une relaxation dont les
contraintes ont été tirées aléatoirement, pour le problème Sevvoth. Le groupe des points
correspondant aux meilleurs temps de résolution quelque soit la taille correspond unique-
ment à des relaxations utilisant la contrainte de somme sur toutes les variables. . . . . . 82
139
TABLE DES MATIÈRES
1 Introduction 1
1.1 La programmation par contraintes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Problématique de la thèse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Nos contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 Organisation du document . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 État de l’art 9
2.1 Modélisation des problèmes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.1.1 En programmation par contraintes . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.1.2 En programmation mathématique . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.3 Généralisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 Techniques de résolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.1 Filtrage des domaines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.2 Exploration de l’espace de recherche . . . . . . . . . . . . . . . . . . . . . . . 25
141
142 TABLE DES MATIÈRES
4.3.2 Résultats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
6 Spécialisation aux domaines d’entiers du filtrage basé sur l’arithmétique des intervalles 87
6.1 Travaux en rapport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.2 Opérations arithmétiques sur intervalles d’entiers . . . . . . . . . . . . . . . . . . . . . 89
6.2.1 Nouvelles définitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
6.2.2 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.3 Gestion des contraintes d’intégralité . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
6.3.1 Troncature aux bornes entières . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
6.3.2 Amélioration du mécanisme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
6.4 Implémentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
6.4.1 Nouvelles opérations arithmétiques . . . . . . . . . . . . . . . . . . . . . . . . 96
6.4.2 Gestion au plus tôt des contraintes d’intégralité . . . . . . . . . . . . . . . . . . 96
6.5 Expériences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
6.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
Bibliographie 129
Abstract
Constraints are a generic way of expressing regularities that rule our world. Given a set of constraints,
an important question is to determine whether there exists a way of simultaneously satisfying them all.
This problematic lies deep inside constraint programming, a powerful paradigm for efficiently solving
problems that appear in many areas of human activity. During the 80’s originally, it was dedicated to
solving discrete problems in artificial intelligence. Then during the 90’s, constraint programming was
extended so as to solve continuous problems. However, mixed problems —with both discrete and conti-
nuous variables— have only been a very small concern so far.
In this thesis, we take the continuous solving framework point of view. We propose and implement
several improvement of this framework :
• Integration of a rigorous global optimum search to the classical, without objective solving framework,
in order to model and solve a conception problem in robotics ;
• Cooperation of two solvers, a discrete one and a continuous one, cooperation which is more efficient
than both solvers at solving problems involving both discrete and continuous constraints ;
• Comparison of several ways for modeling and filtering the alldifferent discrete, global constraint
within a solver dedicated to continuous problems ;
• Specialization of the filtering techniques that make use of interval arithmetic, so as to improve the
filtering power of discrete and mixed discrete/continuous arithmetic constraints.