Differential Algebraic Equations (DAEs)

3.3. Differential Algebraic Equations (DAEs)#

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()

3.3.1. External Documentaiton 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).

3.3.2. Dynamic Optimization Overview#

3.3.2.1. 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 resepct 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 tpyically specified as initial value problems with initial coniditions at zero. These equations are experessed with resepct to an independent variable, \(t\), and are typically autnomous, or \(t\) does not explicilty appear in the equation.

DAEs are tpyically of the form:

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

\(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:

(3.2)#\[\begin{equation} \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 \end{equation}\]

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

(3.3)#\[\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)=0 \end{align}\]

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.

3.3.3. General Breakdown of DAEs#

With resepct to constructing DAEs:

\(z(t)\): Differential equations are typically dervied from conservation laws.
\(y(t)\) : Algebraic equations are typiclaly derived from constitutive and equilbrium 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 parameteres, initial conditions, and other steady state conditions.

3.3.4. 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 bondary conditions. With respect to solving DAE optimization problems, BVP DAEs are of the form:

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

In this form the solutions may be nonunique or not exist over a speficied 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) \(|| 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:

(3.5)#\[\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}\]

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.

3.3.5. 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:

(3.6)#\[\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}\]

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 desciption and further explanation of these applications can be found in Chapter 8 of the Biegler (2010) textbook.

3.3.6. 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 emliminated 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.

3.3.6.1. DAE Index and Index Reduction#

Definition 8.6: Consider the DAE systems of the form:

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

OR

(3.8)#\[\begin{equation} \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 \end{equation}\]

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 no 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:

(3.9)#\[\begin{equation} \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. \end{equation}\]

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.

3.3.6.2. 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

Figure made with TikZ

The following DAE system describes the pendulum system:

(3.10)#\[\begin{align} & x' = u \\ & y' = v \\ & u' = -Tx \\ & v' = g - Ty \\ & x^2 + y^2 = 1 \\ \end{align}\]

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

(3.11)#\[\begin{align} & x^2 + y^2 = 1 \\ & 2xx' + 2yy' = 0 \\ & 2xu + 2yv = 0 \end{align}\]

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

(3.12)#\[\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}\]

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

(3.13)#\[\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}\]

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:

(3.14)#\[\begin{align} & 2xx' + 2yy' = 0\\ \end{align}\]

which can be solved for \(x'\). Substitute for \(x'\) and \(y'\) to yield:

(3.15)#\[\begin{align} & 2xu + 2yv = 0\\ \end{align}\]

The new DAE system is now:

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

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

(3.17)#\[\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}\]

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

(3.18)#\[\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}\]

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))
    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)

    
    # 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)
    
    # 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('$\sqrt{x^2 + y^2}$')
    plt.show()
    
    #return results

3.3.6.3. Formulation 1: Index-3 DAE#

Consider the following model:

\[\begin{split}\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}\end{split}\]

This assumes mass and length of unity.

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')
tsim, profiles = sim.simulate(numpoints=100, integrator='idas',integrator_options=int_ops)

# Plot solution
plot_results(sim, tsim, profiles)

psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
At t = 0 and h = 1.06624e-14, the corrector convergence failed repeatedly or with |h| = hmin.

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/var/folders/8v/lbvy4xb128d687t1pwzl_cfh0000gp/T/ipykernel_66414/739539068.py in ?()
     63 int_ops = {'print_stats':True,"abstol":1E-8,"reltol":1E-6}
     64 
     65 # Solve DAEs
     66 sim = Simulator(index3, package='casadi')
---> 67 tsim, profiles = sim.simulate(numpoints=100, integrator='idas',integrator_options=int_ops)
     68 
     69 # Plot solution
     70 plot_results(sim, tsim, profiles)

~/anaconda3/envs/controls/lib/python3.10/site-packages/pyomo/dae/simulator.py in ?(self, numpoints, tstep, integrator, varying_inputs, initcon, integrator_options)
    928                 tsim, profile = self._simulate_with_casadi_with_inputs(
    929                     initcon, tsim, varying_inputs, integrator, integrator_options
    930                 )
    931             else:
--> 932                 tsim, profile = self._simulate_with_casadi_no_inputs(
    933                     initcon, tsim, integrator, integrator_options
    934                 )
    935 

~/anaconda3/envs/controls/lib/python3.10/site-packages/pyomo/dae/simulator.py in ?(self, initcon, tsim, integrator, integrator_options)
    990 
    991         integrator_options['grid'] = tsim
    992         integrator_options['output_t0'] = True
    993         F = casadi.integrator('F', integrator, dae, integrator_options)
--> 994         sol = F(x0=initcon)
    995         profile = sol['xf'].full().T
    996 
    997         if len(self._algvars) != 0:

~/anaconda3/envs/controls/lib/python3.10/site-packages/casadi/casadi.py in ?(self, *args, **kwargs)
   7768         else:
   7769           return tuple(ret)
   7770       else:
   7771     # Named inputs -> return dictionary
-> 7772         return self.call(kwargs)

~/anaconda3/envs/controls/lib/python3.10/site-packages/casadi/casadi.py in ?(self, *args)
   6941         pre-CasADi 3.2
   6942 
   6943 
   6944         """
-> 6945         return _casadi.Function_call(self, *args)

RuntimeError: .../casadi/interfaces/sundials/idas_interface.cpp:591: 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:

casadi-error-1

casadi-error-2

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!

3.3.6.4. Formulation 2: Pure ODE Model/ Index-0 Example#

\[\begin{split}\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} &= 4 T (x u + y v) + 3 g v \end{align}\end{split}\]

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.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: 167
Number of calls to the user’s f function: 242
Number of calls made to the linear solver setup function: 31
Number of error test failures: 7
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.00246466
Step size to be attempted on the next internal step: 0.00492933
Current internal time reached: 0.00492933
Number of nonlinear iterations performed: 240
Number of nonlinear convergence failures: 0

../../_images/9f839cf4f05e631514896e7768e846b198be8192a67fbf1ae38f5f64e612d036.png

../../_images/e5c8ad54593f7664c7d2efe82b4fcff0d80ae5a0329747bd6bec75a459899e86.png

Discussion

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

3.3.6.5. Formulation 3: Index-1 DAE Model#

\[\begin{split}\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}\end{split}\]

(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')
tsim, profiles = sim.simulate(numpoints=100, integrator='idas',integrator_options=int_ops)
# tsim, profiles = sim.simulate(numpoints=100, integrator='collocation')

Constraint 1:
0

Constraint 2:
0

Constraint 3:
0.0

psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
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)
/var/folders/8v/lbvy4xb128d687t1pwzl_cfh0000gp/T/ipykernel_66414/4113190258.py in ?()
     80 int_ops = {'print_stats':True,"abstol":1E-6,"reltol":1E-4}
     81 
     82 # Solve DAEs
     83 sim = Simulator(index1, package='casadi')
---> 84 tsim, profiles = sim.simulate(numpoints=100, integrator='idas',integrator_options=int_ops)
     85 # tsim, profiles = sim.simulate(numpoints=100, integrator='collocation')

~/anaconda3/envs/controls/lib/python3.10/site-packages/pyomo/dae/simulator.py in ?(self, numpoints, tstep, integrator, varying_inputs, initcon, integrator_options)
    928                 tsim, profile = self._simulate_with_casadi_with_inputs(
    929                     initcon, tsim, varying_inputs, integrator, integrator_options
    930                 )
    931             else:
--> 932                 tsim, profile = self._simulate_with_casadi_no_inputs(
    933                     initcon, tsim, integrator, integrator_options
    934                 )
    935 

~/anaconda3/envs/controls/lib/python3.10/site-packages/pyomo/dae/simulator.py in ?(self, initcon, tsim, integrator, integrator_options)
    990 
    991         integrator_options['grid'] = tsim
    992         integrator_options['output_t0'] = True
    993         F = casadi.integrator('F', integrator, dae, integrator_options)
--> 994         sol = F(x0=initcon)
    995         profile = sol['xf'].full().T
    996 
    997         if len(self._algvars) != 0:

~/anaconda3/envs/controls/lib/python3.10/site-packages/casadi/casadi.py in ?(self, *args, **kwargs)
   7768         else:
   7769           return tuple(ret)
   7770       else:
   7771     # Named inputs -> return dictionary
-> 7772         return self.call(kwargs)

~/anaconda3/envs/controls/lib/python3.10/site-packages/casadi/casadi.py in ?(self, *args)
   6941         pre-CasADi 3.2
   6942 
   6943 
   6944         """
-> 6945         return _casadi.Function_call(self, *args)

RuntimeError: .../casadi/interfaces/sundials/idas_interface.cpp:591: 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
tsim, profiles = sim.simulate(numpoints=100, integrator='idas',integrator_options=int_ops)

Constraint 1:
0.010000000000000009

Constraint 2:
0.1

Constraint 3:
0.0

psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
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)
/var/folders/8v/lbvy4xb128d687t1pwzl_cfh0000gp/T/ipykernel_66414/575022909.py in ?()
     15 # Solve DAEs
     16 sim = Simulator(index1_again, package='casadi')
     17 
     18 # Simulator
---> 19 tsim, profiles = sim.simulate(numpoints=100, integrator='idas',integrator_options=int_ops)

~/anaconda3/envs/controls/lib/python3.10/site-packages/pyomo/dae/simulator.py in ?(self, numpoints, tstep, integrator, varying_inputs, initcon, integrator_options)
    928                 tsim, profile = self._simulate_with_casadi_with_inputs(
    929                     initcon, tsim, varying_inputs, integrator, integrator_options
    930                 )
    931             else:
--> 932                 tsim, profile = self._simulate_with_casadi_no_inputs(
    933                     initcon, tsim, integrator, integrator_options
    934                 )
    935 

~/anaconda3/envs/controls/lib/python3.10/site-packages/pyomo/dae/simulator.py in ?(self, initcon, tsim, integrator, integrator_options)
    990 
    991         integrator_options['grid'] = tsim
    992         integrator_options['output_t0'] = True
    993         F = casadi.integrator('F', integrator, dae, integrator_options)
--> 994         sol = F(x0=initcon)
    995         profile = sol['xf'].full().T
    996 
    997         if len(self._algvars) != 0:

~/anaconda3/envs/controls/lib/python3.10/site-packages/casadi/casadi.py in ?(self, *args, **kwargs)
   7768         else:
   7769           return tuple(ret)
   7770       else:
   7771     # Named inputs -> return dictionary
-> 7772         return self.call(kwargs)

~/anaconda3/envs/controls/lib/python3.10/site-packages/casadi/casadi.py in ?(self, *args)
   6941         pre-CasADi 3.2
   6942 
   6943 
   6944         """
-> 6945         return _casadi.Function_call(self, *args)

RuntimeError: .../casadi/interfaces/sundials/idas_interface.cpp:591: 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}\]

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

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

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} \]

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
tsim, profiles = sim.simulate(numpoints=20, integrator='idas',integrator_options=int_ops2)

Constraint 1:
0.0

Constraint 2:
0.0

Constraint 3:
0.0

psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:852: Linear solve failed
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)
/var/folders/8v/lbvy4xb128d687t1pwzl_cfh0000gp/T/ipykernel_66414/2451325184.py in ?()
     20 int_ops2 = {'print_stats':True,"abstol":1E-6,"reltol":1E-4,
     21             "verbose":False,"calc_ic":True}
     22 
     23 # Simulator
---> 24 tsim, profiles = sim.simulate(numpoints=20, integrator='idas',integrator_options=int_ops2)

~/anaconda3/envs/controls/lib/python3.10/site-packages/pyomo/dae/simulator.py in ?(self, numpoints, tstep, integrator, varying_inputs, initcon, integrator_options)
    928                 tsim, profile = self._simulate_with_casadi_with_inputs(
    929                     initcon, tsim, varying_inputs, integrator, integrator_options
    930                 )
    931             else:
--> 932                 tsim, profile = self._simulate_with_casadi_no_inputs(
    933                     initcon, tsim, integrator, integrator_options
    934                 )
    935 

~/anaconda3/envs/controls/lib/python3.10/site-packages/pyomo/dae/simulator.py in ?(self, initcon, tsim, integrator, integrator_options)
    990 
    991         integrator_options['grid'] = tsim
    992         integrator_options['output_t0'] = True
    993         F = casadi.integrator('F', integrator, dae, integrator_options)
--> 994         sol = F(x0=initcon)
    995         profile = sol['xf'].full().T
    996 
    997         if len(self._algvars) != 0:

~/anaconda3/envs/controls/lib/python3.10/site-packages/casadi/casadi.py in ?(self, *args, **kwargs)
   7768         else:
   7769           return tuple(ret)
   7770       else:
   7771     # Named inputs -> return dictionary
-> 7772         return self.call(kwargs)

~/anaconda3/envs/controls/lib/python3.10/site-packages/casadi/casadi.py in ?(self, *args)
   6941         pre-CasADi 3.2
   6942 
   6943 
   6944         """
-> 6945         return _casadi.Function_call(self, *args)

RuntimeError: .../casadi/interfaces/sundials/idas_interface.cpp:591: 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
solver.solve(index1_take_two,tee=True)

Ipopt 3.13.2: 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 http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

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 3.29e-01 1.89e-04  -1.7 1.14e+00  -4.0 1.00e+00 1.00e+00h  1
   2  1.0000000e+00 1.72e-01 8.57e-04  -1.7 4.16e-01  -3.6 1.00e+00 1.00e+00h  1
   3  1.0000000e+00 1.72e-01 8.43e-04  -2.5 2.05e+02  -4.1 1.00e+00 2.44e-04h 13
   4  1.0000000e+00 8.00e-02 3.19e-04  -2.5 3.03e-01  -3.6 1.00e+00 1.00e+00h  1
   5  1.0000000e+00 1.50e-01 1.98e-04  -2.5 6.09e-01    -  1.00e+00 1.00e+00h  1
   6  1.0000000e+00 3.53e-02 4.96e-05  -2.5 1.88e-01    -  1.00e+00 1.00e+00h  1
   7  1.0000000e+00 8.82e-03 1.24e-05  -2.5 9.41e-02    -  1.00e+00 1.00e+00h  1
   8  1.0000000e+00 2.20e-03 3.10e-06  -3.8 4.70e-02    -  1.00e+00 1.00e+00h  1
   9  1.0000000e+00 5.51e-04 7.74e-07  -3.8 2.35e-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.38e-04 1.94e-07  -5.7 1.18e-02    -  1.00e+00 1.00e+00h  1
  11  1.0000000e+00 3.45e-05 4.84e-08  -5.7 5.88e-03    -  1.00e+00 1.00e+00h  1
  12  1.0000000e+00 8.61e-06 1.21e-08  -5.7 2.94e-03    -  1.00e+00 1.00e+00h  1
  13  1.0000000e+00 2.15e-06 3.02e-09  -8.6 1.47e-03    -  1.00e+00 1.00e+00h  1
  14  1.0000000e+00 5.38e-07 7.56e-10  -8.6 7.35e-04    -  1.00e+00 1.00e+00h  1
  15  1.0000000e+00 1.35e-07 1.89e-10  -8.6 3.67e-04    -  1.00e+00 1.00e+00h  1
  16  1.0000000e+00 3.36e-08 4.73e-11  -8.6 1.84e-04    -  1.00e+00 1.00e+00h  1
  17  1.0000000e+00 8.41e-09 1.18e-11  -8.6 9.19e-05    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 17

                                   (scaled)                 (unscaled)
Objective...............:   1.0000000000000000e+00    1.0000000000000000e+00
Dual infeasibility......:   1.1815139856129635e-11    1.1815139856129635e-11
Constraint violation....:   8.4114670784174450e-09    8.4114670784174450e-09
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   8.4114670784174450e-09    8.4114670784174450e-09


Number of objective function evaluations             = 30
Number of objective gradient evaluations             = 18
Number of equality constraint evaluations            = 31
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 18
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 17
Total CPU secs in IPOPT (w/o function evaluations)   =      0.009
Total CPU secs in NLP function evaluations           =      0.001

EXIT: Optimal Solution Found.

{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 320, 'Number of variables': 322, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.12675833702087402}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}

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.

3.3.6.6. Formulation 4: Index-1 DAE Model#

\[\begin{split}\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}\end{split}\]