Dynamic Optimization with Collocation and Pyomo.DAE

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

Guide to the transcription¶

Use this page for the full derivation behind Lecture 7: polynomial coordinates, integrated derivative basis and Runge-Kutta equations, continuity and counting, time scaling and implementation, and sparse structure.

The existing quadrature implementation below remains useful for experiments. The orthogonality proof is taught in the numerical-integration lecture and developed in Numeric Integration for DAEs. For exact coefficient construction, see Collocation Points and Runge-Kutta Coefficients; for a solved model, see Pyomo.DAE Example: Race Car.

Primary sources: Biegler (2010), Chapter 8, p. 216, and Chapter 10, pp. 289–298, with the Chapter 10 errata; Ascher & Petzold (1998), Section 4.7, pp. 101–104.

We now turn a DAE-constrained dynamic optimization problem into a sparse nonlinear program. We approximate state and control profiles with piecewise polynomials and enforce the DAE at collocation points. This simultaneous, or direct-transcription, approach is closely related to the implicit Runge--Kutta methods developed in the numerical-integration lecture.

Quadrature Methods¶

We begin with the quadrature ideas needed to derive collocation. The next notebooks apply the method to dynamic optimization examples. For more detail, see:

• Chapter 9 in Biegler (2010)

• Chapter 16, Gauss Quadrature and Multi-dimensional Integrals, in McClarren (2018)

Main Idea¶

The quadrature rule approximates the integral $I$ with weighted sum of function evaluations:

where:

• $x_l$ are nodes (also known as abscissas)

• $w_l$ are weights

• $L$ number of nodes and weights

Central questions:

• How to choose nodes and weights?

• How does this choice impact approximation error?

• How does $L$ impact approximation error?

Polynomial degree and coefficient count¶

A polynomial $p(t)=\sum_{j=0}^K a_jt^j$ has $K+1$ coefficients and degree at most $K$, with equality when $a_K\ne0$. Some older texts call that polynomial “order $K+1$.” Here degree denotes the highest power; order denotes the rate at which a numerical error decreases.

Gauss-Legendre quadrature¶

With the Gauss-Legendre nodes and weights, the $L$-point rule is exact for every polynomial of degree at most $2L-1$. No $L$-node rule with real nodes can be exact through degree $2L$: the squared node polynomial $\prod_l(x-x_l)^2$ is zero at all nodes but has a positive integral. This is optimal polynomial exactness, rather than a claim of minimum error for every nonpolynomial integrand.

The table gives nodes and weights for $L\le8$.

Why $2L-1$?¶

Let $f(x) = a_{2L-1} x^{2L-1} + a_{2L-2} x^{2L-2} + ... + a_{1} x^1 + a_{0}$

Exactness through degree $2L-1$ imposes $2L$ moment equations (one for each monomial $1,x,\ldots,x^{2L-1}$). These match the $2L$ unknown node and weight values in an $L$-point rule. This count motivates the attainable degree; orthogonal-polynomial theory supplies the nodes and proves exactness.

How to determine the nodes and weights?¶

One way to derive the Gauss-Legendre quadrature rules is by looking at the integral generic monomials of degree 0 up to $2L-1$ and setting each equal to the $L$ point Gauss-Legendre quadrature rule:

Notice that the $a_i$ constants cancel out of each equation so they do not matter. This system is $2L$ equations with $L$ weights, $w_l$, and $L$ nodes, $x_l$. We could solve these equations to get the weights and nodes, though this is not how it is done in practice generally---this is accomplished by using the theory of orthogonal polynomials.

Code¶

Visualize Weights¶

Polynomial Example¶

As a simple demonstration of the Gauss-Legendre quadrature, let’s show that it integrates polynomials of degree $2L-1$ exactly. Consider the integral

We can use the GLQuad function to evaluate the error.

Generalization¶

We are generally interested in integrals not just over the domain $x\in [-1,1]$. We’ll now make a function that does Gauss-Legendre quadrature over a general range using the formula

Let’s show that this version integrates polynomials of degree $2L-1$ exactly. Consider the integral

Radau Quadrature¶

Sometimes we desire either:

• The first node is at the beginning of the interval, i.e., $x_1 = -1$ or

• The last node is at the end of the interval, i.e., $x_{L} = 1$.

This is known as Gauss-Radau quadrature. The remaining nodes and weights are chosen optimally; thus, Radau quadrature exactly integrates all polynomials of degree $2L-2$ or less. Specifying one of the nodes leaves only $2L-1$ degrees of freedom.

More info: https://mathworld.wolfram.com/RadauQuadrature.html

Differential Equations¶

Consider solving the differential equation

to determine $\mathbf{x}(t_f)$. This can be expressed as an integral:

Thus, Runge-Kutta methods can be interpreted as a quadrature rule.

Collocation Methods¶

See the following resources for more information on collocation:

• Chapter 10 in Biegler (2010)

What is a Collocation Method?¶

https://encyclopediaofmath.org/wiki/Collocation_method

A collocation method approximates the state with a polynomial and enforces the differential and algebraic equations at selected collocation points. The resulting finite system of algebraic equations can be solved directly or embedded in an optimization problem.

Derivation of Collocation Methods¶

Derivations are followed from Section 10.2 of Biegler (2010).

Consider the following multiperiod dynamic optimization problem (8.5).

Each physical period $l$ may contain many finite elements. Period junctions use the coupling function $\mathbf h$; finite-element boundaries within a period use state continuity. The decision variables include the nodal states, control coefficients and any free parameters. A running cost can be included by adding an accumulator state, $\dot q=L(\mathbf z,\mathbf y,\mathbf u,\mathbf p,t)$, and minimizing $q(t_f)$ plus terminal cost.

This source-controlled course schematic was re-authored from the mathematical layout of Biegler (2010), Figure 10.1; it is not a reproduction of the textbook artwork.

State coordinates on one element¶

The collocation method can be developed by solving the differential equation ($dz/dt$) at specific time points. For example, Figure 10.1 shows the state variable approximation for a single finite element:

The coefficients $\boldsymbol\alpha_j$ are vectors when the state is a vector. A power basis is mathematically valid, but nodal Lagrange coefficients are directly interpretable as state values. This profile can also be represented as a Lagrange interpolation polynomial. For Gauss-Legendre or right Radau, choose $K+1$ distinct interpolation points within element $i$, the state can be represented as:

The key property of this polynomial is that $z^K(t_{ij})=z_{ij}$ where $t_{ij}=t_{i-1} + \tau_j h_i$.

At $\tau=\tau_k$, every $\ell_j$ with $j\ne k$ contains a zero numerator; $\ell_k(\tau_k)=1$. Thus $\mathbf z^K(t_{ik})=\mathbf z_{ik}$. The basis also satisfies $\sum_j\ell_j=1$: both sides interpolate the constant polynomial at $K+1$ distinct nodes. These identities justify interpreting coefficients as nodal values.

Why change coordinates? Power coefficients are not nodal states, so componentwise bounds on those coefficients do not directly impose physical state bounds. Monomial Vandermonde systems can also become poorly conditioned as the degree grows or the time interval is badly scaled; nearly dependent columns make coefficients sensitive even when the represented profile changes little. Normalizing each element and using nodal coordinates makes the variables easier to interpret. It does not remove every conditioning problem, and nodal bounds still need the between-node check below.

Bounds are nodal. Imposing $\mathbf z_L\le\mathbf z_{ij}\le\mathbf z_U$ does not bound the whole polynomial. Lagrange basis functions can be negative between nodes. Check additional points, extrema of the interpolation polynomial, or the independently integrated trajectory. For a continuous path constraint, mesh refinement and a residual check are part of validation.

The counterexample below uses $p(\tau)=4.2\tau(1-\tau)$ at the initial point and three right-Radau nodes. Every nodal value lies in $[0,1]$, but $p(1/2)=1.05$.

Integrated derivative basis gives Runge-Kutta equations¶

Let $K$ distinct derivative nodes be $c_j$. Define

The derivative polynomial has degree at most $K-1$; integrating it with the initial value gives

Here $\mathbf v_{ik}$ denotes the physical-time derivative, not the state. Differentiating cancels the factor $h_i$, and $\bar\ell_k(c_j)=\delta_{kj}$ gives $d\mathbf z^K/dt(t_{ij})=\mathbf v_{ij}$. The factor arises from $dt=h_i\,d\tau$.

For a semi-explicit DAE, enforce

Evaluating the state at the stages and endpoint gives

These are the Runge-Kutta stage and update equations. For a state-only ODE they are standard RK; for a DAE the algebraic stage equations are additional. A degree-$K$ state polynomial is fixed by its initial value and its derivative at $K$ distinct points, so this derivative representation is equivalent to the state-coordinate equations when the state interpolation nodes are distinct.

Not every node set is implicit. For one node $c$, $a_{11}=c$ and $b_1=1$. The choices $c=0,1/2,1$ give forward Euler, implicit midpoint and backward Euler. Standard Gauss, Radau IIA and Lobatto IIIA are implicit families (the first Lobatto stage is an explicit evaluation). For Lobatto the derivative nodes include zero; do not append a duplicate zero to the $K+1$-point state formula.

Partition of unity gives stage times. Since $\sum_k\bar\ell_k=1$, $\sum_k\Omega_k(\tau)=\tau$. Hence $\mathbf A\mathbf1=\mathbf c$ and $\mathbf b^T\mathbf1=1$. For $dt/dt=1$ this reproduces $t_{ij}=t_{i-1}+h_ic_j$ and $t_i=t_{i-1}+h_i$. With normalized local coordinate, $dt/d\tau=h_i$, not 1. Always use $b_k$ for the endpoint update; the last row of $\mathbf A$ equals $\mathbf b^T$ only when the last node is 1.

Continuity and degree-of-freedom counting¶

The incoming value on element $i+1$ equals the endpoint reconstructed on element $i$:

The state remains continuous when its derivative has a finite jump. Right Radau has $c_K=1$, so its endpoint is the last stage. Gauss-Legendre has $c_K<1$; it needs the endpoint evaluation, which generally differs from the last stage.

Controls may be piecewise constant and jump at a boundary. Do not impose continuity unless the physical actuator or chosen approximation requires it. Algebraic variables require a separate check. If $\mathbf g(\mathbf z,\mathbf y,\mathbf u,\mathbf p,t)=0$ defines a unique regular local branch $\mathbf y=\mathbf Y(\mathbf z,\mathbf u,\mathbf p,t)$ and its inputs are continuous, then $\mathbf y$ is continuous. Index 1 alone does not ensure this if $\mathbf u$ jumps. In $\dot z=y$, $0=y-u$, switching $u$ from 1 to 2 makes $y$ jump while $z$ stays continuous. At a discontinuity, distinguish the left and right values at the shared time.

For the common degree-$(K-1)$ control and algebraic profile representations, use the $K$-node derivative basis:

These profiles need not match across elements. A lower-order control parameterization can instead set all $\mathbf U_{ij}$ equal on an element. For a regular algebraic branch, an endpoint algebraic value may also be recovered by solving $\mathbf g=0$ from the endpoint differential state and the appropriate left/right control; this does not justify equating both sides when the control jumps. Nodal bounds on either interpolated profile still need a between-node check.

For one differential state, fixed mesh and prescribed controls/parameters, an expanded derivative transcription has:

There are $N(2K+1)+1$ variables and $N(2K+1)$ equations before fixing $z_0$. The initial condition removes the remaining degree of freedom. Each regular algebraic stage variable adds a variable and its equation. Free controls, parameters, durations and terminal constraints change the count. Pyomo may eliminate boundary/stage duplicates or use a different derivative representation; its raw counts need not equal this pedagogical expanded formulation. A square count is necessary for a regular simulation, but does not prove independence or numerical accuracy.

A free final time on a fixed domain¶

Let $\Delta t=t_f-t_0>0$ and $\sigma=(t-t_0)/\Delta t\in[0,1]$. Then

The transformed derivative is with respect to $\sigma$. The duration multiplies every physical-time right-hand side, and any explicit time argument must also be transformed. Minimize $t_0+\Delta t$ for minimum arrival time with fixed $t_0$. In Biegler’s notation $p_f$ is this duration; the car code calls it tf because $t_0=0$.

Normalization fixes the domain on which Pyomo creates points; it does not impose a positive duration. Add a physically meaningful lower bound, initialize within the bounds, and use units to check both sides. The car example uses tf >= 0.001 s. For moving normalized mesh fractions $\eta_i>0$, require $\sum_i\eta_i=1$ and use physical lengths $h_i=\Delta t\eta_i$; do not multiply by the duration twice. Alternatively optimize physical $h_i$ with $\sum_i h_i=\Delta t$.

Pyomo implementation choices¶

Transformation scheme reference¶

Here nfe is the number of finite elements and ncp the stages per element. The finite-difference transformation has no ncp option, and there is no Pyomo Lobatto scheme. These are spatial derivative formulas too; Euler interpretations refer to discretizing time, not to a time integrator created merely by discretizing space.

The official Pyomo.DAE documentation describes these transformations. Read the generated constraints when a derivative appears at an endpoint: declaration alone does not create its defining equation. For boundary-value problems, use the heat-conduction stencil and diagnostics examples. For initialization of a nonlinear process model, see Pyomo.DAE Example: Temperature Control Lab; use a simulated or coarse solution as a starting profile and check structural diagnostics before relying on a solver message.

Node families and endpoint accuracy¶

For Gauss-Legendre, the $K$ nodes are the roots of the shifted degree-$K$ Legendre polynomial. The corresponding quadrature integrates all degrees through $2K-1$. The orthogonality proof and its assumptions are in Numeric Integration for DAEs.

Right Radau adds the endpoint 1 to the roots of a degree-$(K-1)$ polynomial orthogonal with weight $1-\tau$. Lobatto adds 0 and 1 to the roots of a degree-$(K-2)$ polynomial orthogonal with weight $\tau(1-\tau)$. These are free-node Jacobi constructions, with the fixed endpoints added afterward. Exact symbolic tables and verification are in Collocation Points and Runge-Kutta Coefficients.

For a sufficiently smooth ODE, the standard Gauss, Radau IIA and Lobatto IIIA endpoint orders are $2K$, $2K-1$ and $2K-2$. An order-$p$ method has an unnormalized one-step defect $O(h^{p+1})$, or normalized local truncation error $O(h^p)$. Stage errors and DAE algebraic-variable errors need not have that same order. Quadrature moment tests alone do not prove every nonlinear Runge-Kutta order condition.

Gauss and Radau IIA are A-stable and algebraically stable. Algebraic stability implies AN-stability; the converse applies for distinct abscissas. Radau IIA also has homogeneous stiff decay. A-stability bounds amplification for left-half-plane eigenvalues; it does not ensure accurate transient decay or convergence of the nonlinear stage solve. Gauss large-step oscillations can be reduced by refining the mesh until the transient is resolved. For a very stiff problem, that can require far more elements than a method with stiff decay.

See Biegler, Section 10.2.2, pp. 290–293; Ascher & Petzold, pp. 101–104 and Exercise 4.13, p. 119. See the numerical-integration guide for the formal stability definitions, rather than treating their names as interchangeable.

A practical choice. Start with right Radau for a stiff model when damping fast homogeneous transients matters. For a smooth nonstiff problem, compare Gauss-Legendre: it has one higher endpoint order at the same stage count. That is a reason to measure error and work, not a universal speed ranking. Gauss oscillations can be resolved by sufficiently fine elements. Lobatto is useful for understanding endpoint choices but is not supplied by Pyomo.DAE.

Local equations and sparse structure¶

Order the variables by element. A collocation row involves stages in its own element, the incoming boundary state, that element’s control coefficients, and any shared parameters. Continuity connects neighboring elements. The Jacobian therefore has local blocks and a border for shared parameters or a free final time. It need not be a plain banded matrix. A stage derivative can couple every state component, so a block can itself be dense. Sparse NLP solvers exploit the incidence pattern and compute derivatives of the transcribed equations.

A simultaneous NLP iterate can violate the dynamic equations. At a converged point it satisfies the discrete residual tolerance, which still leaves a discretization error. In a sequential architecture, successful integration supplies a dynamic trajectory for a trial control; it does not automatically enforce path constraints, state bounds or terminal conditions. The architecture comparison is taught in Differential Algebraic Equations (DAEs).

Difficult cases and further reading¶

• Free element lengths add products of $h_i$ with stage derivatives. They let a discontinuity align with an element boundary, but require positive lengths and sensible movement bounds. Inactive pieces can make many meshes nearly equivalent, leaving weak curvature in the mesh variables. See the executed moving-element car example, and Biegler, pp. 296–298.

• Scale both variables and residuals. A typical-state scale changes $\mathbf z=\mathbf z_{\rm ref}+\mathbf S_z\hat{\mathbf z}$; a residual scale changes the nonlinear solver’s numerical balance. Preserve units before choosing either. A duration bound is a physical domain constraint, not a substitute for scaling.

• An index-1 semi-explicit DAE has locally nonsingular $\partial\mathbf g/\partial\mathbf y$. Higher-index DAEs and active path constraints can require differentiated consistency conditions or index reduction; general smooth-ODE order claims cannot be transferred unchanged. Biegler, Section 10.4.5, pp. 317–322, develops high-index path constraints.

• Spatial semi-discretization can generate a stiff ODE or DAE. The discretization, boundary equations and constraint index determine which integrator assumptions apply; the fact that the result is finite-dimensional is insufficient. See the heat-conduction guide.

Sources¶

Biegler, Nonlinear Programming (2010), single/multiperiod formulation p. 216; interpolation and integrated derivative representation pp. 289–290, corrected by errata2.pdf; nodes, continuity and RK coefficients pp. 291–293; state/control/algebraic profiles and free-element lengths pp. 296–298; sparse NLP considerations pp. 302–304. Ascher & Petzold, Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations (1998), Section 4.7, pp. 101–104, and method of lines pp. 12–13. The source books and their errata, rather than an unchecked transcription, determine the formulas.

Additional topic map¶