Differential Algebraic Equations (DAEs) - Optimization for Decision Science

Prepared by: Prof. Alexander Dowling, Myia Dickens (mdicken2@nd.edu, 2023), Molly Dougher (mdoughe6@nd.edu, 2023)

import sys

if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    !pip install casadi
    import helper

    helper.easy_install()
else:
    sys.path.insert(0, "../")
    import helper
helper.set_plotting_style()

External Documentation for libraries used¶

Pyomo.dae documentation:

https://pyomo.readthedocs.io/en/latest/modeling_extensions/dae.html

CasADi (need to integrate DAEs):

https://web.casadi.org/get/
For local installation: pip install casadi. Warning: installing CasADi with conda will install an “okay” version of Ipopt. If you really want to install CasADi with conda, you’ll likely need to add import idaes to your notebook to load the “good” version of Ipopt (Linux and Windows users).

Dynamic Optimization Overview¶

Introduction to Dynamic Optimization¶

Dynamic systems can be found in a wide range of engineering fields and subfields. This section is meant to be an introduction to solve differentiable-algebraic equations (DAEs) with respect to optimizing dynamic systems. To allow for uniform understanding, the examples and explanations in this notebook will focus on basic physics examples (i.e. the pendulum example) and explanations.

DAEs are typically specified as initial value problems with initial conditions at zero. These equations are expressed with respect to an independent variable, $t$ , and are typically autonomous, or $t$ does not explicitly appear in the equation.

DAEs are typically of the form:

\begin{align} F\left(x, \frac{d x}{d t}, u(t), p, t\right)=0, \quad h(x(0))=0 \end{align}

(1)

$x(t)$ : state variables, functions of time

$u(t)$ : control variables, functions of time

$p$ : variables, independent of $t$

With respect to optimizing DAEs, the equations are structured to the simpler form:

\frac{d z}{d t}= f(z(t), y(t), u(t), p), \quad z(0)=z_0, \quad g(z(t), y(t), u(t), p)=0

(2)

where the state variables $x(t)$ is partitioned into differentiable variables, $z(t)$ , and algebraic variables $y(t)$ . It is assumed that $y(t)$ can be solved independently from $z(t)$ , $u(t)$ , and $p$ , which allows for the DAE to be treated as an ODE of the form:

\begin{align} \frac{d z}{d t}=f(z(t), y[z(t),u(t),p], u(t), p)= \bar f(z(t), u(t), p), \quad z(0)=z_0 \end{align}

(3)

The section, Handling Path Constraints, will explain DAE indexing and how to reduce the index of a DAE to allow it to be treated as an ODE.

General Breakdown of DAEs¶

With respect to constructing DAEs:

$z(t)$ : Differential equations are typically derived from conservation laws.
$y(t)$ : Algebraic equations are typically derived from constitutive and equilibrium laws.
The decision variables of the system fall under the two following variables:
$u(t)$ : the control variables include manipulated variables that change over time.
$p$ : time-independent variables that correspond to parameters, initial conditions, and other steady state conditions.

DAEs as Boundary Value Problems¶

Because DAEs are typically initial value problems, they can also be boundary value problems (BVP), where the initial condition is replaced by boundary conditions. With respect to solving DAE optimization problems, BVP DAEs are of the form:

\begin{align} \frac{d z}{d t}=\hat f(z(t)), \quad h(z(0),z[t_f;z(0)])=0 \end{align}

(4)

In this form the solutions may be nonunique or not exist over a specified reason; therefore, a key property is finding a locally unique solution, expressed by the theorem:

Theorem 8.2: Consider the BVP DAE with solution $\hat z(t)$ . Also let $\bar f(z(t))$ be Lipschitz continuous for all z(t) with $|| z(t) - \hat z(t) || \leq \epsilon$ for some $\epsilon > 0$ and $t \in [0, t_f]$ . Then the solution $\hat z(t)$ is locally unique if and only if the matrix:
$\begin{align} Q(t) = \frac{\partial h\left(\hat{z}(0), \hat{z}\left(t_f\right)\right)}{\partial z(0)}+Z(t) \frac{\partial h\left(\hat{z}(0), \hat{z}\left(t_f\right)\right)}{\partial z\left(t_f\right)} \end{align}$
(5)
is nonsingular, where the fundamental solution matrix $Z(t) = \frac{dz(t_f)}{dz(0)}$ is evaluated at the solution $\hat z(t)$ .

Theorem 8.2 is important when designing DAE constrained optimization problems.

DAE Optimization¶

Typically dynamic systems needed to be optimized over multiple periods, $l = 1,..., N_T$ , in a time range, $t$ . The model, states, and decisions can change during each period, $t \in (t_{l-1}, t_l]$ , or over multiple periods. The multiperiod, dynamic optimization problem is of the form:

\begin{align} \min & \sum_{l=1}^{N_T} \Phi^l\left(z^l\left(t_l\right), y^l\left(t_l\right), u^l\left(t_l\right), p^l\right) \\ \text { s.t. } & \frac{d z^l}{d t}=f^l\left(z^l(t), y^l(t), u^l(t), p^l\right), \quad z^l\left(t_{l-1}\right)=z_0^l, \\ & g^l\left(z^l(t), y^l(t), u^l(t), p^l\right)=0, \\ & u_L^l \leq u^l(t) \leq u_U^l, \\ & p_L^l \leq p^l \leq p_U^l, \\ & y_L^l \leq y^l(t) \leq y_U^l, \\ & z_L^l \leq z^l(t) \leq z_U^l, \quad t \in\left(t_{l-1}, t_l\right], l=1, \ldots, N_T, \\ & h\left(p, z_0^1, z^1\left(t_1\right), z_0^2, z^2\left(t_2\right), \ldots, z_0^{N_T}, z^{N_T}\left(t_{N_T}\right)\right)=0 \end{align}

(6)

In this form, it is assumed that the state variables are not continuous across periods; therefore, the last line is meant to connect the states of each period. The initial conditions and inequality constraints are represented as simple bounds in this form.

Key applications of DAEs and Dynamic Optimization are as follows:

Chemical Reactor Design
Parameter Estimation of a Dynamic System
Batch Process Estimation
Dynamic Real-Time Optimization

A description and further explanation of these applications can be found in Chapter 8 of the Biegler (2010) textbook.

Handling Path Constraints¶

Information taken from Section 8.4 of Biegler (2010).

The motivation for defining the index of a DAE system starts with considering the general algebraic equality constraint $g_E(z,y,u,p)=0$ . After the algebraic and control variables are established, deriving the Euler-Lagrange equations requires variable and equation nesting. In order to be nested, the algebraic variables must be able to be implicitly eliminated from their paired algebraic equation. If $y(t)$ cannot be implicitly eliminated from this algebraic equality, reformulation is needed. This reformulation begins with establishing an index of a DAE system.

DAE Index and Index Reduction¶

Definition 8.6: Consider the DAE systems of the form:
$\begin{align} F\left(x, \frac{d x}{d t}, u(t), p, t\right)=0, \quad h(x(0))=0 \end{align}$
(7)
OR
$\frac{d z}{d t}= f(z(t), y(t), u(t), p), \quad z(0)=z_0, \quad g(z(t), y(t), u(t), p)=0$
(8)
with decisions $u(t$ ) and $p$ fixed. The index is the integer $s$ that represents the minimum number of differentiations of the DAE system (with respect to time) required to determine an ODE for the variables $z(t)$ and $y(t)$ .

Generally, for semiexplicit systems, the index of a DAE system can be determined by first differentiating the algebraic equations and then substituting resulting differential terms into the respective differential equation. Recall the number of differentiations required to isolate an ODE system corresponds to the index.

As the goal is usually to isolate a DAE of index 1, the following algorithm describes the typical procedure to reduce the index of a DAE system.

Algorithm 8.1: (Reduction of High Index DAEs) Start with a DAE of the form:
$\frac{d z}{d t}= f(z(t), y(t), u(t), p), \quad z(0)=z_0, \quad g(z(t), y(t), u(t), p)=0.$
(9)
Check if the DAE system is index 1. If yes, stop.
Identify a subset of algebraic equations that can be solved for a subset of algebraic variables.
Consider the remaining algebraic equations that contain differential variables $z_j$ . Differentiating these remaining algebraic equations with respect to time leads to terms $\frac{dz_j}{dt}$ in the differentiated equations.
For the differential terms $\frac{dz_j}{dt}$ , substitute the right-hand sides of the corresponding differential equations $f_j(z,y,u,p)$ into the differentiated algebraic equations, and eliminate (some of) these differential equations. This leads to new algebraic equations that replace the same number of existing differential equations.
With this new DAE system, go to step 1.

An example of this algorithm in use can be found within the following pendulum example.

Pendulum model and index reduction¶

Pendulum example:

Python version of example: http://apmonitor.com/wiki/index.php/Apps/PendulumMotion
More details on index reduction for example: https://www.lehigh.edu/~wes1/apci/11may00.pdf

$Simple pendulum: mass m=1 on a rigid link of length L=1, with angle \theta from vertical.$

The following DAE system describes the pendulum system:

\begin{align} & x' = u \\ & y' = v \\ & u' = -Tx \\ & v' = g - Ty \\ & x^2 + y^2 = 1 \\ \end{align}

(10)

To find the index of the DAE system, first take the derivative of the algebraic equation and substitute in for the differential variables:

\begin{align} & x^2 + y^2 = 1 \\ & 2xx' + 2yy' = 0 \\ & 2xu + 2yv = 0 \end{align}

(11)

Next, differentiate the obtained equation using the product rule, make similar substitutions, and rearrange:

\begin{align} & 2xu + 2yv = 0 \\ & 2(xu'+ x'u) + 2(yv' + vy') = 0 \\ & 2(-Tx^2 + u^2) + 2(y(g-Ty) + v^2) = 0 \\ & 2(u^2 + v^2) + 2(-Tx^2 + gy - Ty^2) = 0 \\ & 2(u^2 + v^2) - 2T(x^2 + y^2) + 2gy = 0 \end{align}

(12)

Finally, implicitly differentiate the obtained equation, make similar substitutions, and rearrange to isolate the $T'$ term:

\begin{align} & 2(u^2 + v^2) - 2T(x^2 + y^2) + 2gy = 0 \\ & 4(uu' + vv') - 4T(xx' + yy') - 2T'(x^2 + y^2) + 2gy' = 0 \\ & 4(uu' + vv') - 4T(xx' + yy') + 2gy' = 2T'(x^2 + y^2) \\ & 4(-Tux + v(g-Ty)) - 4T(xu + yv) + 2gv = 2T' \\ & 4(-Tux -Tyv + gv) - 4T(xu + yv) + 2gv = 2T' \end{align}

(13)

Now the equations are a system of ODEs. Three differentiations were completed, indicating that this is an Index 3 DAE.

To reduce this DAE to index 1, follow Algorithm 8.1. The above analysis shows that this is not an index 1 DAE system. Therefore, the following algebraic equation is identified: $x^2 + y^2 = 1$ . Differentiate the algebraic equation to yield:

\begin{align} & 2xx' + 2yy' = 0\\ \end{align}

(14)

which can be solved for $x'$ . Substitute for $x'$ and $y'$ to yield:

\begin{align} & 2xu + 2yv = 0\\ \end{align}

(15)

The new DAE system is now:

\begin{align} & y' = v \\ & u' = -Tx \\ & v' = g - Ty \\ & x^2 + y^2 = 1 \\ & 2xu + 2yv = 0 \\ \end{align}

(16)

where $x$ is now an algebraic variable and the final equation must be solved for $T$ . Differentiating the algebraic (last) equation leads to:

\begin{align} & 2(xu'+ x'u) + 2(yv' + vy') = 0 \\ & 2(xu' + u^2) + 2(yv' + v^2) = 0 \\ & 2(u^2 - Tx^2) + 2(v^2 + gy - Ty^2) = 0 \\ \end{align}

(17)

The equations can replace $u' = -Tx$ , leading to the following DAE system:

\begin{align} & y' = v \\ & v' = g - Ty \\ & x^2 + y^2 = 1 \\ & 2xu + 2yv = 0 \\ & (u^2 + v^2) - T\underbrace{(x^2 + y^2)}_{1} + gy = 0 \\ \end{align}

(18)

which is an Index 1 DAE.

## Load libraries
import pyomo.environ as pyo
import pyomo.dae as dae
from pyomo.dae.simulator import Simulator
import matplotlib.pyplot as plt
import numpy as np


## Define function for plotting results
def plot_results(sim, tsim, profiles):
    """
    inputs:

    sim: pyomo.DAE simulator results from a simulation of ODE or DAE
    tsim: 1D array of time samples of the DAE/ODE simulatiion
    profiles: 2D array of simulated differential and algebraic equations

    outputs:

    Plot 1: A plot of the curated profiles
    Plot 2: A plot of the simulated results

    time = list(m.t)
    x = [value(m.x[t]) for t in m.t]
    y = [value(m.y[t]) for t in m.t]

    plt.plot(time, x, '-b', label='x')
    plt.plot(time, y, '-r', label='y')
    plt.xlabel('Time')
    plt.ylabel('Position')
    plt.legend(loc='best')
    plt.show()
    """

    plt.figure(1, figsize=(4, 4))

    # Every series also carries a linestyle, so the figure survives a
    # black-and-white printout (figures/README.md).
    linestyles = ["-", "--", "-.", ":", (0, (3, 1, 1, 1, 1, 1))]

    varorder = sim.get_variable_order()
    algorder = sim.get_variable_order(vartype="algebraic")

    # Create empty dictionary
    results = {}

    # Collect Different Profiles
    for idx1, v in enumerate(varorder):
        i = idx1
        v_ = str(v)
        results[v_] = profiles[:, i]
        plt.plot(tsim, results[v_], label=v, linestyle=linestyles[i % len(linestyles)])

    # Collect Algebraic Profiles
    for idx2, v in enumerate(algorder):
        i = len(varorder) + idx2
        v_ = str(v)
        results[v_] = profiles[:, i]
        plt.plot(tsim, results[v_], label=v, linestyle=linestyles[i % len(linestyles)])

    # Plot the simulated algebraic and differential profiles as they change in time
    plt.xlabel("t", fontsize=16, fontweight="bold")
    plt.tick_params(direction="in", labelsize=15)
    plt.legend(loc="best")

    plt.show()

    # Plot the results of the simulation
    plt.figure(2, figsize=(4, 4))
    x_ = results["x[{t}]"]
    y_ = results["y[{t}]"]
    plt.plot(tsim, np.sqrt(x_**2 + y_**2), "-b", label="length")
    plt.xlabel("t", fontsize=16, fontweight="bold")
    plt.tick_params(direction="in", labelsize=15)
    plt.ylabel(r"$\sqrt{x^2 + y^2}$")
    plt.show()

    # return results

Formulation 1: Index-3 DAE¶

Consider the following model:

\begin{align} \frac{d x}{dt} &= u \\ \frac{d y}{dt} &= v \\ \frac{d u}{dt} &= -T x \\ \frac{d v}{dt} &= g - Ty \\ & x^2 + y^2 = 1 \end{align}

(19)

This assumes mass and length of unity.

Why this notebook does not carry Pyomo units.

The other notebooks in this section declare units= on every Var and Param and check with assert_units_consistent -- see Units in Pyomo.dae Models. This one deliberately does not, and the reason is worth stating rather than hiding.

The pendulum here is normalized: mass and length are both set to unity, so the constraint is written $x^2 + y^2 = 1$ against a dimensionless literal, and the “tension” $T$ is not a force. Reading $\frac{du}{dt} = -Tx$ dimensionally forces

[T] = \frac{\mathrm{m}/\mathrm{s}^2}{\mathrm{m}} = \mathrm{s}^{-2},

(20)

i.e. $T$ is really a tension divided by $mL$ . Restoring dimensions would mean rewriting $x^2 + y^2 = 1$ as $x^2 + y^2 = L^2$ and carrying $mL$ through the index reduction -- changing the very equations whose algebraic structure is the lesson of this notebook.

So: no units here, on purpose. A forced, arbitrary set of units is worse than none, because it produces a check that looks meaningful and is not.

That said, dimensional analysis by hand still pays -- see the next section.

def create_model_index3():
    """
    Creates a model to define the Index 3 DAE system.

    Output: the model, m
    """

    m = pyo.ConcreteModel()

    # Declare time
    m.t = dae.ContinuousSet(bounds=(0.0, 1))

    # Declare parameter - acceleration due to gravity
    m.g = pyo.Param(initialize=9.81)  # [m/s^2]

    # Declare variables indexed over time
    m.x = pyo.Var(m.t)  # horizontal position
    m.y = pyo.Var(m.t)  # vertical position
    m.u = pyo.Var(m.t)  # horizontal velocity
    m.v = pyo.Var(m.t)  # vertical velocity
    m.T = pyo.Var(m.t)  # tension

    # Declare derivative variables
    m.dx = dae.DerivativeVar(m.x)  # with respect to t is implied
    m.dy = dae.DerivativeVar(m.y)
    m.du = dae.DerivativeVar(m.u)
    m.dv = dae.DerivativeVar(m.v)

    # Declare differential equations
    def _dx_eqn(m, t):
        return m.dx[t] == m.u[t]

    m.dx_eqn = pyo.Constraint(m.t, rule=_dx_eqn)

    def _dy_eqn(m, t):
        return m.dy[t] == m.v[t]

    m.dy_eqn = pyo.Constraint(m.t, rule=_dy_eqn)

    def _du_eqn(m, t):
        return m.du[t] == -m.T[t] * m.x[t]

    m.du_eqn = pyo.Constraint(m.t, rule=_du_eqn)

    def _dv_eqn(m, t):
        return m.dv[t] == m.g - m.T[t] * m.y[t]

    m.dv_eqn = pyo.Constraint(m.t, rule=_dv_eqn)

    # Declare algebraic equation
    def _alg_eqn(m, t):
        return m.x[t] ** 2 + m.y[t] ** 2 == 1

    m.alg_eqn = pyo.Constraint(m.t, rule=_alg_eqn)

    # Specify initial conditions
    m.x[0] = 0
    m.y[0] = 1
    m.u[0] = 1
    m.v[0] = 0
    m.T[0] = 1 + m.g

    return m

index3 = create_model_index3()

# Specify integrator options
int_ops = {"print_stats": True, "abstol": 1e-8, "reltol": 1e-6}

# Solve DAEs
sim = Simulator(index3, package="casadi")
try:
    tsim, profiles = sim.simulate(
        numpoints=100, integrator="idas", integrator_options=int_ops
    )
except RuntimeError as err:
    # Expected: IDAS supports index-0 and index-1 DAEs, not this index-3 model.
    print(f"Expected IDAS failure: {err}")
else:
    plot_results(sim, tsim, profiles)

CasADi - 2026-08-21 12:57:46 WARNING("The options 't0', 'tf', 'grid' and 'output_t0' have been deprecated.
The same functionality is provided by providing additional input arguments to the 'integrator' function, in particular:
 * Call integrator(..., t0, tf, options) for a single output time, or
 * Call integrator(..., t0, grid, options) for multiple grid points.
The legacy 'output_t0' option can be emulated by including or excluding 't0' in 'grid'.
Backwards compatibility is provided in this release only.") [.../casadi/core/integrator.cpp:698]
At t = 0 and h = 1.06624e-14, the corrector convergence failed repeatedly or with |h| = hmin.

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[3], line 73
     69 int_ops = {"print_stats": True, "abstol": 1e-8, "reltol": 1e-6}
     70 
     71 # Solve DAEs
     72 sim = Simulator(index3, package="casadi")
---> 73 tsim, profiles = sim.simulate(
     74     numpoints=100, integrator="idas", integrator_options=int_ops
     75 )
     76 

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/pyomo/dae/simulator.py:935, in Simulator.simulate(self, numpoints, tstep, integrator, varying_inputs, initcon, integrator_options)
    931         tsim, profile = self._simulate_with_casadi_with_inputs(
    932             initcon, tsim, varying_inputs, integrator, integrator_options
    933         )
    934     else:
--> 935         tsim, profile = self._simulate_with_casadi_no_inputs(
    936             initcon, tsim, integrator, integrator_options
    937         )
    939 self._tsim = tsim
    940 self._simsolution = profile

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/pyomo/dae/simulator.py:997, in Simulator._simulate_with_casadi_no_inputs(self, initcon, tsim, integrator, integrator_options)
    995 integrator_options['output_t0'] = True
    996 F = casadi.integrator('F', integrator, dae, integrator_options)
--> 997 sol = F(x0=initcon)
    998 profile = sol['xf'].full().T
   1000 if len(self._algvars) != 0:

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/casadi/casadi.py:24049, in Function.__call__(self, *args, **kwargs)
  24045         else:
  24046           return tuple(ret)
  24047       else:
  24048     # Named inputs -> return dictionary
> 24049         return self.call(kwargs)

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/casadi/casadi.py:20698, in Function.call(self, *args)
  20694 
  20695 
  20696 
  20697         """
> 20698         return _casadi.Function_call(self, *args)

RuntimeError: Error in Function::call for 'F' [IdasInterface] at .../casadi/core/function.cpp:1466:
Error in Function::call for 'F' [IdasInterface] at .../casadi/core/function.cpp:362:
.../casadi/interfaces/sundials/idas_interface.cpp:599: IDASolve returned "IDA_CONV_FAIL". Consult IDAS documentation.

Warning: If you run this notebook in Colab, you may get the following runtime error and your kernel may crash:

Why did the IDAS integrator in SUNDIALS fail? It is only meant for index 0 or 1 DAEs! Integrating high index DAEs is really hard!

Formulation 2: Pure ODE Model/ Index-0 Example¶

\begin{align} \frac{d x}{dt} &= u \\ \frac{d y}{dt} &= v \\ \frac{d u}{dt} &= -T x \\ \frac{d v}{dt} &= g - Ty \\ \frac{d T}{dt} &= \frac{-4T(xu+yv)+3gv}{x^2+y^2} \end{align}

(21)

def create_model_ode():
    """
    Creates a model to define the Index 0 DAE system.

    Output: the model, m
    """

    m = pyo.ConcreteModel()

    # Declare time
    m.t = dae.ContinuousSet(bounds=(0.0, 5.0))

    # Declare parameter - acceleration due to gravit
    m.g = pyo.Param(initialize=9.81)  # [m/s^2]

    # Declare variables indexed over time
    m.x = pyo.Var(m.t)  # horizontal position
    m.y = pyo.Var(m.t)  # vertical position
    m.u = pyo.Var(m.t)  # horizontal velocity
    m.v = pyo.Var(m.t)  # vertical velocity
    m.T = pyo.Var(m.t)  # tension

    # Declare derivative variables
    m.dx = dae.DerivativeVar(m.x)  # with respect to t is implied
    m.dy = dae.DerivativeVar(m.y)
    m.du = dae.DerivativeVar(m.u)
    m.dv = dae.DerivativeVar(m.v)
    m.dT = dae.DerivativeVar(m.T)

    # Declare differential equations
    def _dx_eqn(m, t):
        return m.dx[t] == m.u[t]

    m.dx_eqn = pyo.Constraint(m.t, rule=_dx_eqn)

    def _dy_eqn(m, t):
        return m.dy[t] == m.v[t]

    m.dy_eqn = pyo.Constraint(m.t, rule=_dy_eqn)

    def _du_eqn(m, t):
        return m.du[t] == -m.T[t] * m.x[t]

    m.du_eqn = pyo.Constraint(m.t, rule=_du_eqn)

    def _dv_eqn(m, t):
        return m.dv[t] == m.g - m.T[t] * m.y[t]

    m.dv_eqn = pyo.Constraint(m.t, rule=_dv_eqn)

    def _dT_eqn(m, t):
        return (
            m.dT[t]
            == (
                -4 * m.T[t] * (m.x[t] * m.u[t] + m.y[t] * m.v[t])
                + 3 * m.g * m.v[t]
            ) / (m.x[t] ** 2 + m.y[t] ** 2)
        )

    m.dT_eqn = pyo.Constraint(m.t, rule=_dT_eqn)

    # Specify initial conditions
    m.x[0] = 0
    m.y[0] = 1
    m.u[0] = 1
    m.v[0] = 0
    m.T[0] = 1 + m.g

    return m


ode = create_model_ode()

# Specify integrator options
int_ops = {"print_stats": True, "abstol": 1e-6, "reltol": 1e-4}

# Solve DAEs
sim = Simulator(ode, package="casadi")
tsim, profiles = sim.simulate(
    numpoints=100, integrator="idas", integrator_options=int_ops
)

# Plot solution
results = plot_results(sim, tsim, profiles)

FORWARD INTEGRATION:
Number of steps taken by SUNDIALS: 172
Number of calls to the user's f function: 236
Number of calls made to the linear solver setup function: 25
Number of error test failures: 6
Method order used on the last internal step: 5
Method order to be used on the next internal step: 5
Actual value of initial step size: 7.90569e-07
Step size taken on the last internal step: 0.0158106
Step size to be attempted on the next internal step: 0.0316212
Current internal time reached: 5
Number of nonlinear iterations performed: 234
Number of nonlinear convergence failures: 0

The drift, over a longer horizon¶

The plot above is the same quantity as the figure below, run for $t \in [0, 5]$ with IDAS. The figure below is the canonical version used in the lecture handout: the same index-reduced ODE, integrated with scipy.integrate.solve_ivp (RK45) out to $t = 50$ , at two tolerances, and with a third curve that adds the stabilization term $-\gamma F(x) h(x)$ of Ascher & Petzold (1998), eq. (9.40), p. 251.

Drift off the constraint manifold for the index-reduced pendulum

The takeaway: drift is not a tolerance bug. Tightening rtol buys decades, not a fix, because nothing in the reduced ODE refers to $x^2 + y^2 = 1$ — only to its derivatives.

Source: figures/plots/pendulum-drift.py. The figure and Pyomo model both use $\dot T = [-4T(xu+yv) + 3gv]/(x^2+y^2)$ , obtained by differentiating the acceleration constraint.

Dimensional check¶

Restore units: $x,y$ in m; $u,v$ in m/s; and $T$ in $\mathrm{s}^{-2}$ .

Term	Units
$-4T(xu+yv)+3gv$	$\mathrm{m}^2\,\mathrm{s}^{-3}$
$x^2+y^2$	$\mathrm{m}^2$
quotient, $\dot T$	$\mathrm{s}^{-3}$

The denominator is required dimensionally. Normalizing the rod length to one can hide this factor algebraically, but not physically.

Discussion

Are all of the algebraic constraints in the original formulation satisfied?

Formulation 3: Index-1 DAE Model¶

\begin{align} \frac{d y}{dt} &= v \\ \frac{d v}{dt} &= g - Ty \\ & x^2 + y^2 = 1 \\ & 2 x u + 2 y v = 0 \\ & (u^2 + v^2) - T \underbrace{(x^2 + y^2)}_{1} + g y = 0 \end{align}

(22)

(This reformulation is NOT unique... could have written $\frac{dx}{dt}$ and $\frac{du}{dt}$ instead.)

Consistent initial conditions:

Specify $y(0)$ and $v(0)$ .
Solve for $x(0)$ , $u(0)$ , and $T(0)$

def create_model_index1():
    """
    Creates a model to define the Index 1 DAE system.

    Output: the model, m
    """

    m = pyo.ConcreteModel()

    # Declare time
    m.t = dae.ContinuousSet(bounds=(0, 5))

    # Declare parameter - acceleration due to gravity
    m.g = pyo.Param(initialize=9.81)  # [m/s^2]

    # Declare variables indexed over time
    m.x = pyo.Var(m.t)  # horizontal position
    m.y = pyo.Var(m.t)  # vertical position
    m.u = pyo.Var(m.t)  # horizontal velocity
    m.v = pyo.Var(m.t)  # vertical velocity
    m.T = pyo.Var(m.t)  # tension

    # Declare derivative variables
    m.dy = dae.DerivativeVar(m.y)
    m.dv = dae.DerivativeVar(m.v)

    # Declare differential equations
    def _dy_eqn(m, t):
        return m.dy[t] == m.v[t]

    m.dy_eqn = pyo.Constraint(m.t, rule=_dy_eqn)

    def _dv_eqn(m, t):
        return m.dv[t] == m.g - m.T[t] * m.y[t]

    m.dv_eqn = pyo.Constraint(m.t, rule=_dv_eqn)

    # Declare algebraic equations
    def _alg_eqn1(m, t):
        return m.x[t] ** 2 + m.y[t] ** 2 == 1

    m.alg_eqn1 = pyo.Constraint(m.t, rule=_alg_eqn1)

    def _alg_eqn2(m, t):
        return m.x[t] * m.u[t] + m.y[t] * m.v[t] == 0

    m.alg_eqn2 = pyo.Constraint(m.t, rule=_alg_eqn2)

    def _alg_eqn3(m, t):
        return (
            m.u[t] ** 2
            + m.v[t] ** 2
            - m.T[t] * (m.x[t] ** 2 + m.y[t] ** 2)
            + m.g * m.y[t]
            == 0
        )

    m.alg_eqn3 = pyo.Constraint(m.t, rule=_alg_eqn3)

    # Specify initial conditions
    m.x[0] = 0
    m.y[0] = 1
    m.u[0] = 1
    m.v[0] = 0
    m.T[0] = 1 + m.g

    return m

def index1_check_constraints(m):
    """Check if the three constraints are feasible."""

    print("Constraint 1:")
    r1 = m.x[0]() ** 2 + m.y[0]() ** 2 - 1
    print(r1)

    print("\nConstraint 2:")
    r2 = m.x[0]() * m.u[0]() + m.y[0]() * m.v[0]()
    print(r2)

    print("\nConstraint 3:")
    r3 = m.u[0]() ** 2 + m.v[0]() ** 2 - m.T[0]() + m.g * m.y[0]()
    print(r3)


index1 = create_model_index1()

# Check initial condition
index1_check_constraints(index1)

# Specify integrator options
int_ops = {"print_stats": True, "abstol": 1e-6, "reltol": 1e-4}

# Solve DAEs
sim = Simulator(index1, package="casadi")
try:
    tsim, profiles = sim.simulate(
        numpoints=100, integrator="idas", integrator_options=int_ops
    )
except RuntimeError as err:
    # Expected here: the reformulation is singular at x = 0.
    print(f"Expected point-singularity failure: {err}")
# tsim, profiles = sim.simulate(numpoints=100, integrator='collocation')

Constraint 1:
0

Constraint 2:
0

Constraint 3:
0.0

The residual routine or the linear setup or solve routine had a recoverable error, but IDACalcIC was unable to recover.

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[5], line 96
     92 int_ops = {"print_stats": True, "abstol": 1e-6, "reltol": 1e-4}
     93 
     94 # Solve DAEs
     95 sim = Simulator(index1, package="casadi")
---> 96 tsim, profiles = sim.simulate(
     97     numpoints=100, integrator="idas", integrator_options=int_ops
     98 )
     99 # tsim, profiles = sim.simulate(numpoints=100, integrator='collocation')

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/pyomo/dae/simulator.py:935, in Simulator.simulate(self, numpoints, tstep, integrator, varying_inputs, initcon, integrator_options)
    931         tsim, profile = self._simulate_with_casadi_with_inputs(
    932             initcon, tsim, varying_inputs, integrator, integrator_options
    933         )
    934     else:
--> 935         tsim, profile = self._simulate_with_casadi_no_inputs(
    936             initcon, tsim, integrator, integrator_options
    937         )
    939 self._tsim = tsim
    940 self._simsolution = profile

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/pyomo/dae/simulator.py:997, in Simulator._simulate_with_casadi_no_inputs(self, initcon, tsim, integrator, integrator_options)
    995 integrator_options['output_t0'] = True
    996 F = casadi.integrator('F', integrator, dae, integrator_options)
--> 997 sol = F(x0=initcon)
    998 profile = sol['xf'].full().T
   1000 if len(self._algvars) != 0:

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/casadi/casadi.py:24049, in Function.__call__(self, *args, **kwargs)
  24045         else:
  24046           return tuple(ret)
  24047       else:
  24048     # Named inputs -> return dictionary
> 24049         return self.call(kwargs)

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/casadi/casadi.py:20698, in Function.call(self, *args)
  20694 
  20695 
  20696 
  20697         """
> 20698         return _casadi.Function_call(self, *args)

RuntimeError: Error in Function::call for 'F' [IdasInterface] at .../casadi/core/function.cpp:1466:
Error in Function::call for 'F' [IdasInterface] at .../casadi/core/function.cpp:362:
.../casadi/interfaces/sundials/idas_interface.cpp:599: IDACalcIC returned "IDA_NO_RECOVERY". Consult IDAS documentation.

What happened? Point singularity at $x=0$ .

Let’s try $x=0.1$ as the initial point.

index1_again = create_model_index1()

# Specify alternative initial conditions
small_number = 0.1

index1_again.x[0] = small_number
index1_again.y[0] = 1
index1_again.u[0] = 1
index1_again.v[0] = 0
index1_again.T[0] = 1 + index1_again.g

# Check initial condition
index1_check_constraints(index1_again)

# Solve DAEs
sim = Simulator(index1_again, package="casadi")

# Simulator
try:
    tsim, profiles = sim.simulate(
        numpoints=100, integrator="idas", integrator_options=int_ops
    )
except RuntimeError as err:
    # Expected here: these trial initial values violate the invariants.
    print(f"Expected inconsistent-initialization failure: {err}")

Constraint 1:
0.010000000000000009

Constraint 2:
0.1

Constraint 3:
0.0

The residual routine or the linear setup or solve routine had a recoverable error, but IDACalcIC was unable to recover.

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[6], line 19
     15 # Solve DAEs
     16 sim = Simulator(index1_again, package="casadi")
     17 
     18 # Simulator
---> 19 tsim, profiles = sim.simulate(
     20     numpoints=100, integrator="idas", integrator_options=int_ops
     21 )

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/pyomo/dae/simulator.py:935, in Simulator.simulate(self, numpoints, tstep, integrator, varying_inputs, initcon, integrator_options)
    931         tsim, profile = self._simulate_with_casadi_with_inputs(
    932             initcon, tsim, varying_inputs, integrator, integrator_options
    933         )
    934     else:
--> 935         tsim, profile = self._simulate_with_casadi_no_inputs(
    936             initcon, tsim, integrator, integrator_options
    937         )
    939 self._tsim = tsim
    940 self._simsolution = profile

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/pyomo/dae/simulator.py:997, in Simulator._simulate_with_casadi_no_inputs(self, initcon, tsim, integrator, integrator_options)
    995 integrator_options['output_t0'] = True
    996 F = casadi.integrator('F', integrator, dae, integrator_options)
--> 997 sol = F(x0=initcon)
    998 profile = sol['xf'].full().T
   1000 if len(self._algvars) != 0:

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/casadi/casadi.py:24049, in Function.__call__(self, *args, **kwargs)
  24045         else:
  24046           return tuple(ret)
  24047       else:
  24048     # Named inputs -> return dictionary
> 24049         return self.call(kwargs)

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/casadi/casadi.py:20698, in Function.call(self, *args)
  20694 
  20695 
  20696 
  20697         """
> 20698         return _casadi.Function_call(self, *args)

RuntimeError: Error in Function::call for 'F' [IdasInterface] at .../casadi/core/function.cpp:1466:
Error in Function::call for 'F' [IdasInterface] at .../casadi/core/function.cpp:362:
.../casadi/interfaces/sundials/idas_interface.cpp:599: IDACalcIC returned "IDA_NO_RECOVERY". Consult IDAS documentation.

Our initial point does not satisfy the algebraic constraints! We need a consistent initial point.

Given $x = \epsilon$ , solve $x^2 + y^2 = 1$ for $y$ :

y = \sqrt{1^2 - \epsilon^2}

(23)

Then, assume $u = 1$ and solve $2 x u + 2 y v = 0$ for $v$ :

v = \frac{-x u}{y}

(24)

Finally, we can solve $(u^2 + v^2) - T (x^2 + y^2) + g y = 0$ for $T$ :

T = \frac{(u^2 + v^2) + gy}{x^2 + y^2} = \frac{(u^2 + v^2) + gy}{1}

(25)

index1_take_two = create_model_index1()

# Specify alternative initial conditions
small_number = 0.1

index1_take_two.x[0] = small_number
index1_take_two.y[0] = np.sqrt(1 - small_number**2)
index1_take_two.u[0] = 1
index1_take_two.v[0] = (
    -index1_take_two.x[0]() * index1_take_two.u[0]() / index1_take_two.y[0]()
)
index1_take_two.T[0] = (
    index1_take_two.u[0]() ** 2
    + index1_take_two.v[0]() ** 2
    + index1_take_two.g * index1_take_two.y[0]()
)

# Check initial condition
index1_check_constraints(index1_take_two)

# Solve DAEs
sim = Simulator(index1_take_two, package="casadi")

# Specify integrator options
int_ops2 = {
    "print_stats": True,
    "abstol": 1e-6,
    "reltol": 1e-4,
    "verbose": False,
    "calc_ic": True,
}

# Simulator
try:
    tsim, profiles = sim.simulate(
        numpoints=20, integrator="idas", integrator_options=int_ops2
    )
except RuntimeError as err:
    # This documented failure motivates discretizing the model next.
    print(f"Expected formulation failure: {err}")

Constraint 1:
0.0

Constraint 2:
0.0

Constraint 3:
0.0

The residual routine or the linear setup or solve routine had a recoverable error, but IDACalcIC was unable to recover.

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[7], line 34
     30     "calc_ic": True,
     31 }
     32 
     33 # Simulator
---> 34 tsim, profiles = sim.simulate(
     35     numpoints=20, integrator="idas", integrator_options=int_ops2
     36 )

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/pyomo/dae/simulator.py:935, in Simulator.simulate(self, numpoints, tstep, integrator, varying_inputs, initcon, integrator_options)
    931         tsim, profile = self._simulate_with_casadi_with_inputs(
    932             initcon, tsim, varying_inputs, integrator, integrator_options
    933         )
    934     else:
--> 935         tsim, profile = self._simulate_with_casadi_no_inputs(
    936             initcon, tsim, integrator, integrator_options
    937         )
    939 self._tsim = tsim
    940 self._simsolution = profile

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/pyomo/dae/simulator.py:997, in Simulator._simulate_with_casadi_no_inputs(self, initcon, tsim, integrator, integrator_options)
    995 integrator_options['output_t0'] = True
    996 F = casadi.integrator('F', integrator, dae, integrator_options)
--> 997 sol = F(x0=initcon)
    998 profile = sol['xf'].full().T
   1000 if len(self._algvars) != 0:

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/casadi/casadi.py:24049, in Function.__call__(self, *args, **kwargs)
  24045         else:
  24046           return tuple(ret)
  24047       else:
  24048     # Named inputs -> return dictionary
> 24049         return self.call(kwargs)

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/casadi/casadi.py:20698, in Function.call(self, *args)
  20694 
  20695 
  20696 
  20697         """
> 20698         return _casadi.Function_call(self, *args)

RuntimeError: Error in Function::call for 'F' [IdasInterface] at .../casadi/core/function.cpp:1466:
Error in Function::call for 'F' [IdasInterface] at .../casadi/core/function.cpp:362:
.../casadi/interfaces/sundials/idas_interface.cpp:599: IDACalcIC returned "IDA_NO_RECOVERY". Consult IDAS documentation.

Hmm, this does not make sense. Perhaps there is something subtle wrong with the model.

Let’s try solving the model with Ipopt after discretizing with collocation.

# discretize the model
index1_take_two.Obj = pyo.Objective(expr=1)  # Add a dummy objective
discretizer = pyo.TransformationFactory("dae.collocation")
discretizer.apply_to(index1_take_two, nfe=15, scheme="LAGRANGE-RADAU", ncp=3)

# initialize
for t in index1_take_two.t:
    index1_take_two.x[t] = small_number
    index1_take_two.y[t] = np.sqrt(1 - small_number**2)
    index1_take_two.u[t] = 1
    index1_take_two.v[t] = (
        -index1_take_two.x[t]() * index1_take_two.u[t]() / index1_take_two.y[t]()
    )
    index1_take_two.T[t] = (
        index1_take_two.u[t]() ** 2
        + index1_take_two.v[t]() ** 2
        + index1_take_two.g * index1_take_two.y[t]()
    )

# solve the discretized model
solver = pyo.SolverFactory("ipopt")
solver.options["max_iter"] = 300
results = solver.solve(index1_take_two, tee=True)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

Ipopt 3.14.19: max_iter=300


******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:     1186
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      368

Total number of variables............................:      322
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:      320
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0000000e+00 9.07e-01 0.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.0000000e+00 2.58e+00 1.51e-21  -1.7 4.10e+00    -  1.00e+00 1.00e+00h  1

   2  1.0000000e+00 3.54e+00 9.17e-03  -1.7 2.11e+00  -4.0 1.00e+00 1.00e+00h  1

   3  1.0000000e+00 3.17e+00 2.87e-01  -1.7 1.78e+00  -3.6 1.00e+00 1.00e+00h  1
   4  1.0000000e+00 6.54e-01 7.16e-02  -1.7 2.49e+00    -  1.00e+00 1.00e+00h  1

   5  1.0000000e+00 1.56e-01 1.79e-02  -1.7 3.95e-01  -3.1 1.00e+00 1.00e+00h  1
   6  1.0000000e+00 3.66e-02 4.48e-03  -2.5 1.92e-01  -2.7 1.00e+00 1.00e+00h  1

   7  1.0000000e+00 1.87e-02 5.74e-02  -2.5 2.25e-01  -2.3 1.00e+00 5.00e-01h  2
   8  1.0000000e+00 1.84e-02 1.44e-02  -2.5 1.39e-01    -  1.00e+00 1.00e+00h  1
   9  1.0000000e+00 5.18e-03 3.59e-03  -3.8 7.39e-02    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.0000000e+00 1.30e-03 8.97e-04  -3.8 3.69e-02    -  1.00e+00 1.00e+00h  1
  11  1.0000000e+00 3.24e-04 2.24e-04  -5.7 1.85e-02    -  1.00e+00 1.00e+00h  1

  12  1.0000000e+00 8.10e-05 5.61e-05  -5.7 9.23e-03    -  1.00e+00 1.00e+00h  1
  13  1.0000000e+00 2.02e-05 1.40e-05  -5.7 4.62e-03    -  1.00e+00 1.00e+00h  1
  14  1.0000000e+00 5.06e-06 3.50e-06  -5.7 2.31e-03    -  1.00e+00 1.00e+00h  1
  15  1.0000000e+00 1.26e-06 8.76e-07  -8.6 1.15e-03    -  1.00e+00 1.00e+00h  1
  16  1.0000000e+00 3.16e-07 2.19e-07  -8.6 5.77e-04    -  1.00e+00 1.00e+00h  1

  17  1.0000000e+00 7.91e-08 5.48e-08  -8.6 2.88e-04    -  1.00e+00 1.00e+00h  1
  18  1.0000000e+00 1.98e-08 1.37e-08  -8.6 1.44e-04    -  1.00e+00 1.00e+00h  1
  19  1.0000000e+00 4.94e-09 3.42e-09  -9.0 7.21e-05    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 19

                                   (scaled)                 (unscaled)
Objective...............:   1.0000000000000000e+00    1.0000000000000000e+00
Dual infeasibility......:   3.4227969388363734e-09    3.4227969388363734e-09
Constraint violation....:   4.9410901192814551e-09    4.9410901192814551e-09
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   4.9410901192814551e-09    4.9410901192814551e-09


Number of objective function evaluations             = 22
Number of objective gradient evaluations             = 20
Number of equality constraint evaluations            = 22
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 20
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 19
Total seconds in IPOPT                               = 0.330

EXIT: Optimal Solution Found.

Take away: There is something strange with formulation 3 that is causing the numeric integrator to fail. We can still solve this problem after discretizing.

Formulation 4: Index-1 DAE Model¶

\begin{align} \frac{d x}{dt} &= u \\ \frac{d y}{dt} &= v \\ \frac{d u}{dt} &= -T x \\ \frac{d v}{dt} &= g - Ty \\ 0 &= (u^2 + v^2) - T (x^2 + y^2) + g y \end{align}

(26)

Consistent initial conditions:

Specify $x(0)$ , $u(0)$ , $y(0)$ , and $v(0)$ .
Solve for $T(0)$

def create_model_index1_b():
    """
    Creates a model to define the Index 1 DAE system.

    Output: the model, m
    """

    m = pyo.ConcreteModel()

    # Declare time
    m.t = dae.ContinuousSet(bounds=(0.0, 5))

    # Declare parameter - acceleration due to gravity
    m.g = pyo.Param(initialize=9.81)  # [m/s^2]

    # Declare variables indexed over time
    m.x = pyo.Var(m.t)  # horizontal position
    m.y = pyo.Var(m.t)  # vertical position
    m.u = pyo.Var(m.t)  # horizontal velocity
    m.v = pyo.Var(m.t)  # vertical velocity
    m.T = pyo.Var(m.t)  # tension

    # Declare derivative variables
    m.dx = dae.DerivativeVar(m.x)  # with respect to t is implied
    m.dy = dae.DerivativeVar(m.y)
    m.du = dae.DerivativeVar(m.u)
    m.dv = dae.DerivativeVar(m.v)

    # Declare differential equations
    def _dx_eqn(m, t):
        return m.dx[t] == m.u[t]

    m.dx_eqn = pyo.Constraint(m.t, rule=_dx_eqn)

    def _dy_eqn(m, t):
        return m.dy[t] == m.v[t]

    m.dy_eqn = pyo.Constraint(m.t, rule=_dy_eqn)

    def _du_eqn(m, t):
        return m.du[t] == -m.T[t] * m.x[t]

    m.du_eqn = pyo.Constraint(m.t, rule=_du_eqn)

    def _dv_eqn(m, t):
        return m.dv[t] == m.g - m.T[t] * m.y[t]

    m.dv_eqn = pyo.Constraint(m.t, rule=_dv_eqn)

    def _alg_eqn3(m, t):
        return (
            m.u[t] ** 2
            + m.v[t] ** 2
            - m.T[t] * (m.x[t] ** 2 + m.y[t] ** 2)
            + m.g * m.y[t]
            == 0
        )

    m.alg_eqn3 = pyo.Constraint(m.t, rule=_alg_eqn3)

    # Specify initial conditions
    m.x[0] = 0
    m.y[0] = 1
    m.u[0] = 1
    m.v[0] = 0
    m.T[0] = 1 + m.g

    return m

index1_b = create_model_index1_b()

# Specify integrator options
int_ops = {"print_stats": True, "abstol": 1e-8, "reltol": 1e-6}

# Solve DAEs
sim = Simulator(index1_b, package="casadi")
tsim, profiles = sim.simulate(
    numpoints=100, integrator="idas", integrator_options=int_ops
)

# Plot solution
plot_results(sim, tsim, profiles)

FORWARD INTEGRATION:
Number of steps taken by SUNDIALS: 379
Number of calls to the user's f function: 485
Number of calls made to the linear solver setup function: 31
Number of error test failures: 3
Method order used on the last internal step: 5
Method order to be used on the next internal step: 5
Actual value of initial step size: 7.90569e-09
Step size taken on the last internal step: 0.0120796
Step size to be attempted on the next internal step: 0.0120796
Current internal time reached: 5
Number of nonlinear iterations performed: 483
Number of nonlinear convergence failures: 0

Take Away Messages¶

Differential algebraic equations (DAEs) are really powerful modeling tools.
Integrating DAEs requires special care. Make sure your model is index 1.
Often there are many ways to reformulate the DAE model. But the numeric integrator only enforces error tolerances on the equations that are explicitly modeled. If an algebraic constraint must be satisfied to a specific tolerance, include it in the DAE model (as long as it is not high index!).

Important Topics from the Textbook¶

Chapter 8: Dynamic Optimization Introduction (Biegler, 2010)

Chemical engineering examples
Classical (variational) approaches including Hamiltonian and Euler-Lagrange equations

Classification of DAE optimization strategies

Re-authored as TikZ from Biegler (2010), Figure 8.10, p. 245; source at figures/tikz/dae-optimization-strategies.tex.

Chapter 9: Sequential Methods (Biegler, 2010)

DAE integration
Single Shooting
Multiple Shooting

Sequential dynamic optimization workflow with an optimizer repeatedly calling a DAE integrator

Chapter 10: Simultaneous Methods (Biegler, 2010)

Gauss quadrature
Orthogonal collocation on finite elements
Examples, benchmarks, and large-scale extensions