SciPy (pronounced as sigh pie) is a mathematical package based on numpy, which also includes libraries in C and Fortran. With SciPy, an interactive Python session turns into a fully-fledged data processing environment akin to MATLAB, IDL, Octave, R, or SciLab.
In this article, we will discuss the basic techniques of mathematical programming — solving constrained optimization problems for a scalar function of several variables using the scipy.optimize package. Unconstrained optimization algorithms have already been covered in . More detailed and up-to-date documentation on scipy functions can always be accessed using the help() command, Shift+Tab, or in .
Introduction
The general interface for solving both constrained and unconstrained 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 selection of an appropriate method always falls on the researcher.
The suitable optimization algorithm is specified using the argument of the function minimize(..., method="").
For constrained optimization of a scalar function of several variables, the following methods are implemented:
trust-constr— searching for a local minimum in a trust region. , ;SLSQP— Sequential Quadratic Programming with constraints, Newton's method for solving the Lagrange system. .TNC— Truncated Newton Constrained, limited iterations, good for nonlinear functions with a large number of independent variables. .L-BFGS-B— the method of Broyden–Fletcher–Goldfarb–Shanno, implemented with reduced memory consumption through partial loading of vectors from the Hessian matrix. , .COBYLA— COBYLA Constrained Optimization By Linear Approximation, constrained optimization with linear approximation (without computing gradients). .
Depending on the chosen method, conditions and constraints for solving the problem are set differently:
- an instance of the class
Boundsfor methods L-BFGS-B, TNC, SLSQP, trust-constr; - a list of
(min, max)for the same methods L-BFGS-B, TNC, SLSQP, trust-constr; - an object or a list of objects
LinearConstraint,NonlinearConstraintfor methods COBYLA, SLSQP, trust-constr; - a dictionary or a list of dictionaries
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}for methods COBYLA, SLSQP.
Article Outline:
1) Examine the application of the constrained optimization algorithm in a trust region (method="trust-constr") with constraints specified as objects. Bounds, LinearConstraint, NonlinearConstraint ;
2) Consider sequential programming using the least squares method (method='SLSQP') with constraints defined in the form of a dictionary {'type', 'fun', 'jac', 'args'};
3) Analyze an example of optimizing product output using a web studio.
Conditional optimization method='trust-constr'
Implementation of the method trust-constr is based on for equality constraint problems and for for inequality constraint problems. Both methods are implemented as local minimum search algorithms in a trusted region and are well-suited for large-scale problems.
Mathematical formulation of the task for finding the minimum in general form:
For strict equality constraints, the lower bound is set equal to the upper bound 
.
For one-sided constraints, the upper or lower bound is set np.inf with the corresponding sign.
Suppose we need to find the minimum of the well-known Rosenbrock function of two variables:
At the same time, the following constraints are set for its domain of definition:
In our case, there is a unique solution at the point 
, for which only the first and fourth constraints are valid.
We will review the constraints from the bottom up and consider how to express them in scipy.
Restrictions 
and 
define using the Bounds object.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Restrictions 
and 
express in linear form:
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, a nonlinear constraint in matrix form:
Define the Jacobian matrix for this constraint and a linear combination of the Hessian matrix with an arbitrary vector 
:
Now we can define the 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 a 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 
requires significant resources, one can use the class The following strategies are available: BFGS and SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())The Hessian can also be computed using finite differences:
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac='2-point', hess='2-point')The Jacobian matrix for constraints can also be computed using finite differences. However, in this case, the Hessian cannot be computed via 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 function for computing the Hessian 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 via 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 optimized function can be approximated. For example, the Hessian can be approximated using a SR1 (quasi-Newton approximation). The gradient can be approximated using 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='SLSQP'
The SLSQP method is designed to solve minimization problems of the form:
Where 
and 
— a set of indexes of expressions that describe constraints in the form of equalities or inequalities. 
— a set of lower and upper bounds for the domain of the function.
Linear and nonlinear constraints are described as dictionaries with keys type, fun and 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 ]Example of optimization
In connection with the transition to the fifth technological order, let’s consider optimizing production using the example of a web studio that brings us a small but stable income. Let’s imagine ourselves as the director of a gallery producing three types of products:
- x0 — landing pages, starting from 10,000.
- x1 — corporate websites, starting from 20,000.
- x2 — online stores, starting from 30,000.
Our friendly working team consists of four juniors, two mid-levels, and one senior. Their monthly working time fund is:
- juniors:
4 * 150 = 600 man-hours, - mid-levels:
2 * 150 = 300 man-hours, - senior:
150 man-hours.
Let’s assume a junior must spend (10, 20, 30) hours on developing and deploying each type of website (x0, x1, x2), a mid-level spends (7, 15, 20), and a senior spends (5, 10, 15) hours of their best time.
Like any normal director, we want to maximize monthly profit. The first step to success is recording the objective function value as the sum of the income from the products produced over the 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 a negative sign.
The next step is to prohibit overwork for our employees and introduce constraints 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]])
}The formal constraint is that the production output must be positive:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])And finally, the brightest assumption is that due to low prices and high quality, we constantly have a line of satisfied customers. We can choose our monthly production volumes based on solving an optimization problem 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]Let's round to whole numbers and calculate the monthly workload of the rowers under the optimal production scheme. x = (8, 6, 3) :
- juniors:
8 * 10 + 6 * 20 + 3 * 30 = 290 people * hour; - mid-levels:
8 * 7 + 6 * 15 + 3 * 20 = 206 people * hour; - senior:
8 * 5 + 6 * 10 + 3 * 15 = 145 people * hour.
Conclusion: For the director to receive his well-deserved maximum, it is optimal to produce 8 landing pages, 6 medium websites, and 3 stores per month. The senior should work tirelessly at the machine, while the workload for mid-levels will be about 2/3, and for juniors less than half.
Conclusion
This article outlines the main methods for using the package scipy.optimize, used to solve constrained minimization problems. Personally, I use scipy purely for academic purposes, so the example provided is somewhat humorous in nature.
Much theory and numerous examples can be found in the book by I.L. Akulich "Mathematical Programming in Examples and Problems." A more hardcore application scipy.optimize for building a 3D structure from a set of images () can be seen in .
The main source of information is , those wishing to contribute to the translation of this and other sections scipy are welcome at .
Thank you for participating in the preparation of the publication.
Source: habr.com
