How to solve integral in Python Gekko?

Question

How do I include the integral of a variable x in my optimization model in Python Gekko? I found m.vsum(x) that is a summation over the time horizon but doesn't consider variable time steps. The derivative of a variable x.dt() is built-in to Gekko but what about the integral?

from gekko import GEKKO
import numpy as np
m = GEKKO()
m.time = np.linspace(0,2,5)
m.options.IMODE=4; m.options.NODES=3
x = m.Var(5)
m.Equation(x.dt()==-x)
#y = m.Var(0)
#m.Equation(y=m.integral(x))
m.options.SOLVER=1
m.solve()
print(x.value)

Result

[5.0, 3.0434782609, 1.8525519849, 1.1276403386, 0.68638977133]

Answer 1

Gekko solves differential and algebraic equations. You can include the integral by defining a new variable y that has a derivative equal to x.

from gekko import GEKKO
import numpy as np
m = GEKKO()
m.time = np.linspace(0,2,5)
m.options.IMODE=4; m.options.NODES=6
x = m.Var(5)
m.Equation(x.dt()==-x)
y = m.Var(0)
m.Equation(y.dt()==x)
m.solve()
print(np.round(y.value,5))

The result is [0., 1.96735, 3.1606, 3.88435, 4.32332] that is the integral of x at the time points [0, 0.5, 1.0, 1.5, 2.0]. The initial condition for y is zero but it could be non-zero such as with a Proportional Integral Derivative (PID) controller where the integral term may accumulate each cycle.

In Gekko v0.2.6 (pip install gekko --upgrade), there will be a new integral function (thanks to your question). Your commented lines are correct but don't forget the double equal sign in your equation expression.

m.Equation(y==m.integral(x))

Here is a comparison of the numerical and exact solutions.

from gekko import GEKKO
import numpy as np
m = GEKKO()
m.time = np.linspace(0,2,5)
m.options.IMODE=4; m.options.NODES=6
x = m.Var(5)
m.Equation(x.dt()==-x)
y = m.Var(0)
m.Equation(y==m.integral(x))
z = m.Var(0)
m.Equation(z.dt()==x)
m.options.SOLVER=1
m.solve()

# results
print('Numerical x')
print(np.round(x.value,6))
print('Exact x')
print(np.round(5*np.exp(-m.time),6))
print('Numerical Integral y')
print(np.round(y.value,6))
print('Numerical Integral z End Point')
print(np.round(z.value[-1],6))
print('Exact Integral')
print(np.round(-5*(np.exp(-2)-np.exp(0)),6))

The results have 1e-6 differences because they are numerical solutions that are calculated only at the requested time points. The solution is more accurate with additional time points but then the solution is slower (especially for large problems).

Numerical x
[5.       3.032654 1.839398 1.115651 0.676677]
Exact x
[5.       3.032653 1.839397 1.115651 0.676676]
Numerical Integral y
[0.       1.967346 3.160602 3.884349 4.323323]
Numerical Integral z End Point
4.323323
Exact Integral
4.323324