 |
Introduction à la Recherche OpérationelleM2 Statistique Science des Données (SSD) |
Anatoli Juditsky Franck Iutzeler |
import cvxpy as cp
import pandas as pd
import numpy as np
foods = pd.read_csv("data/food.dat",sep='\t')
n,p = foods.shape
Le tableau foods contient tous les aliments présent à la carte d'un restaurant avec leurs informations nutritionelles.
foods
n
Alimentation idéale en termes de nutriments¶
nutr = pd.read_csv("data/nutr_ideal.dat",sep=' ')
La tableau nutr contient les valeurs pour une alimentation optimale en terme de nutriments.
nutr
2- Trouver un régime optimal¶
Combinaisons pour une alimentation pauvre en calories¶
Pour commencer, nous allons essayer de trouver une combinaison $x\in\mathbb{R}^n$ des aliments ($x_i$ correspondant à la quantité d'aliment $i$ prise) qui minimise le total de calories sous contrainte de rester idéal vis à vis des nutriments (CalFat, Fat, VitA, etc.)
Question: Formaliser et résoudre le problème décrit ci-dessus. Qu'observez-vous sur la solution?
- Fonction objective linéaire à minimiser: $c^\mathrm{T}x$ où $c$ contient les calories par aliment
x = cp.Variable(n)
cal = np.array(foods["Cal"]) # calories
diet1 = cp.Minimize(cal @ x)
- Contrainte à vérifier: $x$ est positif (par de quantités négatives d'aliments) et $Ax=b$ où $A$ contient les valeurs des nutriments par aliment et $b$ contient les valeurs idéales
# Constraints matrix
A = np.array(foods.drop(columns=["Food","Cal"])).T
b = nutr.values.T.reshape((-1,))
A
b
contrainte1 = x>= 0.0
contrainte2 = A @ x - b == 0.0
- Resolution par CVXPY
Attention, utilisez le Solver "ECOS"
prob.solve(verbose=True,solver="ECOS")
le solveur par défaut peut donner de fausses informations...
prob = cp.Problem(diet1 , constraints = [contrainte1, contrainte2])
prob.solve(verbose=True)
prob.status
Relaxation des contraintes¶
Le problème précedent n'avait pas de point faisable avec les contraintes imposées. Afin d'obtenir des points faisable, nous allons remplacer la contrainte $$ Ax = b$$ par les deux contraintes de type boite $$ Ax \geq 0.8 b ~~ \text{ et } ~~ Ax \leq b $$ ainsi les contraintes de nutriments sont satisfaites entre 80 et 100% et le problème devient faisable.
Question: Définir les deux contraintes de boite et resoudre le problème. Vérifier qu'il est faisable.
contrainte1 = x>= 0.0
contrainte2a = A @ x >= 0.8*b
contrainte2b = A @ x <= b
prob = cp.Problem(diet1 , constraints = [contrainte1, contrainte2a, contrainte2b])
prob.solve(verbose=True)
Question: Affichez les valeur des nutriments pour la solution obtenue et la comparer à 80% et 100% des valeurs prescrites.
base = pd.DataFrame()
base["0.8*b"] = 0.8*b
base["solution"] = np.dot(A,x.value)
base["1.0*b"] = b
base
prob.constraints[1].dual_value
prob.constraints[2].dual_value
Question: Donner le nombre de calories totales et les aliments pour lesquelles la quantité est supérieure à 0.3.
xopt = x.value
xopt
for j in range(n):
if xopt[j]>0.3:
print("{} \t ({} cal) \t : \t {}".format(foods.iloc[j,0],foods.iloc[j,1],xopt[j]))
On remarque que la solution n'est pas pratique!
D'autres contraintes¶
On va imposer:
- pas de
EQUAL_0_Calorie_Sweetener 1_pkg_(1.0_g)(item $73$) ni deSPLENDA_No_Calorie_Sweetener 1_pkg_(1.0_g) : 30.9360336697206(item $248$) car c'est pas très bon - pas plus de 10 unités dans le menu, pas plus de 2 unités par item
- au moins un petit déjeuner (items 3 à 22)
Relacher au besoin les contraintes de nutrition...
Question : Implémenter ces contraintes une par une et regarder la valeur optimale de la fonction (le nombre de calories). Cette valeur augmente-t-elle ou diminue-t-elle ? Pourquoi ?
contrainte3a = x[73] == 0
contrainte3b = x[248] == 0
prob2= cp.Problem(diet1 , constraints = [contrainte1, contrainte2a,contrainte2b, contrainte3a, contrainte3b])
prob2.solve(verbose=True,solver="ECOS")
print("Sans les contraintes en plus:", prob.value, "Avec les contraintes en plus:", prob2.value)
x.value[73]
prob2.constraints[3].dual_value
Contrainte suivante: pas plus de 10 unités dans le menu, pas plus de 2 unités par item
contrainte4a = cp.sum(x) <= 10
contrainte4b = x <= 2.0
prob3= cp.Problem(diet1 , constraints = [contrainte1, contrainte2a,contrainte2b, contrainte3a, contrainte3b, contrainte4a, contrainte4b])
prob3.solve(verbose=True,solver="ECOS")
print("Valeur optimales successives:", prob.value, prob2.value, prob3.value)
Contrainte suivante: au moins un petit déjeuner (items 3 à 22)
dej = np.zeros(n)
dej[3:23] = 1.0
contrainte5 = dej@x >=1.0
prob4= cp.Problem(diet1 , constraints = [contrainte1, contrainte2a,contrainte2b, contrainte3a, contrainte3b, contrainte4a, contrainte4b, contrainte5])
prob4.solve(verbose=True)
print("Valeur optimales successives:", prob.value, prob2.value, prob3.value,prob4.value)
xoptNew = x.value
for j in range(n):
if xoptNew[j]>0.3:
print("[{:3d}]{} \t ({} cal) \t : \t {}".format(j,foods.iloc[j,0],foods.iloc[j,1],xoptNew[j]))
3- D'autres problèmes¶
Beaucoup d'autres problèmes de ce genre existent. Ci-dessous vous trouverez un nouveau problème mais vous êtes également invités à inventer vos propres problèmes et chercher si ils sont solvables par CVXPY !
Question : Quels sont les aliments à choisir pour obtenir au moins les nutriments recommendés en un minimum d'unités ?
y = cp.Variable(n)
efficace = cp.Minimize( cp.sum(y))
contrainte = A@y >= b
problemEff = cp.Problem(efficace , constraints = [contrainte, y>=0])
problemEff.solve(verbose=True)
xoptEff = y.value
for j in range(n):
if xoptEff[j]>0.3:
print("[{:3d}]{} \t ({} cal) \t : \t {}".format(j,foods.iloc[j,0],foods.iloc[j,1],xoptEff[j]))
Question : Quels nutriments sont les plus dur à obtenir (on pourra regarder la valeur de la variable duale) ?
for i in range(len(b)):
print("{} \t for {} \t {}".format(A.dot(y.value)[i],b[i],nutr.columns[i]))
problemEff.constraints[0].dual_value
4- Programmation entière [*]¶
Les solutions obtenues précédemment sont des solutions continues, alors que l'on commande généralement un nombre entier d'items. Les problèmes d'optimisation avec variables entières appellés (Mixed) Integer Programs, se résolvent avec des méthodes et des solveurs différents. CVXPY appellant différents solveurs nativement, cette opération est transparente ici.
Les différences pour la programmation en nombre entier¶
le type de la variable qui doit être précisé lors de la déclaration à CVXPY
x = cp.Variable(n,integer=True)le solver (qui doit etre capable de faire des MIP) (mais c'est transparent pour nous ici)
Question: Reprendre le problème relaxé plus haut (minimisation des calories entre 80% et 100% des contraintes de nutrition) en ajoutant que la variable $x$ est maintenant à valeur entières. Le problème est-il toujours faisable?
x = cp.Variable(n,integer=True)
cal = np.array(foods["Cal"]) # calories
diet_int = cp.Minimize(cal@x)
# Constraints matrix
A = np.array(foods.drop(columns=["Food","Cal"])).T
b = nutr.values.T.reshape((-1,))
contrainte1 = x>= 0.0
contrainte2a = A @ x >= 0.8*b
contrainte2b = A @ x <= b
prob = cp.Problem(diet_int , constraints = [contrainte1, contrainte2a, contrainte2b])
prob.solve(verbose=True)
Relaxation des contraintes¶
Les problemes MIP (Mixed Integer Programs) sont en général beaucoup plus dur car trouver des points admissibles est difficile.
Question : Relaxez des contraintes jusqu'a obtenir un problème faisable et observez la solution.
contrainte1a = x>=0
contrainte1b = x<=2
contrainte2a = A @ x >= 0.8*b
contrainte2b = A @ x <= 1.5*b
prob2 = cp.Problem(diet_int , constraints = [contrainte1a, contrainte1b, contrainte2a, contrainte2b])
prob2.solve(verbose=True)
prob2.status
x.value
xoptInt = x.value
for j in range(n):
if xoptInt[j]>0.0001:
print("[{:3d}]{} \t ({} cal) \t : \t {}".format(j,foods.iloc[j,0],foods.iloc[j,1],xoptInt[j]))