BYU-PRISM / GEKKO

GEKKO Python for Machine Learning and Dynamic Optimization
https://machinelearning.byu.edu
Other
573 stars 102 forks source link

Can GEKKO Terminate on Final Value for its Variables #96

Closed QiyaoWei closed 3 years ago

QiyaoWei commented 3 years ago

Hi! For all GEKKO examples I looked at, the system runs for a fixed number of timesteps, and outputs a solution. Is it possible to make the system terminate on a final value (eg the system terminates when one of the variables reach e-8 instead of reaching 50 timesteps)?

Also, something I noticed is that the number of timesteps largely affect whether GEKKO finds a solution to the system. For example, when linspace is over 50 timesteps the system might have a solution, but there's no solution found for 49 timesteps. Is there a reason for that? Thanks!

APMonitor commented 3 years ago

If you need to integrate until a condition is met, you can use a for loop in Python and check termination conditions. Here is an adaptation of Problem 4 in Gekko simulation example problems. The solution speed could be improved by integrating larger time section and then checking in that section if the condition is met.

ODE Solution with Sequential Time Steps

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

m = GEKKO(remote=False)
m.time = [0,0.1]

# create GEKKO parameter
u = m.Param(value=0.0)

# create GEKKO variables
x = m.Var(0.0) 
y = m.Var(0.0) 

# create GEEKO equations
m.Equation(2*x.dt()==-x+u) 
m.Equation(5*y.dt()==-y+x) 

# simulate ODE
m.options.IMODE = 4

# store time, u, x, y
ts = [0]; us = [0]; xs = [0]; ys = [0]
while xs[-1]<1.2:
    # step change in u at time=1
    if ts[-1]>=1.0:
        u.value = 2.0

    m.solve(disp=False)

    # store next step
    ts.append(ts[-1]+0.1)
    us.append(u.value[0])
    xs.append(x.value[1])
    ys.append(y.value[1])

# plot results
plt.plot(ts,us,'g:',label='u(t)')
plt.plot(ts,xs,'b-',label='x(t)')
plt.plot(ts,ys,'r--',label='y(t)')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

If you just need an ODE simulation tool, you may want to consider a Scipy package such as ODEINT (see same tutorials but with ODEINT) or solve_ivp that also includes the concept of events.

It is easier for solvers to find a solution with smaller time steps, especially with stiff systems. You can create non-uniform time-steps for your simulation such as at the initial condition or around large input changes. Here is an example:

m.time = [0,0.01,0.02,0.05,0.1,0.2,0.5,1.0,2.0,5.0,10.0]

An effective strategy is to use m.options.IMODE=7 (sequential simulation) with fine time steps to calculate a solution. There is also m.options.COLDSTART=2 with m.options.IMODE=4 (simultaneous simulation) that performs a Lower Block Triangular decomposition of the model sparsity structure for a more robust approach to the sequential time simulation. It can solve independent parts of your model, even within the time-step. There is more information on this approach in this article:

Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E., Initialization Strategies for Optimization of Dynamic Systems, Computers and Chemical Engineering, Vol. 78, pp. 39-50, DOI: 10.1016/j.compchemeng.2015.04.016.