"""Synthesising the mechanism itself, rather than the dimensions of a given one.
What this is for
----------------
§5.4 measures the EXlink topology against a slider-crank and finds it worth
17.6 %. Two topologies establish a contrast and not a trend, and §6.3 recorded
the missing third -- with the warning that choosing one by hand moves the
arbitrariness rather than removing it. This module removes it instead, by
asking the optimizer for the topology.
The method is the second family of §2.2: the spring-connected model of
:cite:t:`kim2007spring`, in the reduced *link* form of :cite:t:`tran2024slm`,
with the gear-linkage extension of :cite:t:`yim2019gearlinkage`. A design
domain is populated with candidate elements, every candidate is a spring whose
stiffness is a design variable, and the stiffnesses are penalised towards their
extremes. A spring driven stiff is a rigid connection; a spring driven soft is
no connection at all. One continuous parameterisation therefore covers every
mechanism the domain can hold, and the topology falls out of a gradient solve
rather than being assumed.
The synthesis problem
---------------------
The input is the **half-speed shaft**, turning once per cycle, and the output is
the piston height. Over one input revolution the target motion has two equal
maxima and two unequal minima -- :func:`exlink.synthesis.target_motion` builds
the one that realises ``STE = 74`` mm and ``epsilon = 16`` exactly -- so the
synthesis is asked for a mechanism that turns one circular input into two
unequal up-and-downs of a point on the cylinder axis.
That period is the whole difficulty. A mechanism driven by a shaft is periodic
in that shaft's angle, so one driven from the *crankshaft* alone cannot produce
a motion that differs between its two revolutions: extended expansion needs
something that turns once per cycle. §5.6 sharpens this into an identity --
reach the piston from one shaft only and every odd harmonic vanishes exactly,
so the asymmetry, which lives entirely in the first harmonic, is *zero* rather
than small. The piston must be reached from both shafts.
What the domain can hold
------------------------
Three kinds of candidate, each with one presence variable :math:`\\rho`:
**Bars.** A spring between two nodes.
**Bodies.** A rigid triangle on three nodes, carried as three springs sharing a
single presence. §5.6 measured why this matters: assembled from three separate
bars a triangle is unreachable by descent, because each bar alone does nothing
and the objective is flat in its geometry until all three switch on together.
One variable turns the whole three-cornered link on, which is the shape the
identity above says extended expansion needs.
**Gears.** A pair between two shafts, at a ratio from a small catalogue. Every
shaft carries its own angle; the input's is prescribed and the rest are *free
coordinates of the equilibrium*, so a shaft with no gear on it is simply loose
and one with two is over-constrained -- both visible in the strain the answer
runs at. A pair between shafts at centres :math:`c_i, c_j` meshes externally,
so its pitch radii are fixed by the centre distance and the ratio,
.. math:: r_i = \\frac{d\\,r}{1 + r}, \\qquad r_j = \\frac{d}{1 + r},
\\qquad d = \\lVert c_i - c_j \\rVert,
and rolling without slip at the pitch point costs
.. math:: E = \\tfrac12 k(\\rho)\\,
\\bigl(r_i \\vartheta_i + r_j \\vartheta_j - \\varphi\\bigr)^2 .
Writing it as slip rather than as an angular error is what keeps the mesh
geometrically consistent for free -- no separate centre-distance constraint is
needed, because the radii are read off the centres the search chose -- and keeps
the term in the same units as every other spring. Choosing the ratio from a
catalogue of one presence per (pair, ratio) needs no new machinery either: two
ratios switched on at once is over-constraint, which the strain already detects.
The model
---------
Nodes carry planar coordinates. Some are grounded, some are *pins* on shafts,
one is the piston -- abscissa on the cylinder axis, ordinate the output -- and
the rest are free. Every spring :math:`m` joining nodes :math:`i` and :math:`j`
contributes
.. math:: E_m = \\tfrac12 k_m \\bigl(\\lVert x_i - x_j \\rVert - L_m\\bigr)^2,
\\qquad k_m = k_{\\min} + (1 - k_{\\min})\\,\\rho_m^{\\,p}
with :math:`L_m` measured on the initial pose, so no design is prestressed. The
configuration at an input angle minimises the total over the free coordinates --
free node positions *and* free shaft angles -- by a ridged Newton with a
backtracking line search, warm-started from the previous angle. That is the
incremental analysis the published method uses, and the reason the branch a
mechanism runs on is followed rather than re-chosen at every angle.
The penalty exponent :math:`p` is raised on a schedule
(:data:`PENALTY_SCHEDULE`). At :math:`p = 1` the problem is nearly convex in the
stiffnesses and the search moves freely; raising it makes intermediate
stiffnesses uneconomic, so the candidates separate into present and absent.
Why the objective can afford to be only the motion
--------------------------------------------------
Nothing here prices a part, and that is a real restriction rather than an
oversight: §2.6 records that the topology-synthesis literature targets a
kinematic quantity and therefore sizes nothing. What this module produces is a
*candidate*, and a candidate is exactly what the rest of the package is for --
the range chain of §3.3 scores it afterwards. Synthesising against range
directly would need the sizing fixed point inside every equilibrium solve at
every angle, which is a study of its own (§6.3).
"""
from __future__ import annotations
from dataclasses import dataclass, field
from typing import TYPE_CHECKING
import numpy as np
from numpy.typing import NDArray
from .constants import DEFAULT_TARGETS
if TYPE_CHECKING:
from .design import Design
FloatArray = NDArray[np.float64]
IntArray = NDArray[np.int64]
STIFFNESS_FLOOR = 1.0e-6
"""Stiffness of a fully absent candidate, relative to a fully present one.
Not zero, for the reason SIMP keeps a floor: a node reachable only through
absent springs would otherwise leave the Hessian singular, and the damped
Newton would have nothing to say about where it goes.
"""
PENALTY_SCHEDULE: tuple[float, ...] = (1.0, 2.0, 3.0, 4.0)
"""Exponents :math:`p` the continuation walks through."""
GROUND = -1
"""Node kind: coordinates fixed by the design vector."""
PIN = -2
"""Node kind: a point on a shaft's crank, carried round by the shaft's angle."""
SLIDER = -3
"""Node kind: abscissa fixed, ordinate free. The piston."""
FREE = -4
"""Node kind: both coordinates free."""
DRIVEN = PIN
"""Old name for :data:`PIN`, kept so §5.6's vocabulary still reads."""
PRESENCE_THRESHOLD = 0.5
"""Presence above which a candidate is kept when the answer is made discrete."""
[docs]
@dataclass(frozen=True)
class Shaft:
"""A grounded axis with an angle, which may be prescribed or solved for."""
name: str
centre_slot: int | None
"""Design-vector offset of its centre, or ``None`` when pinned at the origin."""
ratio: float | None
"""Its speed as a multiple of the input's, or ``None`` when the angle is free.
``1`` is the input shaft itself. ``-2`` is a shaft geared to it two to one
and turning the other way, which is the relation :mod:`exlink.kinematics`
gives the EXlink -- prescribed here rather than discovered, which is what
§5.6's domain does and §5.7's does not.
"""
angle_slot: int | None = None
"""Design-vector offset of its starting angle, for a free shaft."""
@property
def is_free(self) -> bool:
"""Whether the equilibrium solves for this shaft's angle."""
return self.ratio is None
[docs]
@dataclass(frozen=True)
class Pin:
"""A node carried round by a shaft at a radius and a phase."""
node: int
shaft: int
slot: int
"""Design-vector offset of the radius, with the phase next to it."""
[docs]
@dataclass(frozen=True)
class Element:
"""A candidate the search may switch on: one presence variable, one shape."""
kind: str
"""``"bar"`` for two nodes, ``"body"`` for a rigid triangle on three."""
nodes: tuple[int, ...]
@property
def springs(self) -> tuple[tuple[int, int], ...]:
"""The springs it is carried as: one for a bar, three for a body."""
if self.kind == "bar":
return ((self.nodes[0], self.nodes[1]),)
i, j, k = self.nodes
return ((i, j), (j, k), (k, i))
[docs]
@dataclass(frozen=True)
class Gear:
"""A candidate gear pair between two shafts, at one catalogue ratio."""
first: int
second: int
ratio: float
"""How many times faster the second shaft turns, and in the other sense."""
phase_slot: int
"""Design-vector offset of its assembly phase."""
[docs]
@dataclass(frozen=True)
class GroundStructure:
"""The candidate nodes, shafts and elements a mechanism may be built from.
Instances are pure topology -- what exists and what may connect to what.
Everything dimensional lives in the design vector, so one ground structure
serves every design drawn on it.
"""
names: tuple[str, ...]
"""Node names, in coordinate order."""
kinds: tuple[int, ...]
"""One of :data:`GROUND`, :data:`PIN`, :data:`SLIDER`, :data:`FREE`."""
shafts: tuple[Shaft, ...]
pins: tuple[Pin, ...]
elements: tuple[Element, ...]
gears: tuple[Gear, ...] = ()
@property
def n_nodes(self) -> int:
"""Nodes in the structure."""
return len(self.names)
@property
def n_members(self) -> int:
"""Candidate elements -- bars and bodies, not gears."""
return len(self.elements)
@property
def n_presences(self) -> int:
"""Presence variables in all: one per element and one per gear."""
return len(self.elements) + len(self.gears)
@property
def members(self) -> tuple[tuple[int, int], ...]:
"""The bar elements, as node pairs. §5.6's vocabulary."""
return tuple(e.nodes[:2] for e in self.elements if e.kind == "bar") # type: ignore[misc]
@property
def free_shafts(self) -> tuple[int, ...]:
"""Indices of the shafts whose angle the equilibrium solves for."""
return tuple(i for i, s in enumerate(self.shafts) if s.is_free)
[docs]
def springs(self) -> tuple[tuple[int, int], ...]:
"""Every spring in the domain, expanded from the elements."""
return tuple(s for element in self.elements for s in element.springs)
[docs]
def spring_owner(self) -> IntArray:
"""Which presence variable each spring takes its stiffness from."""
owner = [m for m, element in enumerate(self.elements) for _ in element.springs]
return np.asarray(owner, dtype=np.int64)
[docs]
def free_mask(self) -> NDArray[np.bool_]:
"""Which of the ``2 n_nodes`` coordinates the equilibrium solves for."""
mask = np.zeros(2 * self.n_nodes, dtype=bool)
for i, kind in enumerate(self.kinds):
if kind == FREE:
mask[2 * i : 2 * i + 2] = True
elif kind == SLIDER:
mask[2 * i + 1] = True
return mask
@property
def n_reduced(self) -> int:
"""Unknowns per equilibrium: free coordinates, then free shaft angles."""
return int(self.free_mask().sum()) + len(self.free_shafts)
[docs]
def ends(self) -> tuple[IntArray, IntArray]:
"""The springs' first and second node indices, as arrays."""
pairs = np.asarray(self.springs(), dtype=np.int64).reshape(-1, 2)
return pairs[:, 0], pairs[:, 1]
[docs]
def index(self, name: str) -> int:
"""Index of a node by name."""
return self.names.index(name)
[docs]
@dataclass(frozen=True)
class Layout:
"""Where each quantity sits in the design vector.
The vector holds the *initial* configuration, the gear phases and the
presences, and nothing else. Element lengths are deliberately absent: each
is the distance between its endpoints in the initial pose, so every design is
assembled without prestress and an entire class of structures that could
never be built is never visited.
"""
structure: GroundStructure
ground_slots: tuple[tuple[int, int], ...]
"""``(node, offset)`` for every grounded node whose position is free."""
free_slots: tuple[tuple[int, int], ...]
"""``(node, offset)`` for every node whose initial position is chosen."""
slider_slot: tuple[int, int]
"""``(node, offset)`` of the piston: cylinder abscissa then initial height."""
presence_offset: int
"""First index of the presences."""
driven_slots: tuple[tuple[int, int], ...] = ()
"""``(node, offset)`` of each pin's radius, with its phase next to it."""
centre_of: dict[int, int] = field(default_factory=dict)
"""Grounded node each pin turns about. §5.6's vocabulary."""
@property
def size(self) -> int:
"""Length of the design vector."""
return self.presence_offset + self.structure.n_presences
[docs]
def presences(self, x: FloatArray) -> FloatArray:
"""The presence variables of a design vector."""
return np.asarray(x[self.presence_offset :], dtype=float)
[docs]
def element_presences(self, x: FloatArray) -> FloatArray:
"""The presences of the bars and bodies alone."""
return self.presences(x)[: self.structure.n_members]
[docs]
def gear_presences(self, x: FloatArray) -> FloatArray:
"""The presences of the candidate gear pairs."""
return self.presences(x)[self.structure.n_members :]
[docs]
def shaft_centres(layout: Layout, x: FloatArray) -> FloatArray:
"""Where each shaft's axis sits."""
centres = np.zeros((len(layout.structure.shafts), 2), dtype=float)
for s, shaft in enumerate(layout.structure.shafts):
if shaft.centre_slot is not None:
centres[s] = x[shaft.centre_slot : shaft.centre_slot + 2]
return centres
def _shaft_angles(layout: Layout, x: FloatArray, theta: float, free: FloatArray) -> FloatArray:
"""Every shaft's angle: prescribed from the input, or taken from the state."""
angles = np.zeros(len(layout.structure.shafts), dtype=float)
k = 0
for s, shaft in enumerate(layout.structure.shafts):
if shaft.ratio is None:
angles[s] = free[k]
k += 1
else:
angles[s] = shaft.ratio * theta
return angles
def _base_positions(layout: Layout, x: FloatArray) -> FloatArray:
"""Grounded coordinates, and the piston's abscissa."""
pos = np.zeros((layout.structure.n_nodes, 2), dtype=float)
for node, offset in layout.ground_slots:
pos[node] = x[offset : offset + 2]
node, offset = layout.slider_slot
pos[node, 0] = float(x[offset])
return pos
def _place(
layout: Layout, x: FloatArray, theta: float, reduced: FloatArray
) -> tuple[FloatArray, FloatArray]:
"""Node coordinates at one input angle, and each pin's tangent.
The tangent is how the pin moves when its shaft turns, and it is what makes a
free shaft angle a coordinate the equilibrium can solve for rather than a
parameter it is handed.
"""
structure = layout.structure
mask = structure.free_mask()
n_coords = int(mask.sum())
centres = shaft_centres(layout, x)
angles = _shaft_angles(layout, x, theta, reduced[n_coords:])
pos = _base_positions(layout, x)
tangents = np.zeros((structure.n_nodes, 2), dtype=float)
for pin in structure.pins:
radius, phase = float(x[pin.slot]), float(x[pin.slot + 1])
carried = angles[pin.shaft] + phase
direction = np.array([np.cos(carried), np.sin(carried)])
pos[pin.node] = centres[pin.shaft] + radius * direction
tangents[pin.node] = radius * np.array([-np.sin(carried), np.cos(carried)])
pos.reshape(-1)[mask] = reduced[:n_coords]
return pos, tangents
[docs]
def initial_configuration(layout: Layout, x: FloatArray) -> FloatArray:
"""The full node layout at zero input angle, free nodes included."""
structure = layout.structure
mask = structure.free_mask()
reduced = np.zeros(structure.n_reduced, dtype=float)
seed = np.zeros(2 * structure.n_nodes, dtype=float)
for node, offset in layout.free_slots:
seed[2 * node : 2 * node + 2] = x[offset : offset + 2]
node, offset = layout.slider_slot
seed[2 * node + 1] = float(x[offset + 1])
reduced[: int(mask.sum())] = seed[mask]
for k, s in enumerate(structure.free_shafts):
slot = structure.shafts[s].angle_slot
if slot is not None:
reduced[int(mask.sum()) + k] = float(x[slot])
return _place(layout, x, 0.0, reduced)[0]
[docs]
def member_lengths(layout: Layout, x: FloatArray) -> FloatArray:
"""Rest length of every spring in the domain, from the initial pose."""
pos = initial_configuration(layout, x)
springs = layout.structure.springs()
out = np.empty(len(springs), dtype=float)
for m, (i, j) in enumerate(springs):
out[m] = float(np.hypot(*(pos[i] - pos[j])))
return out
[docs]
def stiffnesses(presence: FloatArray, penalty: float) -> FloatArray:
"""SIMP stiffness of each candidate: a floor plus its penalised presence."""
rho = np.clip(np.asarray(presence, dtype=float), 0.0, 1.0)
return STIFFNESS_FLOOR + (1.0 - STIFFNESS_FLOOR) * rho**penalty
[docs]
def gear_radii(layout: Layout, x: FloatArray) -> FloatArray:
"""Pitch radii of every candidate pair, from its centres and its ratio.
An external mesh puts the pitch point between the two axes, so the radii sum
to the centre distance and their ratio is the speed ratio. Reading them off
the centres the search chose is what makes the mesh consistent for free: no
separate centre-distance constraint is needed, because a pair drawn this way
always fits the shafts it is drawn between.
"""
centres = shaft_centres(layout, x)
out = np.zeros((len(layout.structure.gears), 2), dtype=float)
for g, gear in enumerate(layout.structure.gears):
distance = float(np.hypot(*(centres[gear.first] - centres[gear.second])))
out[g] = (distance * gear.ratio / (1.0 + gear.ratio), distance / (1.0 + gear.ratio))
return out
def _energy_terms(
pos: FloatArray,
ends: tuple[IntArray, IntArray],
k: FloatArray,
rest: FloatArray,
) -> tuple[float, FloatArray, FloatArray]:
"""Total spring energy and its first two derivatives in the coordinates.
The blocks are the textbook ones for an axial spring, assembled over the
full coordinate set with the caller restricting to the free ones. It is
written over all springs at once because it is the innermost thing in the
method -- every Newton step of every angle of every finite difference -- and
a Python loop over springs costs more than the arithmetic does.
Args:
pos: ``(n_nodes, 2)`` coordinates.
ends: The springs' first and second node indices.
k: Spring stiffnesses.
rest: Spring rest lengths.
Returns:
Energy, a ``2 n_nodes`` gradient, and the ``2n x 2n`` Hessian.
"""
i, j = ends
n = pos.shape[0]
d = pos[i] - pos[j]
r = np.sqrt(np.einsum("mi,mi->m", d, d))
r = np.maximum(r, 1.0e-12)
stretch = r - rest
total = float(0.5 * np.dot(k, stretch**2))
unit = d / r[:, None]
f = (k * stretch)[:, None] * unit
grad = np.zeros((n, 2), dtype=float)
np.add.at(grad, i, f)
np.add.at(grad, j, -f)
outer = unit[:, :, None] * unit[:, None, :]
eye = np.eye(2)[None, :, :]
block = k[:, None, None] * (outer + (stretch / r)[:, None, None] * (eye - outer))
hess = np.zeros((2 * n, 2 * n), dtype=float)
rows = np.stack([2 * i, 2 * i + 1], axis=1)
cols = np.stack([2 * j, 2 * j + 1], axis=1)
for first, second, sign in (
(rows, rows, 1.0),
(cols, cols, 1.0),
(rows, cols, -1.0),
(cols, rows, -1.0),
):
np.add.at(hess, (first[:, :, None], second[:, None, :]), sign * block)
return total, grad.reshape(-1), hess
def _jacobian(layout: Layout, tangents: FloatArray) -> FloatArray:
"""How the node coordinates move with each reduced unknown."""
structure = layout.structure
mask = structure.free_mask()
n_coords = int(mask.sum())
jac = np.zeros((2 * structure.n_nodes, structure.n_reduced), dtype=float)
jac[np.flatnonzero(mask), np.arange(n_coords)] = 1.0
for k, s in enumerate(structure.free_shafts):
for pin in structure.pins:
if pin.shaft == s:
jac[2 * pin.node : 2 * pin.node + 2, n_coords + k] = tangents[pin.node]
return jac
def _gear_terms(
layout: Layout,
x: FloatArray,
theta: float,
reduced: FloatArray,
stiffness: FloatArray,
radii: FloatArray,
) -> tuple[float, FloatArray, FloatArray, FloatArray]:
"""Slip energy of every candidate pair, in the reduced coordinates.
Returns the energy, its reduced gradient and Hessian, and the slip of each
pair as a fraction of its centre distance -- the gear's own strain, and the
reading that says whether a shaft has ended up with two pairs on it.
"""
structure = layout.structure
n = structure.n_reduced
n_coords = n - len(structure.free_shafts)
angles = _shaft_angles(layout, x, theta, reduced[n_coords:])
column = {s: n_coords + k for k, s in enumerate(structure.free_shafts)}
total = 0.0
grad = np.zeros(n, dtype=float)
hess = np.zeros((n, n), dtype=float)
slip = np.zeros(len(structure.gears), dtype=float)
for g, gear in enumerate(structure.gears):
first, second = radii[g]
residual = (
first * angles[gear.first]
+ second * angles[gear.second]
- float(x[gear.phase_slot])
)
total += 0.5 * stiffness[g] * residual**2
span = first + second
slip[g] = abs(residual) / span if span > 0.0 else 0.0
coefficients = {}
if gear.first in column:
coefficients[column[gear.first]] = first
if gear.second in column:
coefficients[column[gear.second]] = second
for a, ca in coefficients.items():
grad[a] += stiffness[g] * residual * ca
for b, cb in coefficients.items():
hess[a, b] += stiffness[g] * ca * cb
return total, grad, hess, slip
NEWTON_STEPS = 30
"""Newton iterations allowed per input angle."""
NEWTON_TOLERANCE = 1.0e-9
"""Convergence test on the reduced gradient, relative to its scale.
The gradient of a spring energy carries the stiffness and a length, and both
vary by orders of magnitude across a penalisation schedule, so an absolute
threshold would be met instantly at one end of it and never at the other.
"""
RIDGE = 1.0e-9
"""Fraction of the mean diagonal added to the Hessian before solving.
The Hessian is singular exactly when some node or shaft is held only by
candidates the search has switched off -- which, in a topology optimization, is
most of the time. The ridge leaves such a coordinate where the warm start put
it instead of sending it to infinity, and is small enough not to disturb one
that is properly held.
"""
BACKTRACK_STEPS = 20
"""Halvings allowed before a Newton step is abandoned."""
SLACK_TOLERANCE = 1.0e-3
"""Component in a null direction above which the piston counts as undetermined."""
SINGULAR_TOLERANCE = 1.0e-8
"""Singular values below this fraction of the largest count as a null direction."""
def _output_slack(hess: FloatArray, slider_row: int) -> float:
"""How much of the Hessian's null space the piston's coordinate lives in.
A spring network holds a coordinate only through the candidates the search
switched on. Where it holds none, the reduced Hessian is singular in that
direction and the equilibrium is not unique -- so the motion reported is the
one the ridge and the warm start happened to pick. Projecting the slider's
coordinate onto those directions says whether the *output* is among them,
which is the only part that matters: an unused node floating off on its own
is harmless, and a piston floating off on its own is the whole answer.
"""
if hess.size == 0:
return 1.0
_u, singular, right = np.linalg.svd(hess)
largest = float(singular[0]) if singular.size else 0.0
if largest <= 0.0:
return 1.0
null = right[singular <= SINGULAR_TOLERANCE * largest]
if null.size == 0:
return 0.0
return float(np.max(np.abs(null[:, slider_row])))
def _assemble(
layout: Layout,
x: FloatArray,
theta: float,
reduced: FloatArray,
k: FloatArray,
rest: FloatArray,
gear_stiffness: FloatArray,
radii: FloatArray,
ends: tuple[IntArray, IntArray] | None = None,
) -> tuple[float, FloatArray, FloatArray, FloatArray, FloatArray]:
"""Total energy and its reduced derivatives, springs and gears together.
The reduction is a change of variables rather than a restriction: a free
shaft angle moves several pins at once, so its column of the Jacobian
carries a tangent at each of them, and its diagonal picks up the curvature
of the circle each pin runs on.
"""
structure = layout.structure
pos, tangents = _place(layout, x, theta, reduced)
energy, grad_pos, hess_pos = _energy_terms(
pos, structure.ends() if ends is None else ends, k, rest
)
jac = _jacobian(layout, tangents)
grad = jac.T @ grad_pos
hess = jac.T @ hess_pos @ jac
n_coords = structure.n_reduced - len(structure.free_shafts)
centres = shaft_centres(layout, x)
for column, s in enumerate(structure.free_shafts):
curvature = 0.0
for pin in structure.pins:
if pin.shaft != s:
continue
radial = pos[pin.node] - centres[pin.shaft]
curvature -= float(np.dot(grad_pos[2 * pin.node : 2 * pin.node + 2], radial))
hess[n_coords + column, n_coords + column] += curvature
gear_energy, gear_grad, gear_hess, slip = _gear_terms(
layout, x, theta, reduced, gear_stiffness, radii
)
return energy + gear_energy, grad + gear_grad, hess + gear_hess, pos, slip
def _equilibrate(
layout: Layout,
x: FloatArray,
theta: float,
guess: FloatArray,
k: FloatArray,
rest: FloatArray,
gear_stiffness: FloatArray,
radii: FloatArray,
ends: tuple[IntArray, IntArray] | None = None,
) -> tuple[FloatArray, bool]:
"""Minimise the energy over the free coordinates and shaft angles.
Newton with a ridge and a backtracking line search. The ridge handles the
singular Hessian of a partly switched-off mechanism; the line search handles
the fact that the energy is not convex in the coordinates once a spring is in
compression or a shaft has turned past its pin's tangent, which is where a
bare Newton step overshoots.
Args:
layout: The design vector's layout.
x: Design vector.
theta: Input angle [rad].
guess: Reduced coordinates to start from, normally the previous angle's.
k: Spring stiffnesses.
rest: Spring rest lengths.
gear_stiffness: Stiffness of each candidate pair.
radii: Pitch radii of each candidate pair.
Returns:
The reduced coordinates and whether the gradient test was met.
"""
scale = max(1.0, float(np.max(k, initial=0.0)) * float(np.max(rest, initial=1.0)))
tol = NEWTON_TOLERANCE * scale
eye = np.eye(layout.structure.n_reduced)
u = np.array(guess, dtype=float)
for _ in range(NEWTON_STEPS):
energy, grad, hess, _pos, _slip = _assemble(
layout, x, theta, u, k, rest, gear_stiffness, radii, ends
)
if float(np.max(np.abs(grad))) < tol:
return u, True
ridge = RIDGE * max(float(np.mean(np.abs(np.diag(hess)))), 1.0e-30)
try:
step = np.linalg.solve(hess + ridge * eye, -grad)
except np.linalg.LinAlgError:
return u, False
alpha = 1.0
for _back in range(BACKTRACK_STEPS):
trial = u + alpha * step
if (
_assemble(layout, x, theta, trial, k, rest, gear_stiffness, radii, ends)[0]
<= energy
):
u = trial
break
alpha *= 0.5
else:
return u, False
final = _assemble(layout, x, theta, u, k, rest, gear_stiffness, radii, ends)[1]
return u, float(np.max(np.abs(final))) < tol
[docs]
@dataclass(frozen=True)
class Motion:
"""What a candidate mechanism does over one input revolution."""
lam: FloatArray
"""Piston height at each input angle [mm]."""
converged: bool
"""Whether every angle reached equilibrium."""
strain: float
"""Largest strain over the revolution, as a fraction of a rest length.
The honest measure of whether the springs are behaving as rigid links, and
it covers the gears too: a pair's strain is its slip at the pitch point as a
fraction of the centre distance. A mechanism the penalisation has resolved
runs at a strain of order the stiffness floor; anything larger means the
answer is a deforming structure rather than a mechanism.
"""
gear_ratios: tuple[float, ...] = ()
"""Ratio of every pair the design kept, in catalogue order."""
output_slack: float = 0.0
"""How undetermined the piston is, from 0 (fixed by the input) to 1 (free).
Strain catches an *over*-constrained answer. This catches the opposite, and
the opposite is the one that lies. If the equilibrium leaves a direction the
energy does not resist -- a shaft with no gear on it, a chain with a link
missing -- then where the mechanism goes is chosen by the ridge and the warm
start rather than by the mechanism, and the motion that comes back is an
artefact that can look like anything, including extended expansion. The
reading is the largest component the slider's coordinate has in the null
space of the reduced Hessian, over the revolution. A mechanism whose motion
the input actually determines runs at zero; anything above
:data:`SLACK_TOLERANCE` is not a mechanism.
"""
[docs]
def sweep(
layout: Layout,
x: FloatArray,
samples: int = 180,
penalty: float = 1.0,
) -> Motion:
"""Run a candidate through one revolution of the input shaft.
The angles are walked in order from zero and each equilibrium is warm-started
from the last, which is what makes the result a *mechanism's* motion: the
configuration follows one branch of the loop-closure solution rather than
being re-chosen at every angle.
Args:
layout: The design vector's layout.
x: Design vector.
samples: Input angles over ``[0, 2 pi)``.
penalty: SIMP exponent applied to the presences.
Returns:
The piston motion and its diagnostics.
"""
structure = layout.structure
element_rho = layout.element_presences(x)
gear_rho = layout.gear_presences(x)
owner = structure.spring_owner()
k = stiffnesses(element_rho, penalty)[owner]
spring_rho = element_rho[owner]
gear_stiffness = stiffnesses(gear_rho, penalty)
radii = gear_radii(layout, x)
rest = member_lengths(layout, x)
safe_rest = np.where(rest > 0.0, rest, 1.0)
first, second = structure.ends()
# Built once and threaded through, rather than rebuilt inside every Newton
# step of every angle -- which is where an enumeration over topologies would
# otherwise spend most of its time.
ends = (first, second)
solve_k, solve_rest = k, rest
start = initial_configuration(layout, x)
mask = structure.free_mask()
u = np.zeros(structure.n_reduced, dtype=float)
u[: int(mask.sum())] = start.reshape(-1)[mask]
for column, s in enumerate(structure.free_shafts):
slot = structure.shafts[s].angle_slot
if slot is not None:
u[int(mask.sum()) + column] = float(x[slot])
angles = np.linspace(0.0, 2.0 * np.pi, int(samples), endpoint=False)
lam = np.empty(angles.size, dtype=float)
slider_node, _slider_offset = layout.slider_slot
slider_row = int(np.flatnonzero(mask).tolist().index(2 * slider_node + 1))
ok = True
worst = 0.0
slack = 0.0
for step, theta in enumerate(angles):
u, converged = _equilibrate(
layout, x, float(theta), u, solve_k, solve_rest, gear_stiffness, radii, ends
)
ok = ok and converged
hess = _assemble(
layout, x, float(theta), u, solve_k, solve_rest, gear_stiffness, radii, ends
)[2]
slack = max(slack, _output_slack(hess, slider_row))
pos, _tangents = _place(layout, x, float(theta), u)
lam[step] = pos[slider_node, 1]
gap = pos[first] - pos[second]
length = np.sqrt(np.einsum("mi,mi->m", gap, gap))
strain = np.abs(length - rest) / safe_rest
worst = max(worst, float(np.max(strain * spring_rho, initial=0.0)))
slip = _gear_terms(layout, x, float(theta), u, gear_stiffness, radii)[3]
worst = max(worst, float(np.max(slip * gear_rho, initial=0.0)))
kept = tuple(
gear.ratio for g, gear in enumerate(structure.gears) if gear_rho[g] > PRESENCE_THRESHOLD
)
return Motion(lam=lam, converged=ok, strain=worst, gear_ratios=kept, output_slack=slack)
[docs]
def engine_ground_structure() -> tuple[GroundStructure, Layout]:
"""§5.6's candidate set: bars only, and the gear ratio given.
Eight nodes and twelve candidate bars. The input shaft is grounded at the
origin -- a gauge choice, since translating the whole mechanism changes
nothing -- and carries the pin ``P1``. A second grounded shaft carries
``P2`` at twice the input speed and the opposite sense, which is the relation
a 2:1 gear pair imposes and the one EXlink uses; the search is free to leave
it out by switching off every bar that touches it, but not to choose it.
``C`` is a plain grounded pivot, ``F1`` and ``F2`` are free nodes, and ``S``
is the piston, held on the cylinder axis.
Bars between two prescribed nodes are excluded: their length is a function
of the input angle alone, so such a bar is a constraint on the design rather
than a degree of freedom for the mechanism.
Returns:
The structure and the layout of the design vector drawn on it.
"""
names = ("A", "P1", "B", "P2", "C", "F1", "F2", "S")
kinds = (GROUND, PIN, GROUND, PIN, GROUND, FREE, FREE, SLIDER)
idx = {name: i for i, name in enumerate(names)}
movable = ("P1", "P2", "C", "F1", "F2", "S")
elements = tuple(
Element("bar", (idx[a], idx[b]))
for i, a in enumerate(movable)
for b in movable[i + 1 :]
if kinds[idx[a]] in (FREE, SLIDER) or kinds[idx[b]] in (FREE, SLIDER)
)
structure = GroundStructure(
names=names,
kinds=kinds,
shafts=(
Shaft("input", centre_slot=None, ratio=1.0),
Shaft("geared", centre_slot=0, ratio=-2.0),
),
pins=(Pin(idx["P1"], 0, 4), Pin(idx["P2"], 1, 6)),
elements=elements,
)
layout = Layout(
structure=structure,
ground_slots=((idx["B"], 0), (idx["C"], 2)),
free_slots=((idx["F1"], 8), (idx["F2"], 10)),
slider_slot=(idx["S"], 12),
presence_offset=14,
driven_slots=((idx["P1"], 4), (idx["P2"], 6)),
centre_of={idx["P1"]: idx["A"], idx["P2"]: idx["B"]},
)
return structure, layout
GEAR_CATALOGUE: tuple[float, ...] = (1.0, 2.0, 3.0, 4.0)
"""Speed ratios a candidate pair may be cut for.
Whole numbers, and that is a requirement rather than a convenience. The cycle
has to *close*: after one revolution of the input shaft the mechanism must be
back where it started, or the engine's motion differs between successive cycles
and there is no four-stroke to speak of. A geared shaft turning :math:`r` times
per input revolution returns to its start only when :math:`r` is an integer, so
a three-to-two pair would give a mechanism whose period is two input
revolutions. That is a real machine and not this one.
Small on purpose otherwise. A ratio is one integer tooth pair and one module
away from being buildable, and §4.4 already owns that step; what this catalogue
has to settle is only whether the *synthesis* wants a two-to-one relation, which
is the question §5.6 could not ask because it was handed one.
"""
[docs]
def geared_ground_structure() -> tuple[GroundStructure, Layout]:
"""§5.7's candidate set: bodies as well as bars, and the gear chosen.
Three changes from :func:`engine_ground_structure`, each aimed at something
§5.6 measured.
The second shaft's angle is **free**, and four candidate pairs -- one per
catalogue ratio -- may drive it from the input. So the ratio is an outcome:
a design that keeps the two-to-one pair has rediscovered EXlink's gearing,
one that keeps another has found something else, and one that keeps none has
a loose shaft its linkage must hold.
Ten candidate **bodies** join the twelve bars: rigid triangles carried as
three springs sharing a single presence. §5.6's stall was a triangle that
could only be assembled from three separately-switched bars, each useless
alone; one variable now turns a whole three-cornered link on.
Returns:
The structure and the layout of the design vector drawn on it.
"""
names = ("A", "P1", "B", "P2", "C", "F1", "F2", "S")
kinds = (GROUND, PIN, GROUND, PIN, GROUND, FREE, FREE, SLIDER)
idx = {name: i for i, name in enumerate(names)}
movable = ("P1", "P2", "C", "F1", "F2", "S")
loose = {idx["F1"], idx["F2"], idx["S"]}
bars = tuple(
Element("bar", (idx[a], idx[b]))
for i, a in enumerate(movable)
for b in movable[i + 1 :]
if kinds[idx[a]] in (FREE, SLIDER) or kinds[idx[b]] in (FREE, SLIDER)
)
triples = tuple(
Element("body", (idx[a], idx[b], idx[c]))
for i, a in enumerate(movable)
for j, b in enumerate(movable[i + 1 :], start=i + 1)
for c in movable[j + 1 :]
if len(loose & {idx[a], idx[b], idx[c]}) >= 2
)
gears = tuple(
Gear(first=0, second=1, ratio=ratio, phase_slot=15 + g)
for g, ratio in enumerate(GEAR_CATALOGUE)
)
structure = GroundStructure(
names=names,
kinds=kinds,
shafts=(
Shaft("input", centre_slot=None, ratio=1.0),
Shaft("geared", centre_slot=0, ratio=None, angle_slot=14),
),
pins=(Pin(idx["P1"], 0, 4), Pin(idx["P2"], 1, 6)),
elements=bars + triples,
gears=gears,
)
layout = Layout(
structure=structure,
ground_slots=((idx["B"], 0), (idx["C"], 2)),
free_slots=((idx["F1"], 8), (idx["F2"], 10)),
slider_slot=(idx["S"], 12),
presence_offset=15 + len(GEAR_CATALOGUE),
driven_slots=((idx["P1"], 4), (idx["P2"], 6)),
centre_of={idx["P1"]: idx["A"], idx["P2"]: idx["B"]},
)
return structure, layout
[docs]
def default_bounds(layout: Layout) -> tuple[FloatArray, FloatArray]:
"""Box the search runs in, in millimetres and radians.
Deliberately loose. The point of a topology search is that the answer is
not known in advance, and a box tight enough to be informative would be
assuming it.
"""
n = layout.size
lower = np.empty(n, dtype=float)
upper = np.empty(n, dtype=float)
span = 260.0
for _node, offset in layout.ground_slots:
lower[offset : offset + 2] = -span
upper[offset : offset + 2] = span
for pin in layout.structure.pins:
lower[pin.slot], upper[pin.slot] = 5.0, 120.0
lower[pin.slot + 1], upper[pin.slot + 1] = -np.pi, np.pi
for _node, offset in layout.free_slots:
lower[offset : offset + 2] = -span
upper[offset : offset + 2] = span
_node, offset = layout.slider_slot
lower[offset : offset + 2] = -span
upper[offset : offset + 2] = span
for shaft in layout.structure.shafts:
if shaft.angle_slot is not None:
lower[shaft.angle_slot], upper[shaft.angle_slot] = -np.pi, np.pi
for gear in layout.structure.gears:
lower[gear.phase_slot], upper[gear.phase_slot] = -span * np.pi, span * np.pi
lower[layout.presence_offset :] = 0.0
upper[layout.presence_offset :] = 1.0
return lower, upper
COUNT_WEIGHT = 0.05
"""Millimetres of motion error a whole extra candidate is worth.
A tie-break and nothing more. It began at 1.0 mm, which was a mistake worth
recording: a design with no elements at all leaves the piston free, the piston
then stays where the warm start put it, and a *constant* motion scores the
target's own standard deviation -- 23.7 mm here. At 1.0 mm the empty structure
was therefore cheaper than any partly-connected one, and one start duly
converged to it, fully discrete and completely useless. The travel floor below
is the real repair; this weight is now small enough not to matter either way.
"""
TRAVEL_WEIGHT = 2.0
"""Millimetres charged per millimetre of stroke the piston fails to deliver.
Without it "build nothing" is a local minimum of the motion error, for the
reason above. With it a static piston is charged twice the required stroke,
which no moving mechanism can be worse than, so the search is obliged to
connect the piston to something before it can start improving the shape.
"""
STRAIN_WEIGHT = 800.0
"""Millimetres of motion error a unit of strain is worth.
This is the term that decides whether the answer is a mechanism, and it began
an order of magnitude too small. At 50 a start reached a motion error of
2.0 mm while carrying 3 % strain, and paid 1.5 mm for it -- so the search was
being *paid* to leave candidates half present, because a half-stiff member lets
an overconstrained network move by giving rather than by articulating. The gray
design fitted the target well and rounded to nothing. At 800 the same 3 % costs
25 mm and softness stops being a way to buy motion.
"""
HARMONIC_WEIGHT = 1.0
"""Weight on the *relative* error of each harmonic, against the plain distance.
The repair for the stall §5.6 measured. The target carries 32.33 mm at its
second harmonic and 9.02 mm at its first, and the whole Atkinson asymmetry is
the first; but a plain root-mean-square is dominated by the second, so a
mechanism that reproduces the second and none of the first already scores
7.19 mm against the 23.74 mm of doing nothing. That is a deep local minimum
holding an Otto engine. Charging each harmonic by its error *relative to the
target's own amplitude there* makes missing the small one cost as much,
proportionally, as missing the large one.
"""
BRIDGE_WEIGHT = 30.0
"""Millimetres charged for a *driven* pin the linkage does not carry to the piston.
The condition §5.6 *derives* rather than assumes: a piston reached from one
shaft alone is periodic in that shaft's angle, so every odd harmonic vanishes
exactly and the answer is an Otto engine whatever else it does. Extended
expansion therefore needs both pins joined to the piston, and both searches
spent their budget on designs where one of them was not -- start 2 of §5.7 built
two three-cornered bodies and hung all of them off the geared shaft.
Imposing a proven necessary condition is not assuming the answer: nothing here
says *how* the two chains meet, how many elements it takes, or what shape the
body that joins them has.
The charge is weighted by how far each pin's shaft is *driven*
(:func:`driven_fraction`), which is what keeps it honest in both directions. A
pin on a shaft no gear turns is not a second source, so reaching it earns
nothing -- and not reaching it costs nothing either, because a topology with no
gear at all is a legitimate single-input linkage rather than a broken two-input
one. The enumeration of :func:`enumerate_mechanisms` contains both families.
"""
SLACK_WEIGHT = 200.0
"""Millimetres charged for a piston the equilibrium does not determine.
The other half of §5.7's finding. Strain charges an over-constrained answer;
without this, nothing charges an under-constrained one, and an under-constrained
one can report any motion at all -- start 3 of §5.7 reported 70.9 mm of
asymmetry from a mechanism whose geared shaft nothing was driving.
"""
HARMONIC_ORDERS: tuple[int, ...] = (1, 2)
"""Harmonics the relative term is taken over: the asymmetry, and the strokes."""
HARMONIC_FLOOR = 1.0e-6
"""Amplitude, relative to the target's own scale, below which an order is skipped.
A relative error needs something to be relative *to*. A target with no content
at an order still has a coefficient of order 1e-15 there rather than exactly
zero, and dividing by that turns a bounded objective into 1e16 -- which is what
happened the first time this was run against a pure second-harmonic target.
"""
DISCRETENESS_SCHEDULE: tuple[float, ...] = (0.0, 2.0, 6.0, 18.0)
r"""Millimetres charged for a fully undecided candidate, one per penalty rung.
The SIMP exponent alone did not separate them here. Ramping an explicit
:math:`4\rho(1-\rho)` charge alongside it is the standard repair, and the ramp
matters: charged from the start it decides the topology before the geometry is
any good, and never charged at all it leaves an answer that cannot be rounded.
"""
[docs]
@dataclass(frozen=True)
class Candidate:
"""A mechanism the synthesis produced, and how well it does."""
x: FloatArray
"""Design vector."""
rms: float
"""Motion error against the target, mean removed [mm]."""
strain: float
"""Largest strain in a present candidate (:attr:`Motion.strain`)."""
discreteness: float
"""Mean of :math:`4\\rho(1-\\rho)`: zero when every candidate has resolved."""
present: tuple[tuple[str, ...], ...]
"""Elements that survived, by node name."""
gear_ratios: tuple[float, ...]
"""Ratios of the gear pairs that survived."""
motion: FloatArray
"""The piston motion it realises [mm]."""
converged: bool
"""Whether every equilibrium solve met its gradient test."""
@property
def uses_second_shaft(self) -> bool:
"""Whether any surviving element touches the geared shaft's pin."""
return any("P2" in nodes for nodes in self.present)
def _centred(values: FloatArray) -> FloatArray:
"""Values with their mean removed, so an offset is not an error."""
arr = np.asarray(values, dtype=float)
return arr - float(arr.mean())
def _resample(values: FloatArray, size: int) -> FloatArray:
"""Put a motion on a given number of uniformly spaced angles."""
arr = np.asarray(values, dtype=float)
if arr.size == size:
return arr
grid = np.linspace(0.0, 1.0, size, endpoint=False)
return np.interp(grid, np.linspace(0.0, 1.0, arr.size, endpoint=False), arr, period=1.0)
[docs]
def best_datum(target: FloatArray, lam: FloatArray) -> int:
"""The rotation of the input datum that lines a motion up with the target.
Where :math:`\\theta_1 = 0` sits is a choice, not a property of a mechanism:
turn the whole design about the input axis and every crank phase moves with
it, giving the same engine with its cycle starting somewhere else. A
comparison that fixes the datum therefore charges a design for something it
was never asked to control, and it charges a lot -- the real EXlink sits
170 degrees from the target's datum and was scored at 13.97 mm as posed
against 4.05 mm once aligned, which put it *behind* a degenerate answer that
had tuned its phase.
Minimising the squared difference over the shift is the same as maximising
the circular cross-correlation, so one transform pair settles it.
"""
a, b = _centred(target), _centred(lam)
spectrum = np.fft.rfft(a) * np.conj(np.fft.rfft(b))
return int(np.argmax(np.fft.irfft(spectrum, a.size)))
[docs]
def align(target: FloatArray, lam: FloatArray) -> FloatArray:
"""A motion resampled onto the target's angles and rolled onto its datum."""
resampled = _resample(_centred(lam), np.asarray(target).size)
return np.roll(resampled, best_datum(target, resampled))
[docs]
def motion_error(target: FloatArray, lam: FloatArray) -> float:
"""Root-mean-square distance between two motions, offset and datum removed.
The target's mean height is arbitrary -- it is set by where the cylinder is
bolted, not by the mechanism -- and so is the crank angle its cycle is drawn
from, for the reason :func:`best_datum` gives. Both are removed before
anything is charged.
"""
a = _centred(target)
return float(np.sqrt(np.mean((a - align(a, lam)) ** 2)))
[docs]
def harmonic(lam: FloatArray, order: int) -> tuple[float, float]:
"""The cosine and sine coefficients of one harmonic of a motion.
The diagnostic that matters for extended expansion. A motion's *even*
harmonics give the two up-and-downs per input revolution; its *odd* ones are
what make the two halves differ, so the whole Atkinson asymmetry lives in
the first harmonic and nowhere else. A synthesised mechanism can therefore
be read off directly: a first harmonic of zero is an Otto engine whatever
else it does.
"""
arr = _centred(lam)
theta = np.linspace(0.0, 2.0 * np.pi, arr.size, endpoint=False)
return (
float(2.0 * np.mean(arr * np.cos(order * theta))),
float(2.0 * np.mean(arr * np.sin(order * theta))),
)
[docs]
def harmonic_error(target: FloatArray, lam: FloatArray) -> float:
"""Millimetres charged for missing each harmonic, in proportion to its size.
Each order's coefficient error is divided by the target's own amplitude
there and the mean is scaled back into millimetres by the target's standard
deviation, so a design that misses the first harmonic entirely is charged
the same as one that misses the second entirely -- which a plain distance
does not do, and which is why §5.6 stalled on an Otto engine.
Taken on the *aligned* motion, for the reason :func:`best_datum` gives: a
harmonic coefficient carries a phase, and rotating the input datum by
:math:`s` turns the :math:`n`-th of them by :math:`ns`, so comparing
coefficients at a fixed datum charges a design an order-one error per
harmonic for a quantity that is not its business.
"""
reference = float(np.std(_centred(target)))
resampled = align(target, lam)
errors = []
for order in HARMONIC_ORDERS:
wanted = np.array(harmonic(target, order))
got = np.array(harmonic(resampled, order))
amplitude = float(np.hypot(*wanted))
if amplitude <= HARMONIC_FLOOR * reference:
continue
errors.append(float(np.hypot(*(got - wanted))) / amplitude)
return reference * float(np.mean(errors)) if errors else 0.0
DEGENERATE_ERROR = 1.0e3
"""What a motion scores when it is not a piston motion at all.
Above anything a real candidate reaches -- the worst of §5.7's four starts
scored 236 -- and finite, so that a failed sweep does not poison a whole run.
Not larger than it needs to be, because :func:`requirement_error` adds a real
error on top of it and the sum has to stay well conditioned.
"""
[docs]
def requirement_error(lam: FloatArray) -> float:
"""Root-mean-square requirement residual, in millimetres.
The four conditions of :func:`requirement_residuals`, put back on the scale
of the expansion stroke so the number can be read as a motion error and
compared against :func:`motion_error`.
A motion with no four-stroke phasing to speak of scores
:data:`DEGENERATE_ERROR` *plus its precision-point error*, rather than
raising or returning a constant. The constant was tried and it does not
work: descent from a random start sat at 1000022 for five iterations and
stopped, because a flat penalty is a plateau with no way off it. The
precision-point error is defined for any motion at all -- it reads four
heights and four slopes off the grid and needs no phasing -- so adding it
makes the plateau a slope that runs towards a motion with dead centres in
the right places, which is where the requirement becomes readable again.
"""
from .constants import DEFAULT_TARGETS
from .cycle import PhaseError
try:
residual = requirement_residuals(lam)
except (PhaseError, IndexError, ValueError):
residual = None
if residual is None or not np.all(np.isfinite(residual)):
return DEGENERATE_ERROR + precision_error(lam)
return DEFAULT_TARGETS.expansion_stroke * float(np.sqrt(np.mean(residual**2)))
[docs]
def motion_measure(kind: str, target: FloatArray, lam: FloatArray) -> float:
"""How close a motion is to what is wanted, in millimetres, three ways.
Args:
kind: ``"target"`` for the distance to a sampled target motion --
the measure §5.7's first screen used, and a proxy: a design can
beat the studied mechanism on it while missing the requirement.
``"requirements"`` for :func:`requirement_error`, which prescribes
the four dead centres and their heights but lets the two bottom
ones fall where they fall. ``"precision"`` for
:func:`precision_error`, which additionally fixes them a quarter
revolution apart -- stricter than EXlink itself satisfies, and so
reported alongside rather than searched on.
target: Target piston motion, used by ``"target"`` only.
lam: The motion to measure.
Returns:
Millimetres.
Raises:
ValueError: If ``kind`` is not one of the three.
"""
if kind == "target":
return motion_error(target, lam)
if kind == "requirements":
return requirement_error(lam)
if kind == "precision":
return precision_error(lam)
msg = f"unknown motion measure {kind!r}"
raise ValueError(msg)
[docs]
def travel_shortfall(target: FloatArray, lam: FloatArray) -> float:
"""Stroke the piston fails to deliver, in millimetres, never negative.
Overshoot is not charged here: a design that moves too far is wrong in a way
the motion error already sees, and charging it twice would bias the search
towards under-travel -- which is the failure this term exists to prevent.
"""
required = float(np.ptp(np.asarray(target, dtype=float)))
return max(0.0, required - float(np.ptp(np.asarray(lam, dtype=float))))
[docs]
def objective(
layout: Layout,
x: FloatArray,
target: FloatArray,
samples: int,
penalty: float,
discreteness: float = 0.0,
measure: str = "target",
) -> float:
"""What the search minimises: the motion, and four deterrents.
Args:
layout: The design vector's layout.
x: Design vector.
target: Target piston motion.
samples: Input angles per sweep.
penalty: SIMP exponent on the presences.
discreteness: Millimetres charged for a fully undecided candidate.
measure: How the motion is judged; see :func:`motion_measure`. The
default keeps the sampled-target reading the first screen used.
Returns:
The objective, in millimetres of equivalent motion error.
Note:
The harmonic term is charged only under the sampled-target measure, and
deliberately. It compares a candidate's first two harmonics against a
*particular* motion's, so under a precision-point measure it would put
the proxy straight back: on EXlink at 48 samples it is 4.09 mm against a
requirement error of 0.21, twenty to one, so whatever the four points
asked for, the search would still be fitting one nominated curve. What
happens between the precision points is either free or expressed as more
precision points -- not as a shape charge.
"""
motion = sweep(layout, x, samples=samples, penalty=penalty)
if not np.all(np.isfinite(motion.lam)):
return DEGENERATE_ERROR
rho = layout.presences(x)
driven = driven_fraction(layout, x)
unbridged = float(np.sum(driven * (1.0 - reachability(layout, x))))
shape = HARMONIC_WEIGHT * harmonic_error(target, motion.lam) if measure == "target" else 0.0
return (
motion_measure(measure, target, motion.lam)
+ shape
+ TRAVEL_WEIGHT * travel_shortfall(target, motion.lam)
+ STRAIN_WEIGHT * motion.strain
+ SLACK_WEIGHT * motion.output_slack
+ BRIDGE_WEIGHT * unbridged
+ discreteness * float(np.mean(4.0 * rho * (1.0 - rho)))
+ COUNT_WEIGHT * float(np.mean(rho))
)
START_SPAN = 110.0
"""Half-width, in millimetres, of the box a start is drawn from.
Not the same as the box the search runs in, and the difference matters. The
bounds are deliberately loose because the answer is not known in advance; a
*start* drawn that loosely is a different thing -- nodes scattered over half a
metre, elements three hundred millimetres long, and a piston that has to travel
74. The first runs drew from the full box and spent their budget dragging that
back to scale. This is a stroke and a half, which is the size of the mechanism
the specification implies without saying anything about its shape.
"""
[docs]
def random_start(layout: Layout, rng: np.random.Generator) -> FloatArray:
"""A design drawn on the specification's own scale, every candidate half present.
Starting the presences at one half rather than at random is deliberate: at
:math:`p = 1` the stiffnesses enter almost linearly, so a symmetric start
lets the first rung of the continuation decide which candidates matter
instead of inheriting a decision from the draw.
"""
lower, upper = default_bounds(layout)
span = np.minimum(upper - lower, 2.0 * START_SPAN)
middle = 0.5 * (np.clip(lower, -START_SPAN, None) + np.clip(upper, None, START_SPAN))
x = middle + (rng.random(layout.size) - 0.5) * span
for pin in layout.structure.pins:
x[pin.slot] = lower[pin.slot] + rng.random() * min(
60.0 - lower[pin.slot], upper[pin.slot] - lower[pin.slot]
)
x[pin.slot + 1] = -np.pi + rng.random() * 2.0 * np.pi
for shaft in layout.structure.shafts:
if shaft.angle_slot is not None:
x[shaft.angle_slot] = -np.pi + rng.random() * 2.0 * np.pi
for gear in layout.structure.gears:
x[gear.phase_slot] = (rng.random() - 0.5) * 2.0 * START_SPAN * np.pi
x[layout.presence_offset :] = 0.5
return np.clip(x, lower, upper)
def _minimise(
layout: Layout,
x0: FloatArray,
target: FloatArray,
samples: int,
penalty: float,
iterations: int,
frozen: NDArray[np.bool_] | None = None,
discreteness: float = 0.0,
measure: str = "target",
evaluations: int | None = None,
) -> FloatArray:
"""One rung: bound-constrained descent on :func:`objective`.
L-BFGS-B rather than SLSQP because the only constraints are the bounds, and
because a quasi-Newton method tolerates the finite-difference gradient this
objective is stuck with -- an equilibrium sweep is differentiable in
principle but the branch it follows is chosen by a warm start, so an
analytic derivative would have to differentiate the warm start too.
Frozen coordinates are *removed* rather than pinned by equal bounds. With a
finite-difference gradient the difference is the whole cost of the screen:
the geared domain has 45 variables of which 26 are presences, so a screen
that fixes the topology and pins them still spends 26 of every 46 sweeps
differentiating in directions it has forbidden itself to move. Dropping
them cannot change the answer -- the pinned bounds already held those
coordinates at their starting values -- and it makes each gradient 19 sweeps
instead of 45.
Args:
layout: The design vector's layout.
x0: Starting design.
target: Target motion.
samples: Input angles per sweep.
penalty: SIMP exponent on the presences.
iterations: L-BFGS-B iterations allowed.
frozen: Coordinates to hold at their starting values.
discreteness: Millimetres charged for a fully undecided candidate.
measure: How the motion is judged; see :func:`motion_measure`.
evaluations: Cap on objective evaluations, for a screen that has to
budget. L-BFGS-B spends far more than one evaluation per iteration
on this objective -- around five, on the line search -- so an
iteration count alone does not bound the work.
Returns:
The design reached, on all of the layout's coordinates.
"""
from scipy.optimize import minimize
lower, upper = default_bounds(layout)
start = np.clip(np.asarray(x0, dtype=float), lower, upper)
free = np.ones(layout.size, dtype=bool) if frozen is None else ~np.asarray(frozen)
def at(v: FloatArray) -> float:
whole = start.copy()
whole[free] = v
return objective(layout, whole, target, samples, penalty, discreteness, measure)
options: dict[str, float | int] = {"maxiter": int(iterations), "eps": 1.0e-4}
if evaluations is not None:
options["maxfun"] = int(evaluations)
result = minimize(
at,
start[free],
method="L-BFGS-B",
bounds=list(zip(lower[free], upper[free], strict=True)),
options=options,
)
out = start.copy()
out[free] = np.asarray(result.x, dtype=float)
return out
def _describe(layout: Layout, x: FloatArray, target: FloatArray, samples: int) -> Candidate:
"""Score a design at full penalty and name the candidates it kept."""
motion = sweep(layout, x, samples=samples, penalty=PENALTY_SCHEDULE[-1])
rho = layout.presences(x)
names = layout.structure.names
present = tuple(
tuple(names[n] for n in element.nodes)
for m, element in enumerate(layout.structure.elements)
if rho[m] > PRESENCE_THRESHOLD
)
return Candidate(
x=np.array(x, dtype=float),
rms=motion_error(target, motion.lam),
strain=motion.strain,
discreteness=float(np.mean(4.0 * rho * (1.0 - rho))),
present=present,
gear_ratios=motion.gear_ratios,
motion=motion.lam,
converged=motion.converged,
)
[docs]
def synthesis_residuals(
layout: Layout,
x: FloatArray,
samples: int = 48,
penalty: float = 4.0,
measure: str = "requirements",
) -> FloatArray:
"""The precision-point conditions as a *vector*, for a least-squares fit.
The objective of :func:`objective` aggregates these into one number, which
throws away the structure a Gauss-Newton method lives on. Precision-point
synthesis is classically a root-finding problem -- five conditions in
fourteen free coordinates, underdetermined, with a solution manifold rather
than an isolated minimum -- and that is what this exposes.
Two structural residuals ride along, so that a fit cannot buy the motion by
straining the linkage or by leaving the piston under-determined. They enter
as square roots because a least-squares fit squares whatever it is handed,
and the weights are the ones :func:`objective` charges.
Args:
layout: The design vector's layout.
x: Design vector.
samples: Input angles per sweep.
penalty: SIMP exponent on the presences.
measure: ``"requirements"`` for the five conditions the engine imposes,
``"precision"`` for the eight that additionally fix the angles.
Returns:
Seven residuals, or ten under ``"precision"``.
Raises:
ValueError: If ``measure`` is neither.
"""
from .cycle import PhaseError
if measure == "requirements":
width = 5
elif measure == "precision":
width = 8
else:
msg = f"unknown synthesis measure {measure!r}"
raise ValueError(msg)
motion = sweep(layout, x, samples=samples, penalty=penalty)
if not np.all(np.isfinite(motion.lam)):
return np.full(width + 2, np.sqrt(DEGENERATE_ERROR), dtype=float)
structural = np.array(
[
np.sqrt(STRAIN_WEIGHT * max(motion.strain, 0.0)),
np.sqrt(SLACK_WEIGHT * max(motion.output_slack, 0.0)),
],
dtype=float,
)
if measure == "precision":
wanted = precision_heights()
n = motion.lam.size
shifts = range(n)
best = min(
shifts, key=lambda k: float(np.sum(precision_residuals(motion.lam, wanted, k) ** 2))
)
return np.concatenate([precision_residuals(motion.lam, wanted, best), structural])
try:
residual = requirement_residuals(motion.lam)
except (PhaseError, IndexError, ValueError):
residual = None
if residual is None or not np.all(np.isfinite(residual)):
# Unreadable: charge the floor, but slope it towards a motion with dead
# centres in the right places, exactly as requirement_error does.
from .constants import DEFAULT_TARGETS
escape = precision_error(motion.lam) / DEFAULT_TARGETS.expansion_stroke
residual = np.full(width, np.sqrt(DEGENERATE_ERROR) + escape, dtype=float)
return np.concatenate([residual, structural])
[docs]
def fit_to_the_points(
layout: Layout,
x0: FloatArray,
samples: int = 48,
evaluations: int = 60,
measure: str = "requirements",
frozen: NDArray[np.bool_] | None = None,
) -> FloatArray:
"""Drive the precision-point residuals to zero, by Gauss-Newton.
A trust-region least-squares step rather than a quasi-Newton one on the
aggregated objective. The two cost the same per iteration -- both build
their derivative by finite differences over the free coordinates -- but this
one knows that the thing it is driving to zero is a *vector*, and that the
problem has more coordinates than conditions. §5.8's control is what
motivates the change: descent on the scalar could not recover the studied
mechanism's dimensions even when handed its topology.
Args:
layout: The design vector's layout.
x0: Starting design.
samples: Input angles per sweep.
evaluations: Cap on residual evaluations.
measure: Which conditions to drive; see :func:`synthesis_residuals`.
frozen: Coordinates to hold at their starting values.
Returns:
The design reached, on all of the layout's coordinates.
"""
from scipy.optimize import least_squares
lower, upper = default_bounds(layout)
start = np.clip(np.asarray(x0, dtype=float), lower, upper)
free = np.ones(layout.size, dtype=bool) if frozen is None else ~np.asarray(frozen)
def at(v: FloatArray) -> FloatArray:
whole = start.copy()
whole[free] = v
return synthesis_residuals(layout, whole, samples, PENALTY_SCHEDULE[-1], measure)
result = least_squares(
at,
start[free],
bounds=(lower[free], upper[free]),
method="trf",
jac="2-point",
diff_step=1.0e-4,
max_nfev=int(evaluations),
)
out = start.copy()
out[free] = np.asarray(result.x, dtype=float)
return out
[docs]
def synthesise_one(
layout: Layout,
x0: FloatArray,
target: FloatArray,
samples: int = 90,
iterations: int = 120,
schedule: tuple[float, ...] = PENALTY_SCHEDULE,
measure: str = "target",
) -> Candidate:
"""Walk one start up the penalisation schedule and make the answer discrete.
The last step is the one that turns a topology optimization into a
mechanism. Rounding the presences to zero and one gives a linkage; the
geometry is then re-optimised with that topology frozen, so what is reported
is a real mechanism's best fit rather than a blurred one's.
Args:
layout: The design vector's layout.
x0: Starting design.
target: Target motion, on any number of samples.
samples: Input angles per sweep during the search.
iterations: L-BFGS-B iterations per rung.
schedule: Penalisation exponents.
measure: How the motion is judged; see :func:`motion_measure`.
Returns:
The discrete mechanism and its score.
"""
x = np.array(x0, dtype=float)
for penalty, sharpness in zip(schedule, DISCRETENESS_SCHEDULE, strict=False):
x = _minimise(
layout,
x,
target,
samples,
penalty,
iterations,
discreteness=sharpness,
measure=measure,
)
rho = layout.presences(x)
x[layout.presence_offset :] = (rho > PRESENCE_THRESHOLD).astype(float)
frozen = np.zeros(layout.size, dtype=bool)
frozen[layout.presence_offset :] = True
x = _minimise(
layout, x, target, samples, schedule[-1], iterations, frozen=frozen, measure=measure
)
return _describe(layout, x, target, samples=len(target))
[docs]
@dataclass(frozen=True)
class Assessment:
"""What a synthesised mechanism actually delivers, against the requirements."""
rms: float
"""Motion error against the target [mm]."""
strain: float
"""Largest strain in a present candidate; a mechanism runs at the floor."""
four_phases: bool
"""Whether the motion is a four-stroke one at all."""
expansion_stroke: float
"""``STE`` realised [mm]; ``nan`` when the motion is not a four-stroke."""
compression_stroke: float
"""``STC`` realised [mm]; ``nan`` likewise."""
compression_ratio: float
"""``epsilon`` implied by ``STC`` and the clearance volume."""
first_harmonic: float
"""Amplitude of the first harmonic [mm] -- the asymmetry, and nothing else."""
second_harmonic: float
"""Amplitude of the second harmonic [mm] -- the two up-and-downs."""
present: tuple[tuple[str, ...], ...]
"""Elements that survived."""
gear_ratios: tuple[float, ...] = ()
"""Ratios of the gear pairs that survived."""
output_slack: float = 0.0
"""How undetermined the piston is (:attr:`Motion.output_slack`)."""
travel: float = 0.0
"""Peak-to-peak piston movement [mm]."""
@property
def delivers_the_stroke(self) -> bool:
"""Whether the piston moves far enough for the strokes to mean anything.
A guard the results asked for. A design whose piston wobbles by two
millimetres can still show four monotone phases and an "asymmetry" above
a millimetre, and one did -- start 3 of §5.7's re-run, whose harmonics
are 1.9 and 2.5 mm against a target's 9.0 and 32.3. Calling that
extended expansion is an artefact of testing the *shape* of a motion
without testing its *size*.
"""
return self.travel >= 0.5 * DEFAULT_TARGETS.expansion_stroke
@property
def is_mechanism(self) -> bool:
"""Whether the input actually determines where the piston goes."""
return self.output_slack < SLACK_TOLERANCE
@property
def asymmetry(self) -> float:
"""``STE - STC``: what makes the cycle extended-expansion rather than Otto."""
return self.expansion_stroke - self.compression_stroke
@property
def is_extended_expansion(self) -> bool:
"""Whether this is a mechanism, and its two strokes actually differ.
The mobility and travel tests come first on purpose. An
under-constrained answer can report any motion at all, extended
expansion included, and so can a piston that barely moves. Both
happened.
"""
return bool(
self.is_mechanism
and self.delivers_the_stroke
and self.four_phases
and self.asymmetry > 1.0
)
[docs]
def assess(
layout: Layout,
x: FloatArray,
target: FloatArray,
samples: int = 720,
) -> Assessment:
"""Score a synthesised design against what the engine needs.
Deliberately separate from the objective. The objective is a distance to a
prescribed motion, which is what a gradient method can descend; this is the
requirement, which is what a design has to meet. Keeping them apart is what
lets a run report that it minimised the one without satisfying the other.
"""
from .constants import DEFAULT_SPEC
from .cycle import PhaseError, find_phases
motion = sweep(layout, x, samples=samples, penalty=PENALTY_SCHEDULE[-1])
rho = layout.presences(x)
names = layout.structure.names
present = tuple(
tuple(names[n] for n in element.nodes)
for m, element in enumerate(layout.structure.elements)
if rho[m] > PRESENCE_THRESHOLD
)
try:
phases = find_phases(motion.lam)
ste, stc = phases.expansion_stroke, phases.compression_stroke
four = True
except PhaseError:
ste = stc = float("nan")
four = False
ratio = 1.0 + stc * DEFAULT_SPEC.piston_area / DEFAULT_SPEC.dead_volume
return Assessment(
rms=motion_error(target, motion.lam),
strain=motion.strain,
four_phases=four,
expansion_stroke=ste,
compression_stroke=stc,
compression_ratio=ratio,
first_harmonic=float(np.hypot(*harmonic(motion.lam, 1))),
second_harmonic=float(np.hypot(*harmonic(motion.lam, 2))),
present=present,
gear_ratios=motion.gear_ratios,
output_slack=motion.output_slack,
travel=float(np.ptp(motion.lam)),
)
[docs]
def exlink_in_the_domain(
design: Design | None = None,
samples: int = 720,
) -> tuple[GroundStructure, Layout, FloatArray]:
"""Write the studied mechanism into the synthesis domain, as a design vector.
The check that says whether a negative result is about the *domain* or about
the *search*, and here it says the search. EXlink is three elements and a
gear pair -- the swing rod ``P1-F1``, the trigonal link as a body
``P2-F1-F2``, the piston rod ``F2-S``, and a two-to-one pair -- so it is
expressible in :func:`geared_ground_structure` exactly, and the spring model
reproduces its analytic motion to about 2e-3 mm.
Every corner is a design variable, which is what makes this possible: the
two free nodes are the trigonal link's corners :math:`A` and :math:`E`, the
pins carry a radius and a phase, the piston carries its axis, and each
element's lengths follow from where those corners are put. Nothing about
the shape of the three-cornered body is fixed in advance.
Args:
design: The linkage to write in. Defaults to the study's result.
samples: Angles used to read its pose at zero input angle.
Returns:
The geared ground structure, its layout, and the design vector.
"""
from .kinematics import solve
from .reference import RELIABLE_DESIGN
kinematics = solve(RELIABLE_DESIGN if design is None else design, samples=samples)
origin = kinematics.R1[0]
second = kinematics.R2[0] - origin
pin_1 = kinematics.Q[0] - origin
pin_2 = kinematics.D[0] - kinematics.R2[0]
structure, layout = geared_ground_structure()
x = np.zeros(layout.size, dtype=float)
x[0:2] = second
x[2:4] = [-250.0, -250.0]
x[4], x[5] = float(np.hypot(*pin_1)), float(np.arctan2(pin_1[1], pin_1[0]))
x[6], x[7] = float(np.hypot(*pin_2)), float(np.arctan2(pin_2[1], pin_2[0]))
x[8:10] = kinematics.A[0] - origin
x[10:12] = kinematics.E[0] - origin
x[12:14] = kinematics.P[0] - origin
names = [tuple(structure.names[n] for n in e.nodes) for e in structure.elements]
rho = np.zeros(structure.n_presences, dtype=float)
for element in (("P1", "F1"), ("P2", "F1", "F2"), ("F2", "S")):
rho[names.index(element)] = 1.0
pair = next(g for g, gear in enumerate(structure.gears) if gear.ratio == 2.0)
rho[structure.n_members + pair] = 1.0
x[layout.presence_offset :] = rho
return structure, layout, x
[docs]
def reachability(layout: Layout, x: FloatArray) -> FloatArray:
"""How well each driven pin is joined to the piston through the linkage.
A bottleneck path: the value of a route is the weakest presence along it,
and the value of a pin is the best route it has. Fully present elements all
the way give 1, no route at all gives 0, and a half-built chain gives its
weakest link -- which is what makes it something a gradient can climb rather
than a yes-or-no test.
This is the path alone. Whether a pin is worth reaching is a separate
question -- :func:`driven_fraction` answers it -- and keeping the two apart
matters in both directions. Reaching a pin whose shaft nothing drives is
vacuous, and a run duly reached 1.00 on both while keeping no gear at all;
but *charging* for not reaching it would be wrong too, because a topology
with no gear is a legitimate single-input linkage rather than a broken
two-input one.
Returns:
One value per pin, in the order the structure lists its pins.
"""
structure = layout.structure
rho = layout.element_presences(x)
weight = np.zeros((structure.n_nodes, structure.n_nodes), dtype=float)
for m, element in enumerate(structure.elements):
for i, j in element.springs:
weight[i, j] = weight[j, i] = max(weight[i, j], float(rho[m]))
slider = layout.slider_slot[0]
out = np.zeros(len(structure.pins), dtype=float)
for p, pin in enumerate(structure.pins):
value = np.zeros(structure.n_nodes, dtype=float)
value[pin.node] = 1.0
for _sweep in range(structure.n_nodes):
improved = np.max(np.minimum(value[:, None], weight), axis=0)
better = np.maximum(value, improved)
if np.allclose(better, value):
break
value = better
out[p] = value[slider]
return out
[docs]
def driven_fraction(layout: Layout, x: FloatArray) -> FloatArray:
"""How far each pin's shaft is actually driven, from 0 to 1.
The input shaft is driven by definition. Any other is driven only as far as
a gear presence says, because a shaft carrying no gear is a passive grounded
pivot -- a perfectly good linkage element, and not a second *source* of
motion.
Returns:
One value per pin, in the order the structure lists its pins.
"""
structure = layout.structure
gear_rho = layout.gear_presences(x)
driven = np.ones(len(structure.shafts), dtype=float)
for shaft_index, shaft in enumerate(structure.shafts):
if not shaft.is_free:
continue
on_it = [
float(gear_rho[g])
for g, gear in enumerate(structure.gears)
if shaft_index in (gear.first, gear.second)
]
driven[shaft_index] = max(on_it) if on_it else 0.0
return np.array([driven[pin.shaft] for pin in structure.pins], dtype=float)
[docs]
@dataclass(frozen=True)
class Mechanism:
"""One discrete topology: which elements are present, and which gear."""
elements: tuple[int, ...]
"""Indices into the structure's elements."""
gear: int | None
"""Index of the gear pair fitted, or ``None``."""
[docs]
def presences(self, structure: GroundStructure) -> FloatArray:
"""The presence vector this topology stands for."""
rho = np.zeros(structure.n_presences, dtype=float)
for m in self.elements:
rho[m] = 1.0
if self.gear is not None:
rho[structure.n_members + self.gear] = 1.0
return rho
[docs]
def names(self, structure: GroundStructure) -> tuple[tuple[str, ...], ...]:
"""The elements it keeps, by node name."""
return tuple(
tuple(structure.names[n] for n in structure.elements[m].nodes)
for m in self.elements
)
def _joins_the_piston(structure: GroundStructure, elements: tuple[int, ...]) -> bool:
"""Whether both driven pins reach the piston through the chosen elements."""
adjacency: dict[int, set[int]] = {i: set() for i in range(structure.n_nodes)}
for m in elements:
for i, j in structure.elements[m].springs:
adjacency[i].add(j)
adjacency[j].add(i)
slider = next(i for i, kind in enumerate(structure.kinds) if kind == SLIDER)
for pin in structure.pins:
seen, stack = {pin.node}, [pin.node]
while stack:
for neighbour in adjacency[stack.pop()]:
if neighbour not in seen:
seen.add(neighbour)
stack.append(neighbour)
if slider not in seen:
return False
return True
[docs]
def is_determinate(
layout: Layout,
mechanism: Mechanism,
rng: np.random.Generator,
tries: int = 3,
) -> bool:
"""Whether the input angle fixes every unknown, at a general pose.
The rank of the reduced Hessian is the test, taken at random geometries so
that a topology is not rejected for being singular at one unlucky pose. A
topology failing it is under-constrained however its dimensions are chosen,
which is the property §5.7 found reports whatever motion the solver's ridge
happens to pick.
The absent elements are given a stiffness of *exactly* zero rather than the
usual floor. The floor exists so that the continuation always has a Hessian
to invert, and it leaves that Hessian nominally full rank whatever the
topology, so a rank test taken through it answers a question about the
tolerance rather than about the mechanism.
"""
structure = layout.structure
for _try in range(tries):
x = random_start(layout, rng)
x[layout.presence_offset :] = mechanism.presences(structure)
k = layout.element_presences(x)[structure.spring_owner()]
gear_k = layout.gear_presences(x)
radii = gear_radii(layout, x)
rest = member_lengths(layout, x)
reduced = np.zeros(structure.n_reduced, dtype=float)
mask = structure.free_mask()
reduced[: int(mask.sum())] = initial_configuration(layout, x).reshape(-1)[mask]
hessian = _assemble(layout, x, 0.0, reduced, k, rest, gear_k, radii)[2]
if np.linalg.matrix_rank(hessian, tol=1.0e-6 * float(np.max(np.abs(hessian)))) == (
structure.n_reduced
):
return True
return False
[docs]
def enumerate_mechanisms(
layout: Layout,
seed: int = 0,
max_elements: int = 4,
) -> tuple[Mechanism, ...]:
"""Every discrete topology in the domain that is a mechanism at all.
The globalization this problem actually wants. §5.7 leaves a search
problem: the objective ranks the answer an order of magnitude ahead of
anything six local runs found, and the runs commit to a short path in their
first rung and never leave it. The discrete part of the problem is small
enough not to need a heuristic -- it can be *enumerated*, which is the same
treatment §4.4 gives the gear lattice and the exhaustive baseline of
Appendix C.3.
Three conditions, applied in increasing cost. The elements must supply
exactly as many constraints as there are unknowns, counting a bar as one, a
body as three and a fitted gear as one. Both driven pins must reach the
piston, which is §5.6's identity as a condition on the graph. And the
reduced Hessian must have full rank at a general pose, which rejects the
topologies that are under-constrained however they are dimensioned.
Args:
layout: The design vector's layout.
seed: Draws the random poses the rank test uses.
max_elements: Largest number of elements to consider.
Returns:
The admissible topologies, smallest first.
"""
from itertools import combinations
structure = layout.structure
cost = [len(element.springs) for element in structure.elements]
unknowns = structure.n_reduced
rng = np.random.default_rng(seed)
found: list[Mechanism] = []
for size in range(1, max_elements + 1):
for elements in combinations(range(structure.n_members), size):
supplied = sum(cost[m] for m in elements)
if supplied not in (unknowns - 1, unknowns):
continue
if not _joins_the_piston(structure, elements):
continue
gears: list[int | None] = (
list(range(len(structure.gears))) if supplied == unknowns - 1 else [None]
)
for gear in gears:
candidate = Mechanism(elements=elements, gear=gear)
if is_determinate(layout, candidate, rng):
found.append(candidate)
return tuple(found)
[docs]
def screen_mechanism(
layout: Layout,
mechanism: Mechanism,
target: FloatArray,
rng: np.random.Generator,
draws: int = 4,
iterations: int = 20,
samples: int = 24,
measure: str = "target",
evaluations: int | None = None,
) -> tuple[float, FloatArray]:
"""Give one topology a cheap chance, and score it.
With the presences fixed the problem is no longer a topology optimization at
all -- it is the well-conditioned continuous fit :mod:`exlink.synthesis`
already solves reliably on a known linkage. That is the point of enumerating:
the hard, deceptive part of the search is the discrete part, and enumeration
removes it, leaving each topology a problem that local descent handles.
Args:
layout: The design vector's layout.
mechanism: The topology to try.
target: Target piston motion.
rng: Draws the starting geometries.
draws: Random geometries to try before descending from the best.
iterations: L-BFGS-B iterations allowed.
samples: Input angles per sweep -- coarse, because this is a screen.
measure: How the motion is judged; see :func:`motion_measure`.
evaluations: Total objective evaluations the topology is allowed, split
a third to the warm-up and the rest to the polish. An iteration
count does not bound the work on this objective -- L-BFGS-B spends
about five evaluations an iteration on the line search -- and a
screen over hundreds of topologies has to be bounded.
Returns:
The objective reached and the design vector reaching it.
Note:
There was a warm-up here, a third of the iterations spent on the
sampled-target measure before handing over. It existed because a
precision-point error used to be *undefined* wherever the motion had no
four-stroke phasing to read, which is most random starts, so that
measure alone offered a plateau with no gradient off it. Reading the
requirement off the turning points removed the plateau and with it the
warm-up's reason to exist, and measuring says it had become actively
harmful: on the studied mechanism's own topology, dropping it took the
screened error from 18.84 mm to 11.52 and the time from 49 seconds to
31. It was pulling the geometry towards one nominated curve before the
requirement was allowed to speak.
"""
presences = mechanism.presences(layout.structure)
frozen = np.zeros(layout.size, dtype=bool)
frozen[layout.presence_offset :] = True
best_value, best_x = np.inf, None
for _draw in range(draws):
x = random_start(layout, rng)
x[layout.presence_offset :] = presences
value = objective(layout, x, target, samples, PENALTY_SCHEDULE[-1], measure=measure)
if value < best_value:
best_value, best_x = value, x
assert best_x is not None
polished = _minimise(
layout,
best_x,
target,
samples,
PENALTY_SCHEDULE[-1],
iterations,
frozen=frozen,
measure=measure,
evaluations=None if evaluations is None else evaluations - evaluations // 3,
)
reached = objective(
layout, polished, target, samples, PENALTY_SCHEDULE[-1], measure=measure
)
return reached, polished
PRECISION_ANGLES: tuple[float, ...] = (0.0, 90.0, 180.0, 270.0)
"""Input-shaft angles, in degrees, at which the motion is prescribed exactly.
One input revolution is one cycle, so these are 0, 180, 360 and 540 degrees of
crankshaft: the first top dead centre, the deep bottom dead centre that ends
expansion, the second top dead centre, and the shallow one that ends
compression. Prescribing the motion *at points* rather than matching a sampled
curve is classical precision-point synthesis, and it is the repair for what
§5.7's screen exposed -- a distance to a target motion is a proxy, and a design
can beat the studied mechanism on the proxy while missing the requirement.
Here ``STE`` and the compression ratio are conditions rather than outcomes.
"""
[docs]
def precision_heights(
expansion_stroke: float | None = None,
compression_ratio: float | None = None,
) -> FloatArray:
"""Piston height at each precision angle, measured down from top dead centre.
The absolute height is free -- it is set by where the cylinder is bolted --
so only the differences are prescribed: both top dead centres at the same
height, one bottom dead centre a full expansion stroke below it, the other a
compression stroke below it.
"""
from .constants import DEFAULT_SPEC, DEFAULT_TARGETS
stroke = DEFAULT_TARGETS.expansion_stroke if expansion_stroke is None else expansion_stroke
ratio = (
DEFAULT_TARGETS.compression_ratio if compression_ratio is None else compression_ratio
)
compression = (ratio - 1.0) * DEFAULT_SPEC.dead_volume / DEFAULT_SPEC.piston_area
return np.array([0.0, -stroke, 0.0, -compression], dtype=float)
[docs]
def precision_residuals(
lam: FloatArray,
heights: FloatArray | None = None,
shift: int = 0,
) -> FloatArray:
"""How far a motion is from passing the precision points, as extremes.
Eight conditions, made dimensionless by the expansion stroke: four on the
height, taken relative to the first point so that the cylinder's mounting
height is not charged, and four on the slope, because each precision point
has to be a dead centre rather than merely a passing height.
Args:
lam: Piston height over one input revolution, on a grid whose length is
a multiple of four so the precision angles are grid points.
heights: Prescribed heights; :func:`precision_heights` by default.
shift: Rotation of the input datum, in grid steps -- free, for the
reason :func:`best_datum` gives.
Returns:
Four height residuals then four slope residuals.
"""
wanted = precision_heights() if heights is None else np.asarray(heights, dtype=float)
values = np.roll(np.asarray(lam, dtype=float), -shift)
n = values.size
if n % 4:
msg = f"precision angles need a grid divisible by four, got {n}"
raise ValueError(msg)
step = 2.0 * np.pi / n
marks = [(i * n) // 4 for i in range(4)]
scale = float(np.max(np.abs(wanted))) or 1.0
heights_here = values[marks] - values[marks[0]]
slopes = (np.roll(values, -1) - np.roll(values, 1))[marks] / (2.0 * step)
return np.concatenate([(heights_here - wanted) / scale, slopes / scale])
[docs]
def precision_error(lam: FloatArray, heights: FloatArray | None = None) -> float:
"""Root-mean-square precision residual in millimetres, datum profiled out.
The datum is searched rather than assumed, on the same argument as
:func:`best_datum`: which crank angle is called zero is a choice, so a
mechanism that passes all four points starting somewhere else passes them.
"""
wanted = precision_heights() if heights is None else np.asarray(heights, dtype=float)
scale = float(np.max(np.abs(wanted))) or 1.0
n = np.asarray(lam).size
best = np.inf
for shift in range(n):
residual = precision_residuals(lam, wanted, shift)
best = min(best, float(np.sqrt(np.mean(residual**2))))
return scale * best
def _peak(values: FloatArray, index: int) -> float:
"""A turning point's height, refined between samples by a parabola.
Three samples through an extremum fix a parabola, and its vertex is where
the extremum really is. Worth the four lines: on 24 angles the samples are
15 degrees apart, and reading a dead centre off the nearest one costs about
a millimetre of stroke on a 74 mm motion.
"""
n = values.size
before, here, after = values[index - 1], values[index], values[(index + 1) % n]
curvature = before - 2.0 * here + after
if curvature == 0.0:
return float(here)
offset = 0.5 * (before - after) / curvature
return float(here - 0.25 * (before - after) * offset)
[docs]
def requirement_residuals(lam: FloatArray) -> FloatArray:
"""How far a motion is from the requirements, without fixing their angles.
The looser reading of the precision points, and the one the engine actually
imposes. What a four-stroke needs is four dead centres in the right order
with the right heights, and the two top ones **half an input revolution
apart** so that a cycle is 720 degrees of crankshaft and the valve train can
be geared to it. Where the *bottom* dead centres fall between them is an
outcome, not a requirement.
That distinction is worth measuring rather than assuming: EXlink's own
turning points sit at 10.5, 101.5, 189.5 and 286 degrees of input shaft, so
a target demanding them at 0, 90, 180 and 270 is stricter than the studied
mechanism satisfies -- it scores 3.53 mm against that and 0.38 mm against
this.
The turning points are read off directly rather than through
:func:`exlink.cycle.find_phases`, and that is the difference between a
screen that ranks and one that does not. That function demands *exactly*
four monotone phases, which is right for judging a finished design and wrong
for judging a search iterate: one spurious reversal and a motion that is
plainly a four-stroke on a fine grid reads as no motion at all. It happened
to the studied mechanism itself -- screened from scratch on 24 angles it
scored 1281, the degenerate floor, while the same design read 28.7 on 240.
Here a fifth residual charges the spurious reversals instead of refusing to
look: a clean four-phase motion has a total variation of exactly twice its
two strokes, and anything above that is movement nobody asked for.
Returns:
Expansion stroke, compression stroke, top-dead-centre height difference,
top-dead-centre spacing and spurious travel, each as a fraction of the
expansion stroke.
Raises:
PhaseError: If the motion has no two peaks and two troughs to read.
"""
from .constants import DEFAULT_TARGETS
from .cycle import PhaseError
values = np.asarray(lam, dtype=float)
n = values.size
wanted = precision_heights()
stroke = DEFAULT_TARGETS.expansion_stroke
slope = np.sign(np.diff(np.concatenate([values, values[:1]])))
for index in np.flatnonzero(slope == 0.0):
slope[index] = slope[index - 1]
turning = np.flatnonzero(slope != np.roll(slope, 1))
peaks = [int(i) for i in turning if slope[i] < 0.0]
troughs = [int(i) for i in turning if slope[i] > 0.0]
if len(peaks) < 2 or len(troughs) < 2:
msg = (
f"lambda(theta_1) has {len(peaks)} peaks and {len(troughs)} troughs, "
"and needs two of each"
)
raise PhaseError(msg)
heights = {i: _peak(values, i) for i in peaks}
tops = sorted(peaks, key=lambda i: -heights[i])[:2]
first, second = sorted(tops)
top = max(heights[i] for i in tops)
arcs = (
values[first : second + 1],
np.concatenate([values[second:], values[: first + 1]]),
)
shallow, deep = sorted(top - float(np.min(arc)) for arc in arcs)
spacing = (second - first) / n
variation = float(np.sum(np.abs(np.diff(np.concatenate([values, values[:1]])))))
return np.array(
[
(deep - stroke) / stroke,
(shallow + wanted[3]) / stroke,
(heights[tops[0]] - heights[tops[1]]) / stroke,
(min(spacing, 1.0 - spacing) - 0.5) / 0.5,
(variation - 2.0 * (deep + shallow)) / (2.0 * stroke),
],
dtype=float,
)