BOOK2
BOOK2
BOOK2
Reiter
Optimal Path
and Trajectory
Planning for
Serial Robots
Inverse Kinematics for Redundant
Robots and Fast Solution of
Parametric Problems
Optimal Path and Trajectory Planning for
Serial Robots
Alexander Reiter
Springer Vieweg
© Springer Fachmedien Wiesbaden GmbH, part of Springer Nature 2020
This work is subject to copyright. All rights are reserved by the Publisher, whether the whole
or part of the material is concerned, specifically the rights of translation, reprinting, reuse of
illustrations, recitation, broadcasting, reproduction on microfilms or in any other physical way,
and transmission or information storage and retrieval, electronic adaptation, computer software,
or by similar or dissimilar methodology now known or hereafter developed.
The use of general descriptive names, registered names, trademarks, service marks, etc. in this
publication does not imply, even in the absence of a specific statement, that such names are
exempt from the relevant protective laws and regulations and therefore free for general use.
The publisher, the authors and the editors are safe to assume that the advice and information in
this book are believed to be true and accurate at the date of publication. Neither the publisher
nor the authors or the editors give a warranty, expressed or implied, with respect to the material
contained herein or for any errors or omissions that may have been made. The publisher remains
neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This Springer Vieweg imprint is published by the registered company Springer Fachmedien Wiesbaden
GmbH part of Springer Nature.
The registered company address is: Abraham-Lincoln-Str. 46, 65189 Wiesbaden, Germany
Acknowledgement
This work has been supported by the COMET-K2 “Center for Symbiotic
Mechatronics” of the Linz Center of Mechatronics (LCM) funded by the
Austrian federal government and the federal state of Upper Austria.
Alexander Reiter
Contents
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . XI
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XVII
Kurzfassung . . . . . . . . . . . . . . . . . . . . . . . . . . . XIX
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 1
2 Numerical Basics . . . . . . . . . . . . . . . . . . . . . . 9
2.1 Numerical Optimization . . . . . . . . . . . . . . . . . . . 10
2.1.1 Convexity . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.1.2 Constrained Optimization . . . . . . . . . . . . . . . . . 11
2.1.3 Parametric Sensitivities . . . . . . . . . . . . . . . . . . 15
2.2 B-Spline Curves . . . . . . . . . . . . . . . . . . . . . . . . 32
2.2.1 Evaluation Algorithm . . . . . . . . . . . . . . . . . . . 32
2.2.2 Derivatives w.r.t. the Path Parameter . . . . . . . . . . 35
2.2.3 Derivatives w.r.t. the Control Points . . . . . . . . . . . 37
2.2.4 Constrained Approximation . . . . . . . . . . . . . . . . 37
2.3 Optimal Control Problems and Direct Methods . . . . . . 40
2.3.1 Direct Single Shooting . . . . . . . . . . . . . . . . . . . 42
2.3.2 Direct Multiple Shooting . . . . . . . . . . . . . . . . . . 43
2.3.3 Direct Collocation . . . . . . . . . . . . . . . . . . . . . 45
6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . 155
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
Introduction
Due to the rapid progress in developing solutions for industrial applications
of robotic manipulators, today it is not sufficient to simply fulfill a given
task. Now it is required to solve a given problem in an optimal manner. To
this end, not only foundational knowledge about kinematics and dynamics
of multi-body systems and robotics sub-disciplines emerging therefrom are
required, it also necessitates expertise in numerical optimization.
Outline
In Chapter 2, numerical basics relevant for the rest of the work are presented.
Therein, elementary knowledge regarding numerical optimization is repeated
in Section 2.1. In particular, parametric optimization problems are discussed
wherein an iterative scheme to compute approximate solutions to such
problem is proposed as one contribution of this work. Furthermore, B-
spline curves are reviewed in Section 2.2 which are then used for general
optimal control problems in Section 2.3.
Optimal path and trajectory planning, cf. Chapter 4, are robotic applica-
tions of the above. Therein, the evolution of a robot’s joint positions (and
by virtue of the kinematic chain) the end-effector, is computed to fulfill a
specified task, e.g. point-to-point motion. This is accomplished using dif-
ferent strategies, a single-step strategy wherein the geometric and temporal
aspects are treated simultaneously, and a multi-stage approach where these
subtasks are solved subsequently. Furthermore, in order to allow for rapid
computation of the above, the iterative approximation scheme for optimal
solutions using parametric sensitivities is applied.
resolution, but also for the path tracking problem are obtained. For kine-
matically non-redundant robots, neither the inverse kinematics, nor the
path tracking problem are in the scope of this thesis.
Chapter 6 provides final remarks on this thesis and concludes with a brief
summary of the work presented herein.
Examples
In order to ease understanding, theoretical parts are frequently comple-
mented by numerical examples of typical applications. The robots featured
therein are of various types, planar and spatial, kinematically non-redundant
and redundant, as well as industrial manipulators and lab prototypes. The
following brief overview provides an introduction to the robots mentioned
throughout this work.
Figure 1.1: SCARA4 Lab Prototype. Figure 1.2: SCARA4 Planar Model.
Minor contributions that are a side product of the work leading up to this
thesis are:
Numerical Basics
In this chapter, the mathematical and numerical foundation of the subse-
quent chapters will be established.
It is not an aim and not in the scope of this thesis to provide comprehensive
treatment of these well-researched subjects, but only to provide sufficient
2.1.1 Convexity
Definition 2.1: Convex Set
A set S ⊆ Rn is called convex if a line connecting the points x, y ∈ S
lies entirely in S, i.e. x + k(y − x) = (1 − k)x + ky ∈ S ∀k ∈ [0, 1].
minimize
w
J(w) (2.1.1a)
s.t. gi (w) = 0, i = 1, . . . , Neq (2.1.1b)
gi (w) ≤ 0, i = Neq + 1, . . . , Neq + Nineq (2.1.1c)
The feasible set is convex if the equality constraints are linear and the
inequality constraints are concave.
Note that this not only concerns inequality constraints (2.1.1c) but also
all equality constraints (2.1.1b) as they are considered to be always active.
The above definitions rely on the equality operator = to determine the
activity state of a constraint. In practice it is difficult to correctly identify a
constraint’s state due to numerical issues and limited accuracy in computer
implementations. This topic will be later picked up when parametric
sensitivities are discussed as the correct identification of the active set is of
eminent importance.
Condition (2.1.3a) states, that the gradient at the minimizer is zero, i.e. no
local first-order change occurs. A local solution satisfies the (in)equality
constraints which is reflected in the conditions (2.1.3b) and (2.1.3c). Con-
dition (2.1.3d) requires positive sign of the Lagrange multipliers λ. The
complementarity condition (2.1.3e) implies that either the multipliers λ∗i ,
or the constraints gi (w∗ ) (or both) are zero. Multipliers that correspond to
inactive inequality constraints are zero, thus condition (2.1.3a) yields
∂ ∂ ∂
L(w∗ , λ∗ ) = J(w∗ ) + λ∗i gi (w∗ ) = 0. (2.1.4)
∂w ∂w i∈A(w∗ ) ∂w
If exactly one of λ∗i and gi (w∗ ) is zero (and not both), (2.1.3e) is called
the strict complementarity condition. Note that this stronger condition
only applies to inequality constraints. If λ∗i > 0 holds for an inequality
constraint, it is called strictly active.
2.1 Numerical Optimization 15
0.
2
∗ ∗
2 L(w , λ ) must be positive semi-
∂
This means that the Hessian matrix ∂w
∂
definite on the nullspace of the constraint Jacobian ∂w g(w). The expression
∗ ∂2 ∗ ∗ ∗
A (w ) ∂w2 L(w , λ )A(w ) is also called the reduced Hessian matrix.
Note that these conditions are sufficient but not necessary, i.e. they do not
need to be met by a strict local minimizer.
Consider, for example, one of the robots shown in the introductory Chapter 1.
In an optimal trajectory planning task, a load applied to its end-effector
(EE) could be regarded as a parameter that is subjected to uncertainties.
How does an optimal motion trajectory depend on the EE load and how
does it change for load variations?
16 2 Numerical Basics
d ∗
w∗ (ppert ) ≈ ŵ∗ (ppert ) = w∗ (pnom ) + w (ppert − pnom ) (2.1.11)
dp
d a∗
λa∗ (ppert ) ≈ λ̂a∗ (ppert ) = λa∗ (pnom ) + λ (ppert − pnom ) . (2.1.12)
dp
Δp
Henceforth, the abbreviations (·)∗nom := (·)∗ (pnom ) and (·)∗pert := (·)∗ (ppert )
will be used to denote the optimal solution in the nominal and perturbed
case, respectively. A similar notation will be used for the approximate
quantities.
∗
Note that the solution estimate ŵpert , λ̂∗pert obtained by means of the first-
order approximations (2.1.11) and (2.1.12) is not necessarily admissible w.r.t.
the (in)equality constraints (2.1.5b) and (2.1.5c) imposed onto the original
NLP problem. Therefore the question arises how large the admissible
parameter perturbations Δp may become. The limitations of (2.1.11)
and (2.1.12) are given by the change of the active set of constraints due
to Δp. Note, that the vector of the estimate of optimal multipliers in the
perturbed case λ̂∗pert is found by calculating the active multiplier estimate
λa ∗ as in (2.1.12) and then inserting zeros for all inactive constraints.
dλ∗i (pnom )
λ̂∗i (ppert ) = λ∗i (pnom ) + (pj,pert − pj,nom ). (2.1.15)
dpj
λ∗i (pnom )
p̂j,pert (λ∗i (pnom )) = pj,nom − dλ∗i (pnom )
. (2.1.16)
dpj
20 2 Numerical Basics
As only the impact of single parameters on the active set are considered
with the other parameters remaining at their respective nominal values, it is
not necessarily admissible to superpose ranges of perturbations of multiple
parameters in order to find the true admissible range.
d ∗ d ∗
ŵ∗ (ppert ) = w∗ (pnom ) + w (pnom )Δp1 + ŵ (ppert,adm − p)Δp2
dp dp
(2.1.18)
d ∗ d ∗
λ̂∗ (ppert ) = λ∗ (pnom ) + λ (pnom )Δp1 + λ̂ (ppert,adm − p)Δp2 .
dp dp
(2.1.19)
proximation at the bounds of the admissible parameter range. For the latter
part Δp2 = (ppert − ppert,adm ), i.e. the perturbation exceeding the admissi-
d
ble range, additionally sensitivities dp ŵ∗ (ppert,adm − p), dp
d ∗
λ̂ (ppert,adm − p)
are computed at the boundary point which are then used for a similar
approximation. Note that these sensitivities are directional sensitivities,
i.e. they are tangential to the hyperrectangle of constraint approximations.
They are obtained by modifying the KKT system (2.1.9) in order to reflect
the active set at the bounds of the admissible parameter range. Therefore
knowledge of the correspondence of the admissible range limits and the
responsible constraints is required which can be obtained as a byproduct
by (2.1.14) and (2.1.16).
∂w L(w, p, λ)
∂
⎜ ⎟
⎜ ⎟
⎜ ⎟
⎜
⎜
g1 (w, p) ⎟
⎟
⎜ ⎟
⎜ g2 (w, p) ⎟
⎜ ⎟
⎜ .. ⎟
⎜ . ⎟
⎜ ⎟
⎜ ⎟
ρ(w, λ, p) = ⎜
⎜ gN eq
(w, p) ⎟
⎟ (2.1.23)
⎜ ⎟
⎜
⎜
⎜
min−gNeq +1 (w, p), λNeq +1 ⎟ ⎟
⎟
⎜ ⎟
⎜
⎜ min −gNeq +2 (w, p), λNeq +2 ⎟
⎟
⎜ .. ⎟
⎜ . ⎟
⎜ ⎟
⎝ ⎠
min −gNeq +Nineq (w, p), λNeq +Nineq 1
For the parametric NLP problems considered in this section, the function
from the referenced publications required extension by the parameters p,
i.e. ρ(w, λ, p).
Restoring Admissibility
The first-order approximations (2.1.11) and (2.1.12) yield a result that is
not necessarily admissible w.r.t. the (in)equality constraints. However, for
small perturbations (that do not change the active set) the stability of the
solution can be exploited by the Newton-type iterative Algorithm 1 in
order to restore admissibility, cf. [21].
It should be noted, that, by adding constraints to the active set, the LICQ
can be violated in certain cases. Therefore, possible LICQ violations may
be avoided by either pre-selecting constraints such that their gradients
26 2 Numerical Basics
may never be linearly dependent which may not be suitable for practical
applications. Alternatively, the optimization problem can be relaxed by
removing linearly dependent constraints from the iteration scheme, e.g. by
assigning priorities to constraints an applying a Gram-Schmidt process.
w∗ = arg minimize
w
J(w) = (w1 + 1)2 + (w2 − 2)2 (2.1.28a)
⎛ ⎞ ⎛ ⎞
g1 −w1 + p ⎠
s.t. g(w, p) = ⎝ ⎠ = ⎝ ≤0 (2.1.28b)
g2 2w1 + w2
with the optimization variables w = w1 w2 and the parameter p.
Note that, as the objective function is quadratic and the inequality
constraints are linear in w, the NLP is convex which eases its solution.
2.1 Numerical Optimization 27
For a nominal parameter pnom = 1 the solution wnom = 1 2 with the
multipliers λnom = 4 0 was found using qpOASES . The constraint
residues at the nominal solution are g = 0 −2 . Thus only constraint
g1 is found to be active while g2 is inactive, thus ga = g1 .
The derivatives
⎛ ⎞ ⎡ ⎤
∂ 2(w1 + 1)⎠ ∂ −1 0⎦
J(w, p) =⎝ g(w, p) = ⎣
∂w 2(w2 − 2) ∂w 2 1
w λ g admissible tCPU in ms
exact 2.5 1 11 2 · 12.1
directional λ = 11 8 12.5
QP 11.5
iterative+ λ = 11 9 12.1
From the results in Table 2.1 it can be seen that for ppert = 2.5 the
approach based on the QP approximation (2.1.26) converges to the exact
optimal solution. The first-order approximations (2.1.18) and (2.1.19)
based on directional sensitivities also yields the exact solution for the
optimization variables as well as a deviation of the multipliers which is
also the case for the iterative+ method.
All methods produce results that are admissible w.r.t. the constraints
imposed onto the optimization problem. For this example, the QP
is the best solution approximation method as it yields the optimal
solution (which is also admissible) and requires the least computation
time. However, in this particular case, the original NLP problem is
a QP problem of the same dimension as its QP approximation and
thus requires similar computational effort to obtain an exact solution.
2.1 Numerical Optimization 29
Ultimately, computing the exact solution for this rather small problem
(with or without perturbations) is the best option in this case.
For a perturbed parameter ppert = 2.5 the results of the exact solution
as well as for its approximations are reported in Table 2.2.
Table 2.2: Parametric NLP (16.10) from [39]: Results for ppert = 2.5, tCPU in ms.
d ∗
As the only constraint is inactive, sensitivities dp λ do not exist. Fur-
d ∗
thermore, the second column of dp w is zero as p2 is only included in the
inequality constraint g which is inactive at the nominal solution. The
admissible range of parameters is Ppert = (−∞, 0.3125] × [0.1667, ∞).
In this example, the perturbed parameter is set to ppert = 0.25 0
where only the second element changes that forces the inequality con-
2.1 Numerical Optimization 31
While all methods correctly estimate the active set, the improved it-
erative scheme is the only one whose result is admissible w.r.t. the
constraints as the others work only with first-order constraint approxi-
mations but without the constraint functions themselves. Due to this
feature, in this example the improved iterative scheme is superior to
the other methods and also to the exact solution as it requires only 6 %
of its computation time and yields a similar result.
For the solution of the above NLP examples, the interior-point solver Ipopt
[133] was employed. Its enhancement, sIpopt [106], is also capable of
computing solution approximations for parameter perturbations by means
of parameter sensitivities. However, results obtained with sIpopt were
found to be inadmissible, i.e. violate inequality constraints in the perturbed
case. This behavior was confirmed by the developers [13] as the software
– in constrast to the iterative+ method from Algorithm 2 in the present
thesis – cannot properly handle active set changes.
32 2 Numerical Basics
t − tj tj+p+1 − t
Nj,p (t) = Nj,p−1 (t) + Nj+1,p−1 (t) (2.2.2a)
tj+p − tj tj+p+1 − tj+1
⎧
⎨ 1 t ≤ t < t
j j+1
Nj,0 (t) = ⎩ (2.2.2b)
0 otherwise
which are piecewise polynomials. It follows that Nj,p (t) = 0, t ∈ / [tj , tj+p+1 ),
and consequently, only the functions Nj−p,p (t), . . . , Nj,p (t) are non-zero on
the knot interval [tj , tj+p+1 ). This fact can be exploited in the evaluation
of (2.2.1). The curve q(t) is defined for t ∈ [a, b] on which m + 1 knots
2.2 B-Spline Curves 33
a = t0 = t1 = · · · = tp ≤ tp+1 ≤ . . .
· · · ≤ tm−p−1 ≤ tm−p = · · · = tm−1 = tm = b.
Relevant properties
Matrix-Vector Representation
By replacing the sum in (2.2.1) with a vector product, the curve can be
computed as
wherein the basis function vector is denoted n(t) and the vector of control
points is d with
n (t) = N0,p (t) N1,p (t) · · · Nn−1,p (t) Nn,p (t)
and
d = d0 d1 · · · dn−1 dn .
with the vector of parameter powers τ (t) and the coefficients nj,i , i = 0, . . . , p
which are piecewise constant on each knot interval. Now (2.2.3) can be
rewritten as
that is constant on each knot interval. Note that it is also sparse due to
the limited parameter domain with non-zero values of each basis function
Nj,p (t).
d n−1
(1)
q(t) = q (1) (t) = Nj+1,p−1 (t)dj (2.2.6)
dt j=0
36 2 Numerical Basics
with degree p − 1. The n control points of the resulting spline (in contrast
to the n + 1 control points of the original spline) are
(1) p
dj = (dj+1 − dj ), j = 0, . . . , n − 1 (2.2.7)
tj+p+1 − tj+1
which are linear combinations of the control points of the original B-spline
curve. Furthermore, the boundary knots are dropped, i.e. m − 1 knots
(1)
remain. The control points dj can be collected in the vector d(1) =
(1) (1) (1)
d0 d1 . . . dn−1 .
This scheme can of course be repeated for higher derivatives in a similar
manner, e.g. to compute the control points d(2) of the second-derivative
B-spline curve.
Matrix-Vector Representation
Derivation of (2.2.5) w.r.t. t is straightforward as effectively only the vector
of parameter powers needs to be derived as the coefficient matrix is constant
on each interval. With
d
τ (t) = 0 1 · · · (n − 1)tn−2 ntn−1
dt ⎡ ⎤
⎢0 · · · · · · · · · · · · · · · 0⎥ ⎛ 1 ⎞
⎢ ... .. ⎥⎥ ⎜ ⎟
⎢1 .⎥ ⎜ ⎟
⎢ ⎜ t ⎟
⎢ ⎥⎜ ⎟
⎢ .. ⎥⎥ ⎜ 2 ⎟
⎢0
⎢ 2 ... .⎥ ⎜ ⎜ t ⎟
⎟
⎢ .. ⎥⎥ ⎜ .. ⎟
= ⎢⎢ ... ... 3 ... . ⎥⎥ ⎜ ⎜ . ⎟
⎟
⎢
⎢ ..
⎢.
... ... ... .. ⎥⎥ ⎜ ⎜ n−2 ⎟
t
⎟
⎢ .⎥ ⎜ ⎜
⎟
⎟
⎢. . . . n − 1 . . . ... ⎥⎥ ⎜ ⎟
⎢. ⎜tn−1 ⎟
⎢. ⎥⎝ ⎠
⎣ ⎦
0 ··· ··· ··· 0 n 0 tn
T τ (t)
d2
τ (t) = TTτ (t)
dt2
..
.
2.2 B-Spline Curves 37
d
q(t) = τ (t)T N(t)d (2.2.8a)
dt
d2
q(t) = τ (t)T T N(t)d (2.2.8b)
dt2
..
.
d d
q(t) = τ (t)T N(t). (2.2.10)
dd dt
While the properties and applications described in this section are fairly
commonly known in robotics, a software library that includes additional
functionality was developed. A brief overview of features can be found in
Appendix A.
Pure approximation, i.e. computing a curve that fits best for a given set of
points, can be formulated as a linear programming (LP) problem. Here,
it is formalized as a QP in order to perform quadratic curve fitting. A
given point q a,i should be fitted by the curve that assumes the value q(ti )
at a given parameter ti . Approximation points are collected in the vector
qa should be fitted by the curve that assumes the values qa at the same
respective parameter values. The same holds for the curve derivatives, e.g.
for the first derivative qa(1) ≈ q(1)
a .
In this context, the word best is used synonymously with optimal. Thus not
only a control trajectory u(t) that solves the main challenge is required, it
needs to stand out from all other possible trajectories in a quantitative man-
ner. The latter qualification can be expressed by minimizing an objective
function within an numerical optimization problem, cf. Section 2.1.
&T
minimize L(x(t), u(t))dt + E(x(T )) (2.3.1a)
x(t),u(t)
0
s.t. x0 − x(0) = 0 (2.3.1b)
r(x(T )) ≤ 0 (2.3.1c)
ẋ(t) − f (x(t), u(t)) = 0 (2.3.1d)
h(x(t), u(t)) ≤ 0, ∀t ∈ [0, T ] (2.3.1e)
wherein the state x and the control u are functions of time t ∈ [0, T ]
and serve as optimization variables. The terminal time is denoted T . The
objective consists of the integral cost L(x(t), u(t)), referred to as Lagrange
term, and the terminal cost E(x(T )), named the Mayer term. This
particular combination results in a so-called objective of Bolza type. The
OCP is subjected to initial and terminal constraints, (2.3.1b) and (2.3.1c),
respectively, as well as the system dynamics (2.3.1d) and constraints (2.3.1e)
that apply during the path. While (2.3.1b) is a constraint that is linear w.r.t.
the (inital) state, the other constraints are in general nonlinear w.r.t. state
and control. The terminal time T may be fixed as denoted in the above
OCP, or free through reparametrization of t and with T as the terminal
value of an additional state xT = t with its time derivative ẋT = 1.
Below, the most common direct methods, i.e. single and multiple shooting,
as well as collocation, will be discussed. Therein, the shorthand (·)j := (·)(tj )
is used for expressing points along state and control trajectories.
42 2 Numerical Basics
The considered system is only influenced by the control u(t), therefore the
constraint (2.3.1b) on the initial state is no longer included. The OCP
(2.3.2) is constrained by the system dynamics (2.3.1d), included in (2.3.2c)
as numerical simulation by forward time integration, starting at the initial
state x(0). Therefore the state depends nonlinearly from the input. Thus,
even constraints h linear w.r.t. the state may then depend nonlinearly from
the states and from the input.
1 In the textbook [11] the author reports (without reference), that the name shooting stems
from a control example in which a cannon is required be optimally aimed to shoot at a target.
2.3 Optimal Control Problems and Direct Methods 43
The constraints (2.3.2d) can be typically not enforced over the whole interval
[0, T ] but are checked on Nch discrete checkpoints tk ∈ [0, T ], k = 1, . . . , Nch .
The resulting finitely dimensional OCP thus yields
&T
minimize
w
L(x(t), u(t))dt + E(x(T )) (2.3.3a)
u
0
s.t. r(x(T )) ≤ 0 (2.3.3b)
&T
x(t) = x(0) + f (x(t), u(t))dt (2.3.3c)
0
h(x(tk ), u(tk )) ≤ 0, tk ∈ [0, T ], k = 1, . . . , Nch (2.3.3d)
Both, the variables representing the states and the controls at certain points
in time are collected in the vector
w = x 0 u 0 x 1 u 1 · · · x N u N .
44 2 Numerical Basics
N
&tj
minimize
w
L(x(t), u(t))dt + E(xN ) (2.3.4a)
j=1 tj−1
s.t. x0 − x0 = 0 (2.3.4b)
r(xN ) ≤ 0 (2.3.4c)
&tj
xj − f (x(t), u(t))dt = 0, j = 1, . . . , N (2.3.4d)
tj−1
The additional state variables allow for shorter integration time intervals,
reducing numerical errors occurring for longer times. They also enable
direct access to the system behavior, allowing to handle especially unstable
systems and systems with high sensitivity w.r.t. the initial value well.
While the dimension of the resulting KKT system is significantly increased
compared to DSS, its sparsity is improved as intermediate states have no
explicit dependency on all previous states and controls.
The fact that the integration over an interval does not depend on the
result of any of its predecessors enables simultaneous computation of all
integrations which allows for parallel computing.
2.3 Optimal Control Problems and Direct Methods 45
While a basic DMS method can be implemented in less than 100 lines of code,
an enhanced and particularly efficient tool to solve OCPs is MUSCOD-II
[75].
N
&tj
minimize
w ,w
L(x(t), u(t))dt + E(xN ) (2.3.5a)
x u
j=1 tj−1
parametrizing the time evolution of the state x(t). The control trajectory
u(t) is still parametrized in terms of wu .
N
&tj
minimize
w ,w
L(x(t), u(t))dt + E(xN ) (2.3.6a)
x u
j=1 tj−1
s.t. x0 − x0 = 0 (2.3.6b)
r(xN ) ≤ 0 (2.3.6c)
ẋj − f (xj , uj ) = 0, j = 0, . . . , N − 1 (2.3.6d)
h(xk , uk ) ≤ 0, k = 1, . . . , Nch (2.3.6e)
In this thesis, direct collocation with B-spline curves is used for trajectory
planning in Chapter 4.
Numerical Solution
In order to compute a solution to a practical OCP, the considered system
behavior first needs to be formalized in terms of the equations govern-
ing relevant physical processes in a dynamics model. Furthermore, the
optimization goal needs to be expressed in terms of dynamics quantities.
Then, the above OCP formulations are applied to obtain a formulation
ready for numerical solution. Software frameworks used to implement such
formulations are CasADi [2] or YALMIP [79]. The numerical solution is
obtained with an appropriately selected solver program such as Ipopt [133]
for general NLP problems, or qpOASES [44] for QP problems.
Chapter 3
Kinematics of Serial
Robots
Kinematics describes the motion of a system of one or more bodies by
considering their position, velocity, acceleration and possibly higher deriva-
tives, i.e. disregarding forces and torques acting upon them. As this field
is well-researched, the present chapter does not aim to provide a compre-
hensive discussion of existing methods. Interested readers are referred to
general robotics textbooks such as [127] while advanced readers may prefer
publications focusing on kinematics such as [86]. However, the approaches
most relevant for the successive chapters are discussed in great level of
detail.
(A)
Figure 3.1: Position Vector from Point (A) to Point (B).
wherein R xAB , R yAB and R zAB are the distances between A and B measured
in the x, y and z directions of the R frame, respectively. The index on
the left, here R, denotes the frame in whose coordinates the vector is
represented in.
3.1.2 Orientation
Rotation Matrix
Orientation describes the relative rotation between two coordinate frames
which can be uniquely represented by a rotation matrix RRS . Such a
construct not only characterizes the relative rotation between the two coor-
dinate frames indexed with R and S, it is also a coordinate transformation
between these frames. Therefore, the rotation matrix
' (
RRS = R (S e1 ) R (S e2 ) R (S e3 ) ∈ SO(3) (3.1.2)
comprises three mutually orthogonal column vectors that are the unit vectors
aligned with the principal directions of frame S resolved in coordinates of
the R system. It is an element of the special orthogonal group in three
dimensions SO(3). Note that this representation is redundant, i.e. its nine
parameters are more than necessary to fully describe the rotation.
In the larger part of this thesis this rotation matrix representation will be
used for expressing general coordinate transformations as well as orientations
of coordinate frame, in particular a robot’s EE orientation. However, for
52 3 Kinematics of Serial Robots
Quaternions
Another way of uniquely representing orientation are quaternions which
extend the complex numbers. As this representation will not be used in the
present work, only this brief overview is given and the reader is referred to
textbooks such as [86].
Similar to complex numbers that can represent planar rotations on the unit
circle in the complex plane, quaternions represent spatial rotations on the
unit sphere. Using four numbers, a quaternion
⎛ ⎞
q
⎜ 1⎟
⎜ ⎟
Q = (q0 , q) with q = ⎜q2 ⎟,
⎝ ⎠
qi ∈ R, i = 1, . . . , 3
q3
can be used to represent any point in SO(3). Therein, q0 is called the scalar
compontent, and q is called the vector component of the quaternion Q. It is
called a unit quaternion if q = 1 which was assumed in the introductory
explanation using the unit sphere.
Successive Rotations
In contrast to the above representations that provide a unique orientation
representation, spacial rotations can also be expressed by means of three
successive elementary rotations about alternating axes of the accordingly
3.2 Translation and Rotation Velocity 53
such as Euler angles, Kardan angles about axes of moving frames, or roll-
pitch-yaw angles about axes of fixed frames. There are in total 12 different
rotation orders. For the inverse problem, i.e. finding angles ϕi , i = {1, 2, 3},
all of these representations suffer from singularity issues if rotations about
the first and last axis cannot be distinguished.
d
0 vB = 0r (3.2.1)
dt B
= 0 ṙB .
with the spin tensor1 of the angular velocity R ω0R of the frame R w.r.t. the
inertially fixed frame 0.
Note that the angular velocity ω is in general different from the time
derivative of an orientation representation by means of a vector ϕ comprising
three angles, such as Euler angles, cf. Section 3.1.2.
The translation velocity of a point (E) on a body and the angular velocity
of a body-fixed frame, indicated by the index E, are collected in the twist
vector
⎛ ⎞
ωE ⎠
⎝
VE = . (3.2.3)
vE
⎛
⎞ ⎡ ⎤
ωx 0 −ωz ωy
1 spin tensor: ω̃ = ⎜ω ⎟ = ⎢ ω ⎥
−ωx ⎦
⎝ y⎠ ⎣ z 0
ωz −ωy ωx 0
3.3 Joint Space, Work Space, Task Space 55
The reachable work space is the set of positions the EE can assume with
some orientation, i.e.
Such points may include intermediate points of via paths, where the orien-
tation of the EE is often irrelevant. Determining WR is a simple task as
only the joint coordinates need to be swept in their entire admissible range
while recording the respective resulting EE positions.
Example 3.4: SCARA4 (reachable work space)
In order to find the reachable work space WR of the SCARA4 manip-
ulator (with joint limits neglected), it is sufficient to first consider all
points the EE can assume along one radial line (here assumed as the
0 x axis w.l.o.g.). The rest of the reachable work space is then found by
circularly sweeping these reachable points around the first joint.
0y
WR
rout
(E)
(0) 0x
The dexterous work space is the set of positions the EE can assume with
arbitrary orientation, i.e.
Dexterous points are relevant for, e.g. machining tasks in which one
particular point is required to be accessible from all directions, or robot-
based photography for which arbitrary EE orientation is also of significance.
In contrast to the rather effortless task of finding WR , deriving WD is more
involved.
rout WD
(E)
(0) 0x
In any case, once path planning is completed for a specific task, the operation
of the manipulator is restricted to the task space WT ⊆ W. WT denotes
the EE space necessary to execute a given task.
Note that not only the joint positions are relevant for the result of the
forward kinematics mapping, also constant geometric parameters along the
serial kinematic structure of the manipulator under consideration such as
link lengths are required.
As VE is linear w.r.t. the joint velocity q̇, the geometric Jacobian matrix
∂VE
J(q) = (3.4.4)
∂ q̇
VE = J(q)q̇. (3.4.5)
60 3 Kinematics of Serial Robots
with
⎡ ⎤ ⎛ ⎞
*
4 *
4 *
4 *
i
⎢cos qi − sin qi 0⎥ ⎜
⎜i=1
li cos qj ⎟
⎢ i=1 i=1 ⎥ ⎜ j=1 ⎟ ⎟
⎢ ⎥ ⎜* ⎟
R0E (q) = ⎢ *
4 *
4 ⎥ and ⎜
4 *
i
⎟.
⎢ sin qi cos qi 0⎥⎥ 0 rE (q) = ⎜ l sin q ⎟
⎢ ⎜ i=1 i j ⎟
⎣ i=1 i=1 ⎦ ⎝ j=1 ⎠
0 0 1 0
= J(q)q̇
and
⎛
⎞
*4 *i *i
⎜− li sin q j q̇ j ⎟
d0 rE (q) ⎜ i=1
⎜ j=1 j=1
⎟ ⎟
0 vE (q, q̇) = = ⎜ *4 *
i *i ⎟
⎝
dt li cos qj q̇j ⎠
i=1 j=1 j=1
and therefore
⎡
⎢
1 1 ···
⎢ *
4 *
i *
4 *
i
⎢
J(q) =
∂ 0 VE
= ⎢
⎢
− li sin qj − li sin qj ···
⎢ i=1 j=1 i=2 j=1
∂ q̇ ⎢ *
4 *
i *
4 *
i
⎣
li cos qj li cos qj ···
i=1 j=1 i=2 j=1
⎤
··· 1 1 ⎥
*
4 *
i *
4 ⎥
··· − li sin qj −l4 sin qj ⎥⎥⎥.
i=3 j=1 j=1 ⎥
*
4 *
i *
4 ⎥
⎦
··· li cos qj l4 cos qj
i=3 j=1 j=1
62 3 Kinematics of Serial Robots
(E) q4
l4 q3
Ey
l3 q2
l2
l1 q1
0x
(0)
Figure 3.4: Forward Kinematics of the SCARA4 Manipulator.
Also depending on the robot’s structure, more than one joint configuration
can correspond to a particular admissible EE configuration. Depending
on the manipulator and the desired EE configuration, there may be one
single solution, or multiple or no solutions at all to the IK problem on
position-level. However, for kinematically redundant structures, an infinite
amount of configurations correspond to a given EE pose. Of course, now
the question arises which configuration is the best and should therefore be
desired. And also, what does the best actually mean in this case?
Closed-Form Solutions
In rare cases, such an IK relationship can be found in closed-form, i.e.
in terms of functions and mathematical operations, excluding iterations
and time integration. This is facilitated by means of thorough investiga-
tion of a manipulator’s geometry, which of course is only applicable for
kinematically non-redundant structures. In the following, two particular
types of manipulators are considered that are of great relevance for robotics
applications.
0y
(E) q3
l3 ϕE
q2
yE
l2
l1
yG up
q1
0x
down
(0) xG
xE
q1 = α + β
q2 = arctan(yG − l1 sin q1 , xG − l1 cos q1 ) − q1
q3 = ϕ0E − q1 − q2
⎛ ⎞ ⎛ ⎞
xG cos ϕ0E ⎠
with the auxiliary quantities ⎝ ⎠ = 0 rE − l3 ⎝ ,α =
yG
sin ϕ0E
l12 +l2 −l22
arctan(yG , xG ), l2 = x2G + yG
2
, β = ± arccos 2l1 l .
Note, that there are in general up to two joint configurations that solve
the IK problem. In the above case, exactly two solutions exist, i.e. the
3.5 Inverse Kinematics 67
For a spatial manipulator with 6 DOF sufficient conditions for the existence
of a closed-form IK solution are [127]
A notable structure that fulfills the first of the above criteria is a 6 DOF
standard industrial robot with a spherical wrist, cf. the robot on the
linear axis from Figure 1.3. Not only because of the EE positioning and
orientation capabilities, but also due to the simplicity of the corresponding
IK resolution problem, manipulators with this structure are widely used in
many applications. For a thorough investigation of the IK of such structures
the reader is referred to [62].
The key idea is a separation of the joint coordinates q into two disjoint
coordinate sets: the non-redundant joint coordinates qnr and the redundant
coordinates qr . The relationship can be expressed by means of a selector
matrix S, i.e.
⎛ ⎞ ⎡ ⎤
q
⎝ nr ⎠
S
⎣ nr ⎦q.
=
qr Sr
68 3 Kinematics of Serial Robots
Of course, this decomposition alone does not enable solving the IK problem.
However, a solution for the non-redundant coordinates parametrized in
terms of the redundant coordinates can be found, i.e. qnr = qnr (R0E , rE , qr).
0y
(E) q4
ϕ0E
l4 q3
yE q2
l3
l2
l1 q1
0x
(0) xE
redundant coordinate qr , thus the first three joints are the non-redundant
coordinates qnr which is formalized by
⎛ ⎞ ⎛ ⎞
⎛ ⎞ ⎡ ⎤ q1 q
q 1 0 0 0⎥⎜ ⎟ ⎜ 1⎟
⎜ 1⎟ ⎢ ⎜ ⎟ ' (⎜ ⎟
⎜q2 ⎟ ⎜q 2 ⎟
qnr = ⎜ ⎟ ⎢
⎜q2 ⎟ = ⎢0 1 0
⎝ ⎠ ⎣
0⎥⎥⎦⎜
⎜ ⎟
⎟ and qr = q4 = 0 0 0 1 ⎜
⎜ ⎟.
⎟
⎜q3 ⎟ ⎜ ⎟
q3 0 0 1 0 ⎝ ⎠ ⎝q 3 ⎠
Sr
q4 q4
Snr
Iterative Method
The definition of the velocity Jacobian (3.4.4) gives rise to an iterative IK
scheme that operates on position-level. It can be used for point-wise com-
putation of a joint configuration q that corresponds to an EE configuration
zE (q). Starting at a joint configuration that not necessarily already agrees
with the desired EE configuration, the EE configuration error
⎛ ⎞
er (q) ⎠
e(q) = ⎝ (3.5.1)
eR (q)
occurs. It consists of the EE position error er (q) = rE,d − rE (q) and its
orientation error eR (q). While the former is simply obtained from the EE
70 3 Kinematics of Serial Robots
position vector difference, the latter is more involved. With the desired and
actual EE orientation in rotation matrix representation, i.e.
' (
R0E,d = 0 (E e1,d ) 0 (E e2,d ) 0 (E e3,d ) (3.5.2)
' (
R0E (q) = 0 (E e1 (q)) 0 (E e2 (q)) 0 (E e3 (q)) (3.5.3)
with the spin tensor of the orientation error2 0 ẽR = log Re and its explicit
[80] angle-axis representation, cf. Section 3.1.2,
1
0 eR (q) = (E e1 (q)) × 0 (E e1,d ) + 0 (E e2 (q)) × 0 (E e2,d ) + . . .
2 0
+ 0 (E e3 (q)) × 0 (E e3,d ) (3.5.6a)
= 0 uR (q) sin ϑ(q) (3.5.6b)
with the orientation error axis uR and the orientation error angle ϑ. Note
that this representation is w.r.t. the inertially-fixed frame of reference and
only unique for |ϑ| < π2 .
From the forward kinematics relationship (3.5.7) the iterative scheme Algo-
rithm 3 is obtained. Therein, the joint angles q are updated until the EE
error is sufficiently reduced.
2 Note, that the time derivative of e(q) does in general not correspond to the EE velocity
Note, that (3.5.7) is only a linearized form of the actual forward kinematics
relationship (3.4.1). Therefore even an updated joint configuration q + Δq
with a Δq that solves (3.5.7) only approximates, but not necessarily agrees
with the desired EE configuration zE (q).
Not only the resolution of the differential IK problem is of interest, but also
joint positions that correspond to the desired EE pose are relevant. However,
the methods considered in this section only provide such information via
numerical time integration of the resulting joint velocities q̇. By nature
of computer implementations of numerical integration schemes, errors are
introduced that accumulate over (simulation) time yielding a non-negligible
deviation from the desired EE path. Therefore, a stabilizing feedback loop
is introduced that reduces these undesired effects by means of closed-loop
inverse kinematics (CLIK), originally introduced in [137].
Note, that only the solution of a system of linear equation is required rather
than an explicit representation of the Jacobian inverse J(q)−1 .
the solution
−1
q̇ = J J + αI J VE,d (3.5.10)
is obtained.
It is a minimum
as the second-order sufficient condition
2
H = ∂∂ q̇J2 = J J + αI > 0 is satisfied. The result (3.5.10) is also called
the damped least squares solution of the IK problem.
While the above holds for kinematically non-redundant robots, in the redun-
dant case, the IK problem is underdetermined as there are fewer equations
74 3 Kinematics of Serial Robots
wherein the squared joint velocities, weighted by the diagonal matrix W > 0,
are minimized in (3.5.11a) with the path tracking constraint (3.5.11b). The
Lagrange function corresponding to this QP yields
and thus
the solution
−1
(3.5.13), (3.5.14c) → q̇ = W−1 J (q) J(q)W−1 J (q) VE,d (3.5.15)
J+
W (q)
2
is found. It is a minimum of (3.5.11a) as H = ∂∂ q̇L2 = W > 0. Therein,
J+W a generalized inverse that is weighted by W. With the substitution
W = I it is also called the right Moore-Penrose inverse that is denoted
J+ (q). A detailed overview of other generalized inverses is found in [10].
It should be noted that a closed, i.e. circular, work space path does not
necessarily yield a closed joint space path [68]. Thus, such IK schemes are
3.5 Inverse Kinematics 75
not suitable for the unsupervised execution of repetitive tasks [5]. This
issue will be addressed in the augmented IK schemes below where such
additional criteria can be incorporated.
αq̇ − J λ = 0 (3.5.18a)
ė − λ = 0 (3.5.18b)
such that the augmented Jacobian Jaug is square. As a result, the corre-
sponding IK problem yields a single solution
⎛ ⎞
V
q̇ = J−1
aug
⎝ E,d ⎠ (3.5.22)
0
3.5 Inverse Kinematics 77
A special case is the extended Jacobian technique [5] wherein the condition
∂
∈
/ ker J holds. The extended Jacobian technique is not only used
∂q g
in industrial robotics as considered in this thesis, but also for mobile,
non-holonomic manipulators, e.g. in [131, 129].
In the application [74] of the characteristic length method [3], the Jacobian
is augmented such that its condition number is minimized. Therein, an
SQP method is employed to find point-wise optimal poses for a 6 DOF
robot performing a 5 DOF task.
q̇ = J (q)We (3.5.23)
wherein the error w.r.t. the desired joint velocities q̇0 is minimized as a sec-
ondary goal to tracking the prescribed EE velocities VE,d . The Lagrange
function corresponding to the QP (3.5.26) yields
⎛ ⎞
∂L(q̇, λ) ⎠
⎝ = −q̇0 + q̇ − Jλ = 0 (3.5.28a)
∂ q̇
(3.5.28a) → JJ λ = J(q̇ − q̇0 ) (3.5.28b)
−1
(3.5.28b), (3.5.26b) → λ = JJ (VE,d − Jq̇0 ). (3.5.28c)
With the velocities isolated from (3.5.28a) and the result (3.5.28c) for the
Lagrange multipliers λ, the IK scheme yields
(3.5.28a) → q̇ = J λ + q̇0
−1
(3.5.28c) → = J JJ (VE,d − Jq̇0 ) + q̇0
⎛ ⎞
⎜ ⎟
−1 ⎜ −1
= J JJ
VE,d + ⎜I
⎝
− J JJ
J⎟
⎟q̇0
⎠
J+ J+
q̇ = J+ (q)VE,d + I − J+ (q)J(q) q̇0 . (3.5.29)
N(q)
With the relationship (3.5.29) established now regard Figure 3.7. It shows
an illustration of the domains occuring for transformations using the velocity
Jacobian J. It transforms joint velocities q̇ into its image im J resulting
in corresponding EE velocities. Non-zero joint velocities that lie in the
Jacobian nullspace ker J do not result in EE motion. This inherent property
of kinematical redundacy can now be exploited by purposefully assigning
nullspace joint velocities Nq̇0 that do not change the EE behavior but
improve the joint configuration.
Again, the question comes up, what improve means in this context. If
not only an instantaneous velocity relationship is considered but a whole
path, i.e. a sequence of points, is treated with IK methods, an improved
configuration allows the manipulator to track that path so that additional
goals are fulfilled. These typically range from a minimum traversal time
over motion reduction of the individual joints to energy-efficient operation
if additionally considering the system dynamics.
q̇
VE
ker J
0
im J
of kinematic manipulability in an DIK resolution scheme as a local optimization criteria, i.e. the
configuration is drawn to a local maximum of manipulability. It is a common statement, e.g. in
[126], that it represents a distance from singular configurations. However, in general it does not
provide a metric to do so, i.e. it cannot be used to quantify how singular a configuration is.
82 3 Kinematics of Serial Robots
The idea of the result (3.5.32) that tracks a prescribed EE trajectory as its
primary goal while (locally) pursuing an additional, secondary objective
H(q) gives rise to IK schemes that incorporate priorities of different tasks,
e.g. [59, 81, 89]. Therein, tasks of lower priority are projected onto the
Jacobian nullspace of higher prioritized tasks. Therein, methods exist that
make a trade-off between accuracy and strict task prioritization, e.g. [29].
VE = J(q)q̇
= Jnr (q)q̇nr + Jr (q)q̇r (3.5.33)
Assuming the redundant joint velocities q̇r as another input next to the
desired EE twist VE,d , the DIK scheme (3.5.8) is rearranged as
not limit the EE motion. Such singularities are artificial and can therefore
be avoided by selecting a decomposition that is suitable for the current
joint configuration of the considered manipulator, i.e. with a full rank Jnr .
Note that the joints selected to be in the non-redundant set are required to
be capable of performing the prescribed EE velocity motion without the
redundant set.
Example 3.13: SCARA4 (JSD algorithmic singularities)
In this example, JSD is performed for the inherently kinematically
redundant SCARA4 manipulator that is in the configuration depicted
below in Figure 3.8.
For a task that not only involves the EE position but also its orientation,
i.e. the manipulator has a degree of kinematical redundancy of r = 1,
cf. Example 3.8, thus the selection of each joint yields a possible choice
regarding JSD.
For the particular configuration depicted in Figure 3.8, only two of the
choices, namely qr = q1 and qr = q4 yield a full rank non-redundant
Jacobian Jnr . The other choices will result in an algorithmic singularity
of Jnr due to q3 = 0, i.e. the columns corresponding to q̇2 and q̇3 are
linearly dependent.
0y
q4
l4
q2
l3
l2
l1 q1
0x
Applications of DIK JSD are not only found in robots with fixed base
[41] but also in mobile robotics [37, 83]. One particularly interesting work
therein is [37], where DIK JSD is used in conjunction with the reduced
gradient method that is discussed below.
For higher derivatives, e.g. acceleration, jerk and jounce, first time deriva-
tives of (3.5.33) are considered, i.e.
JSD schemes similar to (3.5.34) are then found from solving these equations
for the desired non-redundant quantities, i.e.
q̈nr = J−1
nr V̇E,d − Jr q̈r − J̇q̇ (3.5.36a)
... ...
qnr = J−1
nr V̈E,d − Jr qr − J̈q̇ − 2J̇q̈ (3.5.36b)
.... ... ... ... ...
qnr = J−1
nr VE,d − Jr qr − Jq̇ − 3J̈q̈ − 3J̇q . (3.5.36c)
Therein, terms were collected and simplified as far as possible. Note that
these derivatives are indeed consistent, i.e. (3.5.34) is the time derivative
of (3.5.36a) and so on which is now exemplarily shown:
d d −1
(3.5.34) → q̇nr = Jnr (VE,d − Jr q̇r )
dt dt
d −1 d
= Jnr (VE,d − Jr q̇r ) + J−1 nr (VE,d − Jr q̇r )
dt dt
d
= −J−1 nr (Jnr ) J−1 −1
nr (VE,d − Jr q̇r ) +Jnr V̇E,d − J̇r q̇r − Jr q̈r
dt
q̇nr
= J−1
nr V̇E,d − J̇nr q̇nr − J̇r q̇r − Jr q̈r
(3.5.36a) : q̈nr = J−1
nr V̇E,d − Jr q̈r − J̇q̇ .
the redundant joint velocities q̇r are found such that the time derivative
of a joint pose-dependent performance index H(q), cf. Section 3.5.3, is
maximized, i.e.
d
maximize J(q̇r ) = H(q). (3.5.38)
q̇r dt
⎡ ⎤
−J−1
⎣ nr Jr ⎦
Note, that is an orthogonal complement to J. With
I
d ∂H(q)
H(q) = q̇
dt ∂q
∂H(q) ∂H(q)
(3.5.34) → = q̇nr + q̇r
∂qnr ∂qr
∂H(q) −1 ∂H(q)
= J (VE,d − Jr q̇r ) + q̇r
∂qnr nr ∂qr
⎡ ⎤
∂H(q) −1 ∂H(q) ⎣−J−1nr Jr ⎦
= Jnr VE,d + q̇r (3.5.39)
∂qnr ∂q I
and thus the first-order necessary condition to (3.5.38) follows from (3.5.39)
as
⎛ ⎞
∂J(q̇r )
∂H(q) ⎠
= −(J−1
nr Jr ) I
⎝ = 0. (3.5.40)
∂ q̇r ∂q
(3.5.42)
in which case the EE motion is only governed by the first term in (3.5.37).
With the result (3.5.41) the RG scheme (3.5.37) is rewritten as
⎡ ⎤ ⎡ ⎤ ⎛ ⎞
J−1 −J−1
nr Jr ⎦ ∂H(q) ⎠
q̇ = ⎣ nr ⎦V
E,d + ⎣ −(J−1 J
nr r ) I ⎝
0 I ∂q
⎡ ⎤ ⎡ ⎤⎛ ⎞
J−1
⎣ nr ⎦V
−1 −1
⎢Jnr Jr (Jnr Jr ) −J−1
nr Jr ⎥⎝ ∂H(q) ⎠
= E,d + ⎣ ⎦ (3.5.43)
0 −(J−1nr Jr ) I ∂q
NRG
This method was not only applied to planar [36] but also to spatial and
even mobile non-holonomic manipulators [37].
It should be noted that the RG method does not only exhibit the same
issues regarding artificial singularities as JSD, cf. Section 3.5.3, but also
requires the challenging choice of a pose-dependent performance index H(q)
as in PG schemes discussed in Section 3.5.3.
Any q̇∗ that solves VE,d = Jq̇∗ can be augmented by means of a linear
combination of basis vectors of the Jacobian nullspace, i.e.
r
q̇ = q̇∗ + γi ai (q) (3.5.44)
i=1
= q̇∗ + A(q)γ (3.5.45)
with the nullspace basis vectors ai ∈ ker J(q) columnwisely collected in the
nullspace basis matrix A(q), and the weights γi , collected in the weights
vector γ.
Note that any nullspace augmentation can be rewritten in that form, e.g. a
PG as
⎛ ⎞
∂H(q) ⎠
I − J(q)J (q)
+ ⎝ = A(q)γ. (3.5.46)
∂q
Compared to JSD, for GNA the explicit selection of a redundant set of joint
coordinates is unnecessary and thus artificial singularities do not occur.
Furthermore, GNA does not require the selection of a (pose-dependent)
performance index as for PG and RG methods.
A similar approach has been proposed in archival literature [88, 107, 25,
98]. In the extended Jacobian approach from [96], the weights γi are used
as an additional task that is explicitly prescribed.
which gives
q̇ = J+ VE,d + Aγ
−1
= J+ VE,d − A A A A J+ VE,d
−1
= I − A A A ∗
A q̇ . (3.5.50)
(3.5.51b)
.... ... ... ...
q = J+ VE,d − Jq̇ − 3J̈q̈ − 3J̇q + 6J̇Ȧγ̇ + +3J̇Äγ + 3J̈Aγ̇+
... ...
+3J̈Ȧγ + JAγ + 3J̇Aγ̈ + 3Ȧγ̈ + 3Äγ̇ + Aγ+
−1 ⎫
...
+Aγ
+ 3NJ̇ JJ V̈E,d +
⎪
⎪
⎪
⎪
−1 −1 −1 ⎪
⎪
⎪
+N 3J̈ − 6J̇ JJ JJ̇ − 6J̇ JJ J̇J JJ V̇E,d + ⎪
⎪
⎪
⎪
⎪
... −1 −1 −1 ⎪
⎪
⎪
⎪
+N J − 3J̇ JJ JJ̈ − 6J̇ JJ J̇J̇ − 3J̈ JJ JJ̇ + ⎬ null
−1 −1 −1 −1 ⎪
+6J̇ JJ JJ̇ JJ J̇J + 6J̇ JJ JJ̇ JJ JJ̇ + ⎪
⎪
⎪
⎪
space
−1 −1 ⎪
⎪
⎪
+6J̇ JJ J̇J+ JJ̇ + 6J̇ JJ J̇J+ J̇J + ⎪
⎪
⎪
⎪
−1 −1 −1 ⎪
⎪
⎪
−3J̇ JJ J̇J − 3J̇ JJ J̈J JJ VE,d ⎪
⎭
(3.5.51c)
Other Approaches
Artificial neural networks (ANN) and its application in robotics have been
undergoing extensive research in the last few decades which was mainly
performed by the computer science community. Although many methods in
this field only operate under closed-set conditions, i.e. the situations faced
in real-world operation are already known from the training set wherein
they are included, ANN can also act as a controller instance. Therein,
environmental perception of relevant parameters to the operation conditions,
e.g. computer vision, or fusion of heterogenous sensoric information, plays
an important role. In such open-set conditions, dictated by real-world
operation, the parameters influencing the employment of robotics devices
undergo permanent change, so the training data does not include occuring
in practise. The survey paper [128] provides a comprehensive overview of
current challenges and limitations of the application of ANN in robotics
in general. Genetic algorithms, i.e. population-based stochastic methods,
such as [54] are also used in this context.
Chapter 4
efficiency of the generated plans in terms of select quality criteria, e.g. the
shortest transportation path, a minimum-time or minimum-energy execution
of a trajectory, is increasingly relevant for robot operation. Numerical
optimization can be used for both1 , ensuring admissibility of the results,
and also optimality regarding a specified goal.
In this section, the efficacy of the proposed methods is shown using various
numerical examples. While therein mostly fully actuated industrial robots
are considered, i.e. each joint is driven by its own actuator, the introductory
example is an underactuated model.
1 Note, that such approaches can also be exploited when only admissible, but not optimal
results are required. In this case, optimization is performed w.r.t. a dummy objective.
4.1 Boundary Value Problem 95
q1
gearbox
motor
u
Implementation
The swing-up procedure is thus formalized as the two-point boundary
value problem
ẋ = f (x, u) (4.1.1a)
q(0) = q0 = 0 π (4.1.1b)
q̇(0) = q̇0 = 0 (4.1.1c)
q2 (T ) = q2,T = 0 (4.1.1d)
q̇(T ) = q̇T = 0 (4.1.1e)
with the state x = q q̇ and the additional boundary conditions
u(0) = u(T ) = 0 due to the torque continuity requirement.
Discussion
Of course the function class for u(t) is not restricted to the above
cosine series (4.1.2). Alternatively, for example, a B-spline curve, cf.
Section 2.2, with 6 control points could be used to represent u(t).
4.2 Optimal Point-to-Point Trajectory Planning 97
wherein q(t) describes the time evolution of the robot’s joint angles and
u(t) represents the input trajectory, i.e. actuator forces and torques. The
objective (4.2.1a) is kept very general in order to account for goals frequently
used in robotics such as minimum energy, minimum time, or minimum
4.2 Optimal Point-to-Point Trajectory Planning 99
component wear. For such a PTP problem, particular interest lies in the
specification of initial and terminal conditions (4.2.1b) and (4.2.1c), respec-
tively. The in general nonlinear function r represents the terminal conditions
on position-level, r̃ is used for formalize a velocity-level constraint. Therein,
typically not only the initial and terminal joint angles are considered, but
also their time derivatives which are here (and in the following) omitted
for brevity. The inequality constraints (4.2.1d) and (4.2.1e) represent limi-
tations imposed onto the joint angles (and their derivatives) and onto the
inputs, respectively. In the original OCP, these were included in the path
constraints (2.3.1e) but are explicitly mentioned here due to their relevance
in this application. While previously the system dynamics was considered
in (2.3.1d) in state-space representation, here it is equivalently denoted as
the equations of motion in (4.2.1f). Therein, M(q) is the generalized mass
matrix, g(q, q̇) is a vector of in general nonlinear terms including general
Coriolis, centrifugal and friction forces. The generalized input forces Q(q, u)
are in practice often given by scaled inputs u. Furthermore, the motion is
subjected to actual path constraints (4.2.1g) that impose restrictions arising
from the task to be executed, e.g. work space constraints in multi-robot
cooperation environments, maximum EE velocity in machining tasks, or
limited EE acceleration to avoid sloshing in liquid transportation.
Implementation
The trajectory planning task is formalized as a direct multiple shooting
OCP, cf. (2.3.4), not only to solve the BVP considered previously,
but also to minimize the torque required to drive the pendulum and
100 4 Path and Trajectory Planning
N
&tj
minimize
w ,w
u2 (t)dt (4.2.2a)
x u
j=1 tj−1
Results
The above NLP problem (4.2.2) was implemented using CasADi and
solved with Ipopt. The resulting joint trajectories, velocities and torques
are shown in Figures 4.2 to 4.4. It can be observed, that the result
accomplishes the swing-up motion, i.e. a rest-to-rest motion that drives
both joint angles to zero, and is admissible w.r.t. the technological
constraints of the system (the plots of UA and P are neglected for the
sake of brevity).
4.2 Optimal Point-to-Point Trajectory Planning 101
q in rad q1 q2
0
t in s
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6
−5
−10 t in s
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6
u in Nm
0.4
0.2
−0.2
t in s
−0.4
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6
Problem Statement
Consider the synchronization task depicted in Figure 4.5. Therein, the
SCARA4 manipulator known from the previous examples is required to
pick items from a conveyor belt that is moving with constant velocity
vB . It is assumed, that the object to be picked has the same velocity, i.e.
vObj = vB . In this example, only the computation of the pick trajectory
is considered. Therein, the lateral placement of an item, denoted with
its distance Δx from the base of the robot, is uncertain and therefore
can only be determined immediately before the start of the picking
process. While the initial rest configuration at q0 is the same for all pick
processes, the terminal joint configuration q(T ) varies depending on the
position of the item on the conveyor at pick time T . The pick motion
needs to be executed in minimum time, i.e. the item needs to be picked
as soon as possible. The motion is subjected to physical constraints
regarding the joint angles q, velocities q̇ and motor torques MMot or
accelerations q̈. Moreover, the work space is constrained such that the
EE is not allowed to enter the area behind the conveyor as well as the
area right of the detection line.
4.2 Optimal Point-to-Point Trajectory Planning 103
x
object detection line
vB
terminal
Δx configuration q(T )
y
initial
configuration q0
Implementation
The general OCP (4.2.1) is rewritten in order to fit the present minimum-
time example and reads now
&T
minimize 1dt (4.2.3a)
q(t),T
0
s.t. q0 − q(0) = 0, q̇0 − q̇(0) = 0, q̈(0) = 0 (4.2.3b)
rE (q(T )) − rObj (T ) = 0,
vE (q(T )) − vObj = 0, ωE (q(T )) = 0, q̈(T ) = 0 (4.2.3c)
(4.2.3d)
MMot,min ≤ MMot ≤ MMot,max (4.2.3e)
⎛ ⎞ ⎛ ⎞
x
⎝ min ⎠
xmax ⎠
≤ rE (q(t)) ≤ ⎝ (4.2.3f)
ymin ymax
With the choice for q(t) to be B-spline curves, OCP (4.2.3) is once again
rewritten in order to exploit various properties of this function class.
Therefore, each qi (t), i = 1, . . . , n = 4 is a B-spline curve of degree p = 3
and NCP = 20 control points. The control points of each B-spline curve
describing the time evolution of the i-th joint angle are collected in the
vector dqi ∈ RNCP . Furthermore, the respective j-th control points of
all curves are collected in the vector dq,j ∈ Rn . Thus, the j-th control
point of the i-th curve is denoted dqi ,j . The same notation holds for
time derviatives of the curves. All control points and also the terminal
time T are collected in the vector of optimization variables
w = dq1 ,0 · · · dqn ,0 dq1 ,1 dq2 ,1 · · · dqn ,1 · · · dqn ,0 dqn ,1 · · · dqn ,NCP −1 T .
minimize
w
T (4.2.4a)
s.t. q0 − dq,0 (w) = 0, q̇0 − dq̇,0 (w) = 0, dq̈,0 (w) = 0
(4.2.4b)
rE (q(T ) = dq,NCP −1 (w)) − rObj (T ) = 0, (4.2.4c)
vE (q(T ), q̇(T ) = dq̇,NCP −2 (w)) − vObj = 0,
ωE (q(T ), q̇(T )) = 0, (4.2.4d)
dq̈,NCP −3 (w) = 0 (4.2.4e)
4.2 Optimal Point-to-Point Trajectory Planning 105
(4.2.4f)
q̇i,min 1 ≤ dq̇,j (w) ≤ q̇i,max 1 ≤ 0, i = 1, . . . , n, j = 2, . . . , NCP − 2
(4.2.4g)
q̈i,min 1 ≤ dq̈,j (w) ≤ q̈i,max 1 ≤ 0, i = 1, . . . , n, j = 1, . . . , NCP − 3
(4.2.4h)
⎛ ⎞ ⎛ ⎞
x
⎝ min ⎠
x
⎝ max ⎠
≤ rE (q(t)) ≤ (4.2.4i)
ymin ymax
Results
In the above Figures 4.6 and 4.7 it can be seen that the velocity and
torque constraints are active during most of the trajectory except for
the initial acceleration and terminal deceleration phases in which the
prescribed boundary points are dominant.
4.2 Optimal Point-to-Point Trajectory Planning 107
q (normalized) q1 q2 q3 q4
1
0.5
0
−0.5
t in s
−1
0 5· 10−2 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Note that the terminal motor torques are not zero despite the terminal
accelerations are restricted to q̈(T ) = 0 due to the non-zero terminal
velocity and the viscous friction effects included in the model for the
system dynamics.
Table 4.1 shows the optimal terminal times T and computation times
tCPU for all considered cases Δx. The selected acceleration bounds
appear to be too low, also investigation of the motor torques resulting
from the trajectories reveals that the torque constraints are violated
during multiple times which supports this indication. However, proper
choice of acceleration bounds as a replacement for torque bounds is
complex due to the aforementioned reasons.
q (normalized) q1 q2 q3 q4
1
0.5
0
−0.5
t in s
−1
0 5· 10−2 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65
MMot (normalized) M1 M2 M3 M4
1
0.5
0
−0.5
t in s
−1
0 5· 10−2 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65
Figure 4.7: SCARA4 Minimum-Time Motion (Constrained Motor Torques MMot (t)),
Δx = 0. Limits: qi,max : 1.65 rad, q̇1,max = q̇2,max : 2.51 rad/s,
q̇3,max = q̇4,max : 5.81 rad/s, M1,max = M2,max : 0.336 Nm,
M3,max = M4,max : 0.142 Nm.
The next example will be a variation of the previous, now minimizing the
joint velocities (and thus the energy dissipated via viscous friction).
Problem Statement
Consider the previous Example 4.3 with the cost function rewritten
such that the minimum-time trajectory planning task is converted into
a minimum-(pseudo-)energy task, i.e. (4.2.4a) becomes
4
minimize
w
JE = di d q̇i dq̇i .
i=1
110 4 Path and Trajectory Planning
Discussion
Figure 4.8 shows a typical scenario in which the robot’s EE is synchro-
nized with an object on the conveyor belt. Note, that in this example
the farmost path in the x direction is limited by the path constraints
(4.2.4i). Furthermore, the minimum-energy paths tend to be longer
than the minimum-time ones due to the significantly higher execution
times which can be seen in Table 4.2. Again, paths further away from
the robot base yield higher energy dissipation due to decreasing execu-
tion times. As the rest of the optimization problems is the same, the
trajectories for both kinds of constraints are the equal.
x
object detection line
vB
terminal
Δx configuration q(T )
y
initial
configuration q0
Δx JE T tCPU
1.00 14.56 3.307 13.83
0.95 13.05 3.431 18.29
0.90 11.94 3.526 31.98
0.85 11.10 3.600 18.30
0.80 10.46 3.656 17.95
0.75 9.96 3.697 18.02
0.70 9.60 3.726 23.93
0.65 9.33 3.746 28.68
0.60 9.16 3.758 31.72
For this example, only the torque-constrained plots are shown in Fig-
ure 4.9, as neither the (artificial) acceleration limits, nor the motor
torque limits are reached due to the low speed of the operation. Again,
note the sudden drops and increases of the motor torques in Figure 4.9.
These are a results of stiction effects included in the dynamics model of
the robot.
q (normalized) q1 q2 q3 q4
1
0.5
0
−0.5
t in s
−1
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
q̇ (normalized) q̇1 q̇2 q̇3 q̇4
1
0.5
0
−0.5
t in s
−1
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
112 4 Path and Trajectory Planning
MMot (normalized) M1 M2 M3 M4
1
0.5
0
−0.5
t in s
−1
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
Problem Statement
In this example that was originally published in the author’s work
[116], a rendezvous problem of two identical Comau Racer3 industrial
robots A and B is considered, cf. Figure 4.10. In contrast to robot A,
robot B is mounted on a linear axis and has therefore an increased work
space. Robot B is used to simulate the synchronization target, i.e. it
remains in the depicted pose during the operation, only the linear axis
moves with constant velocity.
vy + Δ
Δx
B
A
Thus, linear work space motion of the EE does not constitute a special
case.
Implementation
In the collocation scheme
N
minimize
w
J(w) = minimize
w
q̇ k q̇k + αM k Mk (4.2.5a)
k=0
s.t. 1qi,min ≤ dqi ≤ 1qi,max (4.2.5b)
1q̇i,min ≤ dq̇i ≤ 1q̇i,max (4.2.5c)
... ...
1q i,min ≤ d...qi ≤ 1q i,max ∀i = 1 . . . n (4.2.5d)
Mk,min ≤ Mk ≤ Mk,max ∀k = 0 . . . N (4.2.5e)
q0 = q0 , q̇0 = q̇0 , q̈0 = q̈0 (4.2.5f)
qN = pqN , q̇N = pq̇N , q̈N = q̈N (4.2.5g)
114 4 Path and Trajectory Planning
While small perturbations have no effect on the active set, for larger
perturbations the active set of the optimal solution to the perturbed
problem may be different from that of the nominal problem.
4.2 Optimal Point-to-Point Trajectory Planning 115
Results
Table 4.3: Comau Racer R3 Trajectory Planning: Results (perturbations Δx, relative
cost increase JΔJ
pert
, number of iterations).
Collocation
Δx in m ΔJ/Jpert # iter
−0.04 +5.65% 4
−0.02 +2.6% 3
−0.01 +1.48% 2
+0.01 +0.39% 2
+0.02 +5.13% 3
+0.04 +8.35% 4
Due to the particular choice of the cost function and the terminal time,
the joint velocities q̇ as well as the motor torques QM rarely reach their
respective limits in this example. The limiting factor is found to be the
...
joint jerks q. Table 4.3 shows an overview of the results wherein the
relative cost increase of J w.r.t. the exact cost function value in the
perturbed case Jpert as well as the number of iterations of Algorithm 2
are reported. The method was found to converge for perturbations
up to Δx = ±0.04 m in a low number of iterations. The results were
successfully validated in experiments featuring the described setup, i.e.
the planned trajectories were executed with negligible joint tracking
error.
Evaluation of the computation time required to solve the above NLP
problems reveals that such formulations can hardly be used in practical
applications. However, by revisiting the methods to approximate solutions
of parametric NLPs discussed in Section 2.1.3 and especially the academic
examples considered in Section 2.1.3.4, these trajectory planning examples
seem to form a suitable basis for the optimal solution approximations as
demonstrated in the above examples.
116 4 Path and Trajectory Planning
Both of the above stages can be further decomposed in substages, e.g. the
geometric planning stage may include determining the initial and terminal
joint configuration for a given EE pose. That way, the problem complexity
of path planning is lowered, allowing efficient joint space formulations of
path planning strategies to be applied.
minimize
q
J(q) (4.2.6a)
s.t. zE,d − zE (q) = 0 (4.2.6b)
qmin ≤ q ≤ qmax (4.2.6c)
h(q) ≤ 0 (4.2.6d)
Implementation (Part 1)
The NLP problem (4.2.6) thus changes into
maximize
q
H 2 (q) = abs det J(q)J (q) (4.2.7a)
s.t. rE,d − rE (q) = 0 (4.2.7b)
qmin ≤ q ≤ qmax (4.2.7c)
vB
WR
0y
Figure 4.11: SCARA4 Path Planning: Optimal Picking Poses for the SCARA4
Manipulator. Red: Nominal Solutions. Green: Admissible Approximate Solutions.
Results (Part 1)
The above NLP problem (4.2.7) is solved for all points on the uniform
grid shown in Figure 4.11 to determine the respective optimal joint
configuration. As the initial guess, all angles were set to zero (joints
centered). Of course, only points within the manipulator’s reachable
work space WR , cf. Section 3.3.2, will yield a corresponding joint con-
4.2 Optimal Point-to-Point Trajectory Planning 119
Implementation (Part 2)
NLP (4.2.7) is thus rewritten as
maximize
q
J(q) = H 2 (q) = abs det J(q)J (q) (4.2.8a)
s.t. g1 (q, p) = p − rE (q) = 0 (4.2.8b)
g2 (q) = qmin − q ≤ 0 (4.2.8c)
g3 (q) = q − qmax ≤ 0 (4.2.8d)
Results (Part 2)
In order to allow for a comprehensive comparison of the methods dis-
cussed in Section 2.1.3, each of the solution points from the first part
120 4 Path and Trajectory Planning
Next, the actual path planning process for the synchronization task is
considered.
The following examples show how optimal path planning is performed for
pick and place tasks. Therein, not only the often computationally involved
exact solution of such problems but also fast solution approximation using
parametric sensitivities, i.e. the iterative+ method, is considered.
122 4 Path and Trajectory Planning
Results
Figure 4.12 shows the grid of reachable picking positions known from
the previous Example 4.6 with circles for each position for which an
admissible path was found by virtue of (4.2.10). In addition, one
particular optimal EE path is shown in red color. It is used as the
nominal solution to (4.2.10) that serves as the basis for approximations
of optimal paths to other picking position. These are displayed as
blue circles in contrast to the black ones that cannot be reached by
an approximation. It is apparent, that many picking positions can
be reached by an appoximation which is mainly due to the benign
NLP problem (4.2.10), regard therefore the objective and the mostly
linear constraints. The reasons for failures of convergence are twofold:
On the one hand, there is the particular combination of work space
constraints and joint limits, on the other hand some of the terminal joint
configurations determined in the preceding example are not suitable
as some joint angle qi,term may be at its lower or upper bound and the
corresponding joint slope qi,term is non-zero. Although these terminal
configurations are admissible w.r.t. the NLP (4.2.8) of the previous
example, they have to be effectively discarded. Modification of (4.2.8)
4.2 Optimal Point-to-Point Trajectory Planning 125
with narrowed joint angle bounds may mitigate this issue. Alternatively
(or additionally) the objective function may be augmented by a term
that penalizes joint displacement and centers the angles, i.e. (4.2.8a) is
modified to J(q) = H 2 (q) + αq q with the weight α.
0x
vB
WR
0y
Figure 4.12: Admissible Picking Positions (Black) with One Nominal Path (Red) and
All Possible Picking Position Reachable by Path Approximation (Blue).
0x
vB
WR
0y
Figure 4.13: Nominal Path (Red) and All Paths to Possible Picking Position
Reachable by Path Approximation (Green).
For the nominal solution shown in Figures 4.12 and 4.13, the success
rate is 78.8% (152/193 samples) which is the highest of all terminal
configurations. An overview of the computational performance in this
particular case – the average, median and maximum computation times
tCPU,avg , tCPU,med and tCPU,max , respectively, and the average, median
and maximum relative cost increases ΔJrel,avg , ΔJrel,med and ΔJrel,max ,
respectively – can be found in Table 4.4. Most paths can be computed
with no or one iterations, which yields low computation times on average
and also little to no cost increase.
Table 4.4: SCARA4 Path Planning: Computation Time and Cost Increase in
Synchronization Tasks (152 Samples).
For an overview of the total performance, all feasible 37442 samples were
analyzed, returning a success rate of 57.7% (21618/37442 samples). An
overview of the total computational performance is shown in Table 4.5.
4.2 Optimal Point-to-Point Trajectory Planning 127
Increased cost, i.e. curvature, is shown for the paths right to the initial
position while increase computation times were found for the longer
paths ending on the far left of the conveyor.
Table 4.5: SCARA4 Path Planning: Computation Time and Cost Increase in
Synchronization Tasks (21618 Samples).
sync2
sync1
Problem Statement
One subtask therein is to perform optimal path planning under geomet-
ric constraints such that the joint path curvature is minimized. The
joint evolutions are expressed by means of B-spline curves that are pa-
rameterized using a path parameter s ∈ [0, 1], i.e. q = q(s). Successful
synchronization of a manipulator’s EE with a target imposes require-
128 4 Path and Trajectory Planning
ments on the initial and terminal joint configuration, q(0) and q(1),
d d
on the path slope ds q(0) and ds q(1) as well as on the path curvature
d2 d2
ds2 q(0) and ds 2 q(1), respectively. Note, that in [65] the curvature
d2
ds2 q(s)
κ(s) = 2 3/2
d
1+ ds q(s)
2
of a single curve q(s) is approximated with κ(s) ≈ ds
d
2 q(s) which is also
Further constraints are imposed on the joint angles due to the joints’
physical turn limits. For collision avoidance, the EE motion is required
to be confined within an admissible region such that no collision with
the second robot occurs beyond xmax = 0.6 m as well as the floor occurs
at zmin = 0.3 m.
In the cited work, in order to reduce the computation time, the path
planning problem is solved in its unconstrained form as only a QP re-
mains for which a closed-form exists. Admissibility w.r.t. the constraints
is then restored with a two-stage iterative process. Firstly, admissibility
w.r.t. the joint limits is repaired by iteratively scaling initial and final
conditions. Then, the work space limits are considered by iteratively
adjusting the curves’ control points.
Proposed Solution
By revisiting the methods to approximate solutions of parametric NLPs
discussed in Section 2.1.3, this path planning example seems to be
efficiently solvable with the improved iterative approximation approach
from Algorithm 2.
which are set according to the synchronization task at hand. The joint
angle limits are incorporated as (4.2.11h) and (4.2.11i). Note, that only
the intermediate control points dqi ,3 , . . . , dqi ,NCP −4 are not influenced by
the boundary conditions which specify the first and last 3 control points.
The robot’s forward kinematics functions are denoted fx and fz for the
inertial x and z direction in the work space constraints (4.2.11j) and
(4.2.11k), respectively. Only the NCP − d − 5 knot spans between td+3
and tNCP −3 are not impacted by the boundary conditions and thus need
to be checked which is facilitated by means of uniformly distributed
check points. Due to this careful selection of constraints, the NLP
problem is kept as small as possible. Furthermore, LICQ violations
between inequality and equality constraints are avoided but they may
arise within the work space limits if both constraints become active in
certain poses of the manipulator. However, these cases are practically
irrelevant as the EE is unlikely to assume a position along this boundary
line as the target path is far away.
Results
In order to allow for comprehensive evaluation of this example, both,
the initial and the target path of the considered manipulators, are
discretized in 32 points each and all resulting 1024 combinations of
4.2 Optimal Point-to-Point Trajectory Planning 131
Table 4.6 gives an overview of the success rate for both paths. Not all
tasks, only 768 out of 1024, were feasible due to the imposed joint angle
limits (4.2.11h) and (4.2.11i). However, it can be seen that nearly all
of the feasible synchronization tasks can be fulfilled using approximate
solutions, i.e. with the iterative+ approach. Note, that admissibility
w.r.t. the constraints of the optimization problem is a property inherent
to this method.
Tables 4.7 and 4.8 show details of the computation time (average tCPU,avg ,
median tCPU,med , maximum tCPU,max ) of the exact solution of the path
planning tasks, as well as the relative cost increase (average ΔJrel,avg ,
median ΔJrel,med , maximum ΔJrel,max ) when using the improved iterative
approach with one particular set of initial and terminal configuration
used as the nominal solution. The iterative+ approach is significantly
faster than the exact solution and yields admissible, but suboptimal
results. More detailed analysis shows a correlation between computation
132 4 Path and Trajectory Planning
time and relative cost increase, i.e. the longer the methods typically
takes to converge, the more expensive the resulting solution is. Note,
that the high cost increase in this example is mainly due to the broad
variety of requested synchronization paths that are mostly unresembling
to the nominal solution, e.g. consider a terminal desired path slope that
is directly opposed to that of the nominal solution.
Table 4.7: sync1 : Computation Time and Cost Increase in Synchronization Tasks
(768 Samples).
Table 4.8: sync2 : Computation time and cost increase in synchronization tasks (768
samples).
0.9 0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
z
z
0.4 0.4
0.3 0.3
0.2 0.2
0.8 0.8
0.6 0.4 0.6
0.4
0.2 0.6 0.2
0 0 0.6
y −0.2 0.2
0.4
y −0.2 0.4
x 0.2 x
0.8 1
0.7 0.9
0.6 0.8
0.7
0.5
0.6
z
0.4
0.5
0.3
0.4
0.2 0.3
0 0.2
0.6 0.8 0.6
0.4
0.2 0.4
0 0.6 0.2 0.6
y −0.2 0.4 0
0.2 x y −0.2 0.2
0.4
x
Figure 4.15: Synchronization Tasks of Two Comau Racer R3, sync1 : Nominal Path
(Red) and Some Approximations of Optimal Paths (Blue). Measures in m.
134 4 Path and Trajectory Planning
1 1
0.9 0.9
0.8 0.8
0.7 0.7
z
z
0.6 0.6
0.5 0.5
0.4 0.4
0.3 0.3
1 0.8 1 0.8
0.6 0.6
0.4 0.6 0.4 0.6
0.2 0.2
0 0.2 0.4 0 0.2 0.4
y x y x
1 1
0.9 0.9
0.8 0.8
0.7 0.7
z
0.6 0.6
0.5 0.5
0.4 0.4
0.3 0.3
1 0.8 1 0.8
0.6 0.6
0.4 0.6 0.4 0.6
0.2 0.2
0 0.2 0.4 0 0.2 0.4
y x y x
Figure 4.16: Synchronization Tasks of Two Comau Racer R3, sync2 : Nominal Path
(Red) and Some Approximations of Optimal Paths (Blue). Measures in m.
4.2.3 Summary
Summarizing, an OCP such as (4.2.1) includes both, geometric and dynamics
planning subtasks, it represents a trade-off between the quality of results
and its use for real-time applications. Note, that such a staggered strategy,
although computationally efficient, may yield different, possibly degraded
results in terms of the optimization goal compared to those of comprehensive
single-step methods. This is an effect of the optimization goal of the
kinematic part often being different than that the originally desired one
as no time-dependent goals (such as minimum traversal time) or dynamic
objectives can be pursued. The examples in this chapter have shown, that
the iterative+ method based on parametric sensitivities is especially suitable
for fast path and trajectory re-planning.
Chapter 5
In the former case, the sequence of a robot’s joint positions is given and a
corresponding temporal interpolation is required to be found such that a
given criteria is optimal. Kinematic redundancy is irrelevant for joint space
path tracking approaches.
The path tracking problem in joint space is not in the scope of this thesis.
In [82, 88] the path tracking problem for kinematically redundant manipula-
tors is formulated in terms of a two point boundary value problem. Therein,
a PG method locally maximizing the kinematic manipulability is employed
that allows to reduce the dimension of the problem which is then solved
with calculus of variations. Further approaches using PG methods in com-
bination with a directed manipulability index are found in [30, 28, 142]. As
discussed in Section 3.5.3, it is difficult to select an appropriate performance
index that locally corresponds to the global trajectory optimization goal.
JSD was previously applied for time-optimal path tracking in [50, 135].
However, these phase-plane-based methods exhibit severe shortcomings as.
For example, the terminal joint velocities at the end of a path cannot be
prescribed, i.e. internal motion at the terminal point cannot be avoided.
Furthermore, the approaches were only applied to manipulators with degree
of redundancy r = 1 which does not confirm their general applicability.
Note, that the time grid is being scaled when T changes in w. The
choice of the IK resolution method determines the state vector x
and the control vector u, e.g. in the case of C 4 trajectories x k =
...
(sk , ṡk , s̈k , sk , sk , q k , q̇ r,k , q̈ r,k , q r,k , qr,k ), u k = (sk , qr,k ) for JSD, and
(3) (4) (4) (5) (5)
...
x k = (sk , ṡk , s̈k , sk , sk , q k , γ k , γ̇ k , γ̈ k , γ k ), u k = (sk , γk ) for GNA.
(3) (4) (5) (4)
Therein, the IK expressions (3.5.34) and (3.5.36) for JSD and (3.5.47)
and (3.5.51) for GNA are used. The control trajectories are piecewise
constant. Therefore, the NLP problem corresponding to minimum-time
path tracking reads
minimize
w
T (5.2.2a)
s.t. s0 = 0, sN = 1 (5.2.2b)
(i) (i) (i) (i)
q0 = q0 , qN = qf , i = 0, . . . , ν (5.2.2c)
xk+1 − RK4(xk , uk ) = 0 (5.2.2d)
IK: JSD or GNA (5.2.2e)
0 ≤ sk ≤ 1 (5.2.2f)
ṡk ≥ 0 (5.2.2g)
(i) (i)
qmin ≤ qk ≤ qmax
(i)
, i = 0, . . . , ν + 1 (5.2.2h)
τk = M(qk )q̈k + g(qk , q̇k ) (5.2.2i)
τmin ≤ τk ≤ τmax (5.2.2j)
The initial and final conditions regarding the path parameter and the
joint coordinates are given by (5.2.2b) and (5.2.2c). Therein, the overlined
quantities, e.g. qf , represent desired values at the respective bounds.
Choosing an IK method used for path tracking, i.e. JSD or GNA, is reflected
in (5.2.2e).
(5.2.2f) and (5.2.2g) constrain the time evolution of the path parameter
s(t). (5.2.2h) represents box constraints imposed on the joint coordinates
and its time derivatives denoted with their respective orders i = 0 for
q(0) = q(t), i = 1 for q(1) = q̇(t), etc. The system dynamics is considered
as the robot’s equations of motion in inverse dynamics form (5.2.2i) and
the torque constraint (5.2.2j).
Initial Guess
The above NLP problem (5.2.2) is in general non-convex and nonlinear
due to the system dynamics and the effect modeled therein. Therefore, the
selection of an initial guess is important, i.e. the closer it is to a solution
to the problem, the lower the computational load. Admissibility w.r.t. the
constraints also improves the quality of an initial guess.
A possible choice for the initial guess is to compute the joint angles by time
integration of the non-augmented IK solution
(4) (3) (3) (3)
k (VE,d,k − Jk q̇k − 3J̈k q̈k − 3J̇k qk ).
qinit,k = J+
Thereby, the path tracking requirement is fulfilled, and also the box con-
straints regarding the joint quantities and the motor torques (5.2.2h)
to (5.2.2j) are satisfied when initialized with a large enough terminal time
T . This also yields the initial guess for the redundant coordinates qr and
their time derivatives for JSD. For GNA, the scaling factors γ are initialized
with zeros.
5.2 Work Space Approaches 143
Task Description
In this example, which was published in the author’s previous publica-
tions [114, 115], the kinematically redundant SCARA4 lab prototype
from Figure 5.1 is considered. The goal is to compute joint trajectories
q(t) with a continuity requirement up to C 4 such that the EE performs
a minimum-time rest-to-rest motion, tracking the quarter circle path
shown in Figure 5.1 in counterclockwise direction. Therein, JSD as well
as the GNA approach will be used to formulate the IK as a means to
satisfy the path following requirement.
path rE (s)
initial configuration q0
Optimization Problem
The formulation of the optimization problem reflecting the problem
described above follows the direct multiple shooting approach (5.2.2)
using N = 200 uniform time intervals. Therein, the boundary conditions
for the path parameter trajectory are specified in (5.2.2b). As a rest-
to-rest motion is required, the initial and final constraints (5.2.2c) are
(3) (4) (3) (4)
specified as q0 = q0 and q̇0 , q̈0 , q0 , q0 , q̇N , q̈N , qN , qN = 0. The
final joint configuration qN is subject to the optimization. In (5.2.2h),
q and its time derivatives are constrained except for the acceleration
which is indirectly limited by the torque constraints (5.2.2i) and (5.2.2j).
Regarding the IK resolution, for JSD all six possible decompositions of
the joint set are tested, as well as GNA.
Initial Guess
For the present example, a B-spline curve of degree 6 is used for s, the
duration of the trajectory is conservatively initialized with tf,init = 3 s.
Results
The optimization problem (5.2.2) was implemented using CasADi [2]
and solved with Ipopt.
Table 5.1 shows the trajectory planning results for the above example.
Therein, optimal trajectory durations for all joint space decompositions
as well as for GNA are reported. The numbers denote the choice of
redundant joints, i.e. JSD (1,2) means q r = (q1 , q2 ). Additionally,
the corresponding computation times tCPU are listed. It should be
mentioned that the resulting minima may not be global as (5.2.2) is
non-convex.
5.2 Work Space Approaches 145
Figures 5.2 to 5.4 show the optimal C 4 time evolutions of the path
parameter s, the joint angles q and the motor torques τ , respectively.
Note that the axes of all joint quantities are normalized w.r.t. their
respective limits and the time axes are normalized w.r.t. T . The
visible steps in the motor torques in Figure 5.4 are due to Coulomb
friction effects, i.e. they occur at sign changes of the corresponding
joint velocities, c.f. Figure 5.3. The optimal final configuration of the
robot in the C 4 case is depicted in Figure 5.1. The time evolutions are
only reported for GNA as the results for JSD are very similar which is
mainly due to the common initial guess.
q/qmax q1 q2 q3 q4
1
0.5
0
−0.5 t/tf
−1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
q̇/q̇max
1
0.5
0
−0.5 t/tf
−1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
(3)
q(3) /qmax
1
0.5
0
−0.5 t/tf
−1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
(3)
q(3) /qmax
1
0.5
0
−0.5 t/tf
−1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
(5)
q(5) /qmax
1
0.5
0
−0.5 t/tf
−1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
(4) (5)
103 rad/s3 , qi,max = 104 rad/s4 , qi,max = 105 rad/s5 .
τ /τmax τ1 τ2 τ3 τ4
1
0.5
0
−0.5 t/tf
−1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
5.2 Work Space Approaches 147
Problem Statement
Consider the redundant 7 DOF manipulator in Figure 5.5. It will serve
as another use case of higher order optimal path tracking. This example
has been largely adopted from the author’s previous publications [114,
115].
tool axis
qlag,1 q1 q2 q3 qlag,1
q4 q5 q6 q7 qlag,2...7
2 200
1 100
0 0
−1 −100
t in s
0 5 · 10−2 0.1 0.15 0.2 0.25 0.3
τ /τmax τ1 τ2 τ3 τ4 τ5 τ6 τ7
1
0
−1 t in s
0 5 · 10−2 0.1 0.15 0.2 0.25 0.3
The problem becomes obvious from Figure 5.6 that shows the tracking
error qlag = qd −q with qd being the desired joint position. More crucial
than this low tracking accuracy is the violation of the actual motor
torque limits as apparent from Figure 5.7 that causes an emergency
stop of the robot at around 0.07 s. The latter is due to the bounded
but discontinuous jerk and the unbounded jounce of the trajectories.
Task Description
The goal is to determine joint trajectories such that a prescribed elliptic
EE path is tracked in a minimum time rest-to-rest motion while the last
axis of the spherical wrist (tool axis) is required to stay perpendicular
to the ellipse (the EE is free to rotate about this axis). The path is
defined in closed-form in terms of a path parameter s ∈ [0, 1]. The joint
trajectories are required to be continuously differentiable up to four times
150 5 Optimal Path Tracking
Optimization Problem
See the previous Example 5.1.
Initial Guess
The initial guess is computed as in the previous Example 5.1 with the
terminal time conservatively initialized to Tinit = 5 s.
Results
The optimization results for the spatial example are found in Table 5.2.
JSD was performed for all 21 possible combinations of redundant coordi-
nates. Table 5.2 only shows the results for separations with the largest
(worst) and the smallest (best) final time. Not all choices converged,
some encountered algorithmic singularities. The corresponding occur-
rences as well as all computation times tCPU are also listed. Again, note
that the resulting minima may not be global as (5.2.2) is non-convex.
The NLP problem (5.2.2) was implemented using CasADi 3.1 [2] avail-
able at the time of the initial submission of [115] and solved with Ipopt.
Computation times for GNA are higher than for JSD as the evalua-
5.2 Work Space Approaches 151
The C 2 and C 3 trajectories obtained by JSD and GNA are very similar
which is mainly due to the common initial guess. Figures 5.8 to 5.10
show the optimal C 4 time evolutions of the path parameter s, the joint
angles q and the motor torques τ , respectively, that were obtained with
JSD. Note that the axes of all joint quantities are normalized w.r.t. their
respective limits as are the time axes w.r.t. T . The visible steps in the
motor torques in Figure 5.10 are again due to Coulomb friction effects,
i.e. they occur at sign changes of the corresponding joint velocities, c.f.
Figure 5.9.
Figure 5.8: Stäubli TX90L Minimum-Time Path Tracking: Path Parameter s and
Derivatives.
Figure 5.9: Stäubli TX90L Minimum-Time Path Tracking: Joint Angles q and
Derivatives Normalized w.r.t. Their Respective Limits
(q1,min = 0, q2,...,7,min = (−π, −2.3, −2.5, −4.7, −2, −3.5) rad,
q1,max = 2.55 m, q2,...,7,max = (π, 2.6, 2.5, 4.7, 2.4, 3.5) rad. Symmetric derivative limits:
q̇max = (4, 7, 7, 7, 8.7, 7.9, 12.6) rad/s, qi,max = 103 rad/s3 , qi,max = 104 rad/s4 , qi,max =
(3) (4) (5)
5 5
10 rad/s .).
τ /τmax τ1 τ2 τ3 τ4 τ5 τ6 τ7
1
0.5
0
−0.5 t/tf
−1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
qlag,1...2 q1 q2 qqlag,1
3 q4 qlag,2
q5 q6 q7 qlag,3...7
1000 100
0 0
−1000 −100
t/tf
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Figure 5.11: Stäubli TX90L Minimum-Time Path Tracking: Tracking Error of Joint
Positions qlag , qlag,1 in μm, qlag,2...7 in mdeg.
154 5 Optimal Path Tracking
For GNA, due to the problem size in the C 4 case, cf. the IK expressions
(3.5.51), a limited-memory quasi-Newton approximation of the Hessian
matrix was used instead of exact derivatives. Convergence was achieved
for relaxed optimality tolerances.
Chapter 6
Conclusion
As its title suggests, this thesis is concerned with optimal path and trajec-
tory planning for serial manipulators. For the treatment of this subject, the
mathematical background regarding numerical optimization and therein es-
pecially parametric problems, B-spline curves and optimal control problems
is required which is provided in this thesis. In addition, robot kinematics, i.e.
forward and inverse kinematics on position, velocity and higher derivative
level, is also examined rigorously in this thesis as a preparation for more
advanced topics.
On the one hand, point-to-point tasks are treated such that trajectories
are computed that are optimal in some given sense, e.g. so that a robot
consumes minimum energy during a transportation task.
work will investigate whether the methods described in [108] can improve
the convergence of the presented active set update approach.
to compute a solution grows with the derivative order. This can be possible
overcome by a recursive IK formulation. Alternatively, this task can be
formalized with a modified OCP, cf. Section 2.3 based on the system’s
forward dynamics. Note, that this approach requires time derivatives of the
EoM which is computationally involved as well. Comparative investigations
about the computational trade-off between an OCP formulation and an
IK-based approach as proposed in this thesis are an open challenge.
Dynamic programming can also be used to express and solve such path
tracking problems. Such approaches are successfully applied for path track-
ing with non-redundant robots as the dynamics can be uniquely re-written
in terms of a path parameter, reducing the problem to finding a single
trajectory. For kinematically redundant structures this is not straightfor-
wardly possible, as the IK resolution is ambiguous and only a parameterized
solution can be obtained. In both, the JSD and the GNA methods, the
number of states grows linearly with the degree of kinematic redundancy.
However, in dynamics programming, the computational complexity grows
exponentially with the number of states which may limit this idea’s practical
feasibility, especially regarding higher derivative levels.
[41] Rajiv V Dubey, James A Euler, and Scott M Babcock. “An effi-
cient gradient projection optimization scheme for a seven-degree-
of-freedom redundant robot with spherical wrist”. In: Robotics and
Automation, 1988. Proceedings., 1988 IEEE International Confer-
ence on. IEEE. 1988, pp. 28–36.
[42] Olav Egeland. “Task-space tracking with redundant manipulators”.
In: IEEE Journal on Robotics and Automation 3.5 (1987), pp. 471–
475.
[43] Francisco Facchinei, Andreas Fischer, and Christian Kanzow. “On
the accurate identification of active constraints”. In: SIAM Journal
on Optimization 9.1 (1998), pp. 14–32.
[44] H.J. Ferreau et al. “qpOASES: A parametric active-set algorithm
for quadratic programming”. In: Mathematical Programming Com-
putation 6.4 (2014), pp. 327–363.
[45] Enrico Ferrentino and Pasquale Chiacchio. “Redundancy
Parametrization in Globally-Optimal Inverse Kinematics”.
In: International Symposium on Advances in Robot Kinematics.
Springer. 2018, pp. 47–55.
[46] Anthony V Fiacco. Introduction to sensitivity and stability analysis
in nonlinear programming. Elsevier, 1983.
[47] Anthony V Fiacco. “Sensitivity analysis for nonlinear programming
using penalty methods”. In: Mathematical programming 10.1 (1976),
pp. 287–311.
[48] Angel Flores-Abad et al. “Optimal Capture of a Tumbling Object
in Orbit Using a Space Manipulator”. In: Journal of Intelligent &
Robotic Systems 86.2 (2017), pp. 199–211.
[49] M Galassi et al. GNU Scientific Library Reference Manual , ISBN
0954612078. GSL. url: https://www.gnu.org/software/gsl/.
[50] Miroslaw Galicki. “Time-optimal controls of kinematically redundant
manipulators with geometric constraints”. In: IEEE Transactions
on Robotics and Automation 16.1 (2000), pp. 89–93.
164 Bibliography
[71] Vyacheslav Kungurtsev and Moritz Diehl. “SQP Methods for Para-
metric Nonlinear Optimization”. In: Optimization Online (2014),
pp. 1–35.
[72] Huibert Kwakernaak and Raphael Sivan. Linear optimal control
systems. Vol. 1. Wiley-Interscience New York, 1972.
[73] Roberto Lampariello et al. “Trajectory planning for optimal robot
catching in real-time”. In: Robotics and Automation (ICRA), 2011
IEEE International Conference on. IEEE. 2011, pp. 3719–3726.
[74] J Leger and J Angeles. “A Redundancy-resolution Algorithm for Five-
degree-of-freedom Tasks via Sequential Quadratic Programming”. In:
TrC-IFToMM Symposium on Theory of Machines and Mechanisms„
Izmir. 2015.
[75] Daniel B Leineweber. “Efficient reduced SQP methods for the op-
timization of chemical processes described by large sparse DAE
models”. PhD thesis. VDI Verlag Düsseldorf, 1999.
[76] Renhui Li. “Strukturierte Modellbildung und Implementierung von
mechanischen Systemen”. MA thesis. Johannes Kepler University,
2009.
[77] Zhi Li et al. “Synthesizing redundancy resolution criteria of the
human arm posture in reaching movements”. In: Redundancy in
Robot Manipulators and Multi-robot systems. Springer, 2013, pp. 201–
240.
[78] Alain Liégeois. “Automatic Supervisory Control of the Configuration
and Behavior of Multibody Mechanisms”. In: IEEE Transactions on
Systems Man and Cybernetics 12 (1977), pp. 868–871.
[79] J. Löfberg. “YALMIP : A Toolbox for Modeling and Optimization
in MATLAB”. In: In Proceedings of the CACSD Conference. Taipei,
Taiwan, 2004.
[80] J Luh, M Walker, and R Paul. “Resolved-acceleration control of
mechanical manipulators”. In: IEEE Transactions on Automatic
Control 25.3 (1980), pp. 468–474.
Bibliography 167
[100] Martin Pfurner. “Closed form inverse kinematics solution for a redun-
dant anthropomorphic robot arm”. In: Computer Aided Geometric
Design 47 (2016), pp. 163–171.
[101] Hung Pham and Quang-Cuong Pham. “A new approach to time-
optimal path parameterization based on reachability analysis”. In:
IEEE Transactions on Robotics (2018).
[102] Quang-Cuong Pham. “A general, fast, and robust implementation
of the time-optimal path parameterization algorithm”. In: IEEE
Transactions on Robotics 30.6 (2014), pp. 1533–1540.
[103] Les Piegl and Wayne Tiller. The NURBS book. Springer Science &
Business Media, 2012.
[104] Donald Lee Pieper. “The kinematics of manipulators under computer
control”. PhD thesis. Stanford Univ., Dept. Of Computer Science,
1968.
[105] François G Pin et al. “Motion planning for mobile manipulators
with a non-holonomic constraint using the FSP (full space param-
eterization) method”. In: Journal of robotic systems 13.11 (1996),
pp. 723–736.
[106] Hans Pirnay, Rodrigo López-Negrete, and Lorenz T Biegler. “Opti-
mal sensitivity based on IPOPT”. In: Mathematical Programming
Computation 4.4 (2012), pp. 307–331.
[107] Ron P Podhorodeski, Andrew A Goldenberg, and Robert G Fenton.
“Resolving redundant manipulator joint rates and identifying spe-
cial arm configurations using Jacobian null-space bases”. In: IEEE
transactions on robotics and automation 7.5 (1991), pp. 607–618.
[108] Daniel Ralph and Stephan Dempe. “Directional derivatives of the
solution of a parametric nonlinear program”. In: Mathematical pro-
gramming 70.1-3 (1995), pp. 159–172.
[109] Alexander Reiter. “An explicit approach for time-optimal trajectory
planning for kinematically redundant robots”. MA thesis. Johannes
Kepler University Linz, Institute of Robotics, Dec. 2014.
170 Bibliography
B-Spline Curves
In Section 2.2, the numerical basics regarding B-spline curves have been
reviewed. Therein, not only the simple evaluation, but also the optimization-
based construction of such functions have been discussed. In this appendix
chapter, although they are not required for the topics in the main part,
additional features of B-splines are documented that find application in a
software tool that was developed as a by-product of the work leading to
the present thesis.
The total length of the spline curve is thus found by setting t = b in (A.1.1).
For the parameterization t(l), the inverse of the above function needs to
be found which is non-trivial. However, an approximation t̂(l), can be
calculated as shown in the following Algorithm 4.
Algorithm 4 Computation of arc length parameterization as B-spline
curve.
q(t) ∈ Rn : B-spline curve parameterized in terms of parameter t
τ ← t 1 t 2 . . . t n : collect knots of all B-spline curves
remove duplicates in τ τ : vector of unique knots
nτ ← length of τ
l ← l l(τ1 ) l(τ2 ) . . . l(τnτ ) : vector of arc lengths at knots
τ ← ( ): vector of intermediate grid points
l ← ( ): vector of arc lengths at intermediate points
for all knot intervals [τk , τk+1 ] ∈ τ , k = 1, . . . , nτ − 1 do
compute N intermediate points τk < τk,1 < τk,2 < · · · < τk,N < τk+1
τ ← τ τk,1 τk,2 . . . τk,N : collect intermediate grid points
l ← l l(τk,1 ) l(τk,2 ) . . . l(τk,N ) : compute arc lengths at interme-
diate grid points
end for
return τ , l, nτ , τ , l
A.1 Arc Length Parameterization 177
Therein, first the knots of all splines are collected in the vector of unique
knots τ of length nτ that has only unique, monotonically increasing el-
ements. The resulting function t̂(l) will interpolate the arc lengths the
knots contained in τ . Therefore, the arc length value at these points are
computed and stored in the vector l of the same length nτ . Then, on each
interval of τ , N intermediate points are computed and stored in the vector
τ of length N (nτ − 1). The corresponding arc lengths are stored in the
vector l of the same length. These will be used to approximate the evolution
of the arc length between the knots.
The arc length parameterization as a B-spline curve t̂(l) with degree p, NCP
control points (which are collected in vector dt̂ ) and open-uniform knot
distribution is finally found by solving the constrained approximation QP
problem
τ −1)
1 N (n 2
minimize τ i − t̂ li (A.1.2a)
dt̂ 2 i=1
s.t. t̂(li ) = τi , ∀ li ∈ l, ∀ τi ∈ τ (A.1.2b)
that is a special form of (2.2.11) wherein τi , li , τ i and li are the i-th elements
of the vectors τ , l, τ and l, respectively.
The degree p should be chosen as the highest degree (or higher) that occurs
in any element of the original curve q(t) in order to preserve the original
degree of continuous differentiability. The number of intermediate points N
and the number of control points dt̂ need to be chosen in order to sufficiently
capture the arc length evolution.
178 A B-Spline Curves
Overview of functionality
Mechatronic Multibody
Systems
In this chapter, the essential subject of computational modeling of mecha-
tronic systems such as robots as well as the solution of the resulting
differential-algebraic equations (DAE) is discussed.
B.1.1 Constraints
A system such as (B.1.1) can be subjected to kinematic equality constraints
on position level, i.e. geometric constraints, or on velocity level. Note, that
velocity-level constraints exist that do not necessarily have a correspondence
on position-level, i.e. non-holonomic constraints.
g(q) = 0 (B.1.2)
d d
ġ(q, q̇) = g(q) = g(q) q̇ = 0 (B.1.3)
dt dq
Jg (q)
d
ġ(q, ṡ) = ġ(q, ṡ) ṡ = 0 (B.1.4)
dṡ
Jġ (q)
d
(A(q)ṡ) = Ȧ(q, q̇)ṡ + A(q)s̈ = 0, (B.1.6)
dt
(B.1.5) can be re-written as
⎡ ⎤⎛ ⎞ ⎛ ⎞
⎣
M(q) −A(q) ⎦⎝ s̈ ⎠
Q(q, u) − h(q, ṡ)
=⎝ ⎠ (B.1.7a)
A(q) O λ −Ȧ(q, q̇)ṡ
subjected to q̇ = H(q)ṡ. (B.1.7b)
The goal of the next step is to solve (B.1.7) for the accelerations s̈. These of
course can be found by simply solving the above (n + m) × (n + m) system
of linear equations.
Now that the generalized constraint forces λ are known, (B.1.5a) is solved
for the accelerations s̈, i.e.
Note that in (B.1.8) and (B.1.9), three systems of linear equations with the
matrix M = M > 0 need to be solved. Therein, the opportunity to further
reduce the computational burden by applying a re-usable decomposition of
M arises, e.g. Cholesky factorization.
Note, that conditions may also depend on the constraint forces λ, as they
are an algebraic function of the state x and the input u, cf. (B.1.8).
188 B Mechatronic Multibody Systems
While this tool does not take the elaborate task of creating a dynamics model
off of the user, it greatly the reduces burden of efficiently implementing it
for use in numeric simulation and control.
Contact Problems
Contact problems are of great relevance for the mechanical part of many
mechatronic systems, e.g. walking machines, or industrial robots interacting
with their environment. Currently, only a very simple form of contact
problems can be solved, i.e. contact that is toggled by geometric but
not dynamic constraints. Further development will be focused on contact
problems that not only allow conditions depending on accelerations but
also on constraint forces. Part of this is described in [53] wherein the
contact problem is posed as a Linear Complementarity Problem (LCP).
LCPs are most efficiently solved by special solvers, e.g. [130], or by virtue
of re-writing as a QP problem by any QP solver, e.g. [44].