python - Solving PENDULUM2 from Schittkowski DAE test suite?

Question

I just tried to solve one of the DAE problems from Schittkowski DAE test suite (http://klaus-schittkowski.de/mc_dae.htm) but no success (Python 3.7).

m=GEKKO(remote=False)
m.time = np.linspace(0.0,100.0,100)
m1 = m.Param(value=1.0) 
g = m.Param(value=9.81)
s = m.Param(value=1.0)

m.options.IMODE=4
m.options.NODES=3
#m.options.RTOL=1e-3

#Initial values t=0
x = m.Var(value=0.0)
y = m.Var(value=-s)
v = m.Var(value=1.0)
w = m.Var(value=0.0)
l = m.Var(value=(m1*(1.0+s*g)/2.0/s**2))

m.Equation(x.dt() == v) 
m.Equation(y.dt() == w)
m.Equation(m1*v.dt() == -2*x*l)
m.Equation(m1*w.dt() == -m1*g - 2*y*l)
m.Equation(x*v  == -y*w)           

m.solve(disp=False)
plt.plot(m.time,x)

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-7-a059f1ac5393> in <module>
     23 m.Equation(x*v  == -y*w)
     24 
---> 25 m.solve(disp=False)
     26 plt.plot(m.time,x)

C:WPy64-3771python-3.7.7.amd64libsite-packagesgekkogekko.py in solve(self, disp, debug, GUI, **kwargs)
   2172             #print APM error message and die
   2173             if (debug >= 1) and ('@error' in response):
-> 2174                 raise Exception(response)
   2175 
   2176             #load results

Exception:  @error: Solution Not Found

I just tried to increase number of points in m.time and number of nodes in m.options.NODES=3. Changing m.options.RTOL does not help as well.

Folllowing the suggestions about resolving the problem of not finding the solution, I managed to get one. Here is the code.

m=GEKKO(remote=False)
m.time = np.linspace(0.0,100.0,500)
m1 = m.Param(value=1.0) 
g = m.Param(value=9.81)
s = m.Param(value=1.0)

m.options.IMODE=4
m.options.NODES=7

#Initial values t=0
x = m.Var(value=0.0)
y = m.Var(value=-1.0)
v = m.Var(value=1.0)
w = m.Var(value=0.0)
l = m.Var(value=4.405)

m.Equation(x.dt() == v) 
m.Equation(y.dt() == w)
m.Equation(m1*v.dt() == -2*x*l)
m.Equation(m1*w.dt() == -m1*g - 2*y*l)
m.Equation(x*v  == -y*w)           

# initialize to get a feasible solution
m.options.COLDSTART=2
m.solve(disp=False)

# optimize, preserving the initial conditions (TIME_SHIFT=0)
m.options.TIME_SHIFT=0
m.options.COLDSTART=0
m.solve(disp=False)
plt.plot(m.time,x)

The resulting plot is

Oscillating behavior will decrees with increasing the time steps or collocation NODES on behave of the calculation time. On the other hand in Julia solving the same problem with the following code

m1 = 1.0 
g = 9.81
s =1.0
u? = [0.0,  -s,  1.0,   0.0,      m1*(1.0+s*g)/2.0/s^2]
du? = [1.0,0.0,0.0,0.0,-g + 2.0 *s*(m1*(1.0+s*g)/2.0/s^2) ]
tspan = (0.0,100.0)
function f(out,du,u,p,t)
    x,y,v,w,l = u
  out[1] = v - du[1]
  out[2] = w - du[2]
  out[3] = -2*x*l/m1 - du[3]
  out[4] = -g - 2.0*y*l/m1 - du[4]
  out[5] = x*v + y*w  
end
#using DifferentialEquations
using Sundials
differential_vars = [true,true,true,true,false]
prob = DAEProblem(f,du?,u?,tspan,differential_vars=differential_vars)
sol = solve(prob,IDA())
n=5251
u1=zeros(n)
u2=zeros(n)
u3=zeros(n)
u4=zeros(n)
u5=zeros(n)
for i=1:n u1[i],u2[i],u3[i],u4[i],u5[i] = sol.u[i] end
plot(sol.t,u1)

will give the following plot in rather short period of time. On the other hand, oscilating behaviour is almost unrecognisible. In gekko, I suppose there will requre quite a large number of time steps and quite significant calculation time. I did not tried that.

It seems that there is no general advice of how to solve this kind of problems in Gekko. I would like someone to comment on this.

Best Regards, Radovan

question from:https://stackoverflow.com/questions/66059473/solving-pendulum2-from-schittkowski-dae-test-suite

Answer 1

Gekko reports the time steps that you request and does not fill in additional points. You are observing aliasing that comes with a different sampling frequency of the simulation than the natural frequency of the oscillating system. To resolve this, you can either increase the sampling frequency or else match the natural frequency so that you are always plotting the minima and maxima. Here is a solution with 5000 data points. With IMODE=7 Gekko scales linearly with increased horizon time.

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m=GEKKO(remote=False)
m.time = np.linspace(0.0,100.0,5000)
m1 = m.Param(value=1.0) 
g = m.Param(value=9.81)
s = m.Param(value=1.0)

m.options.IMODE=7
m.options.NODES=3

#Initial values t=0
x = m.Var(value=0.0)
y = m.Var(value=-1.0)
v = m.Var(value=1.0)
w = m.Var(value=0.0)
l = m.Var(value=4.405)

m.Equation(x.dt() == v) 
m.Equation(y.dt() == w)
m.Equation(m1*v.dt() == -2*x*l)
m.Equation(m1*w.dt() == -m1*g - 2*y*l)

# Index-3 DAE
#m.Equation(x**2+y**2==1)

# Index-2 DAE
m.Equation(x*v  == -y*w)           

m.solve(disp=False)
plt.plot(m.time,x)
plt.show()

With IMODE=7 you also don't need the initialization step. If you are solving parameter optimization then you will benefit from the initialization step. One other suggestion is to use the index-3 DAE form instead of the index-2 DAE form:

# Index-3 DAE
m.Equation(x**2+y**2==1)

# Index-2 DAE
# m.Equation(x*v  == -y*w) 

There is more information on the DAE index and advantages of using a higher index DAE solver like GEKKO versus DAE solvers that are limited to index-1 or index-2 Hessenberg form in the article Initialization strategies for optimization of dynamic systems.

There is no difference for this case study but numerical drift can be an issue if using index-0 (ODE) form.