# Optimal Manoeuvres of Underactuated Linear Mechanical Systems: The Case of Controlling Gantry Crane Operations

Authors: S. Woods; W. Szyszkowski
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

A method of solving optimal manoeuvre control of linear underactuated mechanical systems is presented. The nonintegrable constraints present in such systems are handled by adding dummy actuators and then by applying Lagrange multipliers to reduce their action to zero. The open- and closed-loop control schemes can be analyzed. The method, referred to as the constrained modal space optimal control (CMSOC), is illustrated in the examples of gantry crane operations.


---

## Body

## 1. Introduction

Underactuated mechanical systems have fewer independent actuators than degrees of freedom (DOFs) to be controlled [1]. Typical nonlinear examples of such systems, usually with only several DOFs, are rigid multilink robotic manipulators with passive joints or any manipulator with flexible links (described by at least one mode of vibration). Linear examples include vibrating structures with continuously distributed mass (i.e., with theoretically infinite number of DOFs to describe them) such as masts, antennas, buildings, brides, and car suspension, controlled by discrete actuators. This paper presents a method of analyzing and simulating optimal manoeuvres between two given configurations (often referred to as point-to-point manoeuvres) for linear underactuated systems. The method combines optimal control theory with computational mechanics and the finite element (FE) technique, in particular.The number of DOFs equal to the number of actuators will be referred to as actuated (after [1]), while all remaining DOFs will be referred to as underactuated (however, all DOFs are in fact controlled). The actuated and unactuated DOFs must satisfy a number of constraints equal to the number of unactuated DOFs and resulting from the equations governing the motion of such systems. For mechanical systems we assume that these constraints may be nonintegrable (nonholonomic), meaning unactuated DOFs cannot be explicitly eliminated. Many of the techniques presented in the literature deal with underactuated problems by applying the constraints to eliminate the unactuated DOFs and then by solving the reduced fully actuated problems [2–4]. These approaches are limited to particular problems where the constraints can be simplified to a form making such mathematical manipulations possible. The method presented here is capable of dealing with any linear system, as it does not require the elimination of unactuated DOFs. Instead, the underactuated system is formulated as if it were fully actuated by adding “dummy” (zero-valued) actuators to all unactuated DOFs. The modal space is used in modelling the system motions. The method can be considered as an extension of the independent modal space control (IMSC; e.g., see [5]) into the underactuated problems, therefore it will be referred to as the constrained modal space optimal control (CMSOC) method. The system constraints resulting from underactuation are then determined by eliminating these dummy actuators. The constraints are algebraic in terms of controls but differential (nonintegrable) in terms of the DOFs. The algebraic form of the constraints is used to generate the so-called matrix of constraints, which is utilized to handle the nonintegrable constraints with the help of time-varying Lagrange multipliers. Pontryagin’s principle is used to optimize the trajectory and actuation forces.This paper presents the CMSOC method in a general form and then explains some details of the corresponding numerical procedure in the examples of standard two or three-DOF gantry crane operations. The method is verified by recreating the closed-loop control of the two-DOF gantry crane problem obtained in [1] via applying the classical technique and the open-loop optimal control considered in [3].


## 2. Problem Formulation

### 2.1. Dynamics of a General Underactuated System

The computational model for the motion of a linear mechanical system is represented by a standard form used in FE analysis:(1)
M
q
¨
+
C
q
˙
+
K
q
=
B
F
a
,
where q and F
a are vectors of DOFs and activation forces, respectively, and M, C, and K are constant mass, damping, and stiffness matrices, respectively. In particular, (1) is suitable to model the dynamics of a range of actively controlled structural members undergoing small amplitude oscillations and finite translations. In underactuated systems n
a independent actuation forces are to control n
>
n
a number of DOFs. Matrix B of dimensions n
×
n
a assigns components of vector  F
a to particular DOFs and obviously is not invertible if  n
≠
n
a. Clearly, the actuators via (1) control all DOFs of the system. For the purpose of analysis the DOFs can be divided into actuated (q
a) and unactuated (q
r) ones by rearranging these equations as follows: (2a)
[
M
a
a
M
a
r
M
r
a
M
r
r
]
[
q
¨
a
q
¨
r
]
+
[
C
a
a
C
a
r
C
r
a
C
r
r
]
[
q
˙
a
q
˙
r
]
+
[
K
a
a
K
a
r
K
r
a
K
r
r
]
[
q
a
q
r
]
=
[
F
a
0
]
.The bottom row represents the equations constraining the actuated and unactuated DOFs in the following form:
(3)
M
r
a
q
¨
a
+
C
r
a
q
˙
a
+
K
r
a
q
a
+
M
r
r
q
¨
r
+
C
r
r
q
˙
r
+
K
r
r
q
r
=
0
.
The system can formally be converted to a fully actuated one by using (3) to explicitly determine vector q
r in terms of q
a (i.e, q
r
=
g
(
q
a
)), and then by substituting this vector to the top row of (2a) to obtain F
a
=
F
a
(
q
a
,
g
(
q
a
)
). Unless some matrices in (3) vanish, it is not generally possible, and therefore these constraints are considered as nonholonomic.The control task for vectorF
a in (1) is to manoeuvre the system from an initial state to a final state described by the following boundary conditions (point-to-point manoeuvre):
(4)
q
(
0
)
=
q
0
,
q
˙
(
0
)
=
q
˙
0
,
q
(
t
f
)
=
q
f
,
q
˙
(
t
f
)
=
q
˙
f
.
It should be emphasised that no trajectory is specified in this task. A particular trajectory satisfying (1) and (4) may be determined if extra conditions are imposed on the system. We will identify such a trajectory by optimizing the performance index as discussed in the next section. Note that this problem is different from a typical trajectory tracking problem in which instead of (4) the task is specified as the system output in the form(5a)
y
=
h
(
q
)
.Several methods have been proposed to solve the inverse problems of finding the input F
a for the output y as defined by (1) and (5a), notably the servo-constraint approach [6–9] and the flatness method [10, 11]. In particular, differentiating (5a) with n
a outputs twice one obtains y
¨
=
H
q
¨
+
h
-, where the size of matrix H is n
a
×
n; then square matrix H
M
-
1
B is required to be nonsingular to solve the problem. This condition does not apply in the method presented here since our output is given only in terms of (4), that is, the system’s initial and final configurations.The set of (1) or (2a) is uncoupled when mapped into modal space, where vector of DOFs q (size n) is transformed to the equally sized vector of modal variables η
=
[
η
1
⋯
η
n
]
T. Similarly, vector F is related to an equally sized vector of modal controls U
=
[
u
1
⋯
u
n
]
T. These transformations are(6a)
q
=
Φ
η
,
(6b)
U
=
(
Φ
T
B
)
F
a
=
B
^
F
a
,where B
^
=
Φ
T
B is the transfer matrix of size n
×
n
a between vectors F
a and U and mode shape matrix Φ
=
[
φ
1
⋯
φ
n
] relates vectors q and η. The M-normalized matrix Φ, consisting of n modal shape vectors φ
i (each with n components), satisfies the following orthogonality conditions: (7a)
Φ
T
M
Φ
=
I
,
(7b)
Φ
T
K
Φ
=
Ω
,where I is the unitary matrix and Ω is the diagonal matrix of ordered frequencies with the terms Ω
i
i
=
ω
i
2. Each mode shape vector φ
i and frequency ω
i are solutions to the eigenvalues problem (
K
-
ω
i
2
M
)
φ
i
=
0 (i
=
1
,
…
,
n). The above modal analysis (or operations defined by (6a)–(7b)) is carried out routinely in the FE approach, even for problems with a very large number of DOFs (large n).The equations of motion (1) become uncoupled when applying transformations (6a) and (6b) subject to orthogonality conditions (7a) and (7b) and take the following form:
(8)
I
η
¨
+
Δ
η
˙
+
Ω
η
=
U
,
where for the Rayleigh damping (i.e., C
=
α
M
+
β
K) the diagonal terms of Δ are Δ
i
i
=
2
ς
i
ω
i
=
φ
i
T
C
φ
i and where ς
i
=
α
/
2
ω
i
+
β
ω
i
/
2 are the modal damping ratios. Note that a rigid body translation, for which ω
i
=
0, is also included in the above equation.A continuous system, or an FE model (1) of the system described by n DOFs (where n may be a large number), can be approximated by (8) with only n
m significant modes considered, where usually n
m
≪
n. The number of significant modes that should be sufficient to represent such a system is generally problem related and depends mainly on its physical characteristics, the spatial distribution, and frequency content of the loading [12].In the system approximated byn
m modes (similarly as for the system's DOFs) one can consider n
a modes as actuated and n
r
=
n
m
-
n
a modes as unactuated. Then matrix Φ will be reduced to only n
m columns, and transfer matrix B
^ in (6a) and (6b) will be of dimensions n
m
×
n
a.In order to control alln
m modes this system can be made artificially fully actuated by adding n
r
=
n
m
-
n
a dummy actuation forces (zero valued) forming subvector F
d. For the purpose of analysis vector F
a in (1) is replaced by the augmented force vector F
a
'
=
[
F
a
T
F
d
T
]
T containing n
a real actuation forces forming subvector F
a and n
r dummy actuators forming subvector F
d. Then, in (1) and (6a) and (6b), matrix B of dimensions n
×
n
a is replaced by matrix B
′ of dimensions n
×
n
m (this matrix assigns the component of F
a
′ to particular nodes). Consequently in (6b) matrix B
^ of dimensions n
×
n
a is replaced by a square matrix B
^
′ of dimensions n
m
×
n
m (n
m modes controlled by n
m actuators). The dummy actuators F
d should be placed in such a way that B
^
′ is nonsingular.In the inverse dynamics based control analysis, each controlU
i can be obtained from (8) by substituting the corresponding prescribed mode η
i. Then, for known vector U, the actuation forces F
a
' should be determined by inverting transformation (6b) in which matrix B
^
′ (instead of B
^) is now square and nonsingular (the dummy actuators were added to the system only to ensure that this inversion is possible). In the next step, after computing the inverse of operation (6b), the dummy actuators will be eliminated by giving them zero values. For that purpose the inverse matrix (
B
^
′
)
-
1, representing the mapping from modal controls U
=
[
U
a
T
∣
U
r
T
]
T to actuation forces F
a
′ for any augmented system of size n
m
×
n
m, is partitioned as follows:
(9)
(
B
^
′
)
-
1
U
=
F
a
′
⟹
[
B
~
a
B
~
r
A
a
A
r
]
[
U
a
U
r
]
=
[
F
a
F
d
]
=
[
F
a
0
]
.
Square submatrix B
~
a is of size n
a
×
n
a and square submatrix A
r is of size n
r
×
n
r. To be consistent with modes classifications (actuated and unactuated), vectors U
a
=
[
u
1
⋯
u
n
a
]
T and U
r
=
[
u
n
a
+
1
⋯
u
n
m
]
T are referred to as actuated and unactuated modal controls, respectively. Given the null-valued dummy force vector F
d
=
[
0
⋯
0
]
T (size n
r
×
1) the bottom n
r rows of operation (9) (lower partition) define constraints on the system in terms of all modal controls, in the following form:
(10)
A
U
=
A
a
U
a
+
A
r
U
r
=
0
.
Matrix A
=
[
A
a
∣
A
r
] (size n
r
×
n
m) defines the system constraints written algebraically in terms of modal controls. Since (10) is homogeneous matrix A can be normalized such that the diagonal terms corresponding to controls U
a are set to unity (i.e., A
i
i
=
1). In this form A becomes independent of the choice of dummy actuators, which reflects the fact that these zero-force actuators were added somewhat arbitrarily only to facilitate the elimination process, that is, to satisfy the constraints in (9). Matrix A is discussed with more details in [13, 14].Real actuation force(s) may be obtained from the top partition of operation (9) in terms of all modal controls in vector U. They can also be obtained in terms of only actuated modal controls in vector U
a by applying n
r constraints (10) to eliminate unactuated modal controls U
r. Thus, n
a components of actuator forces in vector F
a can be obtained in terms of n
a actuated modal controls in vector U
a from the following operation:
(11)
F
a
=
B
-
U
a
.
Square matrix B
-
=
B
~
a
-
B
~
r
A
r
-
1
A
a (size n
a
×
n
a) is referred to as the pseudotransfer matrix, and it relates actuated modal controls to real actuator forces. Similar to the normalized constraint matrix A, the pseudotransfer matrix B
- is independent of the choice of dummy actuators.


### 2.2. Optimal Manoeuvres of Underactuated Systems

In linear optimal control [15], the manoeuvre is optimal if, for a given task, it minimizes the performance index:
(12)
J
=
1
2
∫
0
t
f
(
(
η
T
Q
^
d
η
)
+
(
η
˙
T
Q
^
v
η
˙
)
+
(
U
T
R
^
U
)
)
d
t
⟶
min
,
where Q
^
d, Q
^
v, and R
^ are matrices, with the diagonal terms Q
^
d
i
i, Q
^
v
i
i, and R
^
i
i  (i
=
1
,
…
,
n
m), that are weights for the system’s potential energy, kinetic energy, and actuator work, respectively. Note that n
m modal variables and modal controls are included in (12); however, these modes are not independent because of constraint (10), resulting from underactuation. Such a problem can be solved by applying Pontryagin’s principle. Here we use the procedure described in [15]. Hamiltonian H for the constrained optimization problem involving performance index (12), uncoupled equations of motion (8), and constraints (10) is defined in the following form:
(13)
H
=
-
1
2
(
η
T
Q
^
d
η
+
η
˙
T
Q
^
v
η
˙
+
U
T
R
^
U
)
+
P
d
η
˙
+
P
v
(
-
Δ
η
˙
-
Ω
η
+
U
)
+
v
T
A
U
.
P
d and P
v are standard costate vectors related to modal position and velocity states (η and η
˙) of a system, respectively. Vector v
T
=
[
v
1
⋯
v
n
r
] represents the set of time-dependent Lagrange multipliers introduced to enforce constraints (10). These multipliers play a similar role to, for example, that of the multipliers used in the servo-constraint approach [6–9] mentioned before. According to Pontryagin’s principle the costate equations take the following form:(14a)
P
˙
d
=
-
∂
H
∂
η
=
Q
^
d
η
+
Ω
P
v
,
(14b)
P
˙
v
=
-
∂
H
∂
η
˙
=
Q
^
v
η
˙
-
P
d
+
Δ
P
v
.The Hamiltonian is stationary with respect to modal control if
(15)
∂
H
∂
U
=
-
R
^
U
+
P
v
+
A
T
v
=
0
.
Substituting (8) into (15) gives
(16)
P
v
=
R
^
(
I
η
¨
+
Δ
η
˙
+
Ω
η
)
-
A
T
v
.
Substituting (16) into (14b) yields
(17)
P
d
=
Q
^
v
η
˙
-
R
^
(
I

η
∴
+
Δ
η
¨
+
Ω
η
˙
)
+
Δ
R
^
(
I
η
¨
+
Δ
η
˙
+
Ω
η
)
+
A
T
v
˙
-
Δ
A
T
v
.
Finally, substituting (17) into (14a) generates the following set of optimality equations:
(18)
R
^
η
¨
¨
+
(
2
Ω
R
^
-
Q
^
v
-
R
^
Δ
2
)
η
¨
+
(
R
^
Ω
2
+
Q
^
d
)
η
-
(
A
T
v
¨
-
Δ
A
T
v
˙
+
Ω
A
T
v
)
=
0
.
Note that n
m optimality equations (18) contain n
m unknown components in η and n
r unknown components in v. Therefore, additional n
r constraint equations (10) are required in order to obtain all the unknown modal variable functions in vector η and Lagrange multiplier functions in vector v. However, the constraints must be written in terms of η not in terms of U (note the change in the constraints’ form from algebraic to differential). The uncoupled equations of motion (8) are substituted into algebraic constraints (10) to obtain
(19)
A
(
I
η
¨
+
Δ
η
˙
+
Ω
η
)
=
0
.
The number of n
m
+
n
r ((18) and (19)) is equal the unknown components in vectors η and v.Boundary conditions (4) are mapped into modal space by using the inverse of transformation (6a) or through the relation η
=
Φ
T
M
q (obtained by additional substitution of condition (7a)). These transformed boundary conditions are
(20)
η
(
0
)
=
Φ
T
M
q
0
,
η
˙
(
0
)
=
Φ
T
M
q
˙
0
,
η
(
t
f
)
=
Φ
T
M
q
f
,
η
˙
(
t
f
)
=
Φ
T
M
q
˙
f
.
For fully actuated problems, the last term (A
T
v
¨
-
Δ
A
T
v
˙
+
Ω
A
T
v) in optimality equations (18) vanishes because there are no constraints or Lagrange multipliers needed to enforce them. Therefore, a fully actuated problem involves only n
m optimality equations (18) to be solved in terms of n
m uncoupled modal variables in vector η.The solution to the combined set of (18), (19), and (20) can be efficiently obtained using symbolic differential operator D
n
=
d
n
/
d
t
n. Substituting this operator into (19) and (20) and rewriting in matrix notation give
(21)
[
E
-
E
^
T
E
~
0
]
[
η
ν
]
=
0
or
E
p
Y
=
0
,
where
(22)
E
=
R
^
D
4
+
(
2
R
^
Ω
-
Q
^
v
-
R
^
Δ
2
)
D
2
+
(
R
^
Ω
2
+
Q
^
d
)
,
E
^
=
A
(
I
D
2
+
Δ
D
+
Ω
)
,
E
~
=
A
(
I
D
2
-
Δ
D
+
Ω
)
.
Matrix E
p contains submatrices E, E
~, and -
E
^
T. Vector Y
=
[
η
T
∣
ν
T
]
T contains all unknown modal variables and Lagrange multipliers. Note that in a fully actuated case, matrix E
p in (21) consists only of submatrix E and vector Y
=
η.The solution to a system described in form (21) involves the roots r
l (l
=
1
,
…
,
4
n
m) of the characteristic equation for the determinant of  E
p [16], where operator D is replaced by the auxiliary variable r rendering a 4
n
m
th order polynomial. This operation is written as
(23)
det
E
p
|
D
→
r
=
0
.
Generally, the 4
n
m roots of the characteristic equation (23) take the following form:
(24)
r
l
=
±
α
k
±
i
β
k
(
k
=
1
,
…
,
n
m
,
l
=
1
,
…
,
4
n
m
)
.
The positive real numbers α
k and β
k characterize the response of the k
th mode of motion. For nonzero, unique roots, solution vector Y consists of n
m
+
n
r components Y
j that can be written in terms of 4
n
m independent elementary functions related to the roots (24), in the form [16]:
(25)
Y
j
=
∑
k
=
1
n
m
[
e
α
k
t
(
c
k
j
1
sin
(
β
k
t
)
+
c
k
j
2
cos
(
β
k
t
)
)
+
e
-
α
k
t
(
c
k
j
3
sin
(
β
k
t
)
+
c
k
j
4
cos
(
β
k
t
)
)
e
α
k
t
(
c
k
j
1
sin
(
β
k
t
)
+
c
k
j
2
cos
(
β
k
t
)
)
]
,
where
j
=
1
,
…
,
n
m
+
n
r
.
Obviously, the frequency of k
th mode controlled by the actuators can be interpreted as ω
k
a
=
β
k and its rate of active attenuation (or amplification) as ς
k
a
=
α
k
/
β
k. If multiple roots and zero-valued roots are obtained from (23), then solution functions (25) must be modified to mathematically accommodate these situations. There are 4
n
m
(
n
m
+
n
r
) unknown integration constants c
k
j
1
,
…
,
c
k
j
4 contained in the solution functions (25).Integration constantsc
k
j
1
,
…
,
c
k
j
4 are obtained by substituting the assumed form (25) into differential equations (18) and (19) and using the method of undetermined coefficients to generate n
m
+
n
r sets of 4
n
m linear algebraic equations relating these constants. By replacing one set of 4
n
m equations with the set of 4
n
m boundary conditions (20), the integration constants can be solved simultaneously. All these symbolic operations, including the determination of the roots (24) and constants in (25), can be done automatically using the MAPLE mathematical software.For closed-loop control, asymptotically convergent solution functions are required such that the control task is met over an infinite period of time (t
f
→
∞). The resulting number of integration constants is reduced by half, as terms involving positive exponential e
α
k
t in the solution form (25) disappear (c
k
j
1
=
c
k
j
2
=
0).To quantitatively measure the performance of closed-loop control schemes settling timet
f
3
% is defined as the time needed for various variables to be reduced to within 3% of their initial value (i.e., e
-
α
k
min
t
f
3
%
=
0.03
→
t
f
3
%
=
3.5
/
α
k
min).The above procedure was applied to actively suppress vibrations of a spatial antenna mast in [13] and of plane frames in [14], the cases with the number of DOFs much greater than the number of significant modes included in the analysis (i.e., with n
≫
n
m). In both cases only oscillating modes were controlled. Here the application of the above methodology is focused on various control schemes, which are demonstrated in controlling the translational and oscillating modes of a gantry crane.


## 2.1. Dynamics of a General Underactuated System

The computational model for the motion of a linear mechanical system is represented by a standard form used in FE analysis:(1)
M
q
¨
+
C
q
˙
+
K
q
=
B
F
a
,
where q and F
a are vectors of DOFs and activation forces, respectively, and M, C, and K are constant mass, damping, and stiffness matrices, respectively. In particular, (1) is suitable to model the dynamics of a range of actively controlled structural members undergoing small amplitude oscillations and finite translations. In underactuated systems n
a independent actuation forces are to control n
>
n
a number of DOFs. Matrix B of dimensions n
×
n
a assigns components of vector  F
a to particular DOFs and obviously is not invertible if  n
≠
n
a. Clearly, the actuators via (1) control all DOFs of the system. For the purpose of analysis the DOFs can be divided into actuated (q
a) and unactuated (q
r) ones by rearranging these equations as follows: (2a)
[
M
a
a
M
a
r
M
r
a
M
r
r
]
[
q
¨
a
q
¨
r
]
+
[
C
a
a
C
a
r
C
r
a
C
r
r
]
[
q
˙
a
q
˙
r
]
+
[
K
a
a
K
a
r
K
r
a
K
r
r
]
[
q
a
q
r
]
=
[
F
a
0
]
.The bottom row represents the equations constraining the actuated and unactuated DOFs in the following form:
(3)
M
r
a
q
¨
a
+
C
r
a
q
˙
a
+
K
r
a
q
a
+
M
r
r
q
¨
r
+
C
r
r
q
˙
r
+
K
r
r
q
r
=
0
.
The system can formally be converted to a fully actuated one by using (3) to explicitly determine vector q
r in terms of q
a (i.e, q
r
=
g
(
q
a
)), and then by substituting this vector to the top row of (2a) to obtain F
a
=
F
a
(
q
a
,
g
(
q
a
)
). Unless some matrices in (3) vanish, it is not generally possible, and therefore these constraints are considered as nonholonomic.The control task for vectorF
a in (1) is to manoeuvre the system from an initial state to a final state described by the following boundary conditions (point-to-point manoeuvre):
(4)
q
(
0
)
=
q
0
,
q
˙
(
0
)
=
q
˙
0
,
q
(
t
f
)
=
q
f
,
q
˙
(
t
f
)
=
q
˙
f
.
It should be emphasised that no trajectory is specified in this task. A particular trajectory satisfying (1) and (4) may be determined if extra conditions are imposed on the system. We will identify such a trajectory by optimizing the performance index as discussed in the next section. Note that this problem is different from a typical trajectory tracking problem in which instead of (4) the task is specified as the system output in the form(5a)
y
=
h
(
q
)
.Several methods have been proposed to solve the inverse problems of finding the input F
a for the output y as defined by (1) and (5a), notably the servo-constraint approach [6–9] and the flatness method [10, 11]. In particular, differentiating (5a) with n
a outputs twice one obtains y
¨
=
H
q
¨
+
h
-, where the size of matrix H is n
a
×
n; then square matrix H
M
-
1
B is required to be nonsingular to solve the problem. This condition does not apply in the method presented here since our output is given only in terms of (4), that is, the system’s initial and final configurations.The set of (1) or (2a) is uncoupled when mapped into modal space, where vector of DOFs q (size n) is transformed to the equally sized vector of modal variables η
=
[
η
1
⋯
η
n
]
T. Similarly, vector F is related to an equally sized vector of modal controls U
=
[
u
1
⋯
u
n
]
T. These transformations are(6a)
q
=
Φ
η
,
(6b)
U
=
(
Φ
T
B
)
F
a
=
B
^
F
a
,where B
^
=
Φ
T
B is the transfer matrix of size n
×
n
a between vectors F
a and U and mode shape matrix Φ
=
[
φ
1
⋯
φ
n
] relates vectors q and η. The M-normalized matrix Φ, consisting of n modal shape vectors φ
i (each with n components), satisfies the following orthogonality conditions: (7a)
Φ
T
M
Φ
=
I
,
(7b)
Φ
T
K
Φ
=
Ω
,where I is the unitary matrix and Ω is the diagonal matrix of ordered frequencies with the terms Ω
i
i
=
ω
i
2. Each mode shape vector φ
i and frequency ω
i are solutions to the eigenvalues problem (
K
-
ω
i
2
M
)
φ
i
=
0 (i
=
1
,
…
,
n). The above modal analysis (or operations defined by (6a)–(7b)) is carried out routinely in the FE approach, even for problems with a very large number of DOFs (large n).The equations of motion (1) become uncoupled when applying transformations (6a) and (6b) subject to orthogonality conditions (7a) and (7b) and take the following form:
(8)
I
η
¨
+
Δ
η
˙
+
Ω
η
=
U
,
where for the Rayleigh damping (i.e., C
=
α
M
+
β
K) the diagonal terms of Δ are Δ
i
i
=
2
ς
i
ω
i
=
φ
i
T
C
φ
i and where ς
i
=
α
/
2
ω
i
+
β
ω
i
/
2 are the modal damping ratios. Note that a rigid body translation, for which ω
i
=
0, is also included in the above equation.A continuous system, or an FE model (1) of the system described by n DOFs (where n may be a large number), can be approximated by (8) with only n
m significant modes considered, where usually n
m
≪
n. The number of significant modes that should be sufficient to represent such a system is generally problem related and depends mainly on its physical characteristics, the spatial distribution, and frequency content of the loading [12].In the system approximated byn
m modes (similarly as for the system's DOFs) one can consider n
a modes as actuated and n
r
=
n
m
-
n
a modes as unactuated. Then matrix Φ will be reduced to only n
m columns, and transfer matrix B
^ in (6a) and (6b) will be of dimensions n
m
×
n
a.In order to control alln
m modes this system can be made artificially fully actuated by adding n
r
=
n
m
-
n
a dummy actuation forces (zero valued) forming subvector F
d. For the purpose of analysis vector F
a in (1) is replaced by the augmented force vector F
a
'
=
[
F
a
T
F
d
T
]
T containing n
a real actuation forces forming subvector F
a and n
r dummy actuators forming subvector F
d. Then, in (1) and (6a) and (6b), matrix B of dimensions n
×
n
a is replaced by matrix B
′ of dimensions n
×
n
m (this matrix assigns the component of F
a
′ to particular nodes). Consequently in (6b) matrix B
^ of dimensions n
×
n
a is replaced by a square matrix B
^
′ of dimensions n
m
×
n
m (n
m modes controlled by n
m actuators). The dummy actuators F
d should be placed in such a way that B
^
′ is nonsingular.In the inverse dynamics based control analysis, each controlU
i can be obtained from (8) by substituting the corresponding prescribed mode η
i. Then, for known vector U, the actuation forces F
a
' should be determined by inverting transformation (6b) in which matrix B
^
′ (instead of B
^) is now square and nonsingular (the dummy actuators were added to the system only to ensure that this inversion is possible). In the next step, after computing the inverse of operation (6b), the dummy actuators will be eliminated by giving them zero values. For that purpose the inverse matrix (
B
^
′
)
-
1, representing the mapping from modal controls U
=
[
U
a
T
∣
U
r
T
]
T to actuation forces F
a
′ for any augmented system of size n
m
×
n
m, is partitioned as follows:
(9)
(
B
^
′
)
-
1
U
=
F
a
′
⟹
[
B
~
a
B
~
r
A
a
A
r
]
[
U
a
U
r
]
=
[
F
a
F
d
]
=
[
F
a
0
]
.
Square submatrix B
~
a is of size n
a
×
n
a and square submatrix A
r is of size n
r
×
n
r. To be consistent with modes classifications (actuated and unactuated), vectors U
a
=
[
u
1
⋯
u
n
a
]
T and U
r
=
[
u
n
a
+
1
⋯
u
n
m
]
T are referred to as actuated and unactuated modal controls, respectively. Given the null-valued dummy force vector F
d
=
[
0
⋯
0
]
T (size n
r
×
1) the bottom n
r rows of operation (9) (lower partition) define constraints on the system in terms of all modal controls, in the following form:
(10)
A
U
=
A
a
U
a
+
A
r
U
r
=
0
.
Matrix A
=
[
A
a
∣
A
r
] (size n
r
×
n
m) defines the system constraints written algebraically in terms of modal controls. Since (10) is homogeneous matrix A can be normalized such that the diagonal terms corresponding to controls U
a are set to unity (i.e., A
i
i
=
1). In this form A becomes independent of the choice of dummy actuators, which reflects the fact that these zero-force actuators were added somewhat arbitrarily only to facilitate the elimination process, that is, to satisfy the constraints in (9). Matrix A is discussed with more details in [13, 14].Real actuation force(s) may be obtained from the top partition of operation (9) in terms of all modal controls in vector U. They can also be obtained in terms of only actuated modal controls in vector U
a by applying n
r constraints (10) to eliminate unactuated modal controls U
r. Thus, n
a components of actuator forces in vector F
a can be obtained in terms of n
a actuated modal controls in vector U
a from the following operation:
(11)
F
a
=
B
-
U
a
.
Square matrix B
-
=
B
~
a
-
B
~
r
A
r
-
1
A
a (size n
a
×
n
a) is referred to as the pseudotransfer matrix, and it relates actuated modal controls to real actuator forces. Similar to the normalized constraint matrix A, the pseudotransfer matrix B
- is independent of the choice of dummy actuators.


## 2.2. Optimal Manoeuvres of Underactuated Systems

In linear optimal control [15], the manoeuvre is optimal if, for a given task, it minimizes the performance index:
(12)
J
=
1
2
∫
0
t
f
(
(
η
T
Q
^
d
η
)
+
(
η
˙
T
Q
^
v
η
˙
)
+
(
U
T
R
^
U
)
)
d
t
⟶
min
,
where Q
^
d, Q
^
v, and R
^ are matrices, with the diagonal terms Q
^
d
i
i, Q
^
v
i
i, and R
^
i
i  (i
=
1
,
…
,
n
m), that are weights for the system’s potential energy, kinetic energy, and actuator work, respectively. Note that n
m modal variables and modal controls are included in (12); however, these modes are not independent because of constraint (10), resulting from underactuation. Such a problem can be solved by applying Pontryagin’s principle. Here we use the procedure described in [15]. Hamiltonian H for the constrained optimization problem involving performance index (12), uncoupled equations of motion (8), and constraints (10) is defined in the following form:
(13)
H
=
-
1
2
(
η
T
Q
^
d
η
+
η
˙
T
Q
^
v
η
˙
+
U
T
R
^
U
)
+
P
d
η
˙
+
P
v
(
-
Δ
η
˙
-
Ω
η
+
U
)
+
v
T
A
U
.
P
d and P
v are standard costate vectors related to modal position and velocity states (η and η
˙) of a system, respectively. Vector v
T
=
[
v
1
⋯
v
n
r
] represents the set of time-dependent Lagrange multipliers introduced to enforce constraints (10). These multipliers play a similar role to, for example, that of the multipliers used in the servo-constraint approach [6–9] mentioned before. According to Pontryagin’s principle the costate equations take the following form:(14a)
P
˙
d
=
-
∂
H
∂
η
=
Q
^
d
η
+
Ω
P
v
,
(14b)
P
˙
v
=
-
∂
H
∂
η
˙
=
Q
^
v
η
˙
-
P
d
+
Δ
P
v
.The Hamiltonian is stationary with respect to modal control if
(15)
∂
H
∂
U
=
-
R
^
U
+
P
v
+
A
T
v
=
0
.
Substituting (8) into (15) gives
(16)
P
v
=
R
^
(
I
η
¨
+
Δ
η
˙
+
Ω
η
)
-
A
T
v
.
Substituting (16) into (14b) yields
(17)
P
d
=
Q
^
v
η
˙
-
R
^
(
I

η
∴
+
Δ
η
¨
+
Ω
η
˙
)
+
Δ
R
^
(
I
η
¨
+
Δ
η
˙
+
Ω
η
)
+
A
T
v
˙
-
Δ
A
T
v
.
Finally, substituting (17) into (14a) generates the following set of optimality equations:
(18)
R
^
η
¨
¨
+
(
2
Ω
R
^
-
Q
^
v
-
R
^
Δ
2
)
η
¨
+
(
R
^
Ω
2
+
Q
^
d
)
η
-
(
A
T
v
¨
-
Δ
A
T
v
˙
+
Ω
A
T
v
)
=
0
.
Note that n
m optimality equations (18) contain n
m unknown components in η and n
r unknown components in v. Therefore, additional n
r constraint equations (10) are required in order to obtain all the unknown modal variable functions in vector η and Lagrange multiplier functions in vector v. However, the constraints must be written in terms of η not in terms of U (note the change in the constraints’ form from algebraic to differential). The uncoupled equations of motion (8) are substituted into algebraic constraints (10) to obtain
(19)
A
(
I
η
¨
+
Δ
η
˙
+
Ω
η
)
=
0
.
The number of n
m
+
n
r ((18) and (19)) is equal the unknown components in vectors η and v.Boundary conditions (4) are mapped into modal space by using the inverse of transformation (6a) or through the relation η
=
Φ
T
M
q (obtained by additional substitution of condition (7a)). These transformed boundary conditions are
(20)
η
(
0
)
=
Φ
T
M
q
0
,
η
˙
(
0
)
=
Φ
T
M
q
˙
0
,
η
(
t
f
)
=
Φ
T
M
q
f
,
η
˙
(
t
f
)
=
Φ
T
M
q
˙
f
.
For fully actuated problems, the last term (A
T
v
¨
-
Δ
A
T
v
˙
+
Ω
A
T
v) in optimality equations (18) vanishes because there are no constraints or Lagrange multipliers needed to enforce them. Therefore, a fully actuated problem involves only n
m optimality equations (18) to be solved in terms of n
m uncoupled modal variables in vector η.The solution to the combined set of (18), (19), and (20) can be efficiently obtained using symbolic differential operator D
n
=
d
n
/
d
t
n. Substituting this operator into (19) and (20) and rewriting in matrix notation give
(21)
[
E
-
E
^
T
E
~
0
]
[
η
ν
]
=
0
or
E
p
Y
=
0
,
where
(22)
E
=
R
^
D
4
+
(
2
R
^
Ω
-
Q
^
v
-
R
^
Δ
2
)
D
2
+
(
R
^
Ω
2
+
Q
^
d
)
,
E
^
=
A
(
I
D
2
+
Δ
D
+
Ω
)
,
E
~
=
A
(
I
D
2
-
Δ
D
+
Ω
)
.
Matrix E
p contains submatrices E, E
~, and -
E
^
T. Vector Y
=
[
η
T
∣
ν
T
]
T contains all unknown modal variables and Lagrange multipliers. Note that in a fully actuated case, matrix E
p in (21) consists only of submatrix E and vector Y
=
η.The solution to a system described in form (21) involves the roots r
l (l
=
1
,
…
,
4
n
m) of the characteristic equation for the determinant of  E
p [16], where operator D is replaced by the auxiliary variable r rendering a 4
n
m
th order polynomial. This operation is written as
(23)
det
E
p
|
D
→
r
=
0
.
Generally, the 4
n
m roots of the characteristic equation (23) take the following form:
(24)
r
l
=
±
α
k
±
i
β
k
(
k
=
1
,
…
,
n
m
,
l
=
1
,
…
,
4
n
m
)
.
The positive real numbers α
k and β
k characterize the response of the k
th mode of motion. For nonzero, unique roots, solution vector Y consists of n
m
+
n
r components Y
j that can be written in terms of 4
n
m independent elementary functions related to the roots (24), in the form [16]:
(25)
Y
j
=
∑
k
=
1
n
m
[
e
α
k
t
(
c
k
j
1
sin
(
β
k
t
)
+
c
k
j
2
cos
(
β
k
t
)
)
+
e
-
α
k
t
(
c
k
j
3
sin
(
β
k
t
)
+
c
k
j
4
cos
(
β
k
t
)
)
e
α
k
t
(
c
k
j
1
sin
(
β
k
t
)
+
c
k
j
2
cos
(
β
k
t
)
)
]
,
where
j
=
1
,
…
,
n
m
+
n
r
.
Obviously, the frequency of k
th mode controlled by the actuators can be interpreted as ω
k
a
=
β
k and its rate of active attenuation (or amplification) as ς
k
a
=
α
k
/
β
k. If multiple roots and zero-valued roots are obtained from (23), then solution functions (25) must be modified to mathematically accommodate these situations. There are 4
n
m
(
n
m
+
n
r
) unknown integration constants c
k
j
1
,
…
,
c
k
j
4 contained in the solution functions (25).Integration constantsc
k
j
1
,
…
,
c
k
j
4 are obtained by substituting the assumed form (25) into differential equations (18) and (19) and using the method of undetermined coefficients to generate n
m
+
n
r sets of 4
n
m linear algebraic equations relating these constants. By replacing one set of 4
n
m equations with the set of 4
n
m boundary conditions (20), the integration constants can be solved simultaneously. All these symbolic operations, including the determination of the roots (24) and constants in (25), can be done automatically using the MAPLE mathematical software.For closed-loop control, asymptotically convergent solution functions are required such that the control task is met over an infinite period of time (t
f
→
∞). The resulting number of integration constants is reduced by half, as terms involving positive exponential e
α
k
t in the solution form (25) disappear (c
k
j
1
=
c
k
j
2
=
0).To quantitatively measure the performance of closed-loop control schemes settling timet
f
3
% is defined as the time needed for various variables to be reduced to within 3% of their initial value (i.e., e
-
α
k
min
t
f
3
%
=
0.03
→
t
f
3
%
=
3.5
/
α
k
min).The above procedure was applied to actively suppress vibrations of a spatial antenna mast in [13] and of plane frames in [14], the cases with the number of DOFs much greater than the number of significant modes included in the analysis (i.e., with n
≫
n
m). In both cases only oscillating modes were controlled. Here the application of the above methodology is focused on various control schemes, which are demonstrated in controlling the translational and oscillating modes of a gantry crane.


## 3. Dynamics and Optimal Control of the Gantry Crane System

The gantry crane problem is one of the simplest underactuated mechanical systems involving two DOFs—cart translation and suspended load rotation—and a single actuator—a cart-driving force (n
m
=
2, n
a
=
n
r
=
1).The gantry crane model is shown in Figure1. The model includes the mass of the cart M, the mass of the suspended load m, swing angle θ, gravitational acceleration g, horizontal distance a from the cart’s initial position to the origin, and length L of the massless rigid link connecting the cart and load. The task is to manoeuvre the system from an initial resting state at some nonzero horizontal distance (x
=
a, θ
=
0) to a final resting equilibrium state at the origin (x
=
0, θ
=
0) by applying time-varying force F
a. Any finite cart translations are permitted, but swings of the suspended load are assumed to be sufficiently small for a linearized model to be valid. In modal space rigid body translation for such a manoeuvre is easily separated from the oscillatory motion of the suspended load. Dummy force F
d is added to artificially make the system fully actuated and formulate the augmented gantry crane system.Figure 1
Gantry crane system.This same gantry crane model was used in several papers dealing with control or/and optimization. Notably, a Lyapunov function was used in [1] to obtain an asymptotically stable (closed-loop) control (linear and nonlinear) for attenuating disturbances (nonzero initial positions) in the system, and optimal control by applying Pontryagin's principle was considered in [3]. Results for the linearized system are of interest because they serve as a useful comparison for the controls obtained in this paper. Similar problems of controlling the plane motion of gantry cranes were presented in [10, 11] using the concept of flatness. Various aspects of controlling gantry cranes, 3D operations were considered in [9, 17–19].The gantry crane system shown in Figure1 and its coordinate system are chosen to mimic those used in [1]. Matrices and vectors in the general equation of motion (1) take the following forms:
(26)
M
=
[
M
+
m
-
m
-
m
m
]
,
K
=
[
0
0
0
m
g
L
]
,
B
=
[
1
0
]
,
B
′
=
[
1
1
0
-
1
]
,
q
=
[
x
L
θ
]
,
F
a
'
=
[
F
a
F
d
]
.
To be consistent with the assumptions made in [1, 3] no dissipative effects (i.e., friction, etc.) are considered (C
=
0).The initial and final conditions (consistent with [1]) take the following forms:
(27)
q
(
0
)
=
[
a
0
]
T
,
q
˙
(
0
)
=
q
(
t
f
)
=
q
˙
(
t
f
)
=
[
0
0
]
T
.
The modal analysis (n
m
=
2) gives(28a)
Ω
=
[
0
0
0
(
1
+
m
M
)
g
L
]
,
(28b)
Φ
=
[
1
M
+
m
m
M
(
M
+
m
)
0
M
+
m
M
m
]
.The rigid body translational mode of motion is represented in (28a) by the frequency ω
1
=
0 and the second vibrating mode (load swinging) is represented by the frequency ω
2
=
(
1
+
m
/
M
)
(
g
/
L
).The uncoupled modal equations of motion (8) become:
(29)
η
¨
1
=
u
1
,
η
¨
2
+
ω
2
2
η
2
=
u
2
.
The augmented system transfer matrix B
^
′
=
Φ
T
B
′ is obtained by the appropriate substitutions from (26) and (28b) into the general partitioned form (9):
(30)
[
M
M
+
m
M
m
M
+
m
m
M
+
m
-
M
m
M
+
m
]
[
u
1
u
2
]
=
[
F
a
F
d
]
=
[
F
a
0
]
.
Modal controls u
1 and u
2 are considered actuated and unactuated, respectively. The (
n
r
=
1
) constraint equation is obtained by normalizing the bottom row of matrix (
B
^
′
)
-
1 in (30) to obtain
(31)
A
U
=
[
1
-
M
m
]
[
u
1
u
2
]
=
u
1
-
M
m
u
2
=
0
.
The constraint (31) may be applied to eliminate redundant modal control u
2 from the top row operation of (30) to obtain force F
a as a function of independent modal control u
1, giving
(32)
F
a
=
B
-
U
a
=
(
M
+
m
)
u
1
,
where B
-
=
M
+
m is the pseudotransfer matrix (since n
a
=
1, this matrix has only one term). Cart-driving force F
a may be applied using open-loop control (as a known function of time) or using closed-loop control through a set of gains in full-state feedback. Both schemes will be analyzed and simulated using the CMSOC method.The performance index (12) takes the following form, consisting of the gantry crane system’s four states (η
1
,
η
2
,
η
˙
1
,
η
˙
2) and two modal controls (u
1
,
u
2):
(33)
J
=
1
2
∫
0
t
f
(
Q
^
d
11
η
1
2
+
Q
^
d
22
η
2
2
+
Q
^
v
11
η
˙
1
2
+
Q
^
v
22
η
˙
2
2
+
R
^
11
u
1
2
+
R
^
22
u
2
2
)
d
t
⟶
min
.
The n
m
=
2 coupled optimality equations (18) take the following forms:(34a)
R
^
11
η
1
¨
¨
-
Q
^
v
11
η
¨
1
+
Q
^
d
11
η
1
-
ν
¨
=
0
,
(34b)
R
^
22
η
2
¨
¨
+
(
2
R
^
22
ω
2
2
-
Q
^
v
22
)
η
¨
2
+
(
R
^
22
ω
2
4
+
Q
^
d
22
)
η
2
+
M
m
(
ν
¨
+
ω
2
2
ν
)
=
0
,where ν is the Lagrange multiplier used to meet the n
r
=
1 constraint (31). The differential form (19) of constraint equation (31) is written as(35a)
η
¨
1
-
M
m
(
η
¨
2
+
ω
2
2
η
2
)
=
0
.In modal space, the boundary conditions (27) are(36a)
η
1
(
0
)
=
a
M
+
m
,
η
1
(
t
f
)
=
0
,
η
2
(
0
)
=
η
˙
1
(
0
)
=
η
˙
2
(
0
)
=
0
,
η
2
(
t
f
)
=
η
˙
1
(
t
f
)
=
η
˙
2
(
t
f
)
=
0
.Equations (34a), (34b), and (35a) written according to form (21) (with D
n
=
d
n
/
d
t
n) yield
(37)
E
p
Y
=
[
E
1
0
-
E
^
11
0
E
2
-
E
^
21
E
~
11
E
~
21
0
]
[
η
1
η
2
ν
]
=
0
,
where
(38)
E
1
=
R
^
11
D
4
-
Q
^
v
11
D
2
+
Q
^
d
11
,
E
2
=
R
^
22
D
4
+
(
2
R
^
22
ω
2
2
-
Q
^
v
22
)
D
2
+
(
R
^
22
ω
2
4
+
Q
^
d
22
)
,
E
^
11
=
E
~
11
=
D
2
,
E
^
21
=
E
~
21
=
-
M
m
(
D
2
+
ω
2
2
)
.
The characteristic equation of the system represented in (37) is obtained through operation (23), giving the 8th order polynomial equation:
(39)
det
E
p
|
D
→
r
=
E
1
E
^
21
2
+
E
2
E
^
11
2
|
D
→
r
=
M
m
(
R
^
11
r
4
-
Q
^
v
11
r
2
+
Q
^
d
11
)
(
r
2
+
ω
2
2
)
2
+
r
4
(
(
R
^
22
ω
2
4
+
Q
^
d
22
)
R
^
22
r
4
+
(
2
R
^
22
ω
2
2
-
Q
^
v
22
)
r
2
+
(
R
^
22
ω
2
4
+
Q
^
d
22
)
)
=
0
.
Eight roots may be obtained from the characteristic equation (39), which are then substituted into an appropriate assumed solution form (if the roots take the full complex form (24), then the assumed function takes form (25)) to characterize the three unknown solution functions (η
1
,
η
2
,
ν). This leaves twenty-four unknown integration constants to be determined by substituting the appropriate solution form into the three equations (34a), (34b), and (35a). By relating the coefficients corresponding to each of the eight independent elementary functions (i.e., in (25) each is in the form e
(
±
α
i
±
β
i
)
t), eight algebraic equations are obtained for each differential equation in the set (34a),(34b) and (35a), resulting in a total of twenty-four equations in terms of twenty-four unknown integration constants c
k
j
i. However, these twenty-four equations are linearly dependent. To obtain a unique solution, any one set of eight algebraic equations (obtained from either (34a), (34b), or (35a)) must be replaced with the eight boundary conditions (36a).The optimal actuation forces needed to drive the gantry crane from an initially disturbed position (x
=
a, θ
=
0) to the origin (x
=
0, θ
=
0) will be derived for four cases using the CMSOC method. These cases are as follows(A)
an open-loop control that minimizes actuation forces for a fixed time interval as in [3];
(B)
a closed-loop control that mimics the control presented in [1];
(C)
a closed-loop control with response improved over that presented in [1];
(D)
a closed-loop control of the fully actuated system (two actuators).For each case, the gantry crane’s physical parameters are chosen to match those given in [1]; namely, M
=
m
=
1
kg, L
=
1
m, g
=
9.8
m
/
s
2, a
=
-
5
m, and ω
2
=
4.43
s
-
1.As a final case, the CMSOC method is applied to a modified three-DOF gantry crane, with an additional link-mass hinge attached to the existing model in Figure1 and controlled by one or two actuators. This final case involves two subcases.(E1)
A closed-loop control that manoeuvres a modified gantry crane to the origin using the cart-driving force as well as a torque applied to the first rigid link (two actuators).
(E2)
An open-loop control that manoeuvres the modified gantry crane to the origin using only the cart-driving force (one actuator) over a fixed time interval.(A) Open-Loop Control of Gantry Crane Manoeuvre in a Finite Time Interval. The first control manoeuvres the gantry crane from a known initial position to the origin in a finite time interval t
f in an open-loop scheme. The performance index is chosen to be consistent with that presented in [3], corresponding to the weightings R
^
11
=
R
^
22
=
1 in the general form (33) with all other weightings null valued. Thus, the optimal control minimizes
(40)
J
=
∫
0
t
f
(
u
1
2
+
u
2
2
)
d
t
=
1
M
∫
0
t
f
F
a
2
d
t
⟶
min
.
Performance index (40) minimizes the modal controls or the actuation force over the finite manoeuvre time t
f, which is chosen as t
f
=
4
s to represent again one of the cases considered in [3]. The gantry crane’s characteristic polynomial equation (39) is simplified to
(41)
(
1
+
M
m
)
(
r
2
+
ω
2
2
)
2
r
4
=
0
.
The roots of (41) are r
1
,
…
,
8
=
0,0
,
0,0
,
±
i
ω
2
,
±
i
ω
2. There are four zero roots r
1
,
…
,
4
=
0, two imaginary roots r
5,7
=
i
ω
2, and two imaginary roots r
6,8
=
-
i
ω
2. When written in form (24), these roots correspond to β
2
=
ω
2 and α
1
=
β
1
=
α
2
=
0. Because of the zero roots and repeating roots, the solution functions take the following form:
(42)
Y
j
=
c
1
j
+
c
2
j
t
+
c
3
j
t
2
+
c
4
j
t
3
+
(
c
5
j
+
c
7
j
t
)
sin
(
ω
2
t
)
+
(
c
6
j
+
c
8
j
t
)
cos
(
ω
2
t
)
.
Each solution function (42) (j
=
1
,
…
,
3) contains eight unknown integration constants c
k
j (k
=
1
,
…
,
8), which are determined through substitution and comparison of similar terms in any two differential equations in the set (34a), (34b), and (35a) and by substitution of the eight boundary conditions (36a).  With the integration constants determined, the resulting solution functions are
(43)
η
1
=
-
7.09
-
.
104
t
+
1.41
t
2
-
.
235
t
3
+
.
0235
sin
(
4.43
t
)
+
.
0151
cos
(
4.43
t
)
,
(44)
η
2
=
.
144
-
.
0719
t
+
(
.
00451
-
.
0333
t
)
sin
(
4.43
t
)
+
(
-
.
144
+
.
0520
t
)
cos
(
4.43
t
)
,
(45)
ν
=
-
2.82
+
1.41
t
-
.
460
sin
(
4.43
t
)
-
.
295
cos
(
4.43
t
)
.
Substituting (43) into (29) yields
(46)
u
1
=
2.82
-
1.41
t
-
.
460
sin
(
4.43
t
)
-
.
295
cos
(
4.43
t
)
.
The Lagrange multiplier function v
(
t
) (45), which represents a “modal force” enforcing the modal constraints, is not used in further analysis and is shown here only for completeness of the solution.Mapping modal variablesη
1 and η
2 ((43) and (44)) into DOFs x and θ via transformation (6a), the trajectories shown in Figures 2(a) and 2(b) are obtained. Optimal force F
a, shown in Figure 2(c), is obtained by substituting modal control u
1 (46) into transformation (32).Histograms of (a) cart trajectoryx, (b) swing angle θ, and (c) force F
a (open loop).
(a)
x
=
1
2
(
η
1
+
η
2
)
(b)
θ
=
2
η
2
(c)
F
a
=
2
u
1This phase of the solution was done automatically using MAPLE. The solution procedure accepts any problem withn
m modes (obtained from FE analysis for more complex structures) controlled by n
a
≤
n
m actuators. The modal-to-DOF transformations for the gantry crane are indicated in Figures 2(a), 2(b), and 2(c). As shown, the open-loop control is able to perform the task in exactly four seconds, with a peak force of about 3.6 N and a maximum load swing angle of about 0.28 rad (16°). The optimal force accelerates the gantry crane over the first half of the manoeuvre (2 s) and decelerates the cart over the last half with identical, but opposite and mirrored, forces.Similar plots for the closed-loop control presented in [1] are shown in Figure 3. This control requires an effective manoeuvre time of t
f
3
%
≈
6
s to reach the origin, a maximum load rotation angle of 0.73 rad (42 deg), and a maximum force of 15 N. It should be noted that this relatively large rotation angle is mentioned here (and other angles quoted in the sequel) for the purpose of comparison only.Histograms of (a) cart trajectory, (b) swing angle, and (c) optimal force from [1].
(a)
(b)
(c)From Figure2 and Figure 3, one can conclude that the open-loop control performs the manoeuvre in a shorter period of time (t
f

=

4
s versus t
f
3
%
≈
6
s) with much smaller peak force requirements (3.6 N versus 15 N) and much smaller angles of oscillation (16° versus 42°). Also, the open-loop control brings the system to a complete stop after 4 s, while the closed-loop control produces overshoot and the system takes longer to effectively come to rest.Calculations show that if the finite manoeuvre time for the open-loop control is extended (or shortened), the peak force requirement and maximum swing angle are reduced (or increased)—approximately proportional tot
f
-
2. For example, if the open-loop control is modified to settle over the same effective period of time as that of the closed-loop control (t
f
≈
6), the maximum force is reduced to approximately 1.6 N with a maximum swing of about 7 deg.The open-loop control can always provide a faster and more efficient manoeuvre. However, it is possible only when the initial positions and manoeuvre times are known in advance. Closed-loop control is necessary if any initial configuration (unknown explicitly) is treated as disturbance, and its automatic reduction/removal is desired (the final position is at rest). Case (B) demonstrates how the CMSOC method is applied to analyze and simulate a closed-loop system that approximately produces the same dynamic responses as given in [1].(B) Closed-Loop Control of Gantry Crane: Reproducing Control from [1]. A closed-loop control can perform the same task as that of the open-loop control (case (A)); however it does so automatically, without prior knowledge of initial conditions involved. Any disturbance from its resting configuration at the origin (x
=
0, θ
=
0) is relayed through a set of constant gains to generate the cart-driving force F
a to attenuate this disturbance.In general, to simulate the closed-loop process analytically the manoeuvre timet
f is infinite and all parameters are driven asymptotically to the origin with increasing time. For the gantry crane, this requires that all roots of the characteristic equation (39) be nonzero complex numbers in the left half of the complex plane (unlike the open-loop system of case (A), which contained zero roots and purely imaginary roots). It can be verified that the weightings Q
^
d
11 and Q
^
d
22 in the performance index (33) must be nonzero in order to meet these criteria.The gantry crane control as given in [1] is closely reproduced by choosing the weightings in the performance index (33) equal to Q
^
v
11
=
Q
^
v
22
=
0, Q
^
d
11
=
4.5, Q
^
d
22
=
42, and R
^
11
=
R
^
22
=
1. The resulting characteristic polynomial equation (39) has eight roots that take form (24), with real and complex parts equal to
(47)
α
1
=
0.853
,
β
1
=
0.856
,
α
2
=
0.513
,
β
2
=
4.46
.
Note that the first actively controlled mode of frequency ω
1
a
=
0.856 s−1 (ω
1
=
0 for uncontrolled system) is damped with the ratio ς
1
a
=
0.996, while the second mode of frequency ω
2
a
=
4.46 s−1 (ω
2
=
4.43 s−1 for uncontrolled system) is damped with the ratio ς
1
a
=
0.115.Similar to case (A), modal variablesη
1 and η
2 are determined by substituting the parameters from (47) into the assumed solution function (25) and then solving for the unknown coefficients by comparing similar terms in two of the three optimality/constraint equations (34a), (34b) and (35a), and substituting the boundary conditions (36a). Unlike case (A), the closed-loop problem requires that only half as many integration constants must be solved because the coefficients preceding exponential growth functions (e
α
i
t) are assumed to be null valued. This gives
(48)
η
1
=
e
-
α
1
t
(
-
7.70
sin
(
β
1
t
)
-
7.08
cos
(
β
1
t
)
)
+
e
-
α
2
t
(
.
126
sin
(
β
2
t
)
-
.
00564
cos
(
β
2
t
)
)
,
η
2
=
e
-
α
1
t
(
-
.
565
sin
(
β
1
t
)
+
.
534
cos
(
β
1
t
)
)
+
e
-
α
2
t
(
.
149
sin
(
β
2
t
)
-
.
534
cos
(
β
2
t
)
)
,
u
1
=
e
-
α
1
t
(
-
10.3
sin
(
β
1
t
)
+
11.3
cos
(
β
1
t
)
)
+
e
-
α
2
t
(
-
2.44
sin
(
β
2
t
)
-
.
686
cos
(
β
2
t
)
)
.
Using the appropriate transformations (see Figure 2) the modal space variables (48) are mapped into the DOF space variables. The resulting system trajectories and the optimal force histogram are visually indistinguishable from those shown in Figure 3.The CMSOC method can also generate the closed-loop gains from the assumed weighting coefficients to demonstrate that the gains corresponding to the solution (48) are obtained and compared with the gains used in [1].In full-state feedback control the active force is a function of all system states in the following form:(49)
F
a
=
-
G
d
q
-
G
v
q
˙
.
For the general CMSOC method, gains G
d
=
[
g
1
d
⋯
g
n
m
d
] and G
v
=
[
g
1
v
⋯
g
n
m
v
] correspond to the observed positions and velocities of all n
m DOFs of a system. For the gantry crane, (49) takes the following form:
(50)
F
a
=
-
g
1
d
x
-
g
2
d
L
θ
-
g
1
v
x
˙
-
g
2
v
L
θ
˙
.
By substituting known DOF trajectories (x
=
(
1
/
2
)
(
η
1
+
η
2
) and θ
=
2
η
2) and the known force function (F
a
=
2
u
1) into (50) and grouping the terms related to the four independent elementary solution functions (operations are done in MAPLE automatically), one obtains: (51)
e
-
α
1
t
{
(
-
14.6
-
5.85
g
1
d
+
8.95
g
1
v
-
.
799
g
2
d
+
.
0352
g
2
v
)
sin
(
β
1
t
)
+
(
16.0
-
4.63
g
1
d
-
1.06
g
1
v
+
.
755
g
2
d
-
1.33
g
2
v
)
cos
(
β
1
t
)
}
+
e
-
α
2
t
{
(
-
3.45
+
.
195
g
1
d
+
1.56
g
1
v
+
.
211
g
2
d
+
3.26
g
2
v
)
sin
(
β
2
t
)
+
(
-
.
970
-
.
374
g
1
d
+
1.06
g
1
v
-
.
755
g
2
d
+
1.33
g
2
v
)
cos
(
β
2
t
)
}
=
0
.Each of the bracketed terms in (51) (containing the unknown gains) must equal to zero for the equation to be true at any time. This gives four equations in terms of four unknown gains, which may be solved to obtain
(52)
G
d
=
[
g
1
d
g
2
d
]
=
[
3.00
.
732
]
,
G
v
=
[
g
1
v
g
2
v
]
=
[
3.66
-
.
924
]
.
Though initial conditions were assumed in determining the trajectories x and θ and force F
a, it can be verified that gains (52) remain invariant towards any choice of these assumed conditions.The control gains used in [1] were
(53)
G
d
*
=
[
g
1
d
g
2
d
]
=
[
3.0
.
71
]
,
G
v
*
=
[
g
1
d
g
2
d
]
=
[
3.69
-
.
87
]
.
Comparing the gains (52) and (53) confirms that the CMSOC method is able to closely reproduce the closed-loop control in [1] by careful selection of the weighting parameters in performance index (33). However, as shown next in case (C), the performance of this closed-loop control may be improved through better selection of these weighting parameters to produce faster convergence without an increase in the required peak actuation forces.(C) Closed-Loop Control of Gantry Crane: Improving Performance. Case (B) developed a control that closely reproduced the control given in [1] by minimizing a performance index that gave no weight (Q
^
v
11
=
Q
^
v
22
=
0) to states η
˙
1 and η
˙
2, representing the gantry crane’s velocity. Consequently, the control caused the gantry to gain too much speed and then overshoot its target and produce large persistent load swings. These problems are mitigated by a more careful choice of the performance index weighting parameters in (33). To demonstrate the effect these parameters have on the gantry crane’s dynamics and to illustrate how they might be meaningfully selected, several cases, labelled P1 to P5 (each with different performance indices as listed in Table 1), are considered.Table 1
Weightings for five different performance indices in form (33).
Index #
Q
^
d
11
Q
^
d
22
Q
^
v
11
Q
^
v
22
R
^
11
R
^
22
P0
0.01
0.01
0
0
1
1
P1
6
0.01
0
0
1
1
P2
6
50
0
0
1
1
P3
6
50
4
0
1
1
P4
6
50
4
50
1
1Each case reflects a performance index which gives significant weightings to an incrementally increasing number of system states (of four possible statesη
1, η
2, η
˙
1, η
˙
2), while holding the weighting on both modal controls (u
1
,
u
2) at unity. Case P0 gives none of the states a significant weighting, case P1 gives a significant weighting to a single state (η
1), case P2 gives significant weightings to two states (η
1
,
η
2), and so on until case P4 significantly weights all four states. Table 1 summarizes how these weightings are chosen for each case. Since the gantry crane’s asymptotic convergence mathematically requires that weightings Q
^
d
11 and Q
^
d
22 in the index (33) are nonzero, a small value (0.01) is used instead of zero in cases P0 and P1 to demonstrate how the system behaves when these weightings are negligible. The DOF trajectories (x and θ) and force histogram (F
a) for the manoeuvres minimizing the performance indices for cases P0–P4 are presented in Table 2. The settling times t
f
3
% of the DOFs are also listed for each case. All plots in Table 2 are shown over the first 8 s of the manoeuvre period except for P0 (30 s).Table 2
DOF responses and force histograms that optimize performance indices P0–P4.
Index #
x (m)
θ (rad)
F
a (N)
P0
P1
P2
P3
P4Note that the first modal variableη
1 primarily influences the cart’s rigid body mode of motion, while the second modal variable η
2 influences the suspended load rotation. In fact there is a direct relationship between the angle of the load rotation and the second modal variable (θ
=
2
η
2) such that this angular trajectory is directly affected by varying the weights given to η
2
(
Q
^
d
22
) and its derivative η
˙
2
(
Q
^
v
22
) in the performance index (33). Likewise, the speed at which the cart can be made to reach its target is affected through the weightings given to η
1
(
Q
^
d
11
) and its derivative η
˙
1
(
Q
^
v
11
).The performance index in case P0 heavily weights the modal controlsu
1 and u
2 in comparison to modal variables η
1 and η
2 (100 times more) and neglects the modal velocities η
˙
1 and η
˙
2. The resulting control requires a small peak force (0.7 N), producing small maximum load swing angles (0.06 rad or 3.4 deg), but requires a very long manoeuvre time to converge to the origin (t
f
3
%
≈
440
s). If the weightings R
^
11, and R
^
22 were increased even further relative to the weighting Q
^
d
11, the maximum force requirements and angular rotations would become infinitesimal while the settling times would approach infinity.In case P1 a significant weighting value is given to the first modal variableη
1  (
Q
^
d
11
=
6
), while other weightings remain unchanged from case P0. This control greatly increases the speed at which the cart reaches its target position at x
=
0 (~2 s), but upon reaching this position the load undergoes large swing angles (1.0 rad or 57 deg) that persist for a very long time (t
f
3
%
≈
440
s).The maximum force increases significantly (17.3 N) in comparison to case P0 because the rigid body cart motion requires much larger accelerations during the initial 2 s of the manoeuvre in order to quickly attenuateη
1 due to its significant weighting value.Case P2 improves the load swing attenuation, which was poorly dampened in case P1, by including a large weighting value to the second modal variableη
2  (
Q
^
d
22
=
50
) (other weightings remain the same as in the previous case). The maximum load swing angle is reduced (0.8 rad or 46 deg) and the load swinging motion is damped much more quickly (~6.5 s). The cart translation requires similar accelerations and thus approximately the same maximum force (17.3 N) is needed. By inspection, one can see that case P2 produces similar behaviour to the control given in [1] shown in Figure 3 and likewise shares the problem of target overshoot and large persistent load swings.Case P3 reduces the tendency of the cart to overshoot the target by also giving a significant weighting to the first modal velocityη
˙
1  (
Q
^
v
11
=
4
). However, large persistent load swings are still present, and so convergence is not significantly improved over that produced in case P2. The maximum required force (17.3 N) remains essentially unchanged, while the load swing angles are reduced slightly (0.75 rad or 43 deg).The performance index in case P4 includes a large weight on the second modal velocityη
˙
2 (Q
^
v
22
=
50), while keeping all other weightings unchanged from case P3. This produces a control that reduces the magnitude of load swing angles (0.45 rad or 25.8 deg) while attenuating the swinging motion more quickly (t
f
3
%
≈
4.2
s). The gantry crane performs the manoeuvre in essentially a single load swing cycle, with similar initial cart accelerations and thus maximum forces (17.3 N) as in previous cases. Case P4 produces faster convergence then previous cases because, from an optimal control perspective, it incorporates all of the gantry crane’s states in the minimization by assigning all weightings in the performance index (33) with significant numerical values.It is essential that the second modal velocityη
˙
2 is given a significant weighting value to yield fast convergence, because the energy of the suspended load oscillates equally between potential and kinetic energy. Since potential energy and kinetic energy are proportional to the squares of displacement and velocity, respectively, both of the corresponding states η
2 and η
˙
2 should carry a significant and approximately equal weight in the performance index (33). Without weighting the load swing velocity state, the control focuses on eliminating swing angles but not swing velocities. However, when the load is near the bottom of its swing, its velocity is near maximum (η
˙
2
=
2
θ
˙
→
max), while its displacement is near minimum (η
2
=
2
θ
→
min). Therefore, the optimal force derived without consideration for the load swing velocity is unable to eliminate any significant portion of the load swing energy when the load is near the bottom of its swing (θ
→
0). Cases P2 and P3 failed to adequately weight η
˙
2, resulting in larger, more persistent load swings than in case P4.The control produced in case P4 provides a significant improvement over the control presented in [1], as it converges more quickly to the origin, while reducing load swing magnitudes, without any increase in the required maximum forces. To complete the design of this closed-loop control, the gains are obtained from (50) in a similar fashion as in case (B), giving
(54)
G
d
=
[
g
1
d
g
2
d
]
=
[
3.46
9.10
]
,
G
v
=
[
g
1
v
g
2
v
]
=
[
5.43
1.79
]
.
Note that the gains g
1
d and g
1
v are somewhat close to the gains for the control presented in [1] (53), but gains g
2
d and g
2
v are substantially different.(D) Closed-Loop Control of Gantry Crane: Fully Actuated Control. The CMSOC method can also be applied to fully actuated systems. To illustrate this, consider the same gantry crane system, now with both actuators F
a and F
d acting as real actuators (no dummy actuator). This situation may arise practically when a person is employed to guide the suspended load while the cart performs its translations.Since the problem is fully actuated there, are no additional constraints on the system motion and consequently no Lagrange multipliers needed to enforce them. The optimal forces can be solved by calculating the inverse dynamics directly from (6b), which takes form (30) (except F
d
≠
0), written as
(55)
F
=
B
^
-
1
U
⟹
[
F
a
F
d
]
=
[
M
M
+
m
M
m
M
+
m
m
M
+
m
-
M
m
M
+
m
]
[
u
1
u
2
]
.
The optimality equations in the differential operator form (21) become
(56)
E
p
Y
=
[
E
]
[
η
]
=
0
⟹
[
E
1
0
0
E
2
]
[
η
1
η
2
]
=
0
,
where
(57)
E
1
=
R
^
11
D
4
-
Q
^
v
11
D
2
+
Q
^
d
11
,
E
2
=
R
^
22
D
4
+
(
2
R
^
22
ω
2
2
-
Q
^
v
22
)
D
2
+
(
R
^
22
ω
2
4
+
Q
^
d
22
)
.
With weightings chosen according to the performance index in case P4 (Q
^
d
11
=
6, Q
^
d
22
=
50, Q
^
v
11
=
4, Q
^
v
22
=
50, and R
^
11
=
R
^
22
=
1), the roots of the characteristic equation (39) for the system given by (56) (E
1
E
2
|
D
→
r
=
0) take form (24) with the following real and imaginary parts:
(58)
α
1
=
3.62
,
β
1
=
2.78
,
α
2
=
1.49
,
β
2
=
0.474
.
For any fully actuated system, each modal variable η
i is independently controlled by a single modal control u
i, resulting in uncoupled solution functions of the following form:
(59)
η
i
=
e
α
i
t
(
c
i
1
sin
(
β
i
t
)
+
c
i
2
cos
(
β
i
t
)
)
.
For the gantry crane (i
=
1,2) the four unknown integration constants c
i
1,2 are obtained by substituting the four initial conditions for η
i
(
0
) and η
˙
i
(
0
) given by (36a). As in the previous cases, the solved modal variables in form (59) are mapped into the original coordinates to obtain the DOF trajectories and optimal forces. Figure 4 shows the cart trajectory x and the optimal forces on the cart and suspended load F
a and F
d, respectively. The angular trajectory θ of the load is not shown because it remains zero (θ
=
0) all the time. Practically, this means that for the optimal manoeuvre the person (actuator) guiding the suspended load must simply act to prevent it from swinging. Fast convergence (t
f
3
%
=
0.78
s) to the origin is obtained; however the task requires relatively large maximum forces (104 N) compared to previous cases. The required actions of cart-driving force F
a and suspended load guiding force F
d are identical, as the whole gantry crane system moves as a single rigid body. If smaller forces are desired, then a larger weight may be given to the modal controls (R
^
11
,
R
^
22) in the performance index. Note that only a 30 N maximum force would be required to execute the manoeuvre in 1 s by applying an open-loop control scheme.Histograms of (a) cart trajectory and (b) forces for the fully actuated gantry crane.
(a)
(b)(E1) Closed-Loop Control of Modified Three-DOF Gantry Crane (Two Actuators). To illustrate the application of the CMSOC method to a problem of a higher dimension (n
m
=
3), the gantry crane is modified by adding an additional link with an end load, as shown in Figure 5 (a case of a double-pendulum gantry crane in 3D was presented in [19]). In comparison to the previous cases considered the control task is unchanged except that the oscillations of the additional suspended load must also be damped. Consider the control that uses two actuators (n
a
=
2)—the standard cart-driving force F
a and torque T
a, produced by a motor fixed to the cart and applied to the first rigid link of length L
1 which supports the mass m
1. Dummy torque T
d to be used in formulating the augmented system is applied to the second link of length L
2 which carries mass m
2. All other physical variables are the same as in the original gantry crane model with the exception of θ
1 and θ
2, which denote the angles of links of lengths L
1 and L
2, respectively.Figure 5
Modified three-DOF gantry crane model.The standard matrices in the augmented system’s equation of motion (1) for this new model are
(60)
M
=
[
M
+
m
1
+
m
2
-
m
1
-
m
2
-
m
2
-
m
1
-
m
2
m
1
+
m
2
m
2
-
m
2
m
2
m
2
]
,
K
=
[
0
0
0
0
g
(
m
1
+
m
2
)
L
1
0
0
0
g
m
2
L
2
]
,
B
′
=
[
1
0
0
0
1
L
1
-
1
L
1
0
0
1
L
2
]
.
The augmented system consists of DOF vector q
=
[
x
L
1
θ
1
L
2
θ
2
]
T and force vector F
a
'
=
[
F
a
T
a
T
d
]
T. The following numerical values are adopted: M
=
m
1
=
m
2
=
1
kg, L
1
=
L
2
=
1
m, a
=
-
5
m, and g
=
9.8
m
/
s
2.The set ofn
m
=
3 equations of motion (1) are uncoupled in modal space, with matrices of ordered frequencies Ω and mode shapes φ normalized according to (7a) and (7b), taking the following forms:
(61)
Ω
=
[
0
0
0
0
12.4
0
0
0
46.4
]
,
Φ
=
[
.
577
.
577
.
577
0
.
366
1.37
0
1
-
1
]
.
As before, the first mode represents the rigid body mode of motion (ω
1
2
=
0), while the second and third modes, with the squared frequencies ω
2
2
=
12.43
(
rad
/
s
)
2 and ω
3
2
=
46.37
(
rad
/
s
)
2, represent the swinging modes of the rotating link-masses.Augmented force vectorF
a
' is related to modal control vector U through the inverse of transformation (6b) which is partitioned according to (9) to give
(62)
B
^
-
1
U
=
[
1.73
0
0
-
1.73
1.37
.
366
-
.
577
.
789
-
.
211
]
[
u
1
u
2
u
3
]
=
[
F
a
T
a
T
d
]
.
The constraint equation (n
r
=
1) is obtained from the bottom row of (62) (T
d
=
0) and normalized into the following form:
(63)
A
U
=
[
1
-
1.37
.
366
]
[
u
1
u
2
u
3
]
=
0
.
The n
a
=
2 actuation forces may be obtained directly from the top two rows of (62) in terms of all modal controls, but according to (11) these forces may be expressed in terms of two independent modal controls (chosen as u
1 and u
2) in the following form:
(64)
B
-
U
a
=
[
1.73
0
-
2.73
2.73
]
[
u
1
u
2
]
=
[
F
a
T
a
]
.
Selecting a performance index of form (12) gives three (n
m
=
3) optimality equations in the form (18) that, with the constraint equation (63), may be written according to (21) in the following form:
(65)
E
p
Y
=
[
E
1
0
0
E
^
11
0
E
2
0
E
^
21
0
0
E
3
E
^
31
E
~
11
E
~
21
E
~
31
0
]
[
η
1
η
2
η
3
ν
]
=
0
,
where
(66)
E
i
=
R
^
i
i
D
4
-
(
2
R
^
i
i
ω
i
2
-
Q
^
v
i
i
)
D
2
+
(
R
^
i
i
ω
i
4
+
Q
^
d
i
i
)
,
E
^
i
j
=
A
j
i
(
D
2
+
ω
i
2
)
(
i
=
1
,
…
,
3
,
j
=
1
)
.
The parameters A
j
i in the equation above are the j
th row and i
th column components of the constraint matrix A given by (63). The selected weightings for the performance index are Q
^
d
11
=
6, Q
^
d
22
=
Q
^
d
33
=
50, Q
^
v
11
=
4, Q
^
v
22
=
Q
^
v
33
=
50, and R
11
=
R
22
=
R
33
=
1.The twelve (4
n
m) roots of the characteristic equation (23) are obtained in the complex form (24) and are used to generate an assumed solution of form (25) for each unknown modal variable (n
m
=
3) and Lagrange multiplier (n
r
=
1). Half of these roots generate exponential growth functions that are eliminated by assuming their corresponding integration constants to be zero-valued. Then through the method of undetermined coefficients, four (n
m
+
n
r) sets of six linear algebraic equations are obtained. Replacing one set by the set of six initial conditions, the unknown integration constants are obtained by solving the set of twenty-four equations (the number of equations is 2
n
m
(
n
m
+
n
r
)).The boundary conditions for this problem are the same as those chosen for the original gantry crane, written in (36a), with the additional condition that the initial and final positions and velocities of the third modal variable η
3 are also zero. In other words, the manoeuvre requires a horizontal cart translation from a resting position at x
=
a
=
-
5
m with both links hanging vertically to the same resting position at the origin.Figure6 shows trajectories x, θ
1, and θ
2 of the three-DOF gantry crane as well as required actuation forces F
a and T
a.(a) Cart trajectory, (b) swing angles, and (c) force/torque for the modified gantry crane.
(a)
(b)
(c)The manoeuvre, requiring a maximum force of 29 N and a maximum torque of 18 Nm, is effectively completed aftert
eff
3
%
=
4
s. The maximum load swing angle of the first link is 0.11 rad (6.3°) and that of the second link is 0.48 rad (27.5 deg).(E2) Open-Loop Control of Modified Three-DOF Gantry Crane (One Actuator). In order to show the case where one actuator controls three DOFs, the optimal manoeuvre for the modified gantry crane using only a single actuator—cart-driving force F
a—is investigated for an open-loop scheme. Both of the torque actuators T
a and T
d (Figure 5) are treated as dummy actuators and so the inverse transformation, while identical to (62), is repartitioned in the following form:
(67)
B
^
-
1
U
=
[
1.73
0
0
-
1.73
1.37
.
366
-
.
5774
.
789
-
.
211
]
[
u
1
u
2
u
3
]
=
[
F
a
T
a
T
d
]
=
[
F
a
0
0
]
.
The constraint equations (n
r
=
2) are obtained from the two bottom rows of (67) and normalized into the following form:
(68)
A
U
=
[
1
-
1
0
0
1
1
]
[
u
1
u
2
u
3
]
=
0
.
According to (11) the single cart-driving force may be expressed in terms of the independent modal control (chosen as u
1) in the following form:
(69)
B
-
U
a
=
[
1.73
]
[
u
1
]
=
[
F
a
]
.
Note that matrix B
- is the same for all cases involving one actuator. Choosing a performance index in form (12) gives three (n
m
=
3) optimality equations in form (18) that, with constraint equation (68), may be written according to (21) in the following form:
(70)
E
p
Y
=
[
E
1
0
0
E
^
11
E
^
12
0
E
2
0
E
^
21
E
^
22
0
0
E
3
E
^
31
E
^
32
E
^
11
E
^
21
E
^
31
0
0
E
^
12
E
^
22
E
^
32
0
0
]
[
η
1
η
2
η
3
ν
1
ν
2
]
=
0
,
where
(71)
E
i
=
R
^
i
i
D
4
-
(
2
R
^
i
i
ω
i
2
-
Q
^
v
i
i
)
D
2
+
(
R
^
i
i
ω
i
4
+
Q
^
d
i
i
)
,
E
^
i
j
=
A
j
i
(
D
2
+
ω
i
2
)
(
i
=
1
,
…
,
3
,
j
=
2
)
.
The assumed weightings are R
11
=
R
22
=
R
33
=
1 and Q
^
d
11
=
Q
^
d
22
=
Q
^
d
33
=
Q
^
v
11
=
Q
^
v
22
=
Q
^
v
33
=
0 (only the control effort is to be minimized). Consistent with the open-loop control presented in case (A) the finite manoeuvre time is chosen to be t
f
=
4
s.The solution procedure is similar to previous examples. Figure7 shows trajectories x, θ
1, and θ
2 as well as the required cart-driving force F
a. The manoeuvre requires a peak force of 4.8 N and completes the task in exactly 4 s. The maximum load swing angle of the first link is 0.19 rad (11°) and that of the second link is 0.35 rad (20°).(a) Cart trajectory, (b) swing angles, and (c) force for modified gantry crane (open loop).
(a)
(b)
(c)


## 4. Conclusions

The CMSOC methodology was presented as a means of solving linear underactuated (or fully actuated) control problems. The gantry crane problem was selected to illustrate in detail various operations required for different control methodologies. As demonstrated the method can be applied to open-loop control schemes as well as closed-loop (asymptotically convergent) control schemes. In the latter case the weightings of the performance index can be translated to the gains of the full-state feedback closed-loop controllers. The operations would be identical for any similar problems with larger numbers of modes and actuators.


## Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.



---

*Source: 283565-2014-01-12.xml*