Recovering Rational Motions

This tutorial describes how to recover a rational motion curve from Plücker coordinates (screw axes) of a closed 4R linkage.

The prerequisite is to compute the Plücker coordinates, as described in Obtaining Rational Axes.

The input lines of the four-bar mechanism below \(h_1, h_2, h_3, h_4\) come from Hegedüs et al.[1].

In the same paper the algebraic conditions for closed loop linkages are presented. Mainly, the mechanism is described by linear factors and the condition that the chain of links closes, called the closure condition, is given by

\[(t_1 - h_1) (t_2 - h_2) (t_3 - h_3) (t_4 - h_4) \in \mathbb{R} \setminus \{0\}\]

The set of values for \(t_1, t_2, t_3, t_4\) that satisfy the above relation is called the configuration set because it describes all possible configurations of the mechanism. A mechanism is rational if and only if there exist rational transformations such that all four parameters can be expressed as rational functions of a single parameter \(t\).

One starts by computing a Gröbner basis of this set of equations and solving for \(t_2, t_3, t_4\) in terms of \(t_1\). In the example, this yields the following transformations

\[t_1 \mapsto t, \quad t_2 \mapsto t + 1, \quad t_3 \mapsto t, \quad t_4 \mapsto -t - 1.\]

Now, one can pick a base and a moving frame which usually means picking a 2R subchain of the mechanism, e.g., \(h_1, h_2\). The resulting motion polynomial in the example below comes out to be

\[C(t) = (t - h_1)(t + 1 - h_2) = t (t + 1) - (t + 1) \, \mathbf{i} - t \, \mathbf{j} + \mathbf{k} + \mathbf{e} \, (-9 - 9 t \, \mathbf{i} + 9 \, \mathbf{j} + 9 t \, \mathbf{k}).\]

In theory, the same calculation could be done for 6R linkages, but these are not always rational and the calculation might fail since there are no rational transformations for all \(t_i\).

import sympy as sp
from rational_linkages import RationalCurve, Plotter, RationalMechanism, DualQuaternion


# Input lines (h1,h2,h3,h4) for the 4R mechanism
# constructor as_rational accepts tuple as rational numbers (-1,3) = -1/3
h1 = DualQuaternion.as_rational([0, 1, 0, 0, 0, 0, 0, 0])
h2 = DualQuaternion.as_rational([0, 0, 1, 0, 0, 9, 0, -9])
h3 = DualQuaternion.as_rational([0, (-1,3), (-2,3), (2,3), 0, -4, 4, 2])
h4 = DualQuaternion.as_rational([0, (2,3), (1,3), (2,3), 0, 5, 4, -7])

# create symbolic components of DQs representing revolute joints
# t1, t2, t3, t4 are the rotation parameters of the revolute joints
# u is the auxiliary projective variable that ensures the first eq is real and non-zero
# t is the parameter of the motion curve to be recovered, it corresponds to t1
t1, t2, t3, t4, u, t = sp.symbols('t1 t2 t3 t4 u t')
t_dq = DualQuaternion.as_rational([t, 0, 0, 0, 0, 0, 0, 0])
t1_dq = DualQuaternion.as_rational([t1, 0, 0, 0, 0, 0, 0, 0])
t2_dq = DualQuaternion.as_rational([t2, 0, 0, 0, 0, 0, 0, 0])
t3_dq = DualQuaternion.as_rational([t3, 0, 0, 0, 0, 0, 0, 0])
t4_dq = DualQuaternion.as_rational([t4, 0, 0, 0, 0, 0, 0, 0])

# closure condition of the linkage (must equal to identity)
eqs_closure = (t1_dq - h1) * (t2_dq - h2) * (t3_dq - h3) * (t4_dq - h4)

# copy of the closure to be manipulated
eqs = (t1_dq - h1) * (t2_dq - h2) * (t3_dq - h3) * (t4_dq - h4)
# non-zero projective condition for the first equation
eqs[0] = sp.expand(eqs[0]*u - 1)

# list of equations
eqs_list = [sp.expand(el) for el in eqs.array()]

# Groebner basis of the ideal generated by the closure equations
# using graded reverse lexicographic ordering is suggested
gb = sp.groebner(eqs_list, t1, t2, u, t3, t4, order='grevlex')

# collect list of equations from the Groebner basis that do not contain u
eqs_gb = [T for T in gb if u not in T.free_symbols]

# obtain configuration curve by solving for t2,t3,t4
# in terms of t1, i.e. the motion parameter t
config_curve = list(sp.linsolve(eqs_gb, [t2, t3, t4]))

# check if the configuration curve satisfies the closure equations (equals identity)
check_eqs = [sp.simplify(eq.subs([(t2, config_curve[0][0]),
                                (t3, config_curve[0][1]),
                                (t4, config_curve[0][2])])) for eq in eqs_closure]
# check sympy zeros except first eq of check_eqs
if not all(sp.simplify(eq.subs(t1, t)) == 0 for eq in check_eqs[1:]):
    raise ValueError("Not all equations are zero after substitution")

# express t2 as DQ in terms of t1
t2_res = DualQuaternion.as_rational([
    config_curve[0][0].subs(t1,t), 0, 0, 0, 0, 0, 0, 0
])

# construct motion curve of on branch as c = (t - h1)(t2 - h2)
c_dq = (t_dq - h1) * (t2_res - h2)
c = RationalCurve([sp.Poly(eq, t) for eq in c_dq.array()])

# factorize motion curve and create mechanism
m = RationalMechanism(c.factorize())

# plot the motion curve
p = Plotter(mechanism=m, arrows_length=0.025)
p.show()

References