0

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]

Eric Hedengren
  • 455
  • 1
  • 5
  • 19
  • 2
    Here is an example Gekko problem with an integral objective function: https://apmonitor.com/do/index.php/Main/IntegralObjective – John Hedengren Nov 13 '19 at 05:08
  • 2
    you can use scipy package to integrate pass the time `from scipy.integrate import odeint t=np.linspace(0,2,5) def model(x,t): dxdt = -x return dxdt y=odeint(model,3,t)' – Sebin Sunny Nov 13 '19 at 05:24
  • 1
    @SebinSunny - thanks for the `ODEINT` solution for `x`. Do you know how to implement the integral of `x` after you have the ODEINT solution? This is a different variable `y`. – Eric Hedengren Nov 13 '19 at 05:31

1 Answers1

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
John Hedengren
  • 12,068
  • 1
  • 21
  • 25