SciPy, optimisation sous conditions

SciPy (prononcé sai pie) est un package mathématique basé sur numpy qui comprend également les bibliothèques C et Fortran. Avec SciPy, une session Python interactive se transforme en un environnement de traitement de données complet comme MATLAB, IDL, Octave, R ou SciLab.

Dans cet article, nous considérerons les techniques de base de la programmation mathématique - résoudre les problèmes d'optimisation conditionnelle pour une fonction scalaire de plusieurs variables à l'aide du package scipy.optimize. Les algorithmes d'optimisation inconditionnels sont déjà considérés dans le dernier article . Une aide plus détaillée et à jour sur les fonctions scipy peut toujours être obtenue en utilisant la commande help (), Maj + Tab, ou dans la documentation officielle .

Présentation

L'interface générale pour résoudre les problèmes d'optimisation conditionnelle et 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. Par conséquent, comme toujours, le choix d'une méthode adéquate appartient au chercheur.

Un algorithme d'optimisation approprié est spécifié à l'aide de l'argument de la fonction minimize(..., method="") .

Pour l'optimisation conditionnelle de la fonction de plusieurs variables, les méthodes suivantes sont disponibles:

• trust-constr - recherche un minimum local dans la zone de confiance. Article Wiki, article Habr ;
• SLSQP - programmation quadratique séquentielle avec contraintes, la méthode newtonienne de résolution du système de Lagrange. Article sur le wiki .
• TNC - Newton tronqué contraint, un nombre limité d'itérations, est bon pour les fonctions non linéaires avec un grand nombre de variables indépendantes. Article sur le wiki .
• L-BFGS-B - une méthode des quatre Broyden - Fletcher - Goldfarb - Shanno, mise en œuvre avec une consommation de mémoire réduite en raison du chargement partiel des vecteurs de la matrice de Hesse. Article Wiki, article Habr .
• COBYLA - MARES Optimisation contrainte par approximation linéaire, optimisation limitée avec approximation linéaire (sans calcul de gradient). Article sur le wiki .

Selon la méthode choisie, les conditions et restrictions de résolution du problème sont définies différemment:

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

Le plan de l'article:
1) Envisager l'application de l'algorithme d'optimisation conditionnelle dans la zone de confiance (méthode = "trust-constr") avec les restrictions spécifiées sous la forme d'objets Bounds , LinearConstraint , NonlinearConstraint ;
2) Considérons la programmation séquentielle des moindres carrés (méthode = "SLSQP") avec les contraintes spécifiées sous la forme d'un dictionnaire {'type', 'fun', 'jac', 'args'} ;
3) Démonter un exemple d'optimisation de produits sur l'exemple d'un studio web.

Méthode d'optimisation conditionnelle = "trust-constr"

La mise en œuvre de la trust-constr basée sur EQSQP pour les problèmes de restrictions de forme d'égalité et sur TRIP pour les problèmes de restrictions sous forme d'inégalités. Les deux méthodes sont mises en œuvre par des algorithmes pour trouver un minimum local dans la zone de confiance et conviennent bien aux tâches à grande échelle.

La formulation mathématique du problème de recherche minimum sous la forme générale:

$$

$$\ min_x f (x)$$

$$

$$c ^ l \ leq c (x) \ leq c ^ u,$$

$$

$$x ^ l \ leq x \ leq x ^ u$$

Pour les contraintes d'égalité strictes, la borne inférieure est définie comme étant la limite supérieure $$$c ^ l_j = c ^ u_j$ .
Pour une restriction à sens unique, une limite supérieure ou inférieure est définie par np.inf avec le signe correspondant.
Soit qu'il soit nécessaire de trouver le minimum de la fonction de Rosenbrock connue de deux variables:

$$

$$\ min_ {x_0, x_1} 100 (x_0 - x_1 ^ 2) ^ 2 + (1-x_0) ^ 2$$

Dans ce cas, les restrictions suivantes sont définies sur son domaine de définition:

$$

$$x_0 ^ 2 + x_1 \ leq 1$$

$$

$$x_0 ^ 2 - x_1 \ leq 1$$

$$

$$2x_0 + x_1 = 1$$

$$

$$x_0 + 2x_1 \ leq 1$$

$$

$$0 \ leq x_0 \ leq 1$$

$$

$$-0,5 \ leq x_1 \ leq 2.0$$

Dans notre cas, il existe une solution unique au point $$$[x_0, x_1] = [0.4149, 0.1701]$ pour lesquelles seules les première et quatrième restrictions sont valables.
Passons en revue les restrictions de bas en haut et examinons comment elles peuvent être écrites en scipy.
Limitations $$$0 \ leq x_0 \ leq 1$ et $$$-0,5 \ leq x_1 \ leq 2.0$ défini à l'aide de l'objet Bounds.

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

Limitations $$$x_0 + 2 x_1 \ leq 1$ et $$ on écrit sous forme linéaire:

$$

$$\ begin {bmatrix} - \ infty \\ 1 \ end {bmatrix} \ leq \ begin {bmatrix} 1 & 2 \\ 2 & 1 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \ end {bmatrix } \ leq \ begin {bmatrix} 1 \\ 1 \ end {bmatrix}$$

Définissez ces contraintes en tant qu'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:

$$

$$c (x) = \ begin {bmatrix} x_0 ^ 2 + x_1 \\ x_0 ^ 2 - x_1 \ end {bmatrix} \ leq \ begin {bmatrix} 1 \\ 1 \ end {bmatrix}$$

Nous définissons la matrice de Jacobi pour cette restriction et la combinaison linéaire de la matrice de Hesse avec un vecteur arbitraire $$$v$ :

$$

$$J (x) = \ begin {bmatrix} 2x_0 & 1 \\ 2x_0 & -1 \ end {bmatrix}$$

$$

$$H (x, v) = \ sum \ limits_ {i = 0} ^ 1 v_i \ nabla ^ 2 c_i (x) = v_0 \ begin {bmatrix} 2 & 0 \\ 2 & 0 \ end {bmatrix} + v_1 \ begin {bmatrix} 2 & 0 \\ 2 & 0 \ end {bmatrix}$$

Maintenant, nous pouvons définir une contrainte non linéaire comme un 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 une 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 en tant qu'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 de Hesse $$$H (x, v)$ coûteux, 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()) 

La toile de jute peut également être calculée en utilisant des différences finies:

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

La matrice de Jacobi pour les contraintes peut également être calculée en utilisant des différences finies. Cependant, dans ce cas, la matrice de Hesse ne peut pas être calculée en utilisant des différences finies. Hessian doit être défini en tant que fonction ou en utilisant la classe HessianUpdateStrategy.

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

La solution au 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` 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 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 de la Hesse et un vecteur arbitraire via 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ère et deuxième dérivées de la fonction optimisée peuvent être calculées approximativement. Par exemple, la Hesse peut être approximée en utilisant la fonction SR1 (approximation quasi-newtonienne). 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) 

Méthode d'optimisation conditionnelle = "SLSQP"

La méthode SLSQP est conçue pour résoudre le problème de la minimisation de la fonction sous la forme de:

$$

$$\ min_x f (x)$$

$$

$$c_j (x) = 0, j \ in \ mathcal {E}$$

$$

$$c_j (x) \ geq 0, j \ in \ mathcal {I}$$

$$

$$lb_i \ leq x_i \ leq ub_i, i = 1, ..., N$$

O Where $$$\ mathcal {E}$ et $$$\ mathcal {I}$ - ensembles d'indices d'expression décrivant les contraintes sous forme d'égalités ou d'inégalités. $$$[lb_i, ub_i]$ - ensembles de bornes inférieure et supérieure pour le domaine de définition de fonction.

Les contraintes linéaires et non linéaires sont décrites sous forme de dictionnaires avec le type clés, 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 minimale est effectuée 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) 

 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

Dans le cadre du passage au cinquième mode technologique, considérons l'optimisation de la production à l'aide de l'exemple d'un studio web, ce qui nous apporte un revenu modeste mais stable. Présentons-nous en tant que directeur d'une galerie, qui produit trois types de produits:

• x0 - vente de débarquements, à partir de 10 tr
• x1 - sites corporatifs, à partir de 20 tr
• x2 - boutiques en ligne, à partir de 30 tr

Notre équipe de travail amicale comprend quatre Juns, deux intermédiaires et un senior. Fonds de leur temps de travail pendant un mois:

• Juin: 4 * 150 = 600 * ,
• intermédiaires: 2 * 150 = 300 * ,
• Senior: 150 *

Supposons que le premier junior à consacrer au développement et au déploiement d'un site de type (x0, x1, x2) devrait passer (10, 20, 30) heures, moyen - (7, 15, 20), senior - (5, 10, 15) heures des meilleurs moment de sa vie.

Comme tout administrateur normal, nous voulons maximiser notre bénéfice mensuel. La première étape du succès consiste à noter la value de la fonction objective comme la somme des revenus 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 le signe opposé.

La prochaine étape consiste à interdire le traitement à nos employés et à introduire des restrictions sur le fonds de temps de travail:

$$

$$\ begin {bmatrix} 10 & 20 & 30 \\ 7 & 15 & 20 \\ 5 & 10 & 15 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \\ x_2 \ end {bmatrix} \ leq \ begin {bmatrix} 600 \\ 300 \\ 150 \ end {bmatrix}$$

Ce qui est équivalent:

$$

$$\ begin {bmatrix} 600 \\ 300 \\ 150 \ end {bmatrix} - \ begin {bmatrix} 10 & 20 & 30 \\ 7 & 15 & 20 \\ 5 & 10 & 15 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \\ x_2 \ end {bmatrix} \ geq 0$$

 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]]) } 

Restriction formelle - le résultat ne doit être que positif:

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

Et enfin, l'hypothèse la plus optimiste - en raison du prix bas et de la haute qualité, une ligne de clients satisfaits fait constamment la queue pour nous. Nous pouvons choisir nous-mêmes les volumes de production mensuels, en fonction de la solution du 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] 

Nous arrondissons la sucette à l'entier le plus proche et calculons la charge mensuelle des rameurs avec une distribution optimale de la production x = (8, 6, 3) :

• Juin: 8 * 10 + 6 * 20 + 3 * 30 = 290 * ;
• intermédiaires: 8 * 7 + 6 * 15 + 3 * 20 = 206 * ;
• Senior: 8 * 5 + 6 * 10 + 3 * 15 = 145 * .

Conclusion: pour que le réalisateur obtienne son maximum bien mérité, il est optimal de faire 8 débarquements, 6 sites moyens et 3 magasins par mois. En même temps, le seigneur doit labourer sans se détacher de l'engin, la charge de milieu est d'environ 2/3, moins de la moitié du mois de juin.

Conclusion

L'article décrit les techniques de base pour travailler avec le package scipy.optimize qui sont utilisées pour résoudre les problèmes de minimisation conditionnelle. Personnellement, j'utilise scipy uniquement à des fins académiques, donc cet exemple est tellement comique.

De nombreux exemples de théorie et de winrar peuvent être trouvés, par exemple, dans le livre de I. L. Akulich "Programmation mathématique dans les exemples et les problèmes". Une application plus hardcore de scipy.optimize pour construire une structure 3D pour un ensemble d'images ( un article sur le hub ) peut être trouvée dans le livre de recettes scipy .

La principale source d'informations est docs.scipy.org. Si vous voulez traduire ceci et d'autres sections de scipy en traduction, bienvenue sur GitHub .

Merci aux méphistophées d' avoir contribué à cette publication.