Optimisation mathématique en Julia avec JuMP 0.20

Julia est un langage de programmation très dynamique, facilitant notamment la création de langages spécialisés. JuMP.jl en est un exemple : il s'agit d'un paquet Julia qui facilite la modélisation de programmes d'optimisation mathématique.

Cet article est obsolète et n'est gardé que pour des raisons historiques, car une version plus à jour est disponible.

Commentez Donner une note  l'article (5)

Article lu   fois.

L'auteur

Profil Pro

Liens sociaux

Viadeo Twitter Facebook Share on Google+   

1. Écosystème JuliaOpt

JuMP.jl est la pointe de l'iceberg de l'écosystème JuliaOpt. Il s'agit de l'interface principale pour un utilisateur de Julia, grâce à laquelle il peut modéliser des programmes d'optimisation en restant très proche des notations mathématiques.

Les versions actuelles de JuMP.jl utilisent une couche d'abstraction des solveurs, MathOptInterface.jl (souvent abrégée comme MOI). Les versions de JuMP antérieures à la 0.19 utilisaient plutôt MathProgBase.jl (MPB) : certains solveurs ne disposent toujours que d'une interface MathProgBase, mais la majorité est déjà passée à MathOptInterface.

Cette couche d'abstraction met à disposition toute une flopée de solveurs, notamment les plus utilisés. Du côté commercial : CPLEX, Gurobi, KNITRO, Mosek, Xpress ; du côté libre : Clp et Cbc, GLPK, Ipopt, SCS. Bien d'autres solveurs sont disponibles, y compris des solveurs écrits en Julia.

Plusieurs paquets permettent de travailler par-dessus une formulation, comme Dualization.jl pour en obtenir le dual ou MathOptFormat.jl pour l'export vers des formats normalisés (LP, MPS, notamment).

Tous ces paquets sont disponibles depuis le gestionnaire de paquets de Julia. Il est nécessaire d'au moins installer JuMP et la couche d'accès d'un solveur. Pour Clp, Cbc et GLPK, l'installation de la couche d'abstraction déclenche automatiquement l'installation du solveur ; pour les autres, le solveur doit déjà être installé sur votre machine. L'entrée dans le mode gestion de paquets se fait avec la touche ], puis l'installation comme ceci :

 
Sélectionnez
add JuMP MathOptInterface
add Cbc

2. Création d'un modèle d'optimisation

Avant de travailler sur un modèle, il faut créer une instance de Model. Il s'agit du point central d'un modèle d'optimisation : cet objet regroupe toutes les variables et toutes les contraintes d'un modèle d'optimisation. Les opérations de résolution ou de dualisation sont lancées au niveau du modèle. Il contient aussi une référence vers un solveur.

On peut créer un modèle sans définir de solveur :

 
Sélectionnez
using JuMP
m = Model()

Cependant, en l'état, il n'est pas possible de lancer un solveur sur ce modèle. Pour y arriver, il faut lier une fabrique de modèles à ce modèle, grâce à la fonction with_optimizer et à un solveur. Par exemple, pour Cbc :

 
Sélectionnez
using JuMP, Cbc
m = Model(with_optimizer(Cbc.Optimizer))

2.1. Variables

Ensuite, on peut définir des variables d'optimisation. JuMP permet de définir des variables continues, avec ou sans bornes, mais aussi discrètes ; on peut encore définir des matrices semi-définies positives. À chaque fois, on utilise la macro @variable (avec le symbole arobase). Celle-ci prend en argument le modèle auquel lier les variables à créer, le nom de ces variables, puis des spécifications sur leurs bornes ou leur type.

 
Sélectionnez
@variable(m, x) # Variable continue, sans borne
@variable(m, x >= 0) # Variable continue avec 
# borne inférieure, positive ou nulle
@variable(m, x <= 0) # Variable continue avec 
# borne supérieure, négative ou nulle
@variable(m, 0 <= x <= 2) # Variable continue avec 
# borne supérieure et inférieure
@variable(m, x <= 5, Int) # Variable entière avec 
# borne supérieure
@variable(m, x, Bin) # Variable binaire

JuMP permet aussi de créer des vecteurs, des matrices ou des tenseurs de plus haut ordre de variables d'un seul coup :

 
Sélectionnez
@variable(m, x[1:5] >= 0) # Vecteur de cinq variables continues
@variable(m, x[i in 1:5, j in i:5] >= 0) # Matrice triangulaire

Il n'est pas nécessaire de nommer la variable sur laquelle l'itération se fait (comme dans le premier exemple), mais cela aide souvent à la lisibilité du code. Une variable d'itération peut être utilisée pour définir les dimensions suivantes, comme dans la déclaration d'une matrice triangulaire de variables.

On peut aussi créer des matrices semi-définies positives (ce qui correspond à une contrainte dans le cône semi-défini) :

 
Sélectionnez
@variable(model, x[1:2, 1:2], PSD)

JuMP ne vous permettra pas de créer plusieurs variables avec le même nom !

2.2. Contraintes

Finalement, la dernière étape dans la définition d'un modèle est l'ajout de contraintes. Pour ce faire, JuMP offre la macro @constraint. Sa fonctionnalité la plus utilisée est l'ajout de contraintes linéaires :

 
Sélectionnez
@constraint(m, x + y <= z)
@constraint(m, x - z >= 1 + y)
@constraint(m, x <= y <= z)

JuMP permet notamment d'utiliser la fonction sum pour sommer sur des séries d'expressions :

 
Sélectionnez
@constraint(m, z == sum(i * x[i] for i in 1:5)

Comme pour les variables, on peut créer d'un coup une série de contraintes avec la même syntaxe :

 
Sélectionnez
@constraint(m, [i in 1:5], x[i] <= y[i])

Jusqu'à présent, ces contraintes n'ont aucun nom, aucune référence. Ceci est problématique pour accéder aux valeurs duales, par exemple. La même macro permet de donner des noms aux contraintes, qu'elles soient créées à l'unité ou dans une boucle :

 
Sélectionnez
@constraint(m, c1, x + y <= z)
@constraint(m, c2[i in 1:5], x[i] <= y[i])

Tout comme pour les variables, vous ne pourrez pas créer plusieurs contraintes avec le même nom.

En plus de ces contraintes linéaires, JuMP est prévu pour l'optimisation conique : on peut donc ajouter des contraintes d'appartenance à un cône classique. Pour le moment, JuMP dispose de trois cônes : du second ordre (SOCP et RSOCP), semi-défini positif (SDP), exponentiel (EXP), puissance (POWER), moyenne géométrique, logarithme/racine du déterminant, ainsi que leurs cônes duaux.

Pour le cône du second ordre, on peut utiliser une syntaxe similaire aux contraintes linéaires :

 
Sélectionnez
@constraint(model, x^2 + y^2 <= t^2)

Cependant, de manière générale, JuMP préfère utiliser une autre syntaxe pour les cônes : indiquer qu'une certaine fonction appartient à un ensemble donné (ce qui est parfois écrit F-in-S dans la documentation). Les ensembles sont définis dans MathOptInterface et non dans JuMP. Par exemple, la contrainte précédente peut s'écrire comme :

 
Sélectionnez
@constraint(m, [t, x, y] in MOI.SecondOrderCone(3))

L'argument de MOI.SecondOrderCone est sa dimension, le nombre de variables que l'on contraint. Pour ce cône, il en faut au moins deux. Traditionnellement, la première variable est notée t. Ce cône correspond à l'expression mathématique ∑ᵢ xᵢ² <= t².

Les autres cônes sont décrits en détail dans la documentation de MathOptInterface, avec leur expression mathématique complète.

3. Accès à la solution

Une fois le modèle défini, il est temps de passer à sa résolution. Pour lancer le solveur choisi, on utilise la fonction optimize!. Elle renvoie une exception si aucun solveur n'est défini : pour récupérer le résultat de l'appel au solveur, il faut utiliser termination_status. Cette fonction renvoie un état très précis ; en général, cependant, il s'agit de OPTIMAL ou de INFEASIBLE (les autres sont définis dans la documentation de MathOptInterface).

 
Sélectionnez
optimize!(m)
status = termination_status(m)
if status == MOI.OPTIMAL
    # Parfait !
elseif status == MOI.INFEASIBLE
    # Il n'existe aucune solution
else
    # Voir la documentation
end

Pour récupérer la valeur d'une variable, on utilise la fonction value. Pour plusieurs valeurs, il faut utiliser la version avec l'opérateur de diffusion . : value..

 
Sélectionnez
@variable(m, x, Bin)
# ...
value(x)
 
Sélectionnez
@variable(m, x[1:2], Bin)
# ...
value.(x)

Pour accéder aux valeurs duales de contraintes, on peut utiliser la fonction dual. Cette fonction prend en argument une référence à une contrainte :

 
Sélectionnez
@constraint(m, c, x + y <= z)
# ...
dual(c)
 
Sélectionnez
@constraint(m, c[i in 1:5], x[i] <= y[i])
# ...
dual.(c)

4. Organisation d'une application

Ces quelques rudiments devraient être suffisants pour développer de petites applications d'optimisation. Cependant, dès que le code devient légèrement compliqué (grand nombre de paramètres lus de l'extérieur, intégration d'un sous-problème pour une génération de contraintes ou de colonnes, etc.), il vaut mieux utiliser des techniques légèrement plus avancées.

Tout d'abord, pour le stockage de paramètres, il est pratique de définir une structure pour tous les stocker. Cette structure peut notamment stocker une fabrique de solveur (la sortie de la fonction with_optimizer) : on peut utiliser cette fabrique pour résoudre plusieurs modèles en séquence.

 
Sélectionnez
struct ProblemData
  graph::SimpleDiGraph
  incompatibility::Vector{Pair{SimpleEdge, SimpleEdge}}
  enable_something::Bool

  solver::OptimizerFactory
end

Ensuite, vous pouvez organiser votre construction du modèle en plusieurs fonctions : définition des variables de base, ajout des contraintes d'incompatibilité une à une, parties optionnelles du modèle, etc. Dans ce cas, il vaut mieux stocker toutes les références dans un autre objet, cette fois-ci muable.

 
Sélectionnez
mutable struct ProblemModel
  data::ProblemData
  m 
  x
  y

  c_basic
  c_incompatibility
end

Lors de la construction de tels modèles, il est parfois utile de recourir à la définition de variables anonymes. En effet, on ne peut pas attacher deux fois des variables ou des contraintes qui ont le même nom à un modèle. La syntaxe est très proche de la définition de variables et de contraintes ordinaires ; la grande différence est qu'on ne peut pas indiquer de bornes pour les variables de la même manière, on doit utiliser les arguments mots clés lower_bound et upper_bound ; pour des variables entières, il s'agit de integer et de binary pour des variables binaires. La construction de modèles peut alors se faire comme ceci :

 
Sélectionnez
function basic_model(d::ProblemData)
  m = Model(d.solver)
  x = @variable(m, [1:5], binary=true)
  y = @variable(m, [1:5], lower_bound=0, upper_bound=5)
  c_basic = @constraint(m, [i = 1:5], y[i] <= 5 * x[i])

  return ProblemModel(d, m, x, y, c_basic, nothing)
end

function incompatibility_model(d::ProblemData)
  # Inutile de répéter la définition des variables, c'est déjà fait dans basic_model
  pm = basic_model(d)

  # On peut remplir les derniers champs
  pm.c_incompatibility = @constraint(m, x[1] + x[3] <= 1)
  return pm
end

Avec l'approche des variables et des contraintes anonymes, évidemment, les variables et les contraintes n'ont plus de nom. Cela peut être problématique pour afficher le modèle ou l'exporter, car les noms aident à s'y retrouver. Pour pallier ce problème, on utilise la fonction set_name.

 
Sélectionnez
x = @variable(m, [1:5], binary=true)
for i in 1:5
  set_name(x[1], "x_binary_$(i)")
end

Vous avez aimé ce tutoriel ? Alors partagez-le en cliquant sur les boutons suivants : Viadeo Twitter Facebook Share on Google+   

  

Copyright © 2019 Thibaut Cuvelier. Aucune reproduction, même partielle, ne peut être faite de ce site ni de l'ensemble de son contenu : textes, documents, images, etc. sans l'autorisation expresse de l'auteur. Sinon vous encourez selon la loi jusqu'à trois ans de prison et jusqu'à 300 000 € de dommages et intérêts.