Jtam.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

46, 2, pp. 383-394, Warsaw 2008

MODELING OF UNDERACTUATED MECHANICAL

SYSTEMS IN PARTLY SPECIFIED MOTION

Wojciech Blajer
Krzysztof Kołodziejczyk

Instytute of Applied Mechanics, Technical University of Radom, Poland

e-mail: w.blajer@pr.radom.pl

This paper deals with the inverse dynamics problem of underactuated
mechanical systems subjected to execute partly specified motions. The
modeling methodology focuses on a special class of differentially flat
systems, represented by a group of relevant technical examples. The go-
verning equations are obtained as index three differential-algebraic equ-
ations, and a simple numerical code for their solution is reported. The
solution comprises both the dynamic analysis of the underactuated sys-
tems in partly specifiedmotion and the synthesis of control that assures
realization of suchmotion.

Key words: underactuated systems, inverse dynamics, differential flat-
ness, servo-constraints

1. Introduction

Underactuated mechanical systems are mechanical systems with fewer control
inputs than degrees of freedom (Spong, 1997), m<n, and as such k=n−m
degrees of freedomcannot bedirectly actuated.Thedetermination of a control
input strategy that forces an underactuated system to complete a given set of
specified motion tasks is a challenging problem, and a solution to the inverse
dynamics problem at hand is possible only if the motion is defined by as
many system outputs as the number of control inputs – the motion is partly
specified. The situation is thus substantially different from the fully actuated
case, m=n, which requires also that the motion of the systemmust be fully
specified.
The controllability of underactuated systems in partly specifiedmotions is

closely related to the concept of differential flatness (Fliess et al., 1995; Ro-



384 W. Blajer, K. Kołodziejczyk

uchon, 2005), denoted that all the system state variables and control inputs
can be algebraically expressed in terms of desired outputs and their time deri-
vatives up to a certain order. The relationships between the imposed outputs
and the required control inputs provide a basis for the synthesis of control
laws. Such analytical algebraic relations are usually difficult to obtain, ho-
wever. The equivalent/alternative formulation reported in this paper is much
more straightforward and applicable.

2. Relevant technical examples

Themodelingmethodology isdeveloped for agroupof relevant technical exam-
ples of underactuatedmechanical systems in partly specifiedmotion, creating
a class of differentially flat systems (Blajer and Kołodziejczyk, 2004). A sim-
ple representative of the group is a two-mass system (n=2) shown in Fig.1,
where desiredmotion s(t) of themass m2 is actuated by the force F applied
to themass m1, and as such m= k=1=n/2. Another example of qualitati-
vely the same type is the position control of a manipulator with elastic joints
(Spong, 1987). For robot models with compliance lumped between rigid links
and actuators (Fig.2), the number of degrees of freedom is n=2m, related to
m links and m actuator rotors, q= [q1, . . . ,qm]

⊤ and θ= [θ1, . . . ,θm]
⊤. The

specifiedmotion qd(t) of links is actuated bymotor torques τ = [τ1, . . . ,τm]
⊤

applied to the actuator rotors, and as such m= k=n/2.

Fig. 1. A two-mass system

Fig. 2. A sketch of a flexible joint



Modeling of underactuated mechanical systems... 385

The overhead trolley crane and the rotary crane seen in Figures 3
and 4 are both five-degree-of-freedom systems, q = [s1,s2, l,θ1,θ2]

⊤ and
q = [ϕ,s,l,θ1,θ2]

⊤, respectively, the three desired outputs are specified in
time load coordinates rd(t) = [xd(t),yd(t),zd(t)]

⊤, and the respective three
control inputs are u = [F1,F2,Mw]

⊤ and u = [Mb,F,Mw]
⊤; see also Blajer

and Kołodziejczyk (2007) for more details. We have thus n= 5, m = 3 and
k=n−m=2.

Fig. 3. An overhead trolley cranemodel

Fig. 4. A rotary cranemodel

For an aircraft in a prescribed trajectory flight (Blajer et al., 2001), Fig.5,
motion of the six-degree-of-freedom system is prescribed by four demands:
the requested trajectory (two specifications), the condition on the airframe
attitude with respect to the trajectory, and the desired flight velocity. The
control inputs are u = [δa,δe,δr,T ]

⊤. We have then n = 6, m = 4 and
k=n−m=2.



386 W. Blajer, K. Kołodziejczyk

Fig. 5. An aircraft in a prescribed trajectory flight

Finally, for a surface ship tracking a desired path (Fig.6), the two desired
outputs of the three-degree-of-freedom system (Fosse, 1995) are the requested
trajectory and the desired velocity variation in time. The two control inputs
are u= [δ,T ]⊤. We have: n=3,m=2 and k=n−m=1.

Fig. 6. A surface ship tracking a trajectory

3. Modeling prelimaries

Thedynamical systemsquoted inSection 2 canbe representedbyan n-degree-
of-freedom system described by n generalized coordinates q = [q1, . . . ,qn]

⊤,
and enforced by actuator forces due to m control inputs, u= [u1, . . . ,um]

⊤,
where m<n. For the sake of simplicity, let uswrite dynamic equations of the
system in the following generic matrix form

M(q)q̈+d(q, q̇)=f(q, q̇)−B⊤(q)u (3.1)

where M is the n×n generalized mass matrix related to q̇, the n-vectors d
and f contain the related dynamic and applied generalized forces, and B⊤



Modeling of underactuated mechanical systems... 387

is an n×m matrix of coefficients that relate the control inputs u and the
actuating generalized force vector f

u
=−B⊤u. In fact, thedynamic equations

for aircrafts and ships are usually introduced as Mv̇+d=f−B⊤u, and then
supplemented with the kinematic equation q̇ = Av; see Blajer et al. (2001)
and Fossen (1995) for more details. However, the subsequent methodology
introduced for dynamic formulation (3.1) is valid formoregeneral formulations
as well.

The systemdescribed in (3.1) is subjected to mmotion specifications (per-
formance goals) formulated as m desired system outputs γ = [γ1, . . . ,γm]

⊤,
which, expressed in terms of q, lead to servo-constraints (called also active or
program constraints) on the system (Blajer and Kołodziejczyk, 2004; Gutow-
ski, 1971). The constraint equations, in the original and twice-differentiated
with respect to time forms, are

c(q, t)≡Φ(q)−γ(t)=0
(3.2)

c̈≡C(q)q̈−ξ(q, q̇, t)=0

where C= ∂Φ/∂q is an m×n constraintmatrix, ξ= γ̈−Ċq̇ is an m-vector of
constraint-induced accelerations. Evidently, equation (3.2)1 is mathematically
equivalent to m rheonomic holonomic constraints on the system, c(q, t)=0,
and γ(t) = 0 can be interpreted as drift in time of manifold Φ(q) = 0 in
the linear space related to q̇ (Blajer, 2001). The actuating generalized for-
ce f

u
= −B⊤u can be then interpreted as the generalized reaction force of

servo-constraints, called the control reaction (Blajer andKołodziejczyk, 2004)
in the sequel. The realization of servo-constraints (3.2)1 by applying the ac-
tive control force fu = −B

⊤u (control inputs u) is, however, qualitatively
different from the realization of contact (passive) constraints caused by joints,
hard surfaces, rigid links, slip-less rolling contacts, etc. Namely, assumed equ-
ation (3.2)1 represents contact constraints, the actuating force fu =−B

⊤u in
(3.1) would be replaced by the passive constraint reaction force f

c
=−C⊤λ,

where C is the constraint matrix as defined in (3.2)1. While the reaction fc
of (ideal) passive constraints is orthogonal to the instantaneous manifold of
contact constraints, the actuating force f

u
may be arbitrary directed with

respect to themanifold of servo-constraints, and, in the extreme, may be tan-
gent to the manifold (for more geometrical insight the reader is referred to
Blajer (2001), Blajer and Kołodziejczyk (2004)). More specifically, denoted:
N – linear n-space of the system related to q̇, C – constrained (specified)
m-subspace spanned by constraint gradients represented in N as rows in C,
B – controlled m-subspace spannedby vectors represented in N as rows in B,



388 W. Blajer, K. Kołodziejczyk

for all the dynamical systems reported in Section 2, the controllability (dif-
ferential flatness) of this class of underactuated mechanical systems in partly
specified motion is warranted by

C∪B=N ⇐⇒ rank([C⊤
...B⊤]⊤)=max=n =⇒ 2m­n (3.3)

where Q = [C⊤
...B⊤]⊤ is a 2m×n matrix. According to (3.3), the C and B

subspaces cover thus the whole space N (see Fig.7), and then

C∩B=P =⇒ rank(CM−1B⊤)= p =⇒ p=2m−n (3.4)

where P=CM−1B⊤ is a m×mmatrix. The deficiency in rank of P, p<m,
shows that only p desired outputs can directly be governed by the available
control inputs, while the other k =m−p= n−m outputs are actuated in-
directly through the dynamical couplings, referred respectively to orthogonal
and tangent realization of servo constraints (Blajer andKołodziejczyk, 2004).
For the systems in Figs. 1 and 2 we have 2m = n and p = 0 (pure tan-
gent realization of servo-constraints), and for the other systems in Figs. 3-6,
2m = n + 1 and p = 1 (mixed orthogonal-tangent realization of servo-
constraints).

Fig. 7. Illustration of C and B subspaces in N

4. Differential flatness formulation

The initial governing equations for the inverse dynamics problem of an unde-
ractuated system in partly specified motion are

q̇=v

M(q)v̇+d(q,v)=f(q,v)−B⊤(q)u (4.1)

0=Φ(q)−γ(t)



Modeling of underactuated mechanical systems... 389

which form 2n + m differential-algebraic equations (DAEs) in 2n states
x = [q⊤,v⊤]⊤, and m control inputs u. An important characteristic of a
DAE system is its index, which is ameasure of singularity of theDAEs andde-
termines difficulty in their numerical integration (Ascher and Petzold, 1998).
The index of DAEs (4.1) is equal to five (Blajer and Kołodziejczyk, 2004),
which produces nuisances in their numerical treatment (mostDAE solvers are
designed for index-one DAEs (Ascher and Petzold, 1998)).
In general, the index of aDAE system is the number of times one needs to

differentiate it to get a system of differential equations (Ascher and Petzold,
1998). This definition corresponds to the concept of differential flatness (Fliess
et al., 1995; Rouchon, 2005), i.e. an n-degree-of-freedom underactuated sys-
tem, dim(u) = m < n, is differentially flat if for a set of (flat) outputs γ,
dim(γ) =m, all states and controls of the system can algebraically be deter-
mined in terms of γ and their time derivatives up to a certain order. More
specifically, fromDAEs (4.1), one can theoretically obtain

x=x(t,γ, γ̇, . . . ,γr−1) and u=u(t,γ, γ̇, . . . ,γr) (4.2)

which can also be interpreted as index-oneDAEs. For all the dynamic systems
reported in Section 2we have then r=4 since the value of r is by one smaller
than the index of DAEs (4.1). Given relationships (4.2), both the dynamics of
an underactuated system in themotion specified by γ can be studied, and the
requested control that assures realization of the motion can be determined.
The problem is, however, that analytical solutions (4.2) to DAEs (4.1) are
attainable usually only for simple systems, and are very difficult to derive, if
possible at all, for more complex systems of technical relevance.
In order to illustrate the differential flatness concept, let us consider the

two-mass system defined in Fig.1. The respective DAE formulation, (4.1), is

ẋ1 = v1 m1v̇1 = k(x2−x1−d)+F

ẋ2 = v2 m2v̇2 =−k(x2−x1−d)

0=x2−s(t)

(4.3)

It can be then deduced that: M = diag(m1,m2), C = [0,1], B = [1,0],
P = CM−1B⊤ = [0], and p = 0 (pure tangent realization of the servo-
constraint c≡x2−s(t)= 0 is faced).TheDAEs (4.3) canbe thenmanipulated
to

x1 =
m2
k
s̈+(s−d) x2 = s

v1 =
m2
k
s(3)+ ṡ v2 = ṡ

F =
m1m2
k
s(4)+(m1+m2)s̈

(4.4)



390 W. Blajer, K. Kołodziejczyk

which corresponds to (4.2) with x= [q⊤,v⊤] = [x1,x2,v1,v2]
⊤, u= [F], and

γ = [s]. As seen r =4, and as such the index of DAEs (4.3) is equal to five.
The last equation from (4.4) provides immediately a feedforward control law
for the underactuated system in the partly prescribed motion. In accordance
with (4.2), the third order (r− 1 = 3) time derivative of γ is involved in
the determination of the states x. More strictly, the third time derivative
of γ is involved in the determination of velocities v = [v1,v2]

⊤, while the
determination of coordinates q = [x1,x2]

⊤ requires only the second order
time derivative of γ.

5. Index-three DAE formulation and solution code

As previously said, differential flatness formulation (4.2) is usually very dif-
ficult to attain analytically for more complex systems, if not beyond one’s
grasp. Amuchmore applicable formulation was proposed in Blajer andKoło-
dziejczyk (2004), followed the projection of dynamic equations (3.1) into the
complementary constrained (specified) C and unconstrained (unspecified) D
subspaces of N . The constrained m-subspace C has already been defined as
spanned by the vectors represented as rows in the m×n constraintmatrix C.
The k-subspace D (k=n−m) so that C∪D=N and C∩D=O can be then
defined as spanned by vectors represented as columns in an n× k full-rank
matrix D such that CD = 0 ⇔ D⊤C⊤ = 0 (which stands for C ∩D = O),

and then rank([C⊤
...D⊤]⊤)=n (which stands for C∪D=N); for more geo-

metrical insight, the reader is referred to Blajer (2001). For a given matrix C
of full row-rank, its orthogonal complement D as above can sometimes be
guessed or found by inspection (usually for simple systems only) or determi-
ned numerically following a range of computer oriented codes. A convenient
scheme of the latter code patterns after the coordinate partitioning method

(Wehage andHaug, 1982), in which C is first factorized to C= [U
...W], where

U and W are m× k and m×m matrices, respectively, and det(W) 6= 0.

Then, D= [I
...(−W−1U)⊤]⊤, where I is the k×k identity matrix.

The formula of projection of (3.1) into D and C is (Blajer, 2001)

[

D
⊤

CM
−1

]

(

Mv̇+d−f+B⊤u
)

=0 (5.1)



Modeling of underactuated mechanical systems... 391

The projection into D yields k differential equations H(q)v̇ = h(q,v,u, t),
where H=D⊤M is the k×n matrix and h=D⊤(f −d−B⊤u) is the k-
vector. The projection into C, after using (3.2)2 for substituting Cv̇ = Cq̈
with ξ, leads then to ξ + CM−1(d − f + B⊤u) = 0, which represents
m algebraic equations in the system states x = [q⊤,v⊤]⊤ and control in-
puts u, whose symbolic form is b(q,v,u, t) = 0. Since, according to (3.4),
rank(CM−1B⊤) = p<m, these m algebraic equations impose only p condi-
tions on u, and as such k = m− p restrictions on the system motion, and
in particular on the system coordinates q, supplementary to m original re-
strictions c(q, t) = 0. The total number of the supplementary and original
restrictions on q is thus k+m = n, and as such, in this indirect way, the
motion of the system is ’fully specified’, i.e. q and then v can be determi-
ned in terms of γ and their time derivatives. More strictly, since ξ involved
in b(q,v,u, t) = 0 depends on γ̈, q determined this way depends on at
most γ̈, and v depends on at most γ(3). Then, u determined from p inde-
pendent conditions imposed by b(q,v,u, t) = 0 and k differential equations
H(q)v̇ =h(q,v,u, t), p+k =m, depends on at most γ(4). The conclusions
correspond to differential flatness formulation (4.4) for the case study with
p = 0, though the upper time derivatives γ(3) and γ(4) are not introduced
explicitly in the present formulation. The governing equations are then

q̇=v H(q)v̇=h(q,v,u, t)

0= b(q,v,u, t) 0= c(q, t)
(5.2)

which form n+k+m+m=2n+m index-threeDAEs in n+n+m=2n+m
variables q, v and u, which are just in the middle between initial index-five
DAE formulation (4.1) and differential flatness (index-one DAE) formulation
(4.2).
The solution to DAEs (5.2) are variations in time of state variables of

the system executing the prescribedmotion, q
d
(t) and vd(t), and the control

ud(t) that is required to enforce themotion.Due to a specific structure ofDAE
system (5.2), a simple algorithm based on the Euler backward approximation
method can be recommended to solve the DAEs. Namely, given qn and vn
at time tn (un is not required), the values qn+1, vn+1 and un+1 at time
tn+1 = tn+∆t, can be found as solution to the following nonlinear equations

qn+1−qn
∆t

−vn+1 =0

H(q
n+1)
vn+1−vn
∆t

−h(q
n+1,vn+1,un+1, tn+1)=0 (5.3)

b(q
n+1,vn+1,un+1, tn+1)=0 c(qn+1, tn+1)=0



392 W. Blajer, K. Kołodziejczyk

where ∆t is the integration time step. This effective simple scheme is of ac-
ceptable accuracy for an appropriately small value of ∆t. Then, due to the
original servo-constraints c(q, t) = 0 involved in (5.2)/(5.3), supplemented
then by b(q,v,u, t) = 0, dependent respectively on γ(t) and γ̈(t), q

d
(t) is

determined with a numerical accuracy of solving (5.3), and only vd(t) and
ud(t) are determined with an error followed the rough backward difference
method. The truncation errors do not accumulate in the simulation time, ho-
wever. The proposed simple solution code leads thus to reasonable and stable
solutions. The inverse simulation control ud(t) can be used as feedforward
control law for the system subjected to execute the desiredmotion. The solu-
tion scheme can be thenmodified to include feedback control to provide stable
tracking of the specified outputs γ(t) in the presence of possible external per-
turbations modeling inconsistencies; see Blajer and Kołodziejczyk (2004) for
more details.

6. Discussion and conclusions

The formulation developed in this paper relate a class of underactuated me-
chanical systems in partly specifiedmotion, whose relevant technical examples
were reported in Section 2. All of them are differentially flat systems (Fliess
et al., 1995; Rouchon, 2005), which assures solvability of this specific inverse
dynamics problem (controllability of the system in specified motion). From
this point of view, the role of flat systemswithin the set of underactuated sys-
tems is very similar to the role of integrable systemswithin the set of ordinary
differential equation systems (Rouchon, 2005). On the other hand, there is no
a general algorithm to decide whether the system is flat or not. The definition
that all the states and control inputs of a flat system can be determined in
terms of the desired outputs and their time derivatives may be difficult to
implement for more complex systems.

Differential flatness formulation (4.2), due to its inherent complexity, is
applicable rather to simple systems only. Much more convenient in practi-
cal application is index-three DAE formulation (5.2) motivated in this paper.
While both the formulations lead to exactly the same (numerical) solutions,
the excessive complexity of the differential flatness formulation is completely
removed in the present formulation. Of importance is also that, for the con-
sidered class of differentially flat systems, the differential flatness formulation
involves γ(4), while the present formulation requires γ̈ only. The analytical



Modeling of underactuated mechanical systems... 393

solution of the differential flatness formulation is then replaced by simple (and
stable) numerical solution code (5.3).

The present formulation has been successfully applied to numerical simu-
lations of many underactuated mechanical systems executing partly specified
motions, including the two-mass system (Blajer andKołodziejczyk, 2004) (re-
ported also briefly in this paper), an aircraft in a prescribed trajectory flight
(Blajer et al., 2001), and cranes executing prescribedmotions of payloads (Bla-
jer andKołodziejczyk, 2007). Further developments of this useful approach to
the inverse dynamics problem of underactuated systems subjected to execute
partly specifiedmotions are foreseenwithin the frameof theprojectmentioned
below.

Acknowledgment

The work was partically financed from the government support for scientific re-

search within years 2006-2008, under grant No. 4T12C06230.

References

1. Ascher U.M., Petzold L.R., 1998,Computer Methods for Ordinary Diffe-
rential Equations and Differential-Algebraic Equations, SIAM, Philadelphia

2. Blajer W., 2001, A geometrical interpretation and uniform matrix formula-
tion of multibody system dynamics, ZAMM, 81, 247-259

3. BlajerW.,Graffstein J.,KrawczykM., 2001,Predictionof the dynamic
characteristics and control of aircraft in prescribed trajectory flight, Journal of
Theoretical and Applied Mechanics, 39, 1, 79-103

4. Blajer W., Kołodziejczyk K., 2004, A geometric approach to solving pro-
blems of control constraints: theory and a DAE framework,Multibody System
Dynamics, 11, 343-364

5. Blajer W., Kołodziejczyk K., 2007, A DAE formulation for the dynamic
analysis and control design of cranes executing prescribedmotions of payloads,
In: J.C. Garc̀ıa Orden, J.M. Goicolea and J. Cuadrado: Multibody Dynamics,
Series:Computational Methods and Applied Sciencies, Springer, 91-113

6. Fliess M., Lévine J., Martin P., Rouchon P., 1995, Flatness and defect
of nonlinear systems: introductory theory and examples, International Journal
of Control, 61, 1327-1361

7. Fossen T.I., 1995,Guidance and Control of Ocean Vehicles, Willey

8. Gutowski R., 1971,Analytical Mechanics, PWN [in Polish]



394 W. Blajer, K. Kołodziejczyk

9. Rouchon P., 2005, Flatness based control of oscillators,ZAMM, 85, 411-421

10. Spong M.W., 1987, Modeling and control of elastic joint robot, Journal of
Dynamic Systems, Measurement, and Control, 109, 310-319

11. Spong M.W., 1997, Underactuated mechanical systems, In: B. Siciliano and
K.P. Valavanis (eds): Control Problems in Robotics and Automation, Lecture
Notes in Control and Information Sciences, 230, Springer, London, UK

12. Wehage R.A., Haug E.J., 1982, Generalized coordinate partitioning for di-
mension reduction in analysis of constrained dynamic systems, Journal of Me-
chanical Design, 104, 247-255

Modelowanie sterowanych układów mechanicznych realizujących ruch

programowy niezupełny

Streszczenie

Praca dotyczy zagadnienia symulacji dynamicznej odwrotnej klasy układówme-
chanicznych, w których liczba sygnałów sterowania, równa liczbie zadanych charak-
terystyk ruchu (sygnałów wyjściowych), jest mniejsza od liczby stopni swobody, re-
prezentowanej przez liczne przykłady techniczne. Równania tak zdefiniowanego ru-
chu programowego niezupełnego otrzymywane są w postaci równań różniczkowo-
algebraicznych o indeksie równym trzy. Omawiany jest prosty, efektywny i stabilny
algorytm całkowania tych równań. Rozwiązaniem są przebiegi w czasie zmiennych
stanu układu w zadanym ruchu programowym oraz sterowania wymaganego dla re-
alizacji tego ruchu.

Manuscript received October 24, 2007; accepted for print January 3, 2008