科学，有条件的优化

SciPy（发音为sai pie）是一个基于numpy的数学软件包，其中还包括C和Fortran库。 使用SciPy，交互式Python会话将变成一个完整的数据处理环境，例如MATLAB，IDL，Octave，R或SciLab。

在本文中，我们将考虑数学编程的基本技术-使用scipy.optimize包解决多个变量的标量函数的条件优化问题。 上一篇文章已经考虑了无条件优化算法。 始终可以使用help（）命令，Shift + Tab或在官方文档中获得有关scipy功能的更详细和最新的帮助。

引言

scipy.optimize包中用于解决有条件和无条件优化问题的通用接口由minimal minimize()函数提供。 但是，众所周知，不存在解决所有问题的通用方法，因此，一如既往，研究人员应选择适当的方法。

使用minimize(..., method="")函数的参数指定合适的优化算法。

为了有条件地优化几个变量的功能，可以使用以下方法：

• trust-constr在置信区域中搜索局部最小值。 Wiki 文章 ； Habr文章 ；
• SLSQP具有约束的顺序二次规划，牛顿方法求解拉格朗日系统。 维基文章 。
• TNC截断牛顿约束，迭代次数有限，非常适合具有大量自变量的非线性函数。 维基文章 。
• L-BFGS-B四个Broyden-Fletcher-Goldfarb-Shanno中的一种方法，由于部分加载了Hessian矩阵中的向量，因此实现了减少内存消耗的目的。 Wiki 文章，Habr文章 。
• COBYLA - 母马 通过线性近似进行约束优化，使用线性近似进行有限优化（不计算梯度）。 维基文章 。

根据选择的方法，解决问题的条件和限制设置不同：

• 方法L-BFGS-B，TNC，SLSQP，trust-constr的Bounds类的对象；
• 相同方法L-BFGS-B，TNC，SLSQP，trust-constr的列表(min, max) ；
• 方法COBYLA，SLSQP，trust- LinearConstraint对象或对象列表LinearConstraint ， NonlinearConstraint ；
• 字典或字典{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}用于COBYLA，SLSQP方法。

文章概述：
1）考虑条件优化算法在置信区域（方法=“ trust-constr”）中的应用，该约束以对象Bounds ， LinearConstraint ， NonlinearConstraint的形式指定；
2）考虑以字典{'type', 'fun', 'jac', 'args'}形式指定的约束的顺序最小二乘编程（方法=“ SLSQP”）；
3）在网络工作室的示例上分解产品优化的示例。

条件优化方法=“ trust-constr”

trust-constr constr trust-constr的实现基于EQSQP来解决具有相等形式约束的问题，并基于TRIP来解决具有不平等形式约束的问题。 两种方法均由用于在置信区域内找到局部最小值的算法实现，非常适合于大规模任务。

一般形式的最小搜索问题的数学公式：

$$

$$\ min_x f（x）$$

$$

$$c ^ l \ leq c（x）\ leq c ^ u，$$

$$

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

对于严格的相等约束，将下限设置为等于上限 $$$c ^ l_j = c ^ u_j$ 。
对于单向限制，可通过np.inf设置带有相应符号的上限或下限。
让我们有必要找到两个变量的已知Rosenbrock函数的最小值：

$$

$$\ min_ {x_0，x_1} 100（x_0-x_1 ^ 2）^ 2 +（1-x_0）^ 2$$

在这种情况下，对其定义域设置以下限制：

$$

$$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$$

在我们的案例中，有一个独特的解决方案 $$$[x_0，x_1] = [0.4149，0.1701]$ 仅第一和第四限制有效。
让我们从头开始研究这些限制，并考虑如何用scipy编写它们。
局限性 $$$0 \ leq x_0 \ leq 1$ 和 $$$-0.5 \ leq x_1 \ leq 2.0$ 使用Bounds对象定义。

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

局限性 $$$x_0 + 2 x_1 \ leq 1$ 和 $$$2 x_0 + x_1 = 1$ 我们以线性形式写：

$$

$$\开始{bmatrix}-\ infty \\ 1 \结束{bmatrix} \ leq \开始{bmatrix} 1和2 \\ 2＆1 \结束{bmatrix} \开始{bmatrix} x_0 \\ x_1 \结束{bmatrix } \ leq \开始{bmatrix} 1 \\ 1 \结束{bmatrix}$$

将这些约束定义为LinearConstraint对象：

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

最后，矩阵形式的非线性约束：

$$

$$c（x）= \开始{bmatrix} x_0 ^ 2 + x_1 \\ x_0 ^ 2-x_1 \结束{bmatrix} \ leq \开始{bmatrix} 1 \\ 1 \结束{bmatrix}$$

我们为此限制定义了Jacobi矩阵，并定义了Hessian矩阵与任意矢量的线性组合 $$$v$ ：

$$

$$J（x）= \开始{bmatrix} 2x_0和1 \\ 2x_0＆-1 \结束{bmatrix}$$

$$

$$H（x，v）= \ sum \ limits_ {i = 0} ^ 1 v_i \ nabla ^ 2 c_i（x）= v_0 \开始{bmatrix} 2＆0 \\ 2＆0 \ end {bmatrix} + v_1 \开始{bmatrix} 2＆0 \\ 2＆0 \结束{bmatrix}$$

现在我们可以将非线性约束定义为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) 

如果大小很大，也可以用稀疏形式指定矩阵：

 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) 

或作为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) 

什么时候计算黑森州矩阵 $$$H（x，v）$ 成本HessianUpdateStrategy ，您可以使用HessianUpdateStrategy类。 可以使用以下策略： BFGS和SR1 。

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

粗麻布也可以使用有限差分来计算：

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

约束的雅可比矩阵也可以使用有限差分来计算。 但是，在这种情况下，不能使用有限差分来计算Hessian矩阵。 必须将Hessian定义为函数或使用HessianUpdateStrategy类。

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

优化问题的解决方案如下：

 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] 

如有必要，可以使用LinearOperator类定义Hessian计算函数

 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) 

或Hessian与通过hessp参数的任意向量的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) 

可替代地，可以近似地计算优化函数的一阶和二阶导数。 例如，可以使用函数SR1来近似Hessian（拟牛顿近似）。 梯度可以通过有限差来近似。

 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) 

条件优化方法=“ SLSQP”

SLSQP方法旨在解决以下形式的功能最小化问题：

$$

$$\ 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$$

哪里 $$$\数学{E}$ 和 $$$\数学{I}$ -表达索引集，以相等或不相等的形式描述约束。 $$$[lb_i，ub_i]$ -函数定义域的上下限集。

线性和非线性约束以字典的形式描述，键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]) } 

最小搜索如下：

 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 ] 

优化示例

与过渡到第五种技术模式有关，让我们以网络工作室为例考虑生产优化，这为我们带来了少量但稳定的收入。 让我们自我介绍一下画廊总监，该画廊生产三种类型的产品：

• x0-着陆点，从10 tr起
• x1-公司网站，从20 tr
• x2-网上商店，从30 tr

我们友好的工作团队包括四名Juns，两名中级和一名高级人员。 他们一个月的工作时间经费：

• 6月： 4 * 150 = 600 * ，
• 中间人： 2 * 150 = 300 * ，
• 高级： 150 *

假设第一位花在开发和部署（x0，x1，x2）类型站点上的初级人员应该花（10、20、30）小时，中级-（7、15、20），高级-（5、10、15）小时他一生的时间。

像任何普通董事一样，我们希望最大化我们的每月利润。 成功的第一步是将目标函数value为每月生产的收入之和：

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

这不是错误；搜索最大值时，目标函数以相反的符号最小化。

下一步是禁止对我们的员工进行处理，并对工作时间基金实行限制：

$$

$$\开始{bmatrix} 10＆20＆30 \\ 7＆15＆20 \\ 5＆10＆15 \结束{bmatrix} \开始{bmatrix} x_0 \\ x_1 \\ x_2 \结束{bmatrix} \ leq \开始{bmatrix} 600 \\ 300 \\ 150 \结束{bmatrix}$$

这是等效的：

$$

$$\开始{bmatrix} 600 \\ 300 \\ 150 \结束{bmatrix}-\开始{bmatrix} 10＆20＆30 \\ 7＆15＆20 \\ 5＆10＆15 \结束{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]]) } 

正式限制-输出应仅为正数：

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

最后，最乐观的假设-由于价格低廉，质量高，我们一直在排队满足客户的需求。 我们可以根据使用scipy.optimize的条件优化问题的解决方案， 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] 

我们将棒棒糖四舍五入到最接近的整数，并计算出最佳产量x = (8, 6, 3) 8，6，3）的赛艇运动员的月负荷：

• 6月： 8 * 10 + 6 * 20 + 3 * 30 = 290 * ；
• 中间人： 8 * 7 + 6 * 15 + 3 * 20 = 206 * ；
• 高级： 8 * 5 + 6 * 10 + 3 * 15 = 145 * 。

结论：要使导演获得应有的最高回报，最佳的做法是每月进行8次登陆，6个中型站点和3家商店。 同时，主角必须在不脱离机器的情况下犁耕，中间负载约为2/3，不到六月的一半。

结论

本文概述了用于解决条件最小化问题的scipy.optimize软件包的基本技术。 就我个人而言，我纯粹是出于学术目的使用scipy ，因此此示例是如此可笑。

可以在I. L. Akulich的书“示例和问题中的数学编程”一书中找到许多理论和W​​inrar示例。 可以在scipy-cookbook中找到scipy.optimize的更scipy.optimize应用程序， scipy.optimize一组图像（ 集线器上的文章 ）构建3D结构。

信息的主要来源是docs.scipy.org。如果要将scipy的本部分和其他部分翻译成翻译，欢迎使用GitHub 。

感谢mephistopheies对本出版物的贡献。