使用Mystic差分进化函数时参数值的界限不正确

Question

我目前正在使用化学反应网络理论(CRNT)来研究胰岛素信号的双稳态。使用CRNT，我需要解决一个带有约束的全局优化问题，以获得鞍结点。我已经研究了多种算法来解决我的问题，然而，由于我的问题是非线性和非凸性的，可能是多模态的，我发现只有少数方法是合适的。我发现差分进化(DE)似乎是最合适的开始。因为我没有优化领域的专业知识，所以我在寻找一个python库，它可以充当我的目标函数和约束的黑匣子。经过快速搜索，我发现Mystic为DE提供了一个看起来相当容易使用的函数。但是，当我实现DE函数时，我获得的参数值超出了我规定的范围。

我已经在一个非常简单的问题上实现了DE函数，并获得了很好的结果。除此之外，我还尝试为npop、gtol和maxiter设置较大的值。Npop5000左右的值提供了接近我想要的范围的值，但仍有一些值不在我给出的范围内(可能非常大的npop值会给出我想要的结果)。似乎没有任何东西可以解决参数值超出我指定的范围的问题。下面是我正在运行的代码。

import mystic
from mystic.penalty import quadratic_equality
from mystic.solvers import diffev2
from mystic.constraints import as_constraint

def objective(x):
    k15,k51,k16,k82,k28,k27,k93,k39,k34,a1,a2,c4 = x
    c1 = -a2*(k34 + k39)/(c4*k34*k93)
    c2 = c4*k34*k93*(a1 - a2)*(k15 + k16)/(a2*k16*k51*(k34 + k39))
    c3 = -a1*(k27 + k28)/(c4*k27*k82)
    c5 = -(a1 - a2)/k16
    c6 = -a1/k27
    c7 = -a2/k34
    return 1e10*((1.0*c1**2*c4*k34*k51*k82*k93 + 1.0*c1**2*k27*k34*k51*k93 + 1.0*c1**2*k28*k34*k51*k93 - 2.0*c1*c2*c4*k16*k51*k82*k93 +
     1.0*c1*c2*c4*k34*k51*k82*k93 - 2.0*c1*c2*k16*k27*k51*k93 - 2.0*c1*c2*k16*k28*k51*k93 + 1.0*c1*c2*k27*k34*k51*k93 + 1.0*c1*c2*k28*k34*k51*k93 -
     2.0*c1*c3*c4*k27*k51*k82*k93 - 1.0*c1*c3*c4*k34*k51*k82*k93 - 1.0*c1*c3*k27*k34*k51*k82 - 1.0*c1*c3*k27*k39*k51*k82 - 1.0*c1*c4**2*k34*k51*k82*k93 +
     1.0*c1*c4*k15*k34*k82*k93 + 1.0*c1*c4*k16*k34*k82*k93 - 1.0*c1*c4*k27*k34*k51*k93 - 1.0*c1*c4*k28*k34*k51*k93 + 1.0*c1*k15*k27*k34*k93 + 1.0*c1*k15*k28*k34*k93 +
     1.0*c1*k16*k27*k34*k93 + 1.0*c1*k16*k28*k34*k93 - 1.0*c2*c3*k16*k34*k51*k82 - 1.0*c2*c3*k16*k39*k51*k82 - 1.0*c2*c3*k27*k34*k51*k82 - 1.0*c2*c3*k27*k39*k51*k82 -
     1.0*c2*c4*k16*k34*k51*k82 - 1.0*c2*c4*k16*k39*k51*k82 - 1.0*c2*k16*k27*k34*k51 - 1.0*c2*k16*k27*k39*k51 - 1.0*c2*k16*k28*k34*k51 - 1.0*c2*k16*k28*k39*k51 -
     2.0*c3*c4*k15*k27*k82*k93 - 1.0*c3*c4*k15*k34*k82*k93 - 2.0*c3*c4*k16*k27*k82*k93 - 1.0*c3*c4*k16*k34*k82*k93 - 1.0*c3*k15*k27*k34*k82 - 1.0*c3*k15*k27*k39*k82 -
     1.0*c3*k16*k27*k34*k82 - 1.0*c3*k16*k27*k39*k82 - 1.0*c4**2*k15*k34*k82*k93 - 1.0*c4**2*k16*k34*k82*k93 - 1.0*c4*k15*k27*k34*k93 - 1.0*c4*k15*k28*k34*k93 -
     1.0*c4*k16*k27*k34*k93 - 1.0*c4*k16*k28*k34*k93)**2/(k16**2*k27**2*k34**2*k51**2*k82**2*k93**2))

#c1
def penalty1(x):
    k15,k51,k16,k82,k28,k27,k93,k39,k34,a1,a2,c4 = x
    return -a2*(k34 + k39)/(c4*k34*k93)
#c2
def penalty2(x):
    k15,k51,k16,k82,k28,k27,k93,k39,k34,a1,a2,c4 = x
    return c4*k34*k93*(a1 - a2)*(k15 + k16)/(a2*k16*k51*(k34 + k39))

#c3
def penalty3(x):
    k15,k51,k16,k82,k28,k27,k93,k39,k34,a1,a2,c4 = x
    return -a1*(k27 + k28)/(c4*k27*k82)

#c5
def penalty4(x):
    k15,k51,k16,k82,k28,k27,k93,k39,k34,a1,a2,c4 = x
    return -(a1 - a2)/k16

#c6
def penalty5(x):
    k15,k51,k16,k82,k28,k27,k93,k39,k34,a1,a2,c4 = x
    return -a1/k27

#c7
def penalty6(x):
    k15,k51,k16,k82,k28,k27,k93,k39,k34,a1,a2,c4 = x
    return -a2/k34

@quadratic_equality(penalty1, k=1e12)
@quadratic_equality(penalty2, k=1e12)
@quadratic_equality(penalty3, k=1e12)
@quadratic_equality(penalty4, k=1e12)
@quadratic_equality(penalty5, k=1e12)
@quadratic_equality(penalty6, k=1e12)
def penalty(x):
    return 0.0

solver = as_constraint(penalty)

b1 = (1e-1,1e2)
b2 = (1e-1,1e2)
b3 = (1e-1,1e2)
b4 = (1e-1,1e2)
b5 = (1e-1,1e2)
b6 = (1e-1,1e2)
b7 = (1e-1,1e2)
b8 = (1e-1,1e2)
b9 = (1e-1,1e2)
b10 = (-1e-1,0)
b11 = (-1e-1,0)
b12 = (1e-1,1e2)


bounds = [b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12]

#npop should be at least 10*dim, where dim is the number of parameters
result = diffev2(objective, x0=bounds, bounds=bounds, penalty=penalty, npop=120, disp=False, full_output=True, gtol=100, maxiter=6000)

print("")
print("The minimized objection function value:")
print(result[1])
print("")
print("Parameter values:")
print(result[0])

运行此命令，我将获得以下输出。

最小化目标函数值1.216082506137729

参数值1.07383892e-01 9.99893116e+01 8.88912946e+01 9.99859090e+01 1.09022526e-01 9.99587677e+01 9.70349805e+01 1.23842240e+01 4.72484236e+00 -1.01491728e-08 -1.01491720e-08 1.00002390e-01

如您所见，为参数值指定的向量提供了-1.01491728e-08和-1.01491720e-08的值，它们应该在(-0.1,0)范围内。

我只是简单地在Mystic中实现或误解了一些错误的东西，或者我需要探索一种不同的算法来解决我的优化问题？如果您建议使用不同的算法会更好，您是否建议我使用分散搜索(SS)？另外，Mystic是否提供了可以执行SS的函数，或者我是否需要自己实现它？

任何帮助都将不胜感激。

Answer 1

我是mystic的作者。首先要注意的是，您使用的是penalty函数。它们不限制优化器的搜索空间，而是在目标上增加了惩罚，“鼓励”优化器在惩罚为非零的情况下不在空间中找到解决方案。所以这是第一件要注意的事情。如果要将解决方案仅限于约束有效的空间区域，请改用constraints。

然而，这并不能说明为什么你的bounds不被尊重。他们应该是的。它们可能不是有原因的，比如您在优化过程中没有找到任何好的解决方案。众所周知，DE对于具有严格边界的大量参数的性能很差。我会附加一个监视器，这样您就可以看到参数是如何随着时间的推移而变化的。在这种情况下，如果您看到所有的参数都从“有效”变为某些超出范围的参数，那么您就发现了一个错误，请在GitHub上报告它。

即使在后一种情况下，您也可以通过添加使用一系列不等式的约束对象来更强烈地强制执行约束，以确保参数在您选择的范围内……但是，我怀疑发生了什么事情，可以通过连接VerboseMonitor来识别。

mystic确实提供了一些类似于分散搜索算法的东西。查看mystic.ensemble解算器。

更新:显示结果我始终得到

>>> mon = mystic.monitors.VerboseMonitor()
>>> result = diffev2(objective, x0=bounds, bounds=bounds, penalty=penalty, npop=40, disp=False, full_output=True, gtol=100, maxiter=6000, itermon=mon)
Generation 0 has ChiSquare: 2070417340.327223
Generation 10 has ChiSquare: 2070417340.327223
Generation 20 has ChiSquare: 1504378818.925922
Generation 30 has ChiSquare: 49226009.999661
Generation 40 has ChiSquare: 49226009.999661
Generation 50 has ChiSquare: 49226009.999661
Generation 60 has ChiSquare: 49226009.999661
Generation 70 has ChiSquare: 26549433.119812
Generation 80 has ChiSquare: 15513978.070527
Generation 90 has ChiSquare: 2637293.216637
Generation 100 has ChiSquare: 2466027.220625
Generation 110 has ChiSquare: 857082.416445
Generation 120 has ChiSquare: 409397.900243
Generation 130 has ChiSquare: 61023.902060
Generation 140 has ChiSquare: 61023.902060
Generation 150 has ChiSquare: 34911.899051
Generation 160 has ChiSquare: 22564.321601
Generation 170 has ChiSquare: 3078.667678
Generation 180 has ChiSquare: 3078.667678
Generation 190 has ChiSquare: 1233.029461
Generation 200 has ChiSquare: 1233.029461
Generation 210 has ChiSquare: 161.823695
Generation 220 has ChiSquare: 43.540529
Generation 230 has ChiSquare: 42.662921
Generation 240 has ChiSquare: 16.988486
Generation 250 has ChiSquare: 16.988486
Generation 260 has ChiSquare: 16.988486
Generation 270 has ChiSquare: 8.237803
Generation 280 has ChiSquare: 8.237803
Generation 290 has ChiSquare: 5.994087
Generation 300 has ChiSquare: 5.597002
Generation 310 has ChiSquare: 5.597002
Generation 320 has ChiSquare: 4.998805
Generation 330 has ChiSquare: 4.722383
Generation 340 has ChiSquare: 4.544368
Generation 350 has ChiSquare: 4.544368
Generation 360 has ChiSquare: 4.332436
Generation 370 has ChiSquare: 3.711041
Generation 380 has ChiSquare: 1.618530
Generation 390 has ChiSquare: 1.342488
Generation 400 has ChiSquare: 1.279087
Generation 410 has ChiSquare: 1.266669
Generation 420 has ChiSquare: 1.233121
Generation 430 has ChiSquare: 1.225398
Generation 440 has ChiSquare: 1.225398
Generation 450 has ChiSquare: 1.224930
Generation 460 has ChiSquare: 1.220611
Generation 470 has ChiSquare: 1.220611
Generation 480 has ChiSquare: 1.219702
Generation 490 has ChiSquare: 1.217958
Generation 500 has ChiSquare: 1.217265
Generation 510 has ChiSquare: 1.217095
Generation 520 has ChiSquare: 1.216879
Generation 530 has ChiSquare: 1.216421
Generation 540 has ChiSquare: 1.216096
Generation 550 has ChiSquare: 1.215315
Generation 560 has ChiSquare: 1.215274
Generation 570 has ChiSquare: 1.215125
STOP("ChangeOverGeneration with {'tolerance': 0.005, 'generations': 100}")
>>> lb,ub = zip(*bounds)
>>> result[0] < ub
array([ True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True])
>>> result[0] > lb
array([ True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True])
>>>