High Index DAE Solution

APMonitor and Gekko solve high index Differential and Algebraic Equations (DAEs) without rearrangement or differentiation. A common method for high index DAEs is to use Pantelides algorithm to reduce the index order. The index of a DAE is the number of times the equations must be differentiated to return to index-0 (ODE) form. This method has several challenges:

APMonitor solves high index DAEs without differentiation with orthogonal collocation on finite elements and efficient solvers.

Index-0 to Index-3 Pendulum Example

The pendulum weight is located on the end of a massless cord suspended from a pivot, without friction. When given an initial push, it swings back and forth at a constant amplitude. Real pendulums are subject to friction and air drag, so the amplitude of the swing declines.

The following models are mathematically equivalent but are of different index order. The most natural form is an index-3 differential and algebraic (DAE) equation form, posed in terms of absolute position. The index is the number of times the equations must be differentiated to achieve an ordinary differential equation (ODE) form.


Index-0 DAE (ODE) Model

The differentiation of the equations (3 times) creates numerical drift in the solution, especially noticeable with lambda and y.

# Pendulum - Index 0 DAE = ODE
from gekko import GEKKO
import numpy as np

mass = 1
g = 9.81
s = 1

m = GEKKO()
x = m.Var(0)
y = m.Var(-s)
v = m.Var(1)
w = m.Var(0)
lam = m.Var(mass*(1+s*g)/2*s**2)

m.Equation(lam.dt()==(-4*lam*(x*v+y*w) - 1.5*g*mass*w)/(x**2+y**2))
m.Equation(x.dt()==v)
m.Equation(y.dt()==w)
m.Equation(mass*v.dt()==-2*x*lam)
m.Equation(mass*w.dt()==-mass*g - 2*y*lam)

m.time = np.linspace(0,2*np.pi,100)
m.options.IMODE=4
m.options.NODES=3
m.solve(disp=False)

import matplotlib.pyplot as plt
plt.figure(figsize=(10,5))
plt.subplot(3,1,1)
plt.plot(m.time,x.value,label='x')
plt.plot(m.time,y.value,label='y')
plt.ylabel('Position')
plt.legend(); plt.grid()
plt.subplot(3,1,2)
plt.plot(m.time,v.value,label='v')
plt.plot(m.time,w.value,label='w')
plt.ylabel('Velocity')
plt.legend(); plt.grid()
plt.subplot(3,1,3)
plt.plot(m.time,lam.value,label=r'$\lambda$')
plt.legend(); plt.grid()
plt.xlabel('Time')
plt.savefig('index0.png',dpi=600)
plt.show()

! Pendulum - Index 0 DAE = ODE
Model pend0
  Parameters
    m = 1
    g = 9.81
    s = 1
  End Parameters

  Variables
    x = 0
    y = -s
    v = 1
    w = 0
    lam = m*(1+s*g)/2*s^2
  End Variables

  Equations
    $lam = (-4*lam*(x*v+y*w) - 1.5*g*m*w)/(x^2+y^2)
    $x = v
    $y = w
    m*$v = -2*x*lam
    m*$w = -m*g - 2*y*lam
  End Equations
End Model

Index-1 DAE Model

There is less numerical drift with the index-1 form but the solution is still not sufficiently accurate.

# Pendulum - Index 1 DAE
from gekko import GEKKO
import numpy as np

mass = 1
g = 9.81
s = 1

m = GEKKO()
x = m.Var(0)
y = m.Var(-s)
v = m.Var(1)
w = m.Var(0)
lam = m.Var(mass*(1+s*g)/2*s**2)

m.Equation(mass*(v**2+w**2-g*y)-2*lam*(x**2+y**2)==0)
m.Equation(x.dt()==v)
m.Equation(y.dt()==w)
m.Equation(mass*v.dt()==-2*x*lam)
m.Equation(mass*w.dt()==-mass*g - 2*y*lam)

m.time = np.linspace(0,2*np.pi,100)
m.options.IMODE=4
m.options.NODES=3
m.solve(disp=False)

import matplotlib.pyplot as plt
plt.figure(figsize=(10,5))
plt.subplot(3,1,1)
plt.plot(m.time,x.value,label='x')
plt.plot(m.time,y.value,label='y')
plt.ylabel('Position')
plt.legend(); plt.grid()
plt.subplot(3,1,2)
plt.plot(m.time,v.value,label='v')
plt.plot(m.time,w.value,label='w')
plt.ylabel('Velocity')
plt.legend(); plt.grid()
plt.subplot(3,1,3)
plt.plot(m.time,lam.value,label=r'$\lambda$')
plt.legend(); plt.grid()
plt.xlabel('Time')
plt.savefig('index1.png',dpi=600)
plt.show()

! Pendulum - Index 1 DAE
Model pend1
  Parameters
    m = 1
    g = 9.81
    s = 1
  End Parameters

  Variables
    x = 0
    y = -s
    v = 1
    w = 0
    lam = m*(1+s*g)/2*s^2
  End Variables

  Equations
    m*(v^2+w^2-g*y) - 2*lam*(x^2+y^2) = 0
    $x = v
    $y = w
    m*$v = -2*x*lam
    m*$w = -m*g - 2*y*lam
  End Equations
End Model

Index-2 DAE Model

Some DAE numerical methods can solve up to index-2 DAEs in Hessenberg form. This is not required with Gekko and APMonitor that can solve any index-2 or higher DAE form. The index-2 form has no numerical drift but is still not as desirable as index-3 form that requires no rearrangement.

# Pendulum - Index 1 DAE
from gekko import GEKKO
import numpy as np

mass = 1
g = 9.81
s = 1

m = GEKKO()
x = m.Var(0)
y = m.Var(-s)
v = m.Var(1)
w = m.Var(0)
lam = m.Var(mass*(1+s*g)/2*s**2)

m.Equation(x*v + y*w==0)
m.Equation(x.dt()==v)
m.Equation(y.dt()==w)
m.Equation(mass*v.dt()==-2*x*lam)
m.Equation(mass*w.dt()==-mass*g-2*y*lam)

m.time = np.linspace(0,2*np.pi,100)
m.options.IMODE=4
m.options.NODES=3
m.solve(disp=False)

import matplotlib.pyplot as plt
plt.figure(figsize=(10,5))
plt.subplot(3,1,1)
plt.plot(m.time,x.value,label='x')
plt.plot(m.time,y.value,label='y')
plt.ylabel('Position')
plt.legend(); plt.grid()
plt.subplot(3,1,2)
plt.plot(m.time,v.value,label='v')
plt.plot(m.time,w.value,label='w')
plt.ylabel('Velocity')
plt.legend(); plt.grid()
plt.subplot(3,1,3)
plt.plot(m.time,lam.value,label=r'$\lambda$')
plt.legend(); plt.grid()
plt.xlabel('Time')
plt.savefig('index2.png',dpi=600)
plt.show()

! Pendulum - Index 2 DAE
Model pend2
  Parameters
    m = 1
    g = 9.81
    s = 1
  End Parameters

  Variables
    x = 0
    y = -s
    v = 1
    w = 0
    lam = m*(1+s*g)/2*s^2
  End Variables

  Equations
    x * v + y * w = 0
    $x = v
    $y = w
    m*$v = -2*x*lam
    m*$w = -m*g - 2*y*lam
  End Equations
End Model

Index-3 DAE Model

The index-3 form requires no rearrangement and has no numerical drift associated with index reduction. APMonitor and Gekko solve high-index DAEs with orthogonal collocation on finite elements either in simultaneous or sequential modes. This also allows the DAEs to be solved in optimization problems versus only simulation.

# Pendulum - Index 3 DAE
from gekko import GEKKO
import numpy as np

mass = 1
g = 9.81
s = 1

m = GEKKO()
x = m.Var(0)
y = m.Var(-s)
v = m.Var(1)
w = m.Var(0)
lam = m.Var(mass*(1+s*g)/2*s**2)

m.Equation(x**2+y**2==s**2)
m.Equation(x.dt()==v)
m.Equation(y.dt()==w)
m.Equation(mass*v.dt()==-2*x*lam)
m.Equation(mass*w.dt()==-mass*g-2*y*lam)

m.time = np.linspace(0,2*np.pi,100)
m.options.IMODE=4
m.options.NODES=3
m.solve(disp=False)

import matplotlib.pyplot as plt
plt.figure(figsize=(10,5))
plt.subplot(3,1,1)
plt.plot(m.time,x.value,label='x')
plt.plot(m.time,y.value,label='y')
plt.ylabel('Position')
plt.legend(); plt.grid()
plt.subplot(3,1,2)
plt.plot(m.time,v.value,label='v')
plt.plot(m.time,w.value,label='w')
plt.ylabel('Velocity')
plt.legend(); plt.grid()
plt.subplot(3,1,3)
plt.plot(m.time,lam.value,label=r'$\lambda$')
plt.legend(); plt.grid()
plt.xlabel('Time')
plt.savefig('index3.png',dpi=600)
plt.show()

! Pendulum - Index 3 DAE
Model pend3
  Parameters
    m = 1
    g = 9.81
    s = 1
  End Parameters

  Variables
    x = 0
    y = -s
    v = 1
    w = 0
    lam = m*(1+s*g)/2*s^2
  End Variables

  Equations
    x^2 + y^2 = s^2
    $x = v
    $y = w
    m*$v = -2*x*lam
    m*$w = -m*g - 2*y*lam
  End Equations
End Model

Predictions

The mass and length of a pendulum can be determined by tracking the horizontal position of the pendulum (x). The following is a MATLAB script (pendulum.m) that runs the Index-3 DAE through a series of simulations. As additional data is collected, the model predictions are adjusted to match the observed measurements. The starting values for mass are 1 kg and a length of 1 meter.

The technique for aligning measured and model values is termed Moving Horizon Estimation. This is a technique for parameter estimation with differential and algebraic equation models. In this case there is no steady state data available. The mass and length can be determined by observing the time series of horizontal positions.


! APMonitor Modeling Language
! https://www.apmonitor.com

! Pendulum - Index 3 DAE
Model pend3
  Parameters
    m = 1
    g = 9.81
    s = 1
  End Parameters

  Variables
    x = 0
    y = -s
    v = 1
    w = 0
    lam = m*(1+s*g)/2*s^2
  End Variables

  Equations
    x^2 + y^2 = s^2
    $x = v
    $y = w
    m*$v = -2*x*lam
    m*$w = -m*g - 2*y*lam
  End Equations
End Model

See Inverted Pendulum Control and Crane Pendulum Modeling and Control

Home | High Index DAE Solution