I am currently using Chemical Reaction Network Theory (CRNT) to study bistability in insulin signaling. Using CRNT I come to a point where I need to solve a global optimization problem with constraints to obtain saddle nodes. I have looked at multiple algorithms for the solution of my problem, however, since my problem is nonlinear and non-convex with the possibility of being multi-modal, I found that only a handful of methods are appropriate. I found that Differential Evolution (DE) seems the most appropriate to begin with. Since I don't have expertise in the field of optimization, I was looking for a library in python that could act like a sort of black-box for my objective function and constraints. After a quick search, I found that Mystic provides a function for DE that looks fairly simple to use. However, when I implement the DE function, I obtain values for my parameters that are outside of the bounds that I have prescribed.
I have already implemented the DE function on a very simple problem and I obtain great results. In addition to this, I have also tried large values for npop,gtol, and maxiter. A value of around 5000 for npop provides values that are close to the range I want, but there are some values that are still not in the range I have given (possibly a very large npop value will give me the results I desire). Nothing seems to fix this issue of parameter values being outside the bounds that I specified. Below is the exact code that I am running.
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])
Running this, I obtain the following output.
The minimized objection function value 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.00002390e-01]
as you can see, the vector given for the parameter values provides values of -1.01491728e-08 and -1.01491720e-08, which should be in the range (-0.1,0).
Have I simply implemented or misinterpreted something incorrectly in Mystic, or do I need to explore a different algorithm for the solution of my optimization problem? If you suggest that a different algorithm would be better, do you suggest I use Scatter Search (SS)? Also, does Mystic provide functions that can perform SS, or would I need to implement it myself?
Any help would be appreciated.
