High Index DAE Solution
Apps.PendulumMotion History
Hide minor edits - Show changes to output
Changed line 5 from:
to:
Solve high index Differential and Algebraic Equations (DAEs) without rearrangement or differentiation. A common method for high index DAEs is to use [[https://en.wikipedia.org/wiki/Pantelides_algorithm|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. Most ODE and DAE solvers have several challenges:
Changed lines 13-14 from:
to:
(:html:)
<iframe width="560" height="315" src="https://www.youtube.com/embed/Mj4OSUKHdXY" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
(:htmlend:)
<iframe width="560" height="315" src="https://www.youtube.com/embed/Mj4OSUKHdXY" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
(:htmlend:)
Changed lines 401-403 from:
* Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E., Initialization Strategies for Optimization of Dynamic Systems, Computers and Chemical Engineering, 2015, Vol. 78, pp. 39-50, DOI: 10.1016/j.compchemeng.2015.04.016 [[http://www.sciencedirect.com/science/article/pii/S0098135415001179 | Article]]
to:
* Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E., Initialization Strategies for Optimization of Dynamic Systems, Computers and Chemical Engineering, 2015, Vol. 78, pp. 39-50, DOI: 10.1016/j.compchemeng.2015.04.016 [[http://www.sciencedirect.com/science/article/pii/S0098135415001179 | Article]]
%width=550px%Attach:high_index_DAEs.png
%width=550px%Attach:high_index_DAEs.png
Changed line 228 from:
# Pendulum - Index 1 DAE
to:
# Pendulum - Index 2 DAE
Changed line 399 from:
* Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E., Initialization Strategies for Optimization of Dynamic Systems, Computers and Chemical Engineering, 2015, Vol. 78, pp. 39-50, DOI: 10.1016/j.compchemeng.2015.04.016. [[http://www.sciencedirect.com/science/article/pii/S0098135415001179 | Article]]
to:
* Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E., Initialization Strategies for Optimization of Dynamic Systems, Computers and Chemical Engineering, 2015, Vol. 78, pp. 39-50, DOI: 10.1016/j.compchemeng.2015.04.016 [[http://www.sciencedirect.com/science/article/pii/S0098135415001179 | Article]]
Changed line 27 from:
m \left( v^2 + w^2 - g \, y \right) - 2 \lambda \left( x^2 + y^2 \right) = 0
to:
{$m \left( v^2 + w^2 - g \, y \right) - 2 \lambda \left( x^2 + y^2 \right) = 0$}
Changed lines 17-18 from:
to:
A pendulum is used to investigate the effect of Differential and Algebraic (DAE) index with index-0 to index-3 DAE forms of the same model. 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 pendulum equations of motion calculate ''x''-axis and ''y''-axis positions as ''x'' and ''y'' and velocities as ''v'' and ''w'', respectively. Additional parameters include ''m'' as the mass of pendulum, ''g'' as a gravitational constant, and ''s'' as the length of pendulum. The variable {`\lambda`} is a Lagrange multiplier.
Changed lines 21-22 from:
to:
'''Index-0 DAE'''
{$\frac{d\lambda}{dt} = \frac{-4 \lambda \left( x \, v + y \, w \right) }{x^2+y^2}$}
'''Index-1 DAE'''
m \left( v^2 + w^2 - g \, y \right) - 2 \lambda \left( x^2 + y^2 \right) = 0
'''Index-2 DAE'''
{$x \, v + y \, w = 0$}
'''Index-3 DAE'''
{$x^2 + y^2 = s^2$}
{$\frac{d\lambda}{dt} = \frac{-4 \lambda \left( x \, v + y \, w \right) }{x^2+y^2}$}
'''Index-1 DAE'''
m \left( v^2 + w^2 - g \, y \right) - 2 \lambda \left( x^2 + y^2 \right) = 0
'''Index-2 DAE'''
{$x \, v + y \, w = 0$}
'''Index-3 DAE'''
{$x^2 + y^2 = s^2$}
Changed line 39 from:
%width=550px%Attach:high_index_DAE_pendulum.png
to:
%width=450px%Attach:high_index_DAE_pendulum.png
Changed lines 21-23 from:
to:
A pendulum is used to investigate the effect of Differential and Algebraic (DAE) index with different forms of the same model. Index-0 to index-3 DAE forms are shown below with the simulation source code.
The models are mathematically equivalent but are of different index order. The most natural form is an index-3 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. This differentiation can cause numerical drift as seen in the
The models are mathematically equivalent but are of different index order. The most natural form is an index-3 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. This differentiation can cause numerical drift as seen in the
Changed lines 379-383 from:
See [[https://apmonitor.com/do/index.php/Main/InvertedPendulum|Inverted Pendulum Control]] and [[https://apmonitor.com/do/index.php/Main/ControlTypes|Crane Pendulum Modeling and Control]]
to:
See [[https://apmonitor.com/do/index.php/Main/InvertedPendulum|Inverted Pendulum Control]] and [[https://apmonitor.com/do/index.php/Main/ControlTypes|Crane Pendulum Modeling and Control]]
'''References'''
* Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E., Initialization Strategies for Optimization of Dynamic Systems, Computers and Chemical Engineering, 2015, Vol. 78, pp. 39-50, DOI: 10.1016/j.compchemeng.2015.04.016. [[http://www.sciencedirect.com/science/article/pii/S0098135415001179 | Article]]
'''References'''
* Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E., Initialization Strategies for Optimization of Dynamic Systems, Computers and Chemical Engineering, 2015, Vol. 78, pp. 39-50, DOI: 10.1016/j.compchemeng.2015.04.016. [[http://www.sciencedirect.com/science/article/pii/S0098135415001179 | Article]]
Changed line 19 from:
%width=300px%Attach:pendulum.png
to:
%width=200px%Attach:pendulum.png
Changed lines 19-21 from:
Attach:pendulum.png
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.
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.
to:
%width=300px%Attach:pendulum.png
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. This differentiation can cause numerical drift as seen in the
%width=550px%Attach:high_index_DAE_pendulum.png
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. This differentiation can cause numerical drift as seen in the
%width=550px%Attach:high_index_DAE_pendulum.png
Deleted lines 375-415:
!!!! 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.
* %list list-blogroll% [[Main/MATLAB | MATLAB Pendulum Example]]
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.
----
(:html:)<font size=2><pre>
! 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
</pre></font>
(:htmlend:)
Changed lines 5-9 from:
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 [[https://en.wikipedia.org/wiki/Pantelides_algorithm|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:
*A consistent initial condition is required
*Max index-1 or index-2 (Hessenberg form) is required
* There is numerical drift with repeated differentiation
*
*
* There is numerical
to:
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 [[https://en.wikipedia.org/wiki/Pantelides_algorithm|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. Most ODE and DAE solvers have several challenges:
* Consistent initial condition is required
* Required index-1 or index-2 (Hessenberg form)
* Numerical drift with repeated differentiation
* Consistent initial condition is required
* Required index-1 or index-2 (Hessenberg form)
* Numerical drift with repeated differentiation
Changed line 11 from:
APMonitor solves high index DAEs without differentiation with [[https://apmonitor.com/do/index.php/Main/OrthogonalCollocation|orthogonal collocation on finite elements]] and efficient solvers.
to:
APMonitor and Gekko overcome these challenges by solving high index DAEs without differentiation with [[https://apmonitor.com/do/index.php/Main/OrthogonalCollocation|orthogonal collocation on finite elements]] and efficient sparse solvers.
Changed lines 3-5 from:
(:description Solve high index DAEs with APMonitor and Gekko without algebraic rearrangement.:)
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 [[https://en.wikipedia.org/wiki/Pantelides_algorithm|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. APMonitor solves high index DAEs without differentiation with [[https://apmonitor.com/do/index.php/Main/OrthogonalCollocation|orthogonal collocation on finite elements]] and efficient solvers.
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 [[https://en.wikipedia.org/wiki/Pantelides_algorithm|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:
(:description Solve high index Differential Algebraic Equations (DAEs) in APMonitor and Gekko without consistent initial conditions, without differentiation of high index equations, and solution in the original form that does not have numerical drift.:)
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 [[https://en.wikipedia.org/wiki/Pantelides_algorithm|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:
* A consistent initial condition is required
* Max index-1 or index-2 (Hessenberg form) is required
* There is numerical drift with repeated differentiation
APMonitor solves high index DAEs without differentiation with [[https://apmonitor.com/do/index.php/Main/OrthogonalCollocation|orthogonal collocation on finite elements]] and efficient solvers.
%width=550px%Attach:high_index_DAEs.png
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 [[https://en.wikipedia.org/wiki/Pantelides_algorithm|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:
* A consistent initial condition is required
* Max index-1 or index-2 (Hessenberg form) is required
* There is numerical drift with repeated differentiation
APMonitor solves high index DAEs without differentiation with [[https://apmonitor.com/do/index.php/Main/OrthogonalCollocation|orthogonal collocation on finite elements]] and efficient solvers.
%width=550px%Attach:high_index_DAEs.png
Changed lines 195-220 from:
----
!!! Index-3 DAE 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.
* %list list-blogroll% [[Main/MATLAB | MATLAB Pendulum Example]]
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.
----
(:html:)<font size=2><pre>
! APMonitor Modeling Language
! https://www.apmonitor.com
! Pendulum - Index 3 DAE
to:
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.
%width=550px%Attach:index2_dae.png
(:toggle hide gekko2 button show="Python Gekko Index-2 DAE":)
(:div id=gekko2:)
(:source lang=python:)
# 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()
(:sourceend:)
(:divend:)
(:toggle hide apm2 button show="APMonitor Index-2 DAE":)
(:div id=apm2:)
(:source lang=fortran:)
! 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
(:sourceend:)
(:divend:)
----
!!! 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.
%width=550px%Attach:index3_dae.png
(:toggle hide gekko3 button show="Python Gekko Index-3 DAE":)
(:div id=gekko3:)
(:source lang=python:)
# 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()
(:sourceend:)
(:divend:)
(:toggle hide apm3 button show="APMonitor Index-3 DAE":)
(:div id=apm3:)
(:source lang=fortran:)
! 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
(:sourceend:)
(:divend:)
%width=550px%Attach:index2_dae.png
(:toggle hide gekko2 button show="Python Gekko Index-2 DAE":)
(:div id=gekko2:)
(:source lang=python:)
# 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()
(:sourceend:)
(:divend:)
(:toggle hide apm2 button show="APMonitor Index-2 DAE":)
(:div id=apm2:)
(:source lang=fortran:)
! 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
(:sourceend:)
(:divend:)
----
!!! 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.
%width=550px%Attach:index3_dae.png
(:toggle hide gekko3 button show="Python Gekko Index-3 DAE":)
(:div id=gekko3:)
(:source lang=python:)
# 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()
(:sourceend:)
(:divend:)
(:toggle hide apm3 button show="APMonitor Index-3 DAE":)
(:div id=apm3:)
(:source lang=fortran:)
! 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
(:sourceend:)
(:divend:)
Added line 53:
plt.figure(figsize=(10,5))
Added lines 67-68:
plt.xlabel('Time')
plt.savefig('index0.png',dpi=600)
plt.savefig('index0.png',dpi=600)
Added lines 19-20:
The differentiation of the equations (3 times) creates numerical drift in the solution, especially noticeable with ''lambda'' and ''y''.
Deleted lines 72-74:
! http://www.apmonitor.com
Changed lines 104-111 from:
(:div id=apm1:)
(:source lang
(:divend:)
(:
to:
There is less numerical drift with the index-1 form but the solution is still not sufficiently accurate.
%width=550px%Attach:index1_dae.png
(:toggle hide gekko1 button show="Python Gekko Index-1 DAE":)
%width=550px%Attach:index1_dae.png
(:toggle hide gekko1 button show="Python Gekko Index-1 DAE":)
Changed lines 111-184 from:
to:
# 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()
(:sourceend:)
(:divend:)
(:toggle hide apm1 button show="APMonitor Index-1 DAE":)
(:div id=apm1:)
(:source lang=fortran:)
! 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
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()
(:sourceend:)
(:divend:)
(:toggle hide apm1 button show="APMonitor Index-1 DAE":)
(:div id=apm1:)
(:source lang=fortran:)
! 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
Changed line 70 from:
(:source lang=matlab:)
to:
(:source lang=fortran:)
Changed lines 19-20 from:
to:
%width=550px%Attach:index0_dae.png
(:toggle hide gekko0 button show="Python Gekko Index-0 DAE":)
(:div id=gekko0:)
(:source lang=python:)
# 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.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.show()
(:sourceend:)
(:divend:)
(:toggle hide apm0 button show="APMonitor Index-0 DAE":)
(:div id=apm0:)
(:source lang=matlab:)
! APMonitor Modeling Language
! http://www.apmonitor.com
! 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
(:sourceend:)
(:divend:)
(:toggle hide gekko0 button show="Python Gekko Index-0 DAE":)
(:div id=gekko0:)
(:source lang=python:)
# 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.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.show()
(:sourceend:)
(:divend:)
(:toggle hide apm0 button show="APMonitor Index-0 DAE":)
(:div id=apm0:)
(:source lang=matlab:)
! APMonitor Modeling Language
! http://www.apmonitor.com
! 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
(:sourceend:)
(:divend:)
Changed lines 105-117 from:
to:
(:toggle hide apm1 button show="APMonitor Index-0 DAE":)
(:div id=apm1:)
(:source lang=matlab:)
(:sourceend:)
(:divend:)
(:toggle hide gekko1 button show="Python Gekko Index-0 DAE":)
(:div id=gekko1:)
(:source lang=python:)
(:sourceend:)
(:divend:)
(:div id=apm1:)
(:source lang=matlab:)
(:sourceend:)
(:divend:)
(:toggle hide gekko1 button show="Python Gekko Index-0 DAE":)
(:div id=gekko1:)
(:source lang=python:)
(:sourceend:)
(:divend:)
Added lines 1-8:
(:title High Index DAE Solution:)
(:keywords pendulum, Python, model, predictive control, APMonitor, differential, algebraic, modeling language:)
(:description Solve high index DAEs with APMonitor and Gekko without algebraic rearrangement.:)
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 [[https://en.wikipedia.org/wiki/Pantelides_algorithm|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. APMonitor solves high index DAEs without differentiation with [[https://apmonitor.com/do/index.php/Main/OrthogonalCollocation|orthogonal collocation on finite elements]] and efficient solvers.
!!!! Index-0 to Index-3 Pendulum Example
(:keywords pendulum, Python, model, predictive control, APMonitor, differential, algebraic, modeling language:)
(:description Solve high index DAEs with APMonitor and Gekko without algebraic rearrangement.:)
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 [[https://en.wikipedia.org/wiki/Pantelides_algorithm|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. APMonitor solves high index DAEs without differentiation with [[https://apmonitor.com/do/index.php/Main/OrthogonalCollocation|orthogonal collocation on finite elements]] and efficient solvers.
!!!! Index-0 to Index-3 Pendulum Example
Deleted lines 35-37:
* %list list-page% (:html:)<a href="/online/plot.php?d=pendulum&f=pend3.x.xml">Vertical Motion</a>(:htmlend:)
* %list list-page% (:html:)<a href="/online/plot.php?d=pendulum&f=pend3.y.xml">Horizontal Motion</a>(:htmlend:)
Changed lines 77-79 from:
(:htmlend:)
to:
(:htmlend:)
See [[https://apmonitor.com/do/index.php/Main/InvertedPendulum|Inverted Pendulum Control]] and [[https://apmonitor.com/do/index.php/Main/ControlTypes|Crane Pendulum Modeling and Control]]
See [[https://apmonitor.com/do/index.php/Main/InvertedPendulum|Inverted Pendulum Control]] and [[https://apmonitor.com/do/index.php/Main/ControlTypes|Crane Pendulum Modeling and Control]]