Mini tutoriel pour la résolution de programmes linéaires avec Python - ENS Rennes 2021

Pré-requis pour exécuter ce notebook

  • Soit vous le téléchargez et vous l'utilisez localement, dans ce cas il faut que vous ayez installé le module scipy.

    • Utilisez votre gestionnaire de paquet système (apt-get par exemple) si Python installé par ce gestionnaire ;
    • Utilisez pip install scipy ou pip3, ou avec sudo pip (Linux/Mac) ou pip.exe (Windows) si modules Python gérés avec pip (recommandé) ;
    • Utilisez conda install scipy si modules Python gérés avec conda (si installé avec Anaconda.
  • Soit vous utilisez le lien fourni dans Discord pour exécuter ce notebook dans un environnement en ligne, avec MyBinder (gratuit libre et sans connexion).

In [21]:
try:
    import scipy.optimize
except ImportError:
    print("Vous avez lu le paragraphe précédent...?")
    print("Je t'envoie sur https://scipy.org/install.html et tu auras plus d'informations...")
    import webbrowser
    webbrowser.open_new_tab("https://scipy.org/install.html")

Quelques petits problèmes linéaires

Premier problème linéaire

Il vient du tutoriel susmentionné :

In [44]:
# Objective Function: 50x_1 + 80x_2
# Constraint 1: 5x_1 + 2x_2 <= 20
# Constraint 2: -10x_1 + -12x_2 <= -90

result = scipy.optimize.linprog(
    [50, 80],  # Cost function: 50x_1 + 80x_2
    A_ub=[[5, 2], [-10, -12]],  # Coefficients for inequalities
    b_ub=[20, -90],  # Constraints for inequalities: 20 and -90
    bounds=(0, None), # Bounds on x, 0 <= x_i <= +oo by default
)

On utilise les fonctionnalités de scipy pour les problèmes linéaires (doc), et pour commencer la seule fonction scipy.optimize.linprog :

In [45]:
if result.success:
    print(f"X1: {round(result.x[0], 2)} hours")
    print(f"X2: {round(result.x[1], 2)} hours")
else:
    print("No solution")
X1: 1.5 hours
X2: 6.25 hours

Et voilà, pas plus compliqué !

Vous pourrez observer que le résultat de l'appel à cette fonction (si tout marche bien) est un objet qui encapsule plusieurs choses :

  • la valeur de la solution, result.x ;
  • le nombre d'itérations, result.nit ;
  • l'état des slacks variables (cf cours sur le simplexe), result.slack ;
  • évaluation de l'état de l'optimisation, result.success, et result.message ;
  • et même la valeur de la fonction objectif à ce point solution (c'est utile pour gagner du temps, si cette fonction objectif est très couteuse, par exemple quand on apprend un gros réseau de neurone, avec d'autres solveurs).
In [26]:
print(result)
     con: array([], dtype=float64)
     fun: 574.99999998789
 message: 'Optimization terminated successfully.'
     nit: 4
   slack: array([ 3.34711814e-10, -1.83780458e-09])
  status: 0
 success: True
       x: array([1.5 , 6.25])

Problème exemple du cours sur le simplexe

Attention, généralement les solveurs cherchent à minimiser la fonction de coût, pas à la maximiser !

C'est un piège classique, on rentre le problème, le solveur répond [0,0,..,0] comme solution, et on ne sait pas ce qui n'a pas marché...

In [54]:
# Objective Function: x_1 + 6*x_2 + 13*x_3
# Constraint 1: x_1 <= 200
# Constraint 2: x_2 <= 300
# Constraint 3: x_1 + x_2 + x_3 <= 400
# Constraint 4: x_2 + 3*x_3 <= 600
# les variables sont supposées positives par défaut
# x_1 >= 0
# x_2 >= 0
# x_3 >= 0

result = scipy.optimize.linprog(
    [-1, -6, -13],  # Cost function: -x_1 + -6*x_2 + -13*x_3 to MINIMIZE
    A_ub=[  # Coefficients for inequalities
        [1, 0, 0],  # for C1: 1*x_1 + 0*x_2 + 0*x_3 <= 200
        [0, 1, 0],  # for C2: 0*x_1 + 1*x_2 + 0*x_3 <= 300
        [1, 1, 1],  # for C3: 1*x_1 + 1*x_2 + 1*x_3 <= 400
        [0, 1, 3],  # for C4: 0*x_1 + 1*x_2 + 3*x_3 <= 600
    ],
    b_ub=[200, 300, 400, 600],  # Constraints for inequalities: 200, 300, 400, 600
    bounds=(0, None), # Bounds on x, 0 <= x_i <= +oo by default
    method="simplex",
)
In [55]:
print(result)
     con: array([], dtype=float64)
     fun: -3100.0
 message: 'Optimization terminated successfully.'
     nit: 5
   slack: array([2.00000000e+02, 5.68434189e-14, 0.00000000e+00, 0.00000000e+00])
  status: 0
 success: True
       x: array([  0., 300., 100.])
In [56]:
if result.success:
    print(f"X1: {round(result.x[0], 2)} chocolats simples")
    print(f"X2: {round(result.x[1], 2)} pyramides")
    print(f"X2: {round(result.x[2], 2)} pyramides de luxe")
else:
    print("No solution")
X1: 0.0 chocolats simples
X2: 300.0 pyramides
X2: 100.0 pyramides de luxe

On trouve donc la solution commerciale optimale pour Charlie le chocolatier : 300 pyramides, et 100 pyramides de luxe.

Exercice 1 : sur ce problème, cherchez quelles contraintes (inégalités) sont saturées (devenues des égalités)

Indice : lire le tableau result.slack et l'interpréter.

In [ ]:
# TODO
print(result)

print("Variables slack:")
print(result.slack)

Exercice 2 : résoudre un autre problème vu en cours

Bonus : réfléchir à une situation de votre vie personnelle qui pourrait être mis en forme comme un problème linéaire

Un exemple que j'ai beaucoup est le suivant, qui généralise l'idée de « régime diététique optimal » : https://jeremykun.com/2014/06/02/linear-programming-and-the-most-affordable-healthy-diet-part-1/ (très bon blogue à suivre si ça vous plaît).

Comparer différentes méthodes :

Comme le dit la documentation :

The linprog function supports the following methods:

linprog(method=’simplex’)
linprog(method=’interior-point’)
linprog(method=’revised simplex’)
linprog(method=’highs-ipm’)
linprog(method=’highs-ds’)
linprog(method=’highs’)

Certaines méthodes peuvent ne pas être disponible sur votre installation, mais normalement "simplex" et "interior-point" sont disponibles partout.

Exercice 3 : sur un problème de votre choix, testez au moins deux méthodes.

Ce petit morceau de code peut vous aider :

In [33]:
methods = [
    "simplex",
    "interior-point",
    #"revised-simplex",
    #"highs-ipm",
    #"highs-ds",
    #"highs",
]
In [34]:
def solve_problem_1(method):
    return scipy.optimize.linprog(
        [50, 80],  # Cost function: 50x_1 + 80x_2
        A_ub=[[5, 2], [-10, -12]],  # Coefficients for inequalities
        b_ub=[20, -90],  # Constraints for inequalities: 20 and -90
        method=method
    )
In [35]:
for i, method in enumerate(methods):
    # solve problem with this method
    print(f"\n- Pour la méthode #{i}, {method}...")
    solution = solve_problem_1(method)
    print(f"La solution trouvée est {solution}")
- Pour la méthode #0, simplex...
La solution trouvée est      con: array([], dtype=float64)
     fun: 575.0
 message: 'Optimization terminated successfully.'
     nit: 2
   slack: array([0., 0.])
  status: 0
 success: True
       x: array([1.5 , 6.25])

- Pour la méthode #1, interior-point...
La solution trouvée est      con: array([], dtype=float64)
     fun: 574.99999998789
 message: 'Optimization terminated successfully.'
     nit: 4
   slack: array([ 3.34711814e-10, -1.83780458e-09])
  status: 0
 success: True
       x: array([1.5 , 6.25])

Exercice 4 : Chercher un problème qui donne une réponse différente sur deux méthodes différentes.

Avec le code ci-dessous, cherchez un problème plus compliqué qui donne une solution différente. Même si la différence est faible, commentez là :

  • en terme de nombre d'étape ?
  • valeurs de x* ?
  • valeur de f(x*) ?

Conclusion

Si vous êtes très curieux, regardez un peu les références données en haut de ce document.

N'hésitez pas à nous contacter si vous avez des questions !