我目前正在使用化学反应网络理论 (CRNT) 来研究胰岛素信号传导的双稳态。使用 CRNT，我需要解决具有约束的全局优化问题以获得鞍节点。我已经研究了多种算法来解决我的问题，但是，由于我的问题是非线性和非凸的，并且可能是多模态的，我发现只有少数几种方法是合适的。我发现差分进化（DE）似乎是最合适的开始。由于我没有优化领域的专业知识，我一直在寻找一个 Python 库，它可以像我的目标函数和约束的黑盒一样。经过快速搜索，我发现 Mystic 提供了一个看起来相当简单易用的 DE 功能。但是，当我实现 DE 功能时，
我已经在一个非常简单的问题上实现了 DE 函数，并获得了很好的结果。除此之外，我还尝试了 npop、gtol 和 maxiter 的较大值。npop 大约 5000 的值提供了接近我想要的范围的值，但是有些值仍然不在我给出的范围内（可能非常大的 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
Parameter values [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.00002390 e-01]
如您所见，为参数值提供的向量提供了 -1.01491728e-08 和 -1.01491720e-08 的值，它们应该在 (-0.1,0) 范围内。
我是否只是在 Mystic 中错误地实现或误解了某些东西，或者我是否需要探索不同的算法来解决我的优化问题？如果您建议使用不同的算法会更好，您是否建议我使用散点搜索（SS）？另外，Mystic 是否提供可以执行 SS 的功能，还是我需要自己实现它？
任何帮助，将不胜感激。