Optimisation mathématique en Julia avec JuMP 0.21

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. MathOptInterface.jl propose aussi le module FileFormats pour l'export vers des formats normalisés (LP, MPS, notamment) — ce dernier ayant été développé sous le nom de MathOptFormat.jl.

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 :

add JuMP MathOptInterface
add Cbc

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 :

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 un solveur à ce modèle. Par exemple, pour Cbc :

using JuMP, Cbc
m = Model(Cbc.Optimizer)

2-1. Variables▲

@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

@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

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

2-2. Contraintes▲

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

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

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

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

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

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

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).

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..

@variable(m, x, Bin)
# ...
value(x)

@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 :

@constraint(m, c, x + y <= z)
# ...
dual(c)

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

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 un solveur abstrait : on peut l'utiliser pour résoudre plusieurs modèles en séquence.

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.

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 :

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.

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