Dynamics of Mechanical Systems

Course materials and tools for Dynamics of Mechanical Systems. Includes a Python library (dms) for symbolic and numerical analysis, along with simulation examples using MuJoCo.

Install

Clone the repository and install with:

pip install -r requirements.txt

Structure

• src/dms/ - Core library: symbolic dynamics (SymPy), curves, mechanisms, and bundled MuJoCo tools and models.
• examples/ - Standalone demos: MuJoCo simulations, joystick control, mechanisms, curves.
• tutorials/ - Step-by-step notebooks covering rotations, curves, MuJoCo, and EKF.

Quick Start

python examples/mujoco/pendulum_trajectory_control.py

N-bar mechanisms: rigid-link geometry mode

Why: the length-only model is not rigid for 3+ joint links

NBarMechanism(edges, link_lengths) (the original, length-only constructor) models every link as one scalar length and one angle -- exactly right for a binary link (2 joints), but a link with 3+ joints (a ternary+ plate, ground included) has no single "length" to give it. Its extra joint(s) are placed by _compute_joint_positions's cosmetic 0.35*l*perp offset, not by the solved equations, so that link is not actually rigid: the distance between two of its own joints can drift as the mechanism moves.

Concretely, on the Stephenson six-bar G = [(0,1),(1,2),(2,3),(3,0), (2,4),(4,5),(5,0)] with lengths = [1.96, 0.5, 1.34, 0.5, 1.03, 1.22] (node 2 is a ternary link), the ternary-coupler joint distance |E-C| (node 2's own joints) drifted 1.625 -> 1.001 -> 0.445 at input angles 40deg / 45deg / 50deg under the old model -- a >100% swing in a distance that must be physically constant. (The exact pair that drifts depends on BFS traversal order for a given topology/input link; tests/test_geometry_rigidity_property.py's documenting test pins one concrete drifting pair for this exact fixture, |J1_2 - J2_3| = 0.661, 0.657, 0.646 at the same three angles, as a permanent regression check so this limitation is never silently "fixed" inside the frozen length-only path.)

This is not fixed in length-only mode. NBarMechanism(...) is frozen bit-for-bit (old results, including whole GraphThe datasets generated under it, must keep reproducing exactly). It emits a one-time NBarCollinearWarning the first time it sees a link of degree >= 3:

collinear-link model; use NBarMechanism.from_joint_positions for rigid plates

How to migrate: from_joint_positions

Use NBarMechanism.from_joint_positions for any mechanism with a genuine rigid plate (ternary+ link, ground included). Give it every joint's (x, y) position at one assembled reference pose -- the SAME data as the interchange format's edge pos -- and it derives each link's local joint geometry internally; every joint on a link then keeps its reference-pose inter-distance to every other joint on that link exactly (1e-9 relative or better), at every solved angle:

from dms.mechanisms.nbar import NBarMechanism, GetTrajectory

mech = NBarMechanism.from_joint_positions(
    edges_or_graph=[(0, 1), (1, 2), (2, 3), (3, 0), (2, 4), (4, 5), (5, 0)],
    joint_positions={
        (0, 1): (0.0, 0.0), (0, 3): (2.0, 0.0), (0, 5): (1.0, -1.2),
        (1, 2): (0.5, 1.0), (2, 3): (2.3, 1.0),
        (2, 4): (0.9, 1.8), (4, 5): (0.6, -0.3),
    },
    input_link=1,                        # optional; always honored if given
    markers=[(2, (1.0, 1.5))],           # absolute (x, y) reference points
)

• theta_ref: mech.theta_ref is the absolute world angle of the input link's OWN reference vector (ground pivot -> its next joint) -- the same convention FK/ComputePoints/GetTrajectory already use in length-only mode. Internally the solver works with phi = theta1 - theta_ref (each non-input link's rotation FROM the reference pose), but every public angle stays absolute.
• FK(theta1) keeps its length-only signature and return convention ((solution, FKResult)); geometry mode warm-starts from the reference pose and tracks branches continuously -- there are no random restarts in geometry mode, so results are deterministic and reproducible run to run (a requirement for parity against an external reference solver).
• GetTrajectory(mech, n_points=360) marches both ways from the reference pose on the same absolute 360-point grid as length-only mode, and additionally sets mech.reachable_range (a ReachableRange(kind, theta_min, theta_max) named tuple): kind == "full-rotation" if every grid point solved, "bounded" for a rocker/limited-range mechanism (NaN outside a single contiguous arc containing theta_ref), or "none" only if the reference pose itself fails.

Setters and read-only accessors

Both modes support the same four setters -- thin mutators over the already-compiled solver; none of them recompile (the same mech._solver object stays valid across any number of calls), so callers like GraphThe's sampler (which iterates thousands of geometries) never pay a recompilation cost per sample:

• mech.set_lengths(lengths) -- length-only mode; list or {node: length} dict, exactly like the constructor. Resets the guess.
• mech.set_geometry(joint_positions, markers=None, axes=None) -- geometry mode; same joint_positions shape as from_joint_positions. Without an explicit markers argument, the mechanism's EXISTING markers' absolute (x, y) reference points are kept (only their internal offsets are recomputed against the new geometry). Without an explicit axes argument, the mechanism's EXISTING prismatic axes are kept unchanged (see "Prismatic joints and linear input" below). Resets the guess.
• mech.set_markers(markers) -- length-only: (link, (dx, dy)) link-frame offsets; geometry mode: (link, (x, y)) absolute reference-pose points.
• mech.reset_guess() -- resets the warm-start guess to zeros (the reference pose, in geometry mode).

Read-only accessors (no setters -- these mirror already-public, already-frozen internal structure, just without touching a private attribute directly): mech.input_link, mech.unknown_links (same order as FK's solution columns), mech.loops (independent loops from the cycle basis), mech.reference_positions (geometry mode only; None in length-only mode), mech.marker_points (geometry mode only). mech.mode ("lengths" or "geometry") and mech.theta_ref (None in length-only mode) are plain attributes.

Grashof pre-filter in geometry mode

mech.CheckGrashof() dispatches by mode. In geometry mode it uses dms.mechanisms.grashof.check_loop_feasibility_geometry, which applies the Grashof condition (4-link loops) or the polygon inequality (3- and 5+-link loops) using the distance between the two joints each loop actually traverses on every link -- not one scalar length per node shared across every loop, which is only correct in length-only mode. A link that appears in 2+ loops through different joint pairs (the ternary+ case this whole geometry mode exists for) is exactly where that distinction matters.

Prismatic joints and linear input

from_joint_positions also builds mechanisms with prismatic (sliding) joints, via two additional optional dict arguments -- joint_types (default "revolute", so it only needs prismatic entries) and axes (required per prismatic edge):

mech = NBarMechanism.from_joint_positions(
    edges_or_graph=[(0, 1), (1, 2), (2, 3), (0, 3)],  # slider-crank
    joint_positions={
        (0, 1): (0.0, 0.0),      # O2, ground pivot
        (1, 2): (34.64, 20.0),   # A, crank/coupler pin
        (2, 3): (145.6, 10.0),   # P, coupler/slider pin
        (0, 3): (145.6, 10.0),   # rail, same point as P at reference
    },
    input_link=1,                          # the crank (rotary input)
    joint_types={(0, 3): "prismatic"},      # the ground-to-slider edge
    axes={(0, 3): (1.0, 0.0)},              # rail direction (world, at reference)
)

• Owner and axis convention: for a prismatic edge (u, v) (u < v by construction), link u OWNS the axis -- the axis is a world direction at the reference pose (normalized internally; zero-length or non-finite raises ValueError), expressed in link u's own frame (which coincides with world coordinates at the reference pose, same as every other reference-pose convention in this mode). The slip's SIGN convention is +axis moves the non-owner link (v) away from the owner along that axis; the slip's UNITS are the same model length units as every joint position.
• Point keys: ComputePoints/plot emit the OWNER-side point under the usual G{v}/J{u}_{v} key (unchanged from a revolute edge) and, additionally, the non-owner (slider) side under S{u}_{v}. A link whose only two joints are a revolute edge and a prismatic edge that coincide at the reference pose (the classic single-point "slider block") is fully supported -- its own rigidity is degenerate (distance 0) and not something S{u}_{v}/G/J mix up.
• Rotation lock: a link joined to GROUND by a prismatic edge gets a LITERAL rotation of 0 (never a solved unknown) -- it cannot spin, only slide. Any other (non-ground) prismatic edge instead adds one internal slip unknown plus an angle-lock equation pinning the two links' rotations together (they never twist relative to each other).
• CheckGrashof() skips any loop that traverses a prismatic edge (a variable-length side, not a fixed one) -- it only evaluates loops that are entirely revolute.

Linear input (mech.input_kind, derived automatically -- never a separate constructor argument): when the input link's OWN edge to ground is prismatic, mech.input_kind == "linear" instead of "rotary". A linear input still uses FK/ComputePoints/set_geometry/etc with their existing single-scalar-argument signatures, but the argument's UNITS and MEANING change: for "rotary" it is the absolute input angle in radians (unchanged); for "linear" it is the axial slip from the reference pose, in model length units, along the input edge's own axis -- 0.0 is the reference pose itself (mech.theta_ref is always 0.0 for a linear input; it is simply unused, since "0 slip" already means "at the reference pose" by construction, unlike the rotary case's theta_ref offset):

mech = NBarMechanism.from_joint_positions(
    edges_or_graph=[(0, 1), (1, 2), (2, 3), (0, 3)],
    joint_positions={...},               # same slider-crank as above
    input_link=3,                        # drive the SLIDER, not the crank
    joint_types={(0, 3): "prismatic"},
    axes={(0, 3): (1.0, 0.0)},
)
assert mech.input_kind == "linear"
sol, info = mech.FK(12.5)  # slide the block 12.5 length units from reference
pts = mech.ComputePoints(12.5, sol)

Scope (locked, Phase 7.1): linear input support is at the FK/ComputePoints level ONLY. GetTrajectory(mech) and Animate(mech) both raise NotImplementedError (naming FK/ComputePoints) on a linear-input mechanism -- there is no linear sweep grid (a slider's reachable stroke has no periodic "full-rotation" analogue the way a crank's does) and no linear-input cross-repo curve fixture.

Read-only accessors for all of the above (mirroring the existing accessor style -- copies, {}/[]/"rotary" defaults for a revolute-only or length-only mechanism, never None): mech.joint_types ({(u, v): "prismatic"}), mech.axes ({(u, v): ndarray}, unit), mech.input_kind, mech.unknown_slip_edges (edges of the slip columns, in column order after mech.unknown_links). mech.set_geometry(..., axes=None) also accepts a new prismatic-axis dict now (None keeps the existing axes unchanged; an edge that isn't prismatic, or a zero-length/non-finite axis, raises ValueError).

Running the tests

pip install -e .[dev]
python -m pytest