【问题标题】:Python constrained non-linear optimizationPython 约束非线性优化
【发布时间】:2014-03-13 00:19:17
【问题描述】:

在 python 中进行约束非线性优化的推荐包是什么?

我要解决的具体问题是这样的:

我有一个未知的 X (Nx1),我有 M (Nx1) u 向量和 M (NxN) s 矩阵。

max [5th percentile of (ui_T*X), i in 1 to M]
st 
0<=X<=1 and
[95th percentile of (X_T*si*X), i in 1 to M]<= constant

当我开始解决这个问题时,我只对us 有一个点估计,我能够用cvxpy 解决上述问题。

我意识到,我没有对 us 进行一个估计,而是拥有整个值的分布,因此我想更改我的目标函数,以便可以使用整个分布。上面的问题描述是我尝试以有意义的方式包含该信息。

cvxpy 不能用来解决这个问题,我试过scipy.optimize.anneal,但我似乎无法为未知值设置界限。我也看过pulp,但它不允许非线性约束。

【问题讨论】:

  • 要求我们推荐或查找工具、库或最喜欢的场外资源的问题对于 Stack Overflow 来说是题外话,因为它们往往会吸引固执己见的答案和垃圾邮件。相反,请描述问题以及迄今为止为解决该问题所做的工作。
  • 当然。当我开始解决这个问题时,我对 u 和 s 只有一个点估计,并且我能够使用 cvxpy 解决上述问题。我意识到,我没有对 u 和 s 进行一个估计,而是拥有整个值的分布,所以我想改变我的目标函数,以便我可以使用整个分布。上面的问题描述是我尝试以有意义的方式包含该信息。 cvxpy 不能用来解决这个问题,我试过 scipy.optimize.anneal,但我似乎无法为未知值设置界限。我也看过纸浆,但它不允许非线性约束。

标签: python mathematical-optimization cvxpy


【解决方案1】:

虽然scipy.optimize.minimize 中的SLSQP 算法很好,但它有很多限制。第一个是 QP 求解器,因此它适用于非常适合二次编程范式的方程。但是如果你有功能限制会发生什么?此外,scipy.optimize.minimize 不是全局优化器,因此您通常需要在非常接近最终结果的情况下开始。

有一个约束非线性优化包(称为mystic)几乎与scipy.optimize 本身一样长——我建议它作为处理任何一般约束非线性优化的首选。

例如,你的问题,如果我理解你的伪代码,看起来像这样:

import numpy as np

M = 10
N = 3
Q = 10
C = 10

# let's be lazy, and generate s and u randomly...
s = np.random.randint(-Q,Q, size=(M,N,N))
u = np.random.randint(-Q,Q, size=(M,N))

def percentile(p, x):
    x = np.sort(x)
    p = 0.01 * p * len(x)
    if int(p) != p:
        return x[int(np.floor(p))]
    p = int(p)
    return x[p:p+2].mean()

def objective(x, p=5): # inverted objective, to find the max
    return -1*percentile(p, [np.dot(np.atleast_2d(u[i]), x)[0] for i in range(0,M-1)])


def constraint(x, p=95, v=C): # 95%(xTsx) - v <= 0
    x = np.atleast_2d(x)
    return percentile(p, [np.dot(np.dot(x,s[i]),x.T)[0,0] for i in range(0,M-1)]) - v

bounds = [(0,1) for i in range(0,N)]

因此,要在mystic 中处理您的问题,您只需要指定边界和约束。

from mystic.penalty import quadratic_inequality
@quadratic_inequality(constraint, k=1e4)
def penalty(x):
  return 0.0

from mystic.solvers import diffev2
from mystic.monitors import VerboseMonitor
mon = VerboseMonitor(10)

result = diffev2(objective, x0=bounds, penalty=penalty, npop=10, gtol=200, \
                 disp=False, full_output=True, itermon=mon, maxiter=M*N*100)

print result[0]
print result[1]

结果如下所示:

Generation 0 has Chi-Squared: -0.434718
Generation 10 has Chi-Squared: -1.733787
Generation 20 has Chi-Squared: -1.859787
Generation 30 has Chi-Squared: -1.860533
Generation 40 has Chi-Squared: -1.860533
Generation 50 has Chi-Squared: -1.860533
Generation 60 has Chi-Squared: -1.860533
Generation 70 has Chi-Squared: -1.860533
Generation 80 has Chi-Squared: -1.860533
Generation 90 has Chi-Squared: -1.860533
Generation 100 has Chi-Squared: -1.860533
Generation 110 has Chi-Squared: -1.860533
Generation 120 has Chi-Squared: -1.860533
Generation 130 has Chi-Squared: -1.860533
Generation 140 has Chi-Squared: -1.860533
Generation 150 has Chi-Squared: -1.860533
Generation 160 has Chi-Squared: -1.860533
Generation 170 has Chi-Squared: -1.860533
Generation 180 has Chi-Squared: -1.860533
Generation 190 has Chi-Squared: -1.860533
Generation 200 has Chi-Squared: -1.860533
Generation 210 has Chi-Squared: -1.860533
STOP("ChangeOverGeneration with {'tolerance': 0.005, 'generations': 200}")
[-0.17207128  0.73183465 -0.28218955]
-1.86053344078

mystic 非常灵活,可以处理任何类型的约束(例如等式、不等式),包括符号和功能约束。 我在上面将约束指定为“惩罚”,这是传统方式,因为当违反约束时,它们会对目标施加惩罚。 mystic 还提供非线性核变换,通过减少有效解的空间(即通过空间映射或核变换)来限制解空间。

例如,这里的mystic 解决了一个问题,该问题破坏了许多 QP 求解器,因为约束不是约束矩阵的形式。它正在优化压力容器的设计。

"Pressure Vessel Design"

def objective(x):
    x0,x1,x2,x3 = x
    return 0.6224*x0*x2*x3 + 1.7781*x1*x2**2 + 3.1661*x0**2*x3 + 19.84*x0**2*x2

bounds = [(0,1e6)]*4
# with penalty='penalty' applied, solution is:
xs = [0.72759093, 0.35964857, 37.69901188, 240.0]
ys = 5804.3762083

from mystic.symbolic import generate_constraint, generate_solvers, simplify
from mystic.symbolic import generate_penalty, generate_conditions

equations = """
-x0 + 0.0193*x2 <= 0.0
-x1 + 0.00954*x2 <= 0.0
-pi*x2**2*x3 - (4/3.)*pi*x2**3 + 1296000.0 <= 0.0
x3 - 240.0 <= 0.0
"""
cf = generate_constraint(generate_solvers(simplify(equations)))
pf = generate_penalty(generate_conditions(equations), k=1e12)


if __name__ == '__main__':

    from mystic.solvers import diffev2
    from mystic.math import almostEqual
    from mystic.monitors import VerboseMonitor
    mon = VerboseMonitor(10)

    result = diffev2(objective, x0=bounds, bounds=bounds, constraints=cf, penalty=pf, \ 
                     npop=40, gtol=50, disp=False, full_output=True, itermon=mon)

    assert almostEqual(result[0], xs, rel=1e-2)
    assert almostEqual(result[1], ys, rel=1e-2)

在这里找到这个,以及大约 100 个类似的例子:https://github.com/uqfoundation/mystic

我是作者,所以有点偏见。但是,偏差非常轻微。 mystic 既成熟又得到很好的支持,在解决硬约束非线性优化问题的能力方面无与伦比。

【讨论】:

  • 更多非线性约束示例:stackoverflow.com/a/42299338/2379433
  • 听起来很有希望,但从pypi.org/project/mystic (cacr.caltech.edu/~mmckerns) 链接的项目主页已损坏!
  • @feetwet:该页面已移至trac.mystic.cacr.caltech.edu/project/mystic/wiki.html。 pypi 上的链接将在即将发布的版本中更新。
  • 我错过了什么还是应该是**而不是*在这条线上return -1*percentile(p, [np.dot(np.atleast_2d(u[i]), x)[0] for i in range(0,M-1)])
  • @ozataman:基本上我对k的选择一般都是从一个小值开始,如果看起来惩罚不被“尊重”,那么我增加k。是的,这一切都是为了使规模与目标相匹配。理想情况下,您希望惩罚与目标相比可以忽略不计,除了违反约束的地方......您希望惩罚主导目标。这是我能做的最好的了。
【解决方案2】:

scipy 有一个用于约束非线性优化的壮观包。

您可以从阅读 optimize doc 开始,但这里有一个 SLSQP 示例:

minimize(func, [-1.0,1.0], args=(-1.0,), jac=func_deriv, constraints=cons, method='SLSQP', options={'disp': True})

【讨论】:

  • 感谢斯莱特,但计算上述问题的雅可比似乎并不简单,至少对我而言。
  • @akhil the jacobian 是可选输入。阅读实际文档可以比我更能帮助您启动和运行。
  • 是的,确实如此。我能够解决上述问题。谢谢vm!
  • 嗨,SLSQP 可以解决非凸优化问题吗?
  • @Chinni SLSQP 应该,它不假设凸函数。也就是说,非凸优化很难,无法告诉您 SLSQP 的效果如何。
【解决方案3】:

正如其他人所评论的那样,SciPy 最小化包是一个很好的起点。我们还回顾了Python Gekko paper 中的许多其他优化包(参见第 4 节)。我在下面包含了一个示例(Hock Schittkowski #71 基准测试),其中包括 Scipy.optimize.minimize 中的目标函数、等式约束和不等式约束。

import numpy as np
from scipy.optimize import minimize

def objective(x):
    return x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]

def constraint1(x):
    return x[0]*x[1]*x[2]*x[3]-25.0

def constraint2(x):
    sum_eq = 40.0
    for i in range(4):
        sum_eq = sum_eq - x[i]**2
    return sum_eq

# initial guesses
n = 4
x0 = np.zeros(n)
x0[0] = 1.0
x0[1] = 5.0
x0[2] = 5.0
x0[3] = 1.0

# show initial objective
print('Initial SSE Objective: ' + str(objective(x0)))

# optimize
b = (1.0,5.0)
bnds = (b, b, b, b)
con1 = {'type': 'ineq', 'fun': constraint1} 
con2 = {'type': 'eq', 'fun': constraint2}
cons = ([con1,con2])
solution = minimize(objective,x0,method='SLSQP',\
                    bounds=bnds,constraints=cons)
x = solution.x

# show final objective
print('Final SSE Objective: ' + str(objective(x)))

# print solution
print('Solution')
print('x1 = ' + str(x[0]))
print('x2 = ' + str(x[1]))
print('x3 = ' + str(x[2]))
print('x4 = ' + str(x[3]))

Python Gekko 也存在同样的问题:

from gekko import GEKKO
m = GEKKO()
x1,x2,x3,x4 = m.Array(m.Var,4,lb=1,ub=5)
x1.value = 1; x2.value = 5; x3.value = 5; x4.value = 1

m.Equation(x1*x2*x3*x4>=25)
m.Equation(x1**2+x2**2+x3**2+x4**2==40)
m.Minimize(x1*x4*(x1+x2+x3)+x3)

m.solve(disp=False)
print(x1.value,x2.value,x3.value,x4.value)

如果 SLSQP 无法解决您的问题,nonlinear programming solvers for Python 上还有一个更全面的讨论帖。如果您需要有关求解器方法的更多信息,可以使用我在 Engineering Design Optimization 上的课程资料。

【讨论】:

  • 感谢您的回答,在def constraint2(x) 中,您能否将sum_eq 提供给约束函数而不像def constraint2(x, sum_eq) 那样对其进行硬编码?我打算使用像np.linalg.dot(x,p) == 0 这样的点积约束,其中p 是我在优化例程之外计算的静态向量。谢谢。
  • 你可以通过使用像new_constraint = functools.partial(constraint2, p=p)987654335@这样的部分函数来做到这一点
  • 您可以在 Gekko 中使用x=m.Array(m.Var,(4,3))p=m.Array(m.Var,3) 中的数组。然后,您可以将点积包含为m.Equations(np.linalg.dot(x,p) == 0)。这里有例子:github.com/BYU-PRISM/GEKKO/blob/master/examples/test_arrays.pygithub.com/BYU-PRISM/GEKKO/blob/master/examples/test_matrix.py
  • @JohnHedengren,嘿伙计,GEKKO 看起来非常棒,并且有一个非常直观和简单的 API,这是大多数优化包所缺乏的。您能否解释一下 GEKKO 与 SciPy 等其他软件包相比的主要优势是什么?
  • 两者各有利弊。 GEKKO 具有超越非线性规划的其他求解模式:gekko.readthedocs.io/en/latest/overview.html 它使用自动微分以稀疏形式为 IPOPT 和 APOPT 等求解器提供精确的一阶和二阶导数。对于大规模问题,它可以更快。主要缺点是它不能使用外部模型函数调用。方程式必须用 Python 编写。
【解决方案4】:

通常情况下,您可以使用scipy.optimize 函数或lmfit 来进行拟合,它只是扩展了 scipy.optimize 包,以便更容易传递边界之类的东西。就个人而言,我喜欢使用kmpfit,它是 kapteyn 库的一部分,基于 MPFIT 的 C 实现。

scipy.optimize.minimize() 可能是最容易获得且常用的。

【讨论】:

  • 感谢伪立方。这个包看起来很棒。我现在就试试。只是想知道是否有一种简单的方法来进行矩阵乘法?还是我需要扩展每个术语?
  • 其实我的问题不是合适的问题。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2013-02-15
  • 2020-04-15
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多