Gekko - if3() 语句，用于在 ODE

如何解决Gekko - if3() 语句，用于在 ODE

我尝试在 Gekko 中求解具有逻辑条件的微分方程。我知道 Gekko 不喜欢这些东西，但我认为为了在两个给定表达式 (R1,R2) 之间切换的简单 if3() 函数会得到一些解决方案。这是一个失败的简单示例代码 - 未找到解决方案。

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

A=1
B=1e-5
C=2
D=0.01 
G=1
y0=120

m = GEKKO()    # create GEKKO model
nt = 101
m.time = np.linspace(0,100,nt) # time points

y = m.Var(y0) # create GEKKO variable

R1 = m.Intermediate(-(y+A-A/(y/C+1)**(B/D))/G*(D*y+D*C))
R2 = m.Intermediate(-0.1*y)

z = m.Var()              # This way - Solution not found
z = m.if3(y-60,R1,R2)    #
m.Equation(y.dt()== z)   # 

#m.Equation(y.dt()== R1)
#m.Equation(y.dt()== R2)

# solve ODE
m.options.IMODE = 4
m.options.NODES = 5
m.solve(disp=False)

# plot results
plt.plot(m.time,y)
plt.xlabel('time')
plt.ylabel('y(t)')

不幸的是，通过阅读有关 Logical conditions in Optimization 的文章，我似乎还没有弄清楚如何克服这些问题。

此致，

拉多万

解决方法

发布的脚本给出了错误：

Exception:  @error: Solution Not Found

要查看有关错误的其他详细信息，请使用 disp=True 打开显示屏。这提供了表明问题不可行的额外细节。

 Number of state variables:    2800
 Number of total equations: -  2000
 Number of slack variables: -  800
 ---------------------------------------
 Degrees of freedom       :    0
 
 ----------------------------------------------
 Dynamic Simulation with APOPT Solver
 ----------------------------------------------
Iter:     1 I: -1 Tm:     13.07 NLPi:    9 Dpth:    0 Lvs:    0 Obj:  0.00E+00 Gap:       NaN
 Warning: no more possible trial points and no integer solution
 Maximum iterations
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  13.116200000000001 sec
 Objective      :  0.
 Unsuccessful with error code  0
 ---------------------------------------------------
 
 Creating file: infeasibilities.txt
 Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
 @error: Solution Not Found

模型的几个问题：

z = m.Var() 不是必需的，因为 z = if3() 函数定义了变量。
IMODE=6 是需要的，因为问题中存在松弛变量，需要优化模式来解决问题。
尝试使用 if2() 而不是 if3()。 MPCC 形式似乎比混合整数形式更有效。
求解器 APOPT 对这个问题表现更好。试试m.options.SOLVER=1。

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

A=1; B=1e-5; C=2; D=0.01; G=1
y0=120

m = GEKKO(remote=False)    # create GEKKO model
nt = 101
m.time = np.linspace(0,50,nt) # time points

y = m.Var(y0) # create GEKKO variable

R1 = m.Intermediate(-(y+A-A/(y/C+1)**(B/D))/G*(D*y+D*C))
R2 = m.Intermediate(-0.1*y)

z = m.if2(y-60,1) 
m.Equation(y.dt()== (1-z)*R1+z*R2)

# solve ODE
m.options.IMODE  = 6
m.options.SOLVER = 1
m.options.NODES  = 5
m.solve(disp=True)

# plot results
plt.subplot(2,1,1)
plt.plot(m.time,y)
plt.plot([0,50],[60,60],'r--')
plt.ylabel('y(t)'); plt.grid()
plt.subplot(2,2)
plt.plot(m.time,z)
plt.ylabel('z(t)'); plt.grid()
plt.xlabel('time')
plt.show()

在模型中使用 if 语句的替代方法是在循环中集成并检查条件 y>60 以切换到另一个模型。