Direct Collocation for a Pendulum Swing-Up

This demo applies the direct-collocation recipe of the chapter Continuous-Time Transcription and Collocation to a small nonlinear control problem. The chapter develops the recipe for general dynamics and general choices of nodes. Here every step is carried out once with concrete numbers and then written as a few lines of Python, so that each symbol of the recipe can be matched with a quantity in the code. The text explains why each step is needed as well as what it computes.

The demo runs in about half a minute with NumPy, SciPy, and Matplotlib and needs no other files. Run the cells in order.

The Problem¶

A pendulum hangs at rest below its pivot, where a motor applies a torque $u(t)$. At time $T=8$ the pendulum must stand upright and at rest, and the motor should use as little torque as possible. Let $\theta$ be the angle measured from the hanging position and $\omega$ the angular velocity. Time is measured in units of $\sqrt{\ell/g}$ and torque in units of $mg\ell$, where $m$ is the mass, $\ell$ the length, and $g$ the gravitational acceleration. In these units the state $\mathbf x=(\theta,\omega)$ evolves according to

The term $-\sin\theta$ is gravity, which pulls the pendulum back toward the hanging position. Its largest magnitude is 1, reached when the pendulum is horizontal. The motor is weak: its torque is limited to $|u(t)|\le u_{\max}=0.5$, half of that largest gravity torque. It cannot lift the pendulum directly, and it must swing it back and forth to build up energy, which takes time. Writing the dynamics as $\dot{\mathbf x}=\mathbf f(\mathbf x,u)$, the optimal control problem is

Because of $\sin\theta$, the dynamics are nonlinear, and neither the optimal torque nor the resulting motion has a closed-form expression.

From Functions to Numbers¶

SciPy’s minimize, like every nonlinear programming solver, works with a finite vector $\mathbf z$ of unknowns. It needs a function $F(\mathbf z)$ that returns the objective and a function $H(\mathbf z)$ that returns constraint residuals, which must be zero at a solution. The problem above provides neither. Its unknowns are the functions $\theta(t)$, $\omega(t)$, and $u(t)$, its constraint is a condition on a derivative at every instant, and its objective is an integral. Transcription replaces these objects by finitely many numbers and algebraic equations.

Direct collocation makes three replacements:

1. The horizon is divided into $N$ intervals, and on each interval the state is described by a polynomial. A polynomial of low degree is determined by a few numbers, and its derivative and integral can be computed exactly from those numbers.

2. The differential equation is required only at a few chosen times per interval, the collocation nodes, where the slope of the polynomial must equal the slope given by $\mathbf f$.

3. The cost integral is replaced by a weighted sum of cost rates at the same nodes.

After these replacements, every condition is an equation between numbers, which the optimizer can evaluate for any candidate $\mathbf z$.

Many short intervals are used instead of one polynomial over $[0,T]$. This keeps the degree of each polynomial low, lets the accuracy improve by adding intervals, and makes each equation involve only the numbers of one interval.

The sections below follow the steps of the chapter’s integration-form algorithm. Step 1 computes fixed numbers that describe the polynomials. Steps 2 to 4 define $\mathbf z$, $F$, and $H$. Step 5 solves the problem and checks the answer against an accurate simulation. A final section builds the same problem from derivatives instead of integrals.

The first cell defines the model. The rest of the demo uses the dynamics only through calls f(x, u) and the cost only through calls c(x, u), so another model would change only this cell.

Step 1: Compute the Rule Once¶

Divide $[0,T]$ into $N$ intervals of width $h=T/N$, so that interval $k$ runs from $t_k=kh$ to $t_{k+1}=t_k+h$. Inside interval $k$, measure time by the local time $\tau=(t-t_k)/h$, which runs from 0 at the left end to 1 at the right end. Every interval then looks the same in local time, so the numbers computed in this step on $[0,1]$ serve every interval, whatever its position and width.

Three symbols for local time appear below, and they play different roles. The letter $\tau$ stands for any point of the interval. The collocation nodes $\tau_1$ and $\tau_2$ are two fixed values of $\tau$, chosen once and used in every interval. This demo uses $\tau_1=\tfrac{1}{3}$ and $\tau_2=1$, called the two-stage Radau nodes. The letter $\sigma$ is the variable of integration in the integrals below, which run from 0 to some value of $\tau$. The sketch below shows what the method builds. The state trajectory consists of one polynomial piece per interval. Consecutive pieces meet at the mesh states, although their slopes may differ there. Inside each interval, the stage states sit at the nodes, where the slope of the piece must equal the slope given by $\mathbf f$. The torque is constant on each interval. The right panel shows interval $k$ in local time, with the quantities that Step 2 turns into unknowns. The trajectory in the sketch is made up for illustration; the pendulum’s trajectory appears in Step 5.

The differential equation gives the slope of the state through a known function, $\dot{\mathbf x}=\mathbf f(\mathbf x,u)$. By the fundamental theorem of calculus, the state at local time $\tau$ in interval $k$ equals the state at the left end plus the integral of the slope:

The second form substitutes $t=t_k+h\sigma$, so that $dt=h d\sigma$. This is where the factor $h$ comes from: it converts a stretch of local time into physical time.

The function $\mathbf f$ is known, but the integrand also depends on the trajectory $\mathbf x(t)$, which is what we are trying to find, so this integral cannot be evaluated as it stands. Collocation evaluates the integrand only at the nodes and replaces it by the polynomial that takes those values. Write $\mathbf f_{k,j}=\mathbf f(\mathbf x_{k,j},u_k)$ for its value at node $j$, where $\mathbf x_{k,j}$ is the state at that node. Each value is computed from the known function $\mathbf f$. What is not known in advance is the state $\mathbf x_{k,j}$ at which to evaluate it, and Step 2 makes those states unknowns of the optimization. For now, treat $\mathbf f_{k,1}$ and $\mathbf f_{k,2}$ as given numbers.

The number of nodes is a choice of the method. Exactly one polynomial of degree at most $s-1$ passes through $s$ values at $s$ distinct nodes, so more nodes give a more flexible interpolant and usually a more accurate rule. This demo uses $s=2$ nodes, the smallest number whose interpolant can vary inside the interval, so the interpolating polynomial is the line

Each Lagrange basis function $\ell_j$ equals one at its own node $\tau_j$ and zero at the other nodes. In general it is the product

Each factor is zero at one of the other nodes and equals one at $\tau_j$. The denominators are numbers, so $\ell_j$ is a product of $s-1$ linear factors in $\sigma$, which is a polynomial of degree $s-1$. With two nodes each product has a single factor: $\ell_1(\sigma)=(\sigma-\tau_2)/(\tau_1-\tau_2)$ and $\ell_2(\sigma)=(\sigma-\tau_1)/(\tau_2-\tau_1)$. The experiments at the end try other numbers of nodes. Unlike the true integrand, a polynomial can be integrated exactly, which is the purpose of the replacement. Substituting the interpolant into the integral gives

The last equality uses the linearity of the integral. The integral of a sum is the sum of the integrals, and each $\mathbf f_{k,j}$ is a constant that does not depend on $\sigma$, so it moves outside its integral. After this exchange, the slopes only multiply integrals of the basis functions, and those integrals depend on the nodes alone. Evaluating them with $\tau$ equal to each node and to the right end gives the rule of the method:

The entry $A_{ij}$ is the share of slope $j$ in the state change from the left end to node $i$, per unit of interval width. The weight $b_j$ is the share of slope $j$ in the change over the whole interval.

These integrals can be computed exactly, because each $\ell_j$ is a polynomial. The product formula above is one way to write it. Multiplying out its linear factors writes the same polynomial as a sum of powers, $\ell_j(\sigma)=a_0+a_1\sigma+\cdots+a_d\sigma^d$, and for the nodes of this demo this is done by hand below. In that form, the linearity used above applies again: the integral of a sum is the sum of the integrals of its terms, and each coefficient $a_m$ moves outside its integral. Each remaining integral is a power of $\sigma$, which the power rule integrates:

The antiderivative $L_j(\tau)=\int_0^{\tau}\ell_j(\sigma)d\sigma$ is therefore again a polynomial, of one degree higher. The entries $A_{ij}=L_j(\tau_i)$ and $b_j=L_j(1)$ are evaluations of $L_j$, and no numerical approximation of an integral is involved. In code, NumPy’s Polynomial stores a polynomial as its coefficients $(a_0,a_1,\ldots,a_d)$ in increasing powers, and its .integ() method applies the power rule to them, which gives the coefficients $(0,a_0,\tfrac{a_1}{2},\ldots,\tfrac{a_d}{d+1})$. The leading 0 is the integration constant, chosen so that $L_j(0)=0$.

By hand, for $\tau_1=\tfrac{1}{3}$ and $\tau_2=1$, multiplying out the single factor of each product gives the basis functions in powers of $\tau$:

The power rule turns a constant $a$ into $a\tau$ and a term $a\tau$ into $\tfrac{a}{2}\tau^2$, so the antiderivatives are

Differentiating them returns $\ell_1$ and $\ell_2$, which checks the calculation, and both vanish at $\tau=0$. Each entry of the rule is now one evaluation:

The weights $b_j=L_j(1)$ repeat the second row, because $\tau_2=1$ is the right end of the interval:

The code performs the same operations for any distinct nodes: it multiplies the linear factors of each $\ell_j$, integrates with .integ(), and evaluates. It also prints the coefficient lists, which match the expressions above.

In lagrange_basis, Polynomial([0.0, 1.0]) is the polynomial whose coefficients are 0 and 1 in increasing powers, that is, $0+1\sigma$, or $\sigma$ itself. Arithmetic on Polynomial objects returns new polynomials, so sigma - tau_m is the polynomial $\sigma-\tau_m$, stored as the coefficients [-tau_m, 1.0]. The loop starts from the constant polynomial 1 and multiplies in one factor $(\sigma-\tau_m)/(\tau_j-\tau_m)$ for every other node, which is the product formula for $\ell_j$ written as code. The printed coefficient lists read in increasing powers of $\tau$. For example, [3/2, -3/2] is $\ell_1(\tau)=\tfrac{3}{2}-\tfrac{3}{2}\tau$, and .integ() turns it into [0, 3/2, -3/4], which is $L_1(\tau)=\tfrac{3}{2}\tau-\tfrac{3}{4}\tau^2$. Calling L_j(tau_i) evaluates that polynomial at $\tau_i$, which gives the integral from 0 to $\tau_i$. The printed arrays agree with the hand calculation.

Question. The printed row sums of $A$ equal the nodes, and the weights $b$ sum to one. Why must this hold for any choice of nodes? Hint: consider a motion whose slope is the same at every node.

The figure below draws each basis function and shows the entries of the rule as signed areas under it.

The orange dots mark the nodes, where each basis function equals one or zero. The entry $A_{12}=-\tfrac{1}{12}$ is negative because $\ell_2$ is negative between 0 and $\tfrac{1}{3}$, and a negative weight is not an error. The interpolated slope is a line through $\mathbf f_{k,1}$ at $\tau=\tfrac{1}{3}$ and $\mathbf f_{k,2}$ at $\tau=1$. Increasing $\mathbf f_{k,2}$ while $\mathbf f_{k,1}$ stays fixed tilts this line upward to the right of $\tau=\tfrac{1}{3}$ and downward to its left, so the state gains less before the first node.

Any distinct nodes produce a valid rule, but the choice affects accuracy. The Radau nodes make the weights $b$ integrate every polynomial of degree up to two exactly. For example, $\tfrac{3}{4}(\tfrac{1}{3})^2+\tfrac{1}{4}(1)^2=\tfrac{1}{3}$, which equals $\int_0^1\tau^2d\tau$. Two nodes at the ends of the interval, which give the trapezoidal rule, integrate only polynomials of degree up to one exactly. The experiments at the end compare several node choices on the pendulum.

The arrays $A$ and $b$ are computed once, before any optimization, and stay fixed while the optimizer searches.

Step 2: Allocate the Variables¶

The rule needs the slopes at the nodes, $\mathbf f_{k,j}=\mathbf f(\mathbf x_{k,j},u_k)$, and those depend on the states at the nodes, which are unknown before the problem is solved. Direct collocation does not compute these states by simulation. It makes them unknowns of the optimization and adds equations, in Step 3, that make them consistent. On interval $k$ the unknowns are a mesh state, the stage states, and a torque. The diagram after the next code cell shows them with their names in the code and their places in the vector $\mathbf z$.

The mesh state $\mathbf x_k$ is the state at the left end $t_k$. The preceding interval ends at the same time and uses the same variable, and this sharing keeps the trajectory continuous. The stage states $\mathbf x_{k,1}$ and $\mathbf x_{k,2}$ are the states at the two nodes. The torque $u_k$ is held constant over the interval, a choice called a zero-order hold and the simplest control representation. In the chapter’s notation, the control matrix $B$ is then a column of ones. A linear or quadratic torque would add numbers per interval without changing the construction of the state.

All states and torques are unknowns at the same time, and no simulation runs inside the optimizer. This is why direct collocation is called a simultaneous method. In single shooting, by contrast, the optimizer chooses only the torques, and a time integrator computes the states.

Python counts from zero, so the stage state $\mathbf x_{k,j}$ is stored in x_stage[k, j-1]. The solver needs all unknowns in one flat vector $\mathbf z$, and the function unpack cuts that vector into the three arrays.

For the exam. Count the unknowns of a collocation transcription from the number of intervals $N$, the number of nodes $s$, the state dimension $n$, and the control representation. Here the count is $(N+1)n+Nsn+N$, which is $41\cdot2+40\cdot2\cdot2+40=282$: the mesh states, the stage states, and one torque per interval. Step 3 counts the equations in the same way.

Step 3: Assemble Each Interval¶

The slopes $\mathbf f_{k,j}=\mathbf f(\mathbf x_{k,j},u_k)$ are now computed from the unknowns, so they change whenever the optimizer changes $\mathbf z$. The rule $A$, $b$ does not change. On each interval, the unknowns must satisfy the chapter’s stage equations and endpoint equation:

Question. Stop here before reading on. Where do these two equations come from? Derive them from the results of Step 1.

The stage equation for node $i$ requires the stored stage state $\mathbf x_{k,i}$ to lie on the polynomial built from the slopes. It is needed because the slopes were evaluated at the stored stage states. If those states could differ from the polynomial’s values, the slopes would describe some other trajectory. When the stage equation holds, the polynomial has the value $\mathbf x_{k,i}$ at node $i$, and its slope there is $\mathbf f_{k,i}$ by construction of the interpolant. The polynomial therefore satisfies the differential equation at node $i$. This is where the name comes from: to collocate is to place together, from the Latin com- (together) and locare (to place). A collocation method places the approximation and the differential equation together at the chosen nodes, where they must agree exactly. Between the nodes, nothing forces them to agree.

The endpoint equation requires the next mesh state to equal the value of the polynomial at the right end of the interval. The difference between its two sides is called the defect. Because $\mathbf x_{k+1}$ is also the starting state of interval $k+1$, a zero defect joins consecutive pieces into one continuous trajectory.

With the numbers from Step 1, the three equations read

The first two use the rows of $A$ and the third uses $b$. The second and third equations have the same right side because $\tau_2=1$. The second row of $A$ integrates each basis function from 0 to $\tau_2=1$, which is the same integral that defines $b$, so $A_{2j}=b_j$: both rows are $(\tfrac{3}{4},\tfrac{1}{4})$. Two quantities equal to the same expression are equal to each other, so the equations imply $\mathbf x_{k,2}=\mathbf x_{k+1}$. This matches the picture: node 2 sits at the right end of the interval, so the state there is the next mesh state. The code keeps both variables and both equations anyway, so that it stays valid for nodes that do not include the right end.

The objective is built the same way. The cost accumulated since time 0 behaves like an extra state whose slope is the cost rate $c(\mathbf x,u)$, and the same rule integrates it. Its stage equations would use the rows of $A$ and give the cost accumulated up to each node. The objective needs only the total over each interval, which is the value at the right end, so only the weights $b$ appear. Interval $k$ adds $h(\tfrac{3}{4}c_{k,1}+\tfrac{1}{4}c_{k,2})$ to the objective, where $c_{k,j}=c(\mathbf x_{k,j},u_k)$. Here $c=u_k^2$ is constant on the interval and the weights sum to one, so this sum equals $hu_k^2$, the exact integral of the squared torque over the interval. The boundary conditions $\mathbf x_0=(0,0)$ and $\mathbf x_N=(\pi,0)$ and the torque limits $-u_{\max}\le u_k\le u_{\max}$ complete the finite problem. Because the torque is constant on each interval, bounding each $u_k$ bounds the torque at every instant, not only at the nodes.

In equalities, the array slopes has one row per stage and one column per state component. The product A @ slopes computes all the sums $\sum_jA_{ij}\mathbf f_{k,j}$ at once, one row per node $i$, and b @ slopes computes the endpoint sum. Each residual is written as left side minus right side, so a solution makes it zero. The last line flattens the residuals into the single vector $H(\mathbf z)$ that the solver expects.

The count of equations explains what the optimizer is free to choose. The 240 stage and endpoint equations match the 240 states other than $\mathbf x_0$, so for a given initial state and given torques they determine the whole trajectory, as a simulation would. The torque is a zero-order hold with one value on each of the 40 intervals, so the optimizer chooses 40 torques. The four boundary equations fix $\mathbf x_0$ and impose two terminal conditions, which remove two of those choices and leave $40-2=38$ degrees of freedom. The torque limits are inequalities, not equations: they do not change this count, but they rule out part of that family. Among the torque sequences that bring the pendulum to the goal within the limits, the objective selects the one with the least effort.

The table collects the correspondence between the chapter’s symbols and the code.

Step 4: Connect the Controls¶

The state must be continuous, because it is obtained by integrating its rate of change, but the control need not be. A motor command can jump from one value to the next, as it does here between intervals, so no equation links $u_k$ to $u_{k+1}$. State continuity needs no extra equation either, because consecutive intervals share the mesh state $\mathbf x_{k+1}$. A smooth torque could be obtained by making it linear on each interval and sharing its endpoint values between intervals.

Step 5: Solve and Check¶

The solver needs a starting point. The initial guess moves the states along a straight line from the start to the goal and sets the torque to zero. This is not a motion of the pendulum: without torque the pendulum would stay hanging. The largest residual at the guess, printed below, measures how far the guess is from satisfying the equations. A simultaneous method can start from such an infeasible guess, which only has to describe roughly where the motion should go. Such a guess is often easier to provide than a good torque sequence.

SciPy’s SLSQP method then improves $\mathbf z$ step by step. At each iteration it approximates $F$ by a quadratic function and $H$ by a linear one, solves that approximation for a step, and repeats. Intermediate iterates need not satisfy $H(\mathbf z)=\mathbf 0$; at convergence the residuals are zero up to the solver’s tolerance. SLSQP estimates the derivatives of objective and equalities by finite differences, which is adequate for 282 unknowns. Larger problems supply exact derivatives, for example by automatic differentiation, and use solvers that exploit the sparsity of the equations.

The torque limits reach the solver as bounds. SLSQP accepts one (lower, upper) pair per entry of $\mathbf z$, with None for no bound: the states are free, and each of the last $N$ entries, the torques, must lie in $[-u_{\max},u_{\max}]$. The option ftol=1e-10 asks for a tighter tolerance than SLSQP’s default, so that the two formulations compared at the end of the demo agree closely.

For the exam. A small residual does not mean that the trajectory is accurate. Explain why, and say how to check accuracy.

The residual measures one thing: how far the stored numbers are from satisfying the collocation equations. The solver drives it down to its tolerance. Those equations are not the differential equation. They require the ODE only at the nodes, and they replace the integrand by its polynomial interpolant, so even an exact solution of the equations differs from the true motion under the same torques. This discretization error depends on the interval width and the nodes, not on the solver, so tightening the solver’s tolerance does not reduce it. To measure it, the next cell applies the optimized torques to the pendulum with an accurate ODE solver and compares the result with the stored states. This replay is one of the checks described in the chapter’s section on checking the continuous trajectory.

The difference from the replay, about $4\times10^{-3}$, is the discretization error of 40 intervals with the Radau rule. The second experiment at the end shows how quickly it decreases as the number of intervals grows.

In the figure, the grey curves are the replayed pendulum, the dots are the stored states, and the dashed lines mark the torque limits. The dots lie on the curves, and the pendulum arrives upright and at rest. The torque shows how a weak motor swings a pendulum up. For the first two seconds it pushes backward at the limit -0.5, and the pendulum swings back to about $-54^\circ$, which it reaches at $t\approx2.8$. The torque reverses just before the pendulum turns around, and from $t=3$ to $t\approx5.6$ it pushes forward at the limit +0.5 while the pendulum swings through the bottom and climbs. Near the top the motor eases off and then brakes, so that the pendulum arrives at rest instead of swinging past the upright position.

The animation below shows the same replay as a motion: the optimized torques drive the pendulum from hanging at rest to upright at rest. It is computed by the accurate ODE solver, not read from the stored values, so it shows what the torques actually do to the pendulum. Press play, or drag the slider to step through time.

Question. Why does the pendulum swing backward first, instead of going straight up? Compare the torque limit with the largest torque that gravity can exert, and think about the pendulum’s energy.

The Same Construction with Derivatives¶

In the integration form above, the slopes are matched by design: the interpolated slope takes the value $\mathbf f_{k,j}$ at each node. What must be enforced is the consistency of the states. Integrating the interpolated slope from the mesh state $\mathbf x_k$ to node $i$ reconstructs a state value there, and the stage equation requires that reconstructed value to equal the stage state $\mathbf x_{k,i}$ that the optimizer holds as a decision variable. The chapter’s differentiation form reverses the roles. It builds the polynomial through the stored state values, so the states are matched by design, and it enforces the slopes: the derivative of that polynomial must equal the ODE slope at each node. At a solution, both forms express the same two facts: the polynomial passes through the stored states, and its slope equals $\mathbf f$ at each node. They differ in which fact holds by construction and which becomes a constraint.

With $s$ nodes, the interpolated slope has degree $s-1$, so the state polynomial, its integral, has degree $s$, and $s+1$ values determine it. In this demo $s=2$: the slope is a line, the state is a quadratic, and three values determine it. Store its values at $\tau=0$ and at the two collocation nodes. These are $\mathbf x_k$, $\mathbf x_{k,1}$, and $\mathbf x_{k,2}$, the same unknowns as before, so $\mathbf z$, unpack, and objective do not change.

Let $\ell_0,\ell_1,\ell_2$ now denote the Lagrange basis for the three points $(0,\tfrac{1}{3},1)$, and write $\mathbf x_{k,0}=\mathbf x_k$ for the stored value at $\tau=0$. The quadratic through the three stored values is their Lagrange interpolant:

The stored values are constants, so differentiating this sum differentiates only the basis functions, by the same linearity as in Step 1. At node $i$,

The numbers $D_{ir}$ depend only on the points, like $A$ and $b$, and they turn the three stored values into the exact slope of the quadratic at node $i$. That slope is measured per unit of local time. Since $t=t_k+h\tau$, the chain rule gives the physical-time slope $\tfrac{1}{h}\mathbf p_k'(\tau)$. Requiring it to equal the ODE slope at node $i$, and multiplying by $h$, gives the slope equations

By hand, the basis functions are

with derivatives $6\tau-4$, $-9\tau+\tfrac{9}{2}$, and $3\tau-\tfrac{1}{2}$. Evaluating these derivatives at $\tfrac{1}{3}$ and at 1 gives one row of $D$ per collocation node, and hence the two slope equations:

Each row of $D$ sums to zero, because a constant state has zero slope. The endpoint equation evaluates the quadratic at $\tau=1$ with the weights $e_r=\ell_r(1)$. Here $e=(0,0,1)$ because $\tau=1$ is one of the stored points, so the equation reduces to $\mathbf x_{k+1}=\mathbf x_{k,2}$. Compared with the integration form, only the rule and two lines of the constraint function change.

Both programs have the same unknowns and the same objective, and their constraints vanish at the same points, so the solver returns the same trajectory up to its tolerance. Away from a solution the two residual vectors differ, which can change the solver’s iterations but not the set of solutions.

Experiments¶

Each experiment changes a few lines. Run all cells again after a change.

The first experiment replaces the nodes in Step 1. The same code then produces the other schemes of the chapter, because only the numbers in $A$ and $b$ change. The table lists the largest difference from the replay obtained with $N=40$.

Every row that converges satisfies its own collocation equations to the solver’s tolerance. The rows differ in how well those equations represent the differential equation. Explicit Euler is the extreme case: the solver satisfies its equations, and its cost is the lowest of all rows, yet its torques drive the real pendulum nowhere near the top. Implicit Euler stops at SLSQP’s iteration limit with a residual near 0.04. A likely reason is that implicit Euler damps oscillations numerically, removing energy at every step, so that its discrete pendulum cannot reach the top with the weak motor. A node at 0 runs through Step 5 and stops at the first line of the last section, which stores the state at 0 and at the collocation nodes and needs these to be distinct. For the Gauss nodes, show prints rounded fractions because the exact entries are irrational.

The second experiment refines the mesh. With the Radau nodes, set N to 10, 20, 40, and 80; the last takes about half a minute to solve. The difference from the replay falls by a factor of about eight each time N doubles.

The third experiment makes the motor strong. Set T = 4.0, N = 20, and u_max = 1.5, which is above the largest gravity torque. The solution then uses torques up to about 1.35 and never reaches the limit, so the bound is inactive. The motor overpowers gravity at every angle, and the animation in Step 5 shows the pendulum lifted straight up, without a back swing. The collocation construction is the same; only the optimal motion changes.

Summary¶

The demo turned a continuous-time problem into 282 numbers, 244 equations, and limits on 40 of the numbers. The rule $A$, $b$ depends only on the nodes and is computed once by exact polynomial calculus. The stage equations make the polynomial on each interval satisfy the ODE at the nodes, the endpoint equations join the pieces, and the same weights integrate the cost. The solver’s residual shows that the finite equations hold, and the replay measures how well they represent the pendulum. Changing the nodes, the number of intervals, or the model changes a line or a cell, and the construction stays the same. The chapter Continuous-Time Transcription and Collocation develops each step in general, including controls that vary within an interval and constraints checked between nodes.