В настоящее время я использую Сетевую Теорию Химических Реакций (CRNT) для изучения бистабильности при передаче сигналов инсулина. Используя CRNT, я пришел к тому, что мне нужно решить глобальную задачу оптимизации с ограничениями для получения седловых узлов. Я рассмотрел несколько алгоритмов для решения моей проблемы, однако, поскольку моя задача нелинейная и невыпуклая с возможностью быть мультимодальной, я обнаружил, что подходит лишь несколько методов. Я обнаружил, что дифференциальная эволюция (DE) кажется наиболее подходящей для начала. Поскольку у меня нет опыта в области оптимизации, я искал библиотеку на python, которая могла бы действовать как своего рода черный ящик для моих целевых функций и ограничений. После быстрого поиска я обнаружил, что Mystic предоставляет функцию для DE, которая выглядит довольно простой в использовании. Однако, когда я реализую функцию DE, я получаю значения для моих параметров, которые выходят за пределы, которые я предписал.
Я уже реализовал функцию DE для очень простой задачи и получаю отличные результаты. В дополнение к этому я также пробовал большие значения для npop, gtol и maxiter. Значение около 5000 для npop предоставляет значения, близкие к желаемому диапазону, но есть некоторые значения, которые все еще не находятся в указанном диапазоне (возможно, очень большое значение 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.0000 01]
Как видно, вектор, заданный для значений параметров, обеспечивает значения -1.01491728e-08 и -1.01491720e-08, которые должны находиться в диапазоне (-0,1,0).
Я просто неправильно реализовал или неверно истолковал что-то в Mystic или мне нужно изучить другой алгоритм для решения моей проблемы оптимизации? Если вы предполагаете, что другой алгоритм будет лучше, вы предлагаете использовать Scatter Search (SS)? Кроме того, Mystic предоставляет функции, которые могут выполнять SS, или мне нужно было бы реализовать это самостоятельно?
Любая помощь будет оценена.