SciPy (pronounced sai pai) is a numpy-based math package that also includes C and Fortran libraries. With SciPy, an interactive Python session becomes a complete data processing environment like MATLAB, IDL, Octave, R, or SciLab.
In this article, we will look at the basic techniques of mathematical programming - solving conditional optimization problems for a scalar function of several variables using the scipy.optimize package. Unconstrained optimization algorithms have already been discussed in . More detailed and up-to-date help on scipy functions can always be obtained using the help() command, Shift+Tab, or in .
Introduction
A common interface for solving both conditional and unconditional optimization problems in the scipy.optimize package is provided by the function minimize(). However, it is known that there is no universal way to solve all problems, so the choice of an adequate method, as always, falls on the shoulders of the researcher.
A suitable optimization algorithm is specified using the function argument minimize(..., method="").
For conditional optimization of a function of several variables, implementations of the following methods are available:
trust-constrβ search for a local minimum in the confidence region. , ;SLSQPβ Sequential quadratic programming with restrictions, Newtonian method for solving the Lagrange system. .TNC- Truncated Newton Constrained, limited number of iterations, good for non-linear functions with a large number of independent variables. .L-BFGS-B- a method from the Broyden-Fletcher-Goldfarb-Shanno quartet, implemented with reduced memory consumption due to partial loading of vectors from the Hessian matrix. , .COBYLA- MARE Constrained Optimization By Linear Approximation, limited optimization with linear approximation (no gradient calculation). .
Depending on the chosen method, the conditions and restrictions for solving the problem are set differently:
- class object
Boundsfor methods L-BFGS-B, TNC, SLSQP, trust-constr; - list
(min, max)for the same methods L-BFGS-B, TNC, SLSQP, trust-constr; - object or list of objects
LinearConstraint,NonlinearConstraintfor COBYLA, SLSQP, trust-constr methods; - dictionary or list of dictionaries
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}for COBYLA, SLQP methods.
Outline of the article:
1) Consider applying the conditional optimization algorithm in the trust region (method="trust-constr") with constraints specified as objects Bounds, LinearConstraint, NonlinearConstraint ;
2) Consider sequential least squares programming (method="SLSQP") with constraints specified as a dictionary {'type', 'fun', 'jac', 'args'};
3) Analyze an example of product optimization on the example of a web studio.
Conditional optimization method="trust-constr"
Method Implementation trust-constr based on for problems with constraints of the form of equality and on for problems with constraints in the form of inequalities. Both methods are implemented by algorithms for finding a local minimum in the confidence region and are well suited for large-scale problems.
Mathematical formulation of the problem of finding a minimum in general form:
For strict equality constraints, the lower bound is set equal to the upper .
For a one-sided constraint, the upper or lower bound is set np.inf with the corresponding sign.
Let it be necessary to find the minimum of the known Rosenbrock function in two variables:
In this case, the following restrictions on its domain of definition are set:
In our case, there is a unique solution at the point , for which only the first and fourth constraints are valid.
Let's go through the restrictions from the bottom up and see how they can be written in scipy.
Restrictions ΠΈ define using the Bounds object.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Restrictions ΠΈ we write in linear form:
Let's define these constraints as a LinearConstraint object:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])And finally, the nonlinear constraint in matrix form:
We define the Jacobian matrix for this constraint and a linear combination of the Hessian matrix with an arbitrary vector :
Now we can define a nonlinear constraint as an object 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)If the size is large, matrices can also be specified in sparse form:
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)or as an object 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)When calculating the Hessian matrix is expensive, you can use the class . The following strategies are available: BFGS ΠΈ SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())The Hessian can also be calculated using finite differences:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')The Jacobi matrix for constraints can also be computed using finite differences. However, in this case, the Hessian matrix cannot be calculated using finite differences. The Hessian must be defined as a function or using the HessianUpdateStrategy class.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())The solution to the optimization problem is as follows:
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]If necessary, the Hessian calculation function can be defined using the LinearOperator class
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)or the product of the Hessian and an arbitrary vector through the parameter 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)Alternatively, the first and second derivatives of the function being optimized can be calculated approximately. For example, the Hessian can be approximated using the function SR1 (quasi-Newtonian approximation). The gradient can be approximated by finite differences.
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)Conditional optimization method="SLQP"
The SLSQP method is designed to solve problems of function minimization in the form:
Where ΠΈ β sets of indices of expressions describing constraints in the form of equalities or inequalities. are the sets of lower and upper bounds for the domain of the function.
Linear and non-linear constraints are described as dictionaries with keys 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])
}The search for the minimum is carried out as follows:
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 ]Optimization example
In connection with the transition to the fifth technological order, let's consider the optimization of production using the example of a web studio, which brings us a small but stable income. Imagine yourself as the director of a galley that produces three types of products:
- x0 β selling landings, from 10 tr.
- x1 - corporate sites, from 20 tr.
- x2 - online stores, from 30 tr.
Our friendly working team includes four juniors, two middles and one senior. Monthly working time fund:
- Junes:
4 * 150 = 600 ΡΠ΅Π» * ΡΠ°Ρ, - middles:
2 * 150 = 300 ΡΠ΅Π» * ΡΠ°Ρ, - senior:
150 ΡΠ΅Π» * ΡΠ°Ρ.
Suppose that the first junior should spend (0, 1, 2) hours on the development and deployment of one site of type (x10, x20, x30), middle - (7, 15, 20), senior - (5, 10, 15) hours of the best the time of your life.
Like any normal director, we want to maximize our monthly profit. The first step to success is writing the objective function value as the sum of income from products produced per month:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]This is not an error; when searching for the maximum, the objective function is minimized with the opposite sign.
The next step is to prohibit our employees from overworking and introduce restrictions on the working time fund:
Which is equivalent to:
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]])
}Formal restriction - output must be only positive:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])And finally, the most optimistic assumption - because of the low price and high quality, we are constantly lined up with satisfied customers. We can choose the monthly production volumes ourselves, based on the solution of the problem of conditional optimization with 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]Non-strictly round up to integers and calculate the monthly load of rowers with the optimal distribution of products x = (8, 6, 3) :
- Junes:
8 * 10 + 6 * 20 + 3 * 30 = 290 ΡΠ΅Π» * ΡΠ°Ρ; - middles:
8 * 7 + 6 * 15 + 3 * 20 = 206 ΡΠ΅Π» * ΡΠ°Ρ; - senior:
8 * 5 + 6 * 10 + 3 * 15 = 145 ΡΠ΅Π» * ΡΠ°Ρ.
Conclusion: in order for the director to receive his well-deserved maximum, it is optimal to make 8 landing pages, 6 medium sites and 3 stores per month. At the same time, the senior should plow without looking up from the machine, the loading of the middles will be approximately 2/3, less than half of the juniors.
Conclusion
The article outlines the basic methods of working with the package scipy.optimizeused to solve conditional minimization problems. Personally I use scipy purely for academic purposes, which is why the example given is such a joke.
A lot of theory and winrar examples can be found, for example, in the book by I. L. Akulich "Mathematical Programming in Examples and Problems". More hardcore application scipy.optimize to build a 3D structure from a set of images () can be viewed in .
The main source of information is those wishing to contribute to the translation of this and other sections scipy Welcome to .
Thank you for contributing to the publication.
Source: habr.com