Simulation of Lorenz System crashes

user_4884235 0
Viewed 29

I only just started exploring Gekko and tried simulating the Lorenz ODE system. Unfortunately, I get an error ("no solution found") for a simple std case the runs fine using scipy.

The problem solves fine if I only integrate up to say time=0.5 instead of 1.0

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

m = GEKKO()
m.time = np.arange(0.0, 1.0, 0.01)

sigma = 10.; rho   = 28.0; beta  = 8./3.

x = m.Var(value=10); y = m.Var(value=10); z = m.Var(value=10)
t = m.Param(value=m.time)

m.Equation(x.dt()== sigma*(y - x))
m.Equation(y.dt()== x*(rho -z) - y)
m.Equation(z.dt()== x*y - beta*z)
m.options.IMODE = 4
m.options.nodes = 4
m.solve(disp=False)

plt.plot(x.value, y.value)
plt.show()
ode gekko
1 Answers

Sequential simulation m.options.IMODE=7 solves successfully. The simultaneous simulation is a larger problem (1782 Variables/Equations vs. 18 Variables/Equations). It also solves successfully if you reduce m.options.NODES=3 (1188 Variables) or increase the maximum iterations m.options.MAX_ITER=300. It failed previously because it needed 267 iterations to find the solution and the maximum iteration limit for IPOPT is 200 by default.

Lorenz

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

m = GEKKO()
m.time = np.arange(0.0, 1.0, 0.01)

sigma = 10.; rho   = 28.0; beta  = 8./3.

x = m.Var(value=10); y = m.Var(value=10); z = m.Var(value=10)
t = m.Param(value=m.time)

m.Equation(x.dt()== sigma*(y - x))
m.Equation(y.dt()== x*(rho -z) - y)
m.Equation(z.dt()== x*y - beta*z)
m.options.IMODE = 7
m.options.nodes = 4
m.solve(disp=True)

plt.plot(x.value, y.value)
plt.show()
user_23669410
Related