SciPy (prononcé sai pie) est un package mathématique basé sur numpy qui comprend également les bibliothèques C et Fortran. SciPy transforme votre session Python interactive en un environnement complet de science des données comme MATLAB, IDL, Octave, R ou SciLab.
Dans cet article, nous examinerons les techniques de base de la programmation mathématique - résoudre des problèmes d'optimisation conditionnelle pour une fonction scalaire de plusieurs variables à l'aide du package scipy.optimize. Les algorithmes d'optimisation sans contraintes ont déjà été discutés dans . Une aide plus détaillée et à jour sur les fonctions scipy peut toujours être obtenue en utilisant la commande help(), Shift+Tab ou dans .
introduction
Une interface commune pour résoudre les problèmes d'optimisation conditionnelle et sans contrainte dans le package scipy.optimize est fournie par la fonction minimize(). Cependant, on sait qu'il n'existe pas de méthode universelle pour résoudre tous les problèmes, de sorte que le choix d'une méthode adéquate incombe, comme toujours, au chercheur.
L'algorithme d'optimisation approprié est spécifié à l'aide de l'argument de fonction minimize(..., method="").
Pour l'optimisation conditionnelle d'une fonction de plusieurs variables, des implémentations des méthodes suivantes sont disponibles :
trust-constr— recherche d'un minimum local dans la région de confiance. , ;SLSQP— programmation quadratique séquentielle avec contraintes, méthode newtonienne de résolution du système de Lagrange. .TNC- Newton tronqué contraint, nombre d'itérations limité, idéal pour les fonctions non linéaires avec un grand nombre de variables indépendantes. .L-BFGS-B— une méthode de l'équipe Broyden – Fletcher – Goldfarb – Shanno, implémentée avec une consommation de mémoire réduite en raison du chargement partiel des vecteurs de la matrice hessienne. , .COBYLA— MARE Constrained Optimization By Linear Approximation, optimisation contrainte avec approximation linéaire (sans calcul de gradient). .
Selon la méthode choisie, les conditions et restrictions de résolution du problème sont fixées différemment :
- objet de classe
Boundspour les méthodes L-BFGS-B, TNC, SLSQP, trust-constr ; - liste
(min, max)pour les mêmes méthodes L-BFGS-B, TNC, SLSQP, trust-constr ; - un objet ou une liste d'objets
LinearConstraint,NonlinearConstraintpour COBYLA, SLSQP, méthodes trust-constr ; - dictionnaire ou liste de dictionnaires
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}pour les méthodes COBYLA, SLSQP.
Aperçu de l'article :
1) Envisagez l'utilisation d'un algorithme d'optimisation conditionnelle dans la région de confiance (method=”trust-constr”) avec des contraintes spécifiées comme objets Bounds, LinearConstraint, NonlinearConstraint ;
2) Envisager une programmation séquentielle par la méthode des moindres carrés (méthode = "SLSQP") avec des restrictions précisées sous forme de dictionnaire {'type', 'fun', 'jac', 'args'};
3) Analyser un exemple d'optimisation de produits fabriqués en utilisant l'exemple d'un studio web.
Méthode d'optimisation conditionnelle="trust-constr"
Mise en œuvre de la méthode trust-constr basé sur pour les problèmes avec des contraintes de forme d'égalité et sur pour les problèmes avec des contraintes sous forme d’inégalités. Les deux méthodes sont mises en œuvre par des algorithmes permettant de trouver un minimum local dans la région de confiance et conviennent bien aux problèmes à grande échelle.
Formulation mathématique du problème de la recherche d'un minimum sous 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 unidirectionnelle, la limite supérieure ou inférieure est fixée np.inf avec le signe correspondant.
Soit il faut trouver le minimum d'une fonction de Rosenbrock connue de deux variables :
Dans ce cas, les restrictions suivantes sont définies sur son domaine de définition :
Dans notre cas, il existe une solution unique au point , pour lequel seules les première et quatrième restrictions sont valables.
Passons en revue les restrictions de bas en haut et voyons comment nous pouvons les écrire dans scipy.
Restrictions и Définissons-le à l'aide de l'objet Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Restrictions и Écrivons-le sous forme linéaire :
Définissons ces contraintes comme un 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 la contrainte non linéaire sous forme matricielle :
Nous définissons la matrice jacobienne pour cette contrainte et une combinaison linéaire de la matrice hessienne avec un vecteur arbitraire :
Nous pouvons maintenant définir une contrainte non linéaire en tant qu'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 clairsemée :
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)Lors du calcul de la matrice de Hesse demande beaucoup d'efforts, vous pouvez utiliser un cours . Les stratégies suivantes sont disponibles : BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Le Hessian peut également être calculé en utilisant des différences finies :
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')La matrice jacobienne des contraintes peut également être calculée en utilisant des différences finies. Cependant, dans ce cas, la matrice hessienne ne peut pas être calculée en utilisant des différences finies. Le Hessian doit être défini comme une fonction ou à l’aide de la classe HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())La solution au problème d’optimisation ressemble à ceci :
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` termination condition is satisfied.
Number of iterations: 12, function evaluations: 8, CG iterations: 7, optimality: 2.99e-09, constraint violation: 1.11e-16, execution time: 0.033 s.
[0.41494531 0.17010937]Si nécessaire, la fonction de calcul du Hessian peut être définie à l'aide de 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 Hessian 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 dérivées première et seconde de la fonction optimisée peuvent être approchées. Par exemple, le Hessian peut être approximé à l'aide de la fonction SR1 (approximation quasi-newtonienne). Le gradient peut être approché 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)Méthode d'optimisation conditionnelle="SLSQP"
La méthode SLSQP est conçue pour résoudre des problèmes de minimisation d'une fonction sous la forme :
Où и — des ensembles d'indices d'expressions décrivant des restrictions 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 clés type, fun и 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 s'effectue de la manière suivante :
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)Optimization terminated successfully. (Exit mode 0)
Current function value: 0.34271757499419825
Iterations: 4
Function evaluations: 5
Gradient evaluations: 4
[0.41494475 0.1701105 ]Exemple d'optimisation
A propos du passage à la cinquième structure technologique, regardons l'optimisation de la production à l'aide de l'exemple d'un studio web, qui nous rapporte un revenu modeste mais stable. Imaginons-nous en tant que directeur d'une galère qui produit trois types de produits :
- x0 - vente de pages de destination, à partir de 10 tr.
- x1 - sites Web d'entreprise, à partir de 20 tr.
- x2 - boutiques en ligne, à partir de 30 tr.
Notre équipe de travail amicale comprend quatre juniors, deux intermédiaires et un senior. Leur fonds mensuel de temps de travail :
- Juin :
4 * 150 = 600 чел * час, - milieu:
2 * 150 = 300 чел * час, - sénateur :
150 чел * час.
Laissez le premier junior disponible consacrer (0, 1, 2) heures au développement et au déploiement d'un site de type (x10, x20, x30), intermédiaire - (7, 15, 20), senior - (5, 10, 15 ) heures du meilleur moment de votre vie.
Comme tout administrateur normal, nous souhaitons maximiser les profits mensuels. La première étape vers le succès est d'écrire la fonction objectif value comme le montant des revenus des produits fabriqués par mois :
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 le signe opposé.
La prochaine étape consiste à interdire à nos salariés le surmenage et à introduire des restrictions sur les horaires de travail :
Ce qui est équivalent :
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]])
}Une restriction formelle est que la production du produit ne doit être que positive :
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])Et enfin, l’hypothèse la plus optimiste est qu’en raison du prix bas et de la haute qualité, une file de clients satisfaits fait constamment la queue pour nous. Nous pouvons choisir nous-mêmes les volumes de production mensuels, en nous basant sur la résolution du problème d'optimisation contraint 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]Arrondons vaguement aux nombres entiers et calculons la charge mensuelle des rameurs avec une répartition optimale des produits x = (8, 6, 3) :
- Juin :
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - milieu:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - sénateur :
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Conclusion : pour que le réalisateur reçoive son maximum bien mérité, il est optimal de créer 8 landing pages, 6 sites de taille moyenne et 3 boutiques par mois. Dans ce cas, le senior doit labourer sans lever les yeux de la machine, la charge des intermédiaires sera d'environ 2/3, les juniors moins de la moitié.
Conclusion
L'article décrit les techniques de base pour travailler avec le package scipy.optimize, utilisé pour résoudre des problèmes de minimisation conditionnelle. Personnellement j'utilise scipy purement académique, c’est pourquoi l’exemple donné est si comique.
De nombreux exemples théoriques et virtuels peuvent être trouvés, par exemple, dans le livre de I.L. Akulich « Programmation mathématique dans les exemples et les problèmes ». Application plus hardcore scipy.optimize construire une structure 3D à partir d'un ensemble d'images () peut être consulté dans .
La principale source d'information est ceux qui souhaitent contribuer à la traduction de cette section et d'autres scipy Bienvenue à .
merci pour participer à la préparation de la publication.
Source: habr.com