SciPy, optimisation avec contraintes

SciPy (se prononce comme "sci-py") est un package mathématique basé sur numpy, incluant également des bibliothèques en C et Fortran. Avec SciPy, une session Python interactive devient une véritable plateforme de traitement de données, comparable à MATLAB, IDL, Octave, R ou SciLab.

Dans cet article, nous examinerons les principales techniques de programmation mathématique — la résolution de problèmes d'optimisation conditionnelle pour une fonction scalaire de plusieurs variables à l'aide du package scipy.optimize. Les algorithmes d'optimisation sans contrainte ont déjà été examinés dans Dans notre précédent article,. Une documentation plus détaillée et à jour sur les fonctions scipy est toujours accessible via la commande help(), Shift+Tab ou dans la documentation officielle.

Introduction

L'interface générale pour résoudre les problèmes d'optimisation tant conditionnelle qu'inconditionnelle dans le package scipy.optimize est fournie par la fonction minimize(). Cependant, il est connu qu'il n'existe pas de méthode universelle pour résoudre tous les problèmes, donc le choix de la méthode appropriée repose toujours sur les épaules du chercheur.
L'algorithme d'optimisation adapté est spécifié à l'aide de l'argument de la fonction minimize(..., method="").
Pour l'optimisation conditionnelle de fonctions de plusieurs variables, les implémentations des méthodes suivantes sont disponibles :

• trust-constr — recherche d'un minimum local dans une région de confiance. Article sur wiki, un article sur Habré;
• SLSQP — programmation quadratique séquentielle avec contraintes, méthode de Newton pour résoudre le système de Lagrange. Article sur le wiki.
• TNC — Truncated Newton Constrained, nombre d'itérations limité, bon pour les fonctions non linéaires avec un grand nombre de variables indépendantes. Article sur wiki.
• L-BFGS-B — méthode de la paire Broyden–Fletcher–Goldfarb–Shanno, mise en œuvre avec une consommation de mémoire réduite grâce au chargement partiel des vecteurs à partir de la matrice Hessienne. Article sur wiki, un article sur Habré.
• COBYLA — COBYLA Constrained Optimization By Linear Approximation, optimisation contrainte avec approximation linéaire (sans calcul de gradient). Article sur wiki.

En fonction de la méthode choisie, différentes conditions et contraintes sont définies pour résoudre le problème :

• un objet de la classe Bounds pour les méthodes L-BFGS-B, TNC, SLSQP, trust-constr ;
• une liste (min, max) pour ces mêmes méthodes L-BFGS-B, TNC, SLSQP, trust-constr ;
• un objet ou une liste d'objets LinearConstraint, NonlinearConstraint pour les méthodes COBYLA, SLSQP, trust-constr ;
• un dictionnaire ou une liste de dictionnaires {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} pour les méthodes COBYLA, SLSQP.

Plan de l'article :
1) Examiner l'application de l'algorithme d'optimisation conditionnelle dans la région de confiance (method="trust-constr") avec les contraintes spécifiées sous forme d'objets. Bounds, LinearConstraint, NonlinearConstraint ;
2) Considérer la programmation séquentielle par la méthode des moindres carrés (method=»SLSQP») avec des contraintes définies sous forme de dictionnaire {'type', 'fun', 'jac', 'args'};
3) Analyser un exemple d'optimisation de la production à travers un exemple d'une agence web.

Optimisation conditionnelle method=»trust-constr»

Mise en œuvre de la méthode trust-constr est basée sur EQSQP pour les problèmes avec des contraintes de type égalité et sur TRIP pour les problèmes avec des contraintes sous forme d'inégalités. Les deux méthodes sont mises en œuvre par des algorithmes de recherche d'un minimum local dans une région de confiance et conviennent bien aux problèmes à grande échelle.

Formulation mathématique du problème de recherche du minimum dans sa forme générale :

Pour les contraintes d'égalité stricte, la limite inférieure est fixée égale à la limite supérieure .
Pour une contrainte unilatérale, la limite supérieure ou inférieure est fixée np.inf avec le signe correspondant.
Supposons qu'il faut trouver le minimum d'une fonction de Rosenbrock à deux variables :

Dans ce cas, les contraintes suivantes sont données sur son domaine de définition :

Dans notre cas, il n'y a qu'une seule solution au point , pour lequel seules la première et la quatrième contraintes sont valides.
Examinons les contraintes de bas en haut et voyons comment les écrire dans scipy.
Restrictions et définissons avec l'aide de l'objet Bounds.

from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])

Restrictions et écrivons en forme linéaire :

Définissons ces contraintes sous forme d'objet LinearConstraint :

import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])

Et enfin, une contrainte non linéaire sous forme matricielle :

Définissons la matrice Jacobienne pour cette contrainte et la combinaison linéaire de la matrice Hessienne avec un vecteur arbitraire :

Nous pouvons désormais définir la contrainte non linéaire comme objet NonlinearConstraint:

from scipy.optimize import NonlinearConstraint

def cons_f(x):
     return [x[0]**2 + x[1], x[0]**2 - x[1]]

def cons_J(x):
     return [[2*x[0], 1], [2*x[0], -1]]

def cons_H(x, v):
     return v[0]*np.array([[2, 0], [0, 0]]) + v[1]*np.array([[2, 0], [0, 0]])

nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H)

Si la taille est grande, les matrices peuvent également être spécifiées sous forme sparse :

from scipy.sparse import csc_matrix

def cons_H_sparse(x, v):
     return v[0]*csc_matrix([[2, 0], [0, 0]]) + v[1]*csc_matrix([[2, 0], [0, 0]])

nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
                                            jac=cons_J, hess=cons_H_sparse)

ou comme objet LinearOperator:

from scipy.sparse.linalg import LinearOperator

def cons_H_linear_operator(x, v):
    def matvec(p):
        return np.array([p[0]*2*(v[0]+v[1]), 0])
    return LinearOperator((2, 2), matvec=matvec)

nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
                                jac=cons_J, hess=cons_H_linear_operator)

Lorsque le calcul de la matrice Hessienne exige des ressources importantes, vous pouvez utiliser la classe HessianUpdateStrategy. Les stratégies suivantes sont disponibles : BFGS et SR1.

from scipy.optimize import BFGS

nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())

Le Hessien peut également être calculé par des différences finies :

nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')

La matrice Jacobienne pour les contraintes peut également être calculée par des différences finies. Cependant, dans ce cas, la matrice Hessienne ne peut pas être calculée par des différences finies. Le Hessien doit être défini sous forme de fonction ou à l'aide de la classe HessianUpdateStrategy.

nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())

La solution du problème d'optimisation est la suivante :

from scipy.optimize import minimize
from scipy.optimize import rosen, rosen_der, rosen_hess, rosen_hess_prod

x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess,
                constraints=[linear_constraint, nonlinear_constraint],
                options={'verbose': 1}, bounds=bounds)
print(res.x)

`gtol` la condition de terminaison est satisfaite.
Nombre d'itérations : 12, évaluations de fonction : 8, itérations CG : 7, optimalité : 2.99e-09, violation de contrainte : 1.11e-16, temps d'exécution : 0.033 s.
[0.41494531 0.17010937]

Si nécessaire, la fonction de calcul du Hessien peut être définie par la classe LinearOperator.

def rosen_hess_linop(x):
    def matvec(p):
        return rosen_hess_prod(x, p)
    return LinearOperator((2, 2), matvec=matvec)

res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess_linop,
                 constraints=[linear_constraint, nonlinear_constraint],
                 options={'verbose': 1}, bounds=bounds)

print(res.x)

ou le produit du Hessien et d'un vecteur arbitraire à travers le paramètre hessp:

res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hessp=rosen_hess_prod,
                constraints=[linear_constraint, nonlinear_constraint],
                options={'verbose': 1}, bounds=bounds)
print(res.x)

Alternativement, les premières et secondes dérivées de la fonction optimisée peuvent être calculées de manière approchée. Par exemple, le Hessien peut être approximé par une fonction SR1 (approximation quasi-Newton). Le gradient peut être approximé par des différences finies.

from scipy.optimize import SR1
res = minimize(rosen, x0, method='trust-constr',  jac="2-point", hess=SR1(),
               constraints=[linear_constraint, nonlinear_constraint],
               options={'verbose': 1}, bounds=bounds)
print(res.x)

Optimisation conditionnelle method=»SLSQP»

La méthode SLSQP est destinée à résoudre des problèmes de minimisation de fonction sous la forme :

Où et — des ensembles d'indices d'expressions décrivant des contraintes sous forme d'égalités ou d'inégalités. — des ensembles de bornes inférieures et supérieures pour le domaine de définition de la fonction.

Les contraintes linéaires et non linéaires sont décrites sous forme de dictionnaires avec des clés type, fun et jac.

ineq_cons = {'type': 'ineq',
             'fun': lambda x: np.array ([1 - x [0] - 2 * x [1],
                                          1 - x [0] ** 2 - x [1],
                                          1 - x [0] ** 2 + x [1]]),
             'jac': lambda x: np.array ([[- 1.0, -2.0],
                                          [-2 * x [0], -1.0],
                                          [-2 * x [0], 1.0]])
            }

eq_cons = {'type': 'eq',
           'fun': lambda x: np.array ([2 * x [0] + x [1] - 1]),
           'jac': lambda x: np.array ([2.0, 1.0])
          }

La recherche du minimum se déroule comme suit :

x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='SLSQP', jac=rosen_der,
               constraints=[eq_cons, ineq_cons], options={'ftol': 1e-9, 'disp': True},
               bounds=bounds)

print(res.x)

L'optimisation s'est terminée avec succès.    (Mode de sortie 0)
            Valeur actuelle de la fonction : 0.34271757499419825
            Itérations : 4
            Évaluations de la fonction : 5
            Évaluations du gradient : 4
[0.41494475 0.1701105 ]

Exemple d'optimisation

En raison du passage au cinquième ordre technologique, examinons l'optimisation de la production à travers l'exemple d'une agence web, qui nous génère un revenu modeste mais stable. Imaginons que nous soyons le directeur d'une galerie produisant trois types de produits :

• x0 — pages de destination génératrices de ventes, à partir de 10 k.
• x1 — sites web d'entreprise, à partir de 20 k.
• x2 — boutiques en ligne, à partir de 30 k.

Notre équipe de travail soudée comprend quatre juniors, deux intermédiaires et un senior. Le temps de travail qu'ils peuvent consacrer en un mois :

• juniors : 4 * 150 = 600 homme-heures,
• intermédiaires : 2 * 150 = 300 homme-heures,
• senior : 150 homme-heures.

Supposons que pour le développement et le déploiement d'un site de type (x0, x1, x2), le premier junior doive passer (10, 20, 30) heures, l'intermédiaire — (7, 15, 20), le senior — (5, 10, 15) heures du meilleur temps de sa vie.

Comme tout directeur normal, nous souhaitons maximiser le bénéfice mensuel. La première étape vers le succès est d'écrire la fonction objectif value comme la somme des revenus provenant de la production mensuelle :

def value(x):
    return - 10*x[0] - 20*x[1] - 30*x[2]

Ce n'est pas une erreur, lors de la recherche du maximum, la fonction objectif est minimisée avec un signe inversé.

La prochaine étape consiste à interdire aux employés de faire des heures supplémentaires et à établir des limites sur le temps de travail disponible :

Ce qui équivaut à :

ineq_cons = {'type': 'ineq',
             'fun': lambda x: np.array ([600 - 10 * x [0] - 20 * x [1] - 30 * x[2],
                                         300 - 7  * x [0] - 15 * x [1] - 20 * x[2],
                                         150 - 5  * x [0] - 10 * x [1] - 15 * x[2]])
            }

La contrainte formelle est que la production doit être uniquement positive :

bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])

Et enfin, la supposition la plus optimiste est que, grâce à nos prix bas et à notre haute qualité, nous avons constamment une file d'attente de clients satisfaits. Nous pouvons choisir nous-mêmes les volumes mensuels de production en fonction de la résolution d'un problème d'optimisation conditionnelle avec scipy.optimize:

x0 = np.array([10, 10, 10])
res = minimize(value, x0, method='SLSQP', constraints=ineq_cons, bounds=bnds)
print(res.x)

[7.85714286 5.71428571 3.57142857]

Arrondissons modérément à des entiers et calculons la charge mensuelle des rameurs pour le meilleur scénario de production. x = (8, 6, 3) :

• juniors : 8 * 10 + 6 * 20 + 3 * 30 = 290 pers * heure;
• intermédiaires : 8 * 7 + 6 * 15 + 3 * 20 = 206 pers * heure;
• senior : 8 * 5 + 6 * 10 + 3 * 15 = 145 pers * heure.

Conclusion : pour que le directeur reçoive son maximum mérité, il est optimal de produire 8 landing pages, 6 sites moyens et 3 magasins par mois. Le senior doit travailler sans relâche, la charge des mid-levels sera d'environ 2/3, et celle des juniors moins de la moitié.

Conclusion

Cet article présente les principales méthodes de travail avec le paquet scipy.optimize, utilisées pour résoudre des problèmes de minimisation conditionnelle. Personnellement, j'utilise scipy purement à des fins académiques, c'est pourquoi l'exemple présenté a un caractère humoristique.

Beaucoup de théorie et des exemples rarifiés peuvent être trouvés, par exemple, dans le livre d'I. L. Akulich « Programmation mathématique en exemples et en problèmes ». Une application plus hardcore scipy.optimize pour construire une structure 3D à partir d'un ensemble d'images (un article sur Habré) peut être vue dans scipy-cookbook.

La principale source d'information est docs.scipy.org, ceux qui souhaitent contribuer à la traduction de cette section et d'autres sections scipy sont les bienvenus sur GitHub.

Merci mephistopheies pour avoir participé à la préparation de la publication.

Source : habr.com