Dans ce notebook (utilisant Python 3), je souhaite montrer des animations de l'algorithme du Simplexe, un peu comme dans la vidéo suivante :
J'aimerai bien écrire un petit morceau de code Python qui fait les étapes suivantes :
scipy.optimize.linprog(method="simplex"), et s'arrête s'il n'y a pas de solution trouvée par le simplexe ;callback de cette fonction pour afficher des équations en LaTeX représentant l'évolution du système et des variables de la base et hors base ;interact ;Ce document ne sera pas :
from IPython.display import YouTubeVideo
# https://www.youtube.com/watch?v=W_U8ozVsh8s
YouTubeVideo("W_U8ozVsh8s", width=944, height=531)
On a sûrement besoin de Numpy et Matplotlib :
import numpy as np
import matplotlib.pyplot as plt
On a besoin de la fonction scipy.optimize.linprog(method="simplex") du module scipy.optimize :
from scipy.optimize import linprog
On a aussi besoin de la fonction IPython.display.Latex pour facilement afficher du code LaTeX généré depuis nos cellules Python :
from IPython.display import Latex, display
Par exemple :
def display_cos_power(power=1):
return display(Latex(fr"$$\cos(x)^{power} = 0$$"))
for power in range(1, 5):
display_cos_power(power)
On va avoir besoin des widgets IPywidgets, plus tard :
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
Et enfin, de l'extension itikz
%load_ext itikz
Déjà , je vais écrire le problème étudié comme un dictionnaire, que l'on pourra passer à scipy.optimize.linprog(method="simplex") :
# Objective Function: 50x_1 + 80x_2
# Constraint 1: 5x_1 + 2x_2 <= 20
# Constraint 2: -10x_1 + -12x_2 <= -90
problem1 = {
# Cost function: 50x_1 + 80x_2
"cost": [50, 80],
# Coefficients for inequalities
"A_ub": [[5, 2], [-10, -12]],
# Constraints for inequalities: 20 and -90
"b_ub": [20, -90],
# Bounds on x, 0 <= x_i <= +oo by default
"bounds": (0, None),
}
# Objective Function: maximize x_1 + 6*x_2 + 13*x_3
# => so cost will be opposite
# Constraint 1: x1 <= 200
# Constraint 2: x2 <= 300
# Constraint 3: x1+x2+x3 <= 400
# Constraint 2: x2+3x3 <= 600
problem2 = {
# Cost function: minimize -1*x_1 + -6*x_2 + -13*x_3
"cost": [-1, -6, -13],
# Coefficients for inequalities
"A_ub": [
[1, 0, 0],
[0, 1, 0],
[1, 1, 1],
[0, 1, 3],
],
# Constraints for inequalities:
"b_ub": [200, 300, 400, 600],
# Bounds on x, 0 <= x_i <= +oo by default
"bounds": (0, None),
}
Puis une petite fonction qui s'occupe de prendre ce dictionnaire et le donner à scipy.optimize.linprog(method="simplex") :
def linprog_wrapper(problem, **kwargs):
result = linprog(
problem["cost"],
A_ub=problem["A_ub"],
b_ub=problem["b_ub"],
bounds=problem["bounds"],
method="simplex",
**kwargs
)
return result
On va déjà vérifier que l'on peut résoudre ces deux exemples de problème de programmation linéaire :
linprog_wrapper(problem1)
linprog_wrapper(problem2)
C'est bien la solution $x^* = [0, 300, 100]$, avec un objectif valant $+3100$, qui était trouvée dans la vidéo !
Et si on ajoute un callback ?
def round(np_array):
res = np.array(np.round(np_array), dtype=int)
if res.size > 1:
return list(res)
else:
return res
def dummy_callback(r):
print(f"\n- Itération #{r['nit']}, phase {r['phase']} :")
fun = round(r['fun'])
print(f" Valeur objectif = {fun}")
slack = round(r['slack'])
print(f" Variables d'écart = {slack}")
x = round(r['x'])
print(f" Variables objectif = {x}")
# print(r)
linprog_wrapper(problem2, callback=dummy_callback)
def equation_latex_from_step(result):
return r"""
\text{Minimiser} %s
\begin{}
"""
def interactive_latex_exploration(problem):
problem_solved = make_show_latex(problem1)
if problem_solved.status != 0:
print("Error: problem was not solve correctly, stopping this...")
interactive_function = make_show_latex(problem1)
max_step = problem_solved.nitint
return interact(, step=(0, max_step))
Allez on essaie :
interactive_latex_exploration(problem)
%load_ext itikz
Par exemple on peut afficher un premier exemple, avant de chercher à les faire bouger :
%%itikz --temp-dir --file-prefix simplex-example-
\documentclass[tikz]{standalone}
\usepackage{amsfonts}
\begin{document}
% from http://people.irisa.fr/Francois.Schwarzentruber/algo2/ notes
\usetikzlibrary{arrows,patterns,topaths,shadows,shapes,positioning}
\begin{tikzpicture}[scale=0.012, opacity=0.7]
\tikzstyle{point} = [fill=red, circle, inner sep=0.8mm];
\draw[->] (0, 0, 0) -- (300, 0, 0) node[right] {a};
\draw[->] (0, 0, 0) -- (0, 350, 0) node[above] {b};
\draw[->] (0, 0, 0) -- (0, 0, 300) node[below] {c};
\coordinate (O) at (0,0,0);
\coordinate (D) at (200,0,0);
\coordinate (E) at (200, 0, 200);
\coordinate (F) at (0, 0, 200);
\coordinate (G) at (0, 300,0);
\coordinate (C) at (200,200,0);
\coordinate (A) at (100,300, 0);
\coordinate (B) at (0,300, 100);
\draw[fill=blue!20] (O) -- (D) -- (E) -- (F) -- (O) -- cycle;
\draw[fill=blue!20] (D) -- (C) -- (E) -- cycle;
\draw[fill=blue!20] (G) -- (B) -- (F) -- (O) -- cycle;
\draw[fill=blue!20] (B) -- (A) -- (C) --(E) -- cycle;
\draw[fill=blue!20] (B) -- (F) -- (E) -- cycle;
\draw[fill=blue!20] (B) -- (A) -- (G) -- cycle;
\node[point] at (0,0,0) {}; % TODO make this argument of function
\end{tikzpicture}
\end{document}
Maintenant on peut chercher à contrôler la position du point objectif actuel : a,b,c sera $x_1, x_2, x_3$.
simplex_example_str = ""
def default_cost(a, b, c):
"""1*{a} + 6*{b} + 13*{c}"""
return 1*a + 6*b + 13*c
def show_tikz_figure_with_point(a=0, b=0, c=0, cost=default_cost):
# TODO generate nice LaTeX equations
if cost:
current_cost = cost(a, b, c)
cost_doc = cost.__doc__.format(a=a, b=b, c=c)
print(f"Coût = {cost_doc} = {current_cost}")
equation_latex = f"""\
Cout $f(a,b,c) = {cost_doc} = {current_cost}$.\
"""
display(Latex(equation_latex))
# now tikz
global simplex_example_str
simplex_example_str = r"""
\documentclass[tikz]{standalone}
\begin{document}
% from http://people.irisa.fr/Francois.Schwarzentruber/algo2/ notes
\usetikzlibrary{arrows,patterns,topaths,shadows,shapes,positioning}
\begin{tikzpicture}[scale=0.016, opacity=0.7]
\tikzstyle{point} = [fill=red, circle, inner sep=0.8mm];
\draw[->] (0, 0, 0) -- (300, 0, 0) node[right] {a};
\draw[->] (0, 0, 0) -- (0, 350, 0) node[above] {b};
\draw[->] (0, 0, 0) -- (0, 0, 300) node[below] {c};
\coordinate (O) at (0,0,0);
\coordinate (D) at (200,0,0);
\coordinate (E) at (200, 0, 200);
\coordinate (F) at (0, 0, 200);
\coordinate (G) at (0, 300,0);
\coordinate (C) at (200,200,0);
\coordinate (A) at (100,300, 0);
\coordinate (B) at (0,300, 100);
\draw[fill=blue!20] (O) -- (D) -- (E) -- (F) -- (O) -- cycle;
\draw[fill=blue!20] (D) -- (C) -- (E) -- cycle;
\draw[fill=blue!20] (G) -- (B) -- (F) -- (O) -- cycle;
\draw[fill=blue!20] (B) -- (A) -- (C) --(E) -- cycle;
\draw[fill=blue!20] (B) -- (F) -- (E) -- cycle;
\draw[fill=blue!20] (B) -- (A) -- (G) -- cycle;
\node[point] at (""" + f"{a}, {b}, {c}" + """) {};
\end{tikzpicture}
\end{document}
"""
#print(simplex_example_str)
# TODO: run this from this function?
#%itikz --temp-dir --file-prefix simplex-example- simplex_example_str
return get_ipython().run_line_magic(
"itikz", "--temp-dir --file-prefix simplex-example- simplex_example_str"
)
show_tikz_figure_with_point(0, 0, 0)
Et en rendant cela interactif :
interact(
show_tikz_figure_with_point,
a = (-100, 300, 10),
b = (-100, 300, 10),
c = (-100, 300, 10),
cost = fixed(default_cost)
)
linprog_wrapper(problem2, callback=dummy_callback)