Using parameters to solve differential equations in GEKKO python

user_17455072 0
Viewed 128

I need to integrate over a system of differential equations in GEKKO and want to use parameters to alter a gradient.

Conceptually this works, but the output is not what I expected. I added a code snippet for a sample problem below to illustrate the problem.

The solver will integrate over a specified time horizon, in this case t = [0, 1, 2, 3]

A parameter is defined that represents the gradient at each value in time called p = [0, 1, 2, 3].

My expectation is that the gradient at t=0 is 0, at t=1 is 1 and so forth. Instead GEKKO interprets it as the gradient at t=0 is 1, at t=1 is 2 etc.

Is there a reason why GEKKO does not use the gradient information at t=0?

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

n = 3
m = GEKKO(remote=False)
m.time = np.linspace(0,n,n+1)       # [0, 1, 2 .. n]
y = m.Var(value=5)
p = m.Param(value=list(range(n+1))) # [0, 1, 2 .. n]

m.Equation(y.dt()==1*p)
m.options.IMODE=4
m.solve(disp=False)

plt.plot(m.time, y.value, '-x')
plt.xlabel('time'); plt.ylabel('y')
plt.show()

enter image description here

gekko python
1 Answers

Set NODES=3 to get the desired output. The default in Gekko is NODES=2 that is fast but not include interior calculation points for each step. Increasing the Nodes has the effect of increasing solution accuracy but also more variables to solve. For large problems with a long time horizon, improve the speed of simulation by using IMODE=7 (sequential solution) instead of IMODE=4 (simultaneous solution).

Nodes

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

n = 3
t = np.linspace(0,n,n+1)
for nodes in [2,3]:

    m = GEKKO(remote=False)
    m.time = t
    y = m.Var(value=5)
    p = m.Param(t)

    m.Equation(y.dt()==1*p)
    m.options.IMODE=4
    m.options.NODES=nodes
    m.solve(disp=False)
    plt.plot(m.time, y.value, '-x')

plt.legend(['NODES=2','NODES=3'])
plt.xlabel('time'); plt.ylabel('y')
plt.grid()
plt.show()

Using t = np.linspace(0,n,301) and p = m.Param(np.floor(t)) with more time points shows that the solutions converge with either NODES=2 or NODES=3.

Solutions Converge

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

n = 3
t = np.linspace(0,n,301)
for nodes in [2,3]:

    m = GEKKO(remote=False)
    m.time = t
    y = m.Var(value=5)
    p = m.Param(np.floor(t))

    m.Equation(y.dt()==1*p)
    m.options.IMODE=4
    m.options.NODES=nodes
    m.solve(disp=False)
    plt.plot(m.time, y.value, '-.')

plt.legend(['NODES=2','NODES=3'])
plt.xlabel('time'); plt.ylabel('y')
plt.grid()
plt.show()

Additional information on NODES is in the Dynamic Optimization course in the section on Orthogonal Collocation on Finite Elements.

user_23669410
Related