Для использования функции if3 в Gekko требуется только одно исправление (удалите <0).
R1_1 = m.if3(R1, 0, R1)
Функция m.if3 использует условие, которое переключает то, что используется (аргумент 2или 3) в зависимости от того, является ли условие меньше нуля или больше или равно нулю.Вот результат помощи (m.if3) с некоторыми дополнительными пояснениями:
if3(condition, x1, x2) method of gekko.gekko.GEKKO instance
IF conditional with a binary switch variable.
The traditional method for IF statements is not continuously
differentiable and can cause a gradient-based optimizer to fail
to converge.
Usage: y = m.if3(condition,x1,x2)
Inputs:
condition: GEKKO variable, parameter, or expression
x1 and x2: GEKKO variable, parameter, or expression
Output: GEKKO variable y = x1 when condition<0
y = x2 when condition>=0
Еще одна вещь, которую следует учитывать, - это то, что функция if3 использует двоичные переменные, и вам нужно будет использовать Mixed Integerрешатель, чтобы найти целочисленное решение.Вы можете удалить строку, которая переключает решатель на IPOPT (if3 переключается на решатель APOPT по умолчанию), или вы можете переключить опцию на APOPT вручную.
m.options.SOLVER = 1
Некоторые из констант отсутствовали втвой сценарийЯ добавил несколько фиктивных переменных, чтобы запустить его.
from gekko import GEKKO
import numpy as np
m = GEKKO()
help(m.if3)
tstep = 1.0
cs0=[1,1,1,1]
r0 = 1.0
cFeMax = 1.0
kd = [1,1,1]
Keq = [1,1,1]
Dif = [1,1,1,1]
cg0 = [1,1,1]
nus = np.array([[1,1,1],[1,1,1],[1,1,1],[1,1,1]])
m.time = np.array([0,tstep])
cH = m.CV(value=cs0[0])
cM = m.CV(value=cs0[1])
cW = m.CV(value=cs0[2])
cF = m.CV(value=cs0[3])
R1_1 = m.Var()
r3 = m.Intermediate(r0*(1-cF/cFeMax)**(1/3))
r2 = m.Intermediate(r0*((2*cH + 3*cM)/cFeMax)**(1/3))
r1 = m.Intermediate(r0*(2*cH/cFeMax)**(1/3))
x = m.Intermediate(r1/r0)
y = m.Intermediate(r2/r0)
z = m.Intermediate(r3/r0)
A1 = m.Intermediate(1/x**2/(kd[0]*(1+1/Keq[0])))
A2 = m.Intermediate(1/y**2/(kd[1]*(1+1/Keq[1])))
A3 = m.Intermediate(1/z**2/(kd[2]*(1+1/Keq[2])))
B1 = m.Intermediate((y-x)/x/y*r0/Dif[1])
B2 = m.Intermediate((z-y)/y/z*r0/Dif[2])
B3 = m.Intermediate((1-z)/z*r0/Dif[3])
F = 0
W = m.Intermediate((A1+B1)*(A3*(A2+B2+B3+F)+(A2+B2)*(B3+F))+A2*(A3*(B2+B3+F))+B2*(B3+F))
ceq1 = m.Intermediate((cg0[0]+cg0[1])/(1+Keq[0]))
ceq2 = m.Intermediate((cg0[0]+cg0[1])/(1+Keq[1]))
ceq3 = m.Intermediate((cg0[0]+cg0[1])/(1+Keq[2]))
R1 = m.Intermediate(3/r0/W*((A3*(A2+B2+B3+F)+(A2+B2)*(B3+F))*(cg0[0]-ceq1) \
-(A3*(B2+B3+F)+B2*(B3+F))*(cg0[0]-ceq2) \
-A2*(B3+F)*(cg0[0]-ceq3)))
R2 = m.Intermediate(3/r0/W*(-(B2*(A3+B3+F)+A3*(B3+F))*(cg0[0]-ceq1) \
+((A1+B1+B2)*(A3+B3+F)+A3*(B3+F))*(cg0[0]-ceq2) \
-(A1+B1)*(B3+F)*(cg0[0]-ceq3)))
R3 = m.Intermediate(3/r0/W*(-A2*(B3+F)*(cg0[0]-ceq1) \
-(A1+B1)*(B3+F)*(cg0[0]-ceq2) \
+((A1+B1)*(A2+B2+B3+F)+A2*(B2+B1+F))*(cg0[0]-ceq3)))
R1_1 = m.if3(R1, 0, R1)
m.Equation(cH.dt() == nus[0].dot([R1_1, R2, R3]))
m.Equation(cM.dt() == nus[1].dot([R1_1, R2, R3]))
m.Equation(cW.dt() == nus[2].dot([R1_1, R2, R3]))
m.Equation(cF.dt() == nus[3].dot([R1_1, R2, R3]))
m.options.IMODE = 4
m.options.SOLVER = 1
m.options.nodes = 2
m.solve()