Вопрос об условном выражении ('m.if3') в GEKKO

Вам понадобится всего одно исправление (удалите <0), чтобы использовать функцию if3 в Gekko.

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 использует двоичные переменные, и вам нужно будет использовать решатель смешанных целых чисел, чтобы найти целочисленное решение. Вы можете либо удалить строку, которая переключает решатель на 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()