🥇SciPy，条件优化

SciPy（发音为 sai pie）是一个基于 numpy 的数学包，还包括 C 和 Fortran 库。 SciPy 将您的交互式 Python 会话转变为完整的数据科学环境，例如 MATLAB、IDL、Octave、R 或 SciLab。

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

介绍

该函数提供了用于解决 scipy.optimize 包中的条件和无约束优化问题的通用接口 minimize()。然而，众所周知，没有解决所有问题的通用方法，因此选择适当的方法一如既往地落在研究人员的肩上。
使用函数参数指定适当的优化算法 minimize(..., method="").
对于多个变量的函数的条件优化，可以使用以下方法的实现：

trust-constr — 在置信区域中搜索局部最小值。维基文章, 关于哈布雷的文章;
SLSQP — 带约束的顺序二次规划，求解拉格朗日系统的牛顿法。维基文章.
TNC - 截断牛顿约束，迭代次数有限，适合具有大量自变量的非线性函数。维基文章.
L-BFGS-B — Broyden–Fletcher–Goldfarb–Shanno 团队的一种方法，由于从 Hessian 矩阵中部分加载向量，因此减少了内存消耗。维基文章, 关于哈布雷的文章.
COBYLA — MARE Constrained Optimization By Linear Approximation，使用线性近似的约束优化（无梯度计算）。维基文章.

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

类对象 Bounds 对于方法 L-BFGS-B、TNC、SLSQP、trust-constr；
名单 (min, max) 对于相同的方法 L-BFGS-B、TNC、SLSQP、trust-constr；
一个对象或对象列表 LinearConstraint, NonlinearConstraint 对于 COBYLA、SLSQP、信任构造方法；
字典或字典列表 {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} 对于 COBYLA、SLSQP 方法。

文章大纲：
1) 考虑在信任区域（method=”trust-constr”）中使用条件优化算法，并将约束指定为对象 Bounds, LinearConstraint, NonlinearConstraint ;
2) 考虑使用最小二乘法（method =“SLSQP”）进行顺序编程，并以字典的形式指定限制 {'type', 'fun', 'jac', 'args'};
3) 使用网络工作室的示例来分析制造产品优化的示例。

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

方法的实施 trust-constr 基于 EQSQP 对于具有平等形式和约束条件的问题 TRIP 针对不平等形式的约束问题。这两种方法都是通过在置信区域中寻找局部最小值的算法来实现的，并且非常适合大规模问题。

求最小值问题的一般形式的数学表述：

对于严格的相等约束，下限设置为等于上限 .
对于单向约束，设置上限或下限 np.inf 并带有相应的标志。
需要找到两个变量的已知 Rosenbrock 函数的最小值：

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

在我们的例子中，此时有一个独特的解决方案，仅第一和第四限制有效。
让我们从下到上看一下这些限制，看看如何在 scipy 中编写它们。
限制 и 让我们使用 Bounds 对象来定义它。

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

限制 и 我们把它写成线性形式：

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

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

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

我们为此约束定义雅可比矩阵以及 Hessian 矩阵与任意向量的线性组合 :

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

计算Hessian矩阵时需要很多努力，你可以使用一个类 HessianUpdateStrategy。可以采用以下策略： BFGS и SR1.

from scipy.optimize import BFGS

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

Hessian 矩阵也可以使用有限差分来计算：

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:

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)

或者，可以近似正在优化的函数的一阶和二阶导数。例如，Hessian 矩阵可以使用以下函数进行近似 SR1 （准牛顿近似）。梯度可以通过有限差分来近似。

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 方法旨在解决以下形式的函数最小化问题：

哪里 и ——以等式或不等式形式描述限制的表达式索引集。 — 函数定义域的下限和上限集。

线性和非线性约束以带有键的字典的形式描述 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 起。

我们友好的工作团队包括四名初级人员、两名中级人员和一名高级人员。他们每月的工作时间基金：

六月： 4 * 150 = 600 чел * час,
中音： 2 * 150 = 300 чел * час,
先生： 150 чел * час.

让第一个可用的初级人员花费 (0, 1, 2) 小时来开发和部署类型为 (x10, x20, x30) 的一个站点，中级人员 - (7, 15, 20)，高级人员 - (5, 10, 15) ) 一生中最美好的时光。

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

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

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

下一步是禁止我们的员工过度工作并限制工作时间：

什么是等价的：

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:

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 * 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 所著的《示例和问题中的数学编程》一书中可以找到大量理论和虚拟示例。更多硬核应用 scipy.optimize 从一组图像构建 3D 结构（关于哈布雷的文章) 可以在以下位置查看 scipy 食谱.

主要信息来源是 docs.scipy.org希望为本节和其他部分的翻译做出贡献的人 scipy 欢迎来到 GitHub上.

谢谢梅菲斯托菲斯参与出版物的编写。

来源： habr.com

SciPy，条件优化

介绍

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

条件优化方法=“SLSQP”

优化示例

结论

添加评论 取消回复

添加评论取消回复