Multibody and Nonlinear Dynamics

Overview

In the MSc course Multibody and Nonlinear Dynamics (4DM10, Q2 2025–2026) at TU Eindhoven, I completed two substantial take-home assignments. The first assignment covers multibody dynamics — deriving a constrained DAE model of a floating telescopic crane and simulating its forced response. The second covers nonlinear dynamics — passivity-based energy-tracking control of an oscillating mass on a rolling ring, plus a Lyapunov stability analysis of a COVID-25 SIR epidemic model with vaccination.

Assignment 1 — Floating Telescopic Crane

The first assignment models an offshore floating telescopic crane — the kind used for heavy-lift construction — as a planar constrained multibody system. The schematic on the left of the thumbnail shows the three rigid bodies (boat, rotary arm, telescopic rod with a point-mass load) together with the buoyancy spring-damper pair, the cable spring $k_1$, the rotary motor $M$, and the prismatic actuator force $F$.

System and coordinates

The system uses the full redundant coordinate set

With 9 generalized coordinates and 5 holonomic scleronomic constraints (one translational for the boat, a revolute joint coupling body 1 to body 2, and a prismatic joint coupling body 2 to body 3), the system has 4 degrees of freedom. The coordinates are therefore dependent and call for the Lagrange-multiplier formulation.

Problem 1 — Model derivation

Following the standard Lagrangian workflow, I derived:

1. Constraint equations $\underline{h}(\underline{q}) = \underline{0}$ for the three joints (2 equations from the revolute pair, 2 from the prismatic pair, 1 from the boat translation lock).
2. Kinetic energy $T(\underline{q}, \dot{\underline{q}})$ — including the cross-coupling term from the point-mass load $m_L$ rigidly fixed at the rod tip.
3. Potential energy $V(\underline{q})$ — gravitational plus the three elastic terms (buoyancy $k_w$ at each hull point, cable $k_1$, and internal telescopic spring $k_2$).
4. Non-conservative generalized forces $\underline{Q}^{\text{nc}}$ — wave forces $\vec{F}_{w_1}, \vec{F}_{w_2}$, motor torque $M$, prismatic force $F$, and the two viscous dampers $d_w, d_t, d$.
5. Equations of motion in the standard DAE form

With four independent motions and two applied inputs ($M$ and $F$), the system is underactuated — the wave-induced sway and heave of the boat cannot be driven directly, only compensated for through the rotary joint.

Problem 1.7 — Forward dynamic simulation

Given the parameter table, the piecewise inputs

plus sinusoidal wave forces $F_{w_1}(t) = 10^8 \sin^2(0.1 t + \pi/4)$ and $F_{w_2}(t) = 10^8 \sin^2(0.1 t)$, I integrated the index-3 DAE with MATLAB's ode15s (stiff solver) over $t \in [0, 300]$ s. Baumgarte stabilization (with manually tuned $\alpha$, $\beta$) added to the second-derivative form of the constraint equations kept constraint drift within $\|\underline{h}(\underline{q})\|_\infty < 10^{-6}$ throughout the full 300-second simulation — essential at these actuator magnitudes ($\sim 10^7$–$10^8$), where naive integration drifts in seconds.

Problem 2 — Simplified linear dynamics & FRF

For Problem 2 the rod and arm were rigidly fused (no telescoping), reducing the system to the minimal coordinate set $\underline{q}_{\min} = [y_1, \theta_1, \theta_2]^\top$. With $m_L = 0$ and the wave forces set to zero, the nonlinear equations of motion reduce to three coupled ODEs (equations 7a–7c of the assignment), which I linearized around the equilibrium

to obtain $\underline{M}_l, \underline{D}_l, \underline{K}_l, \underline{S}_l$. The eigenvalues of the linearized second-order system all had negative real parts, confirming local asymptotic stability of the operating point. The three undamped eigenmodes correspond to the heave $y_1$ mode (buoyancy spring), the boat-tilt $\theta_1$ mode, and the arm $\theta_2$ mode, and their eigenvectors clearly separate vertical from rotational motion at the low-frequency end of the spectrum.

The frequency response function from actuator torque $M$ to boat heave $y_1$ was then computed analytically from the linearized matrices:

with the low-frequency slope explained by compliance of the buoyancy spring and the two resonance peaks lining up with the computed eigenfrequencies.

Problem 2.3 — Inverse dynamics for wave compensation

Given measured boat motion $y_1(t) = 1.16 \sin(0.2t - \pi/4) + 1.55$ and $\theta_1(t) = 0.09 \sin(0.2t + \pi/6) - 0.06$, the task was to reconstruct the torque $M(t)$ that would generate exactly this motion on the physical system (an inverse-dynamics problem). Substituting the measured trajectory and its analytical derivatives into equation 7b directly yields $M(t)$ with no integration required — a clean demonstration of the inverse-dynamic property of a reference trajectory for a smooth multibody system.

Assignment 2 — Nonlinear Dynamics

Assignment 2 split into two independent problems, both showcasing classical nonlinear-dynamics tools.

Problem 1 — Passivity-based energy control of an oscillating mass on a ring

A uniform ring (mass $M$, inner radius $R_i$, outer radius $R_o$) rolls without slipping on a flat surface. A point mass $m$ is mounted inside the ring, connected to the center via two identical springs of stiffness $k$ and rest length $R_i$, free to slide along a body-fixed radial axis. An actuator applies a force $F$ on the point mass along that same axis, with reaction forces $-\tfrac{1}{2}F$ on each ring-contact point. The independent coordinates are

where $\theta$ is the ring-rotation angle and $\ell$ is the signed mass displacement from the ring center.

After deriving the Euler–Lagrange equations and expressing the total mechanical energy

the key observation is that differentiating $H$ along the flow yields

so the system is passive from input $F$ to output $y = \dot{\ell}$. That one line makes an energy-shaping control law available immediately.

The Lyapunov-based energy-tracking controller

The assignment asks us to design a feedback law that drives $H$ to a prescribed target energy $H_* > 0$. Using the Lyapunov candidate

the derivative along the flow is $\dot{V} = \tfrac{1}{\gamma}(H - H_*)\dot{H} = \tfrac{1}{\gamma}(H - H_*) F\dot{\ell}$. Choosing the control law

yields $\dot{V} = -(\dot{\ell})^2 (H - H_*)^2 / \gamma \le 0$, giving convergence to the invariant set $\{\dot{\ell} = 0\} \cup \{H = H_*\}$. By LaSalle's invariance principle, trajectories either converge to the desired energy level or to an equilibrium with zero mass oscillation — no linearization, no feedback linearization, just passivity + Lyapunov.

Numerical exploration (Problem 1.5)

Simulating the closed loop with $\gamma = 0.5$, $M = 4$ kg, $m = 1$ kg, $R_i = 0.8$ m, $R_o = 1$ m, $k = 20$ N/m showed three distinct asymptotic regimes depending on $H_*$:

The threshold $H_p$ is the potential barrier (energy needed to lift the point mass to the top of the ring) — a clean physical interpretation of a bifurcation in an otherwise abstract Lyapunov analysis.

Problem 2 — COVID-25: SIR dynamics, stability, and vaccination

The second problem treats a population split into vulnerable ($x_1$), contaminated ($x_2$), and immune ($x_3$) compartments of constant total $N$. After normalization to fractions $v = x_1/N$, $c = x_2/N$, $i = x_3/N$, the reduced state-space model becomes

with $i = 1 - v - c$. The set $\mathcal{I} = \{(v, c): v \ge 0, c \ge 0, v + c \le 1\}$ is forward-invariant (each boundary face is crossed inward by the vector field), so the reduced model suffices.

Equilibria and the epidemic threshold $\theta$

Introducing the dimensionless group

the equilibria are the disease-free point $(v^*, c^*) = (1, 0)$ and — when $\theta > 1$ — the endemic point $\bigl(v^*, c^*\bigr) = \bigl(\tfrac{1}{\theta}, \tfrac{\beta}{\beta + \mu}\bigl(1 - \tfrac{1}{\theta}\bigr)\bigr)$. Linearization classifies the disease-free equilibrium as a stable node for $\theta < 1$ and a saddle for $\theta > 1$. The endemic equilibrium's global asymptotic stability for $\theta > 1$ is proven using the radially unbounded Lyapunov candidate

whose time derivative reduces to $-c(\beta + \mu)(2 - v\theta - 1/(v\theta))/2 \le 0$ after the algebraic identity $(\theta v - 1)^2 \ge 0$.

Vaccination control & bifurcation

Adding a vaccination rate $u$ that sends a fraction of newborns directly from vulnerable to immune modifies the vulnerable equation to $\dot{v} = \beta - \beta v - \gamma c v - \beta u$. With a 90%-efficient vaccine, the required vaccination rate to eliminate the endemic equilibrium works out to $u \ge \tfrac{1}{0.9}\bigl(1 - \tfrac{1}{\theta}\bigr)$; with a 100%-efficient vaccine the bifurcation occurs cleanly at $u = 1 - 1/\theta$. Plotting $c^*$ vs. $u$ gives a transcritical bifurcation diagram where the endemic branch collapses into the disease-free one exactly at the predicted $u$.

Results & Reflections

These two assignments were the most mathematically demanding course work of Q2. The multibody assignment taught me how Baumgarte stabilization quietly saves a naive DAE solver from constraint drift — without it, the crane simulation diverges within seconds under $10^8$-magnitude loads. The nonlinear assignment was a reminder that passivity arguments are often shorter than the linearized analysis they replace: one line of $\dot{H} = F\dot{\ell}$ completely determines a globally stabilizing control law, where a linearized design would leave large regions of state space uncovered. The epidemiological half of Assignment 2 was a welcome change of domain — same Lyapunov toolkit, different physical interpretation.

Technologies Used

MATLAB (ode15s, symbolic toolbox, solve, vpa), Lagrangian mechanics with multipliers, Baumgarte constraint stabilization, passivity-based control, LaSalle's invariance principle, bifurcation analysis, frequency response analysis