Invariant-Based Inverse Engineering for Fast and Robust Load Transport in a Double Pendulum Bridge Crane

Abstract

We set a shortcut-to-adiabaticity strategy to design the trolley motion in a double-pendulum bridge crane. The trajectories found guarantee payload transport without residual excitation regardless of the initial conditions within the small oscillations regime. The results are compared with exact dynamics to set the working domain of the approach. The method is free from instabilities due to boundary effects or to resonances with the two natural frequencies.

1. Introduction

The concept of adiabaticity is ubiquitous in physics, but it is not fully exploited in mechanical engineering and control applications. Adiabatic theorems set the existence of approximate adiabatic invariants, such as the action integral in classical mechanics, when the control parameters of a given physical system vary slowly enough in time [1].

Adiabaticity is often used to drive systems in a robust manner. An example is a load hanging as a simple pendulum from a moving trolley on a bridge crane. If the trolley travels slowly enough between two points, the energy of the pendulum is an adiabatic invariant and stays constant along different smooth trolley trajectories for the same initial and final points. In particular, the minimum energy configuration, in which the oscillating mass stays at relative rest with respect to the suspension point, is preserved. More generally, for other initial states the final energy will not suffer excitations. However, the intrinsic slowness of such processes may be problematic, either because long operation times are impractical, or because during a long process time the ideal dynamics can be affected by the accumulation of random and/or uncontrollable perturbations that spoil the desired result.

To overcome these problems, “Shortcuts To Adiabaticity” (STA) methods have been developed in the last decade. The idea is to reach the same results of an adiabatic protocol in short times [2,3]. In STA, the adiabatic invariant is not kept constant throughout the process, but the initial value is recovered at final time. For the simple example of the load hanging from a moving trolley, the shortcuts are certain special and fast driving trajectories of the trolley that induce transitory excitations, but leave the load at final time with the same energy it had initially.

STA methods have been succesfully applied to many different fields and processes in quantum systems, such as quantum computation [4,5,6,7], cooling [8], quantum transport [9,10], quantum state preparation [11,12,13,14], manipulation of cold atoms [15,16,17,18,19,20] or control of polyatomic molecules [21]. They have been also applied to design optical devices [22,23], and recently in mechanical engineering to design fast and robust protocols to control overhead cranes [24,25]. Perhaps surprisingly, because of the differing orders of magnitude involved, the physics of crane control are much related, in some formal aspects and domains even identical, to the physics of microscopic particle transport in moving traps [9]. In both domains the linearized models imply a moving harmonic oscillator. When the setting is more realistic though, beyond the simplest scenarios, the models become specific and require specialized treatment as in the system addressed in this paper, see Figure 1, a planar, double-pendulum, hook ($m_1$) and load ($m_2$) system suspended from a moving trolley. This is a relevant model as cranes behave like moving double pendulums due to different reasons, for example the large scale of the payload, or weighty hooks [26,27]. The control of a moving double pendulum dynamics is significantly more difficult than the single pendulum, with two unactuated degrees of freedom (angles ${\theta }_1$ and ${\theta }_2$ in Figure 1) and only one actuator (the trolley position x). Compared to studies on single pendulum cranes, this system is much less explored, for a recent brief review on recent papers and approaches applied see Reference [27]. Control approaches with and without feedback have been worked out and their pros and cons have been well discussed [26,27,28]. Our STA approach is presented here in an elementary way without feedback but it may be adapted and incorporated into methods with feedback as well.

Among the different STA approaches, dynamical invariant based inverse engineering is one of the most successful and is the one followed here. The essence of the method is to identify exact (rather than adiabatic) dynamical invariants, set boundary conditions to cancel final excitation, design the dynamics compatible with these conditions, and deduce the necessary controls from the dynamics thanks to the relations between dynamical invariants and Hamiltonian. For the moving double pendulum the STA consists on designing trolley trajectories $x\left(t\right)$ from $x=0$ to $x=d$ so that the system ends up at final process time $t_f$ without excitations. From the point of view of STA process design, this system poses interesting, non-trivial challenges with respect to the single pendulum, as we shall see.

The article is organized as follows. The physical model and Hamiltonian of the system are set in Section 2, both in exact form and in the small oscillation regime. Dynamically decoupled normal modes are found in Section 3, and then the STA protocol is designed in Section 4. Numerical results are presented in Section 5 and, finally, in Section 6 we end with the conclusions and discuss some open questions.

2. Physical Model

The physical model and relevant parameters are shown in Figure 1. The model assumes several conditions and idealizations: (i) the mass of the wires and friction are neglected; (ii) point masses; (iii) constant wire lengths $l_1$ and $l_2$; (iv) the trolley position is treated as a control parameter rather than a dynamical variable. This last assumption is a common and simplifying assumption [29] that requires a good controller, but a more fundamental approach considering the trolley as a dynamical variable is also possible as in Reference [24].

2.1. Lagrangian

In terms of the angles ${\theta }_1$ and ${\theta }_2$, see Figure 1, the Cartesian coordinates of each mass in a rest frame are given by

$$\begin{array}{ccc}\hfill x_1=x+l_{1}sin{\theta }_1,& & y_1=-l_{1}cos{\theta }_1,\hfill \\ \hfill x_2=x_1+l_{2}sin{\theta }_2,& & y_2=y_1-l_{2}cos{\theta }_2,\hfill \end{array}$$

(1)

so kinetic $\left(T\right)$ and potential $\left(V\right)$ energies are given by

$$\begin{array}{ccc}\hfill T& =& \frac{1}{2}{m}_1({\dot{x}}_1^2+{\dot{y}}_1^2)+\frac{1}{2}{m}_2({\dot{x}}_2^2+{\dot{y}}_2^2),V=m_{1}g{y}_1+m_{2}g{y}_2,\hfill \end{array}$$

(2)

where the dots represent time derivatives. The Lagrangian of the system using ${\theta }_1$ and ${\theta }_2$ as generalized coordinates and dynamical variables and $x\left(t\right)$ as a control parameter will be given by

$$\begin{array}{ccc}\hfill L& =& L({\theta }_i,{\dot{\theta }}_i;t)=T-V.\hfill \end{array}$$

(3)

To avoid deformations or excessive tensions, cranes usually work in the small oscillations regime, in which ${\theta }_i$ are small so that we may approximate $sin{\theta }_i\approx {\theta }_i$ and $cos{\theta }_i\approx 1-{\theta }_i^2/2$. Angular velocities $\dot{\theta }$ will be considered small as well. This approximation linearizes the dynamical equations of motion of the system. Results found with exact and approximate dynamics will be compared later to check the validity of the approximation and its limits.

In this small oscillation regime and keeping up to second order quadratic terms in ${\theta }_i$ and ${\dot{\theta }}_i$, kinetic and potential energies are given by

$$\begin{array}{ccc}\hfill T& \approx & \frac{1}{2}M{\dot{x}}^2+\frac{1}{2}({\dot{\theta }}_1,{\dot{\theta }}_2)\left(\begin{array}{cc}M{l}_1^2& m_2{l}_1{l}_2\\ m_2{l}_1{l}_2& m_2{l}_2^2\end{array}\right)\left(\begin{array}{c}{\dot{\theta }}_1\\ {\dot{\theta }}_2\end{array}\right)+\dot{x}\left(M{l}_1,m_2{l}_2\right)\left(\begin{array}{c}{\dot{\theta }}_1\\ {\dot{\theta }}_2\end{array}\right),\hfill \\ \hfill V& \approx & -Mg{l}_1\left(1-\frac{{\theta }_1^2}{2}\right)-m_{2}g{l}_2\left(1-\frac{{\theta }_2^2}{2}\right),\hfill \end{array}$$

(4)

where M denotes the total mass $M=m_1+m_2$. The Lagrangian becomes

$$\begin{array}{ccc}\hfill L\approx T-V& =& \frac{1}{2}({\dot{\theta }}_1,{\dot{\theta }}_2)\left(\begin{array}{cc}M{l}_1^2& m_2{l}_1{l}_2\\ m_2{l}_1{l}_2& m_2{l}_2^2\end{array}\right)\left(\begin{array}{c}{\dot{\theta }}_1\\ {\dot{\theta }}_2\end{array}\right)-\frac{1}{2}({\theta }_1,{\theta }_2)\left(\begin{array}{cc}Mg{l}_1& 0\\ 0& m_{2}g{l}_2\end{array}\right)\left(\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right)\hfill \\ & +& \dot{x}\left(M{l}_1,m_2{l}_2\right)\left(\begin{array}{c}{\dot{\theta }}_1\\ {\dot{\theta }}_2\end{array}\right),\hfill \end{array}$$

(5)

where purely time-dependent and constant terms have been omitted since they do not affect the dynamics.

2.2. Hamiltonian

To implement a Hamiltonian formulation, which is more convenient to treat the invariants and inverse engineering of trolley trajectories, we need the conjugate momentum of each ${\theta }_i$,

$$\begin{array}{ccc}\hfill p_{{\theta }_1}& =& \frac{\partial L_{\theta }}{\partial {\dot{\theta }}_1}=M{l}_1\dot{x}+M{l}_1^2{\dot{\theta }}_1+m_2{l}_1{l}_2{\dot{\theta }}_2,\hfill \\ \hfill p_{{\theta }_2}& =& \frac{\partial L_{\theta }}{\partial {\dot{\theta }}_2}=m_2{l}_2\dot{x}+m_2{l}_1{l}_2{\dot{\theta }}_1+m_2{l}_2^2{\dot{\theta }}_2.\hfill \end{array}$$

(6)

These relations can be inverted to have the generalized velocities ${\dot{\theta }}_i$ in terms of the generalized momenta $p_{{\theta }_i}$. The Hamiltonian is found from the Lagrangian as

$$\begin{array}{ccc}\hfill H_{\theta }& =& \sum _{i=1}^2{\dot{\theta }}_i{p}_{{\theta }_i}-L=\frac{p_{{\theta }_1}^2}{2m_1{l}_1^2}+\frac{p_{{\theta }_2}^2}{2\mu l_2^2}-\frac{p_{{\theta }_1}{p}_{{\theta }_2}}{m_1{l}_1{l}_2}+\frac{1}{2}Mg{l}_1{\theta }_1^2+\frac{1}{2}{m}_{2}g{l}_2{\theta }_2^2-\dot{x}\frac{p_{{\theta }_1}}{l_1},\hfill \end{array}$$

(7)

where $\mu =\frac{m_1{m}_2}{m_1+m_2}$ is the reduced mass and where constant terms that do not affect the dynamics have been neglected. In matrix representation, this Hamiltonian can be written as

$$\begin{array}{ccc}\hfill H_{\theta }& =& \frac{1}{2}(p_{{\theta }_1},p_{{\theta }_2})T\left(\begin{array}{c}{p}_{{\theta }_1}\\ p_{{\theta }_2}\end{array}\right)+\frac{1}{2}({\theta }_1,{\theta }_2)K\left(\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right)-\left(\frac{\dot{x}}{l_1},0\right)\left(\begin{array}{c}{p}_{{\theta }_1}\\ p_{{\theta }_2}\end{array}\right),\hfill \end{array}$$

(8)

where

$$T=\left(\begin{array}{cc}\frac{1}{m_1{l}_1^2}& \frac{-1}{m_1{l}_1{l}_2}\\ \frac{-1}{m_1{l}_1{l}_2}& \frac{m_1+m_2}{m_1{m}_2{l}_2^2}\end{array}\right);K=\left(\begin{array}{cc}Mg{l}_1& 0\\ 0& m_{2}g{l}_2\end{array}\right).$$

Whereas the potential matrix K is diagonal, the kinetic matrix T is not, i. e., $p_{{\theta }_1}$ and $p_{{\theta }_2}$ momenta are coupled. We want to find a coordinate transformation, i. e., normal modes, where both the coordinates and momenta are uncoupled so that we can easily get the dynamical invariants to inverse engineer $x\left(t\right)$. In the following section, these normal modes will be calculated following Reference [30]. Normal modes for the double pendulum with fixed suspension point are known [31], but our treatment takes the motion of the trolley into account. Finding dynamical normal modes for quadratic time-dependent Hamiltonians is generically non-trivial [30], but in this system the task is facilitated by the fact that the time-dependence appears in linear terms via $\dot{x}\left(t\right)$.

3. Normal Modes

3.1. Diagonalization of $H_{\theta }$

Let us first define a new set of coordinates $u_1$ and $u_2$ by the linear transformation

$$\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=A\left(\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right),$$

(9)

where the A matrix is yet to be determined. The corresponding momenta transform according to Reference [30]

$$\left(\begin{array}{c}{p}_{u_1}\\ p_{u_2}\end{array}\right)=A^{-T}\left(\begin{array}{c}{p}_{{\theta }_1}\\ p_{{\theta }_2}\end{array}\right),$$

(10)

where $A^{-T}={\left(A^{-1}\right)}^T$ stands for the transpose of the inverse matrix. The Hamiltonian in these variables reads

$$\begin{array}{ccc}\hfill H_u& =& \frac{1}{2}(p_{u_1},p_{u_2})\left(AT{A}^T\right)\left(\begin{array}{c}{p}_{u_1}\\ p_{u_2}\end{array}\right)+\frac{1}{2}(u_1,u_2)\left(A^{-T}K{A}^{-1}\right)\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)-\left(\frac{\dot{x}}{l_1},0\right)A^T\left(\begin{array}{c}{p}_{u_1}\\ p_{u_2}\end{array}\right)\hfill \end{array}$$

(11)

We now look for a transformation matrix A that diagonalizes simultaneously both the $AT{A}^T$ and $A^{-T}K{A}^{-1}$ matrices in the expression above. To do so it is useful to define the following matrix

$$\tilde{T}=K^{1/2}T{K}^{1/2},$$

(12)

which is symmetric and positive definite and therefore can be diagonalized by an orthogonal matrix $\mathcal{O}$. Without loss of generality, this orthogonal matrix $\mathcal{O}$ can be parametrized as

$$\mathcal{O}=\left(\begin{array}{cc}cos\theta & -sin\theta \\ sin\theta & cos\theta \end{array}\right),$$

(13)

and choosing the parameter $\theta$ (not to be confused with the angles ${\theta }_i$) by

$$tan2\theta =\frac{2l_1}{(l_1-l_2)}\sqrt{\frac{m_2{l}_2}{M{l}_1}},$$

(14)

we have that

$${\mathcal{O}}^T\tilde{T}\mathcal{O}=\mathrm{diag}({\omega }_1^2,{\omega }_2^2)=T_d.$$

(15)

The eigenvalues ${\omega }_i^2$ are positive since $\tilde{T}$ is a positive definite matrix, and have dimensions of (angular) frequency square. The explicit expressions are

$${\omega }_1^2=\frac{g}{m_1{l}_1{l}_2}\left(-\sqrt{M{m}_2{l}_1{l}_2}sin2\theta +M{l}_1{sin}^2\theta +M{l}_2{cos}^2\theta \right),$$

$${\omega }_2^2=\frac{g}{m_1{l}_1{l}_2}\left(\sqrt{M{m}_2{l}_1{l}_2}sin2\theta +M{l}_1{cos}^2\theta +M{l}_2{sin}^2\theta \right),$$

in agreement with the eigenfrequencies given in Reference [31]. Now, by writing the transformation matrix as

$$\begin{array}{ccc}\hfill A& =& {\mathcal{O}}^T{K}^{1/2}=\left(\begin{array}{cc}\sqrt{Mg{l}_1}cos\theta & \sqrt{m_{2}g{l}_2}sin\theta \\ -\sqrt{Mg{l}_1}sin\theta & \sqrt{m_{2}g{l}_2}cos\theta \end{array}\right),\hfill \end{array}$$

(16)

both quadratic terms in the transformed Hamiltonian (11) are diagonal since

$$\begin{array}{ccc}\hfill AT{A}^T& =& {\mathcal{O}}^T{K}^{1/2}T{K}^{1/2}\mathcal{O}={\mathcal{O}}^T\tilde{T}\mathcal{O}=T_d,\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill A^{-T}K{A}^{-1}& =& {\mathcal{O}}^T{K}^{-1/2}K{K}^{-1/2}\mathcal{O}=1.\hfill \end{array}$$

(18)

Finally, the Hamiltonian (11) takes the uncoupled form

$$\begin{array}{ccc}\hfill H_u& =& \frac{1}{2}\sum _{i=1}^2\left({\omega }_i^2{p}_{u_i}^2+u_i^2\right)+\dot{x}\sqrt{\frac{Mg}{l_1}}(-p_{u_1}cos\theta +p_{u_2}sin\theta ).\hfill \end{array}$$

(19)

3.2. Lewis-Leach Family of Hamiltonians and Second Canonical Transformation

The Lewis-Leach (LL) family of Hamiltonians are of the form [32]

$$H_{LL}=\frac{1}{2}\left[p^2+\mathsf{\Omega }\left(t\right)q^2\right]-F\left(t\right)q,$$

(20)

i.e., quadratic Hamiltonians with linear in position terms. For them quadratic invariants are explicitly known. By a suitable canonical transformation to some generalized coordinates $\{q_i,p_i\}$, we shall transform $H_u$ into this form. This can be easily achieved just by exchanging momentum and coordinate [33]. The transformation is generated by $F_1=u_1{q}_1+u_2{q}_2$ which gives the new coordinates and momenta in terms of the old ones as follows,

$$\begin{array}{c}\hfill p_{u_i}=\frac{\partial F_1}{\partial u_i}=q_i,p_i=-\frac{\partial F_1}{\partial q_i}=-u_i.\end{array}$$

(21)

By using this canonical transformation, the new Hamiltonian is

$$\begin{array}{ccc}\hfill H_q& =& H_1+H_2,\hfill \end{array}$$

(22)

a sum of two independent forced harmonic oscillators that belong to the LL family,

$$\begin{array}{ccc}\hfill H_1& =& \frac{1}{2}\left(p_1^2+{\omega }_1^2{q}_1^2\right)-q_1\dot{x}\sqrt{\frac{Mg}{l_1}}cos\theta ,\hfill \\ \hfill H_2& =& \frac{1}{2}\left(p_2^2+{\omega }_2^2{q}_2^2\right)+q_2\dot{x}\sqrt{\frac{Mg}{l_1}}sin\theta .\hfill \end{array}$$

(23)

3.3. Explicit Expression of Normal Mode Coordinates

Taking into account the two canonical transformations, the explicit expression of the normal mode coordinates and momenta $\{q_i,p_i\}$ in terms of the original variables $\{{\theta }_i,p_{{\theta }_i}\}$ is

$$\begin{array}{ccc}\hfill \left(\begin{array}{c}{q}_1\\ q_2\\ p_1\\ p_2\end{array}\right)& =& \left(\begin{array}{cc}0& {\mathcal{I}}_2\\ -{\mathcal{I}}_2& 0\end{array}\right)\left(\begin{array}{cc}A& 0\\ 0& A^{-T}\end{array}\right)\left(\begin{array}{c}{\theta }_1\\ {\theta }_2\\ p_{{\theta }_1}\\ p_{{\theta }_2}\end{array}\right),\hfill \end{array}$$

(24)

where ${\mathcal{I}}_2$ is the $2\times 2$ identity matrix and, using the explicit expression of A in (16), we have

$$\begin{array}{ccc}\hfill q_1& =& \frac{cos\theta }{\sqrt{Mg{l}_1}}{p}_{{\theta }_1}+\frac{sin\theta }{\sqrt{m_{2}g{l}_2}}{p}_{{\theta }_2},\hfill \\ \hfill q_2& =& -\frac{sin\theta }{\sqrt{Mg{l}_1}}{p}_{{\theta }_1}+\frac{cos\theta }{\sqrt{m_{2}g{l}_2}}{p}_{{\theta }_2},\hfill \\ \hfill p_1& =& -cos\theta \sqrt{Mg{l}_1}{\theta }_1-sin\theta \sqrt{m_{2}g{l}_2}{\theta }_2,\hfill \\ \hfill p_2& =& sin\theta \sqrt{Mg{l}_1}{\theta }_1-cos\theta \sqrt{m_{2}g{l}_2}{\theta }_2.\hfill \end{array}$$

(25)

4. Designing the STA Protocol

We are now ready to define the invariants and design the driving function $x\left(t\right)$.

4.1. Dynamical Invariants

A dynamical invariant of a Hamiltonian system remains constant during the time evolution [34]. Labelling the dynamical invariant of the Hamiltonian $H_i$ as $I_i$ we have that

$$\begin{array}{ccc}\hfill \frac{d{I}_i}{dt}& =& {\partial }_t{I}_i+\{I_i,H_i\}=0,\hfill \end{array}$$

(26)

with $\{I_i,H_i\}$ being the Poisson bracket. The sum of invariants $I=I_1+I_2$ is invariant with respect to the sum of Hamiltonians $H_q=H_1+H_2$ since

$$\begin{array}{ccc}\hfill \frac{dI}{dt}& =& \{I,H_q\}+{\partial }_{t}I=\{I_1+I_2,H_1+H_2\}+{\partial }_t(I_1+I_2)\hfill \\ & =& \left(\{I_1,H_1\}+{\partial }_t{I}_1\right)+\left(\{I_2,H_2\}+{\partial }_t{I}_2\right)+\{I_1,H_2\}+\{I_2,H_1\}=0.\hfill \end{array}$$

The invariants for (23) have the explicit form [32]

$$\begin{array}{ccc}\hfill I_i& =& \frac{1}{2}{\left(p_i-{\dot{\alpha }}_i\right)}^2+\frac{{\omega }_i^2}{2}{\left(q_i-{\alpha }_i\right)}^2,\hfill \end{array}$$

(27)

provided the functions ${\alpha }_i$ satisfy the following Newton equations,

$$\begin{array}{ccc}\hfill {\ddot{\alpha }}_1+{\omega }_1^2{\alpha }_1& =& \dot{x}\sqrt{\frac{Mg}{l_1}}cos\theta ,\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill {\ddot{\alpha }}_2+{\omega }_2^2{\alpha }_2& =& -\dot{x}\sqrt{\frac{Mg}{l_1}}sin\theta .\hfill \end{array}$$

(29)

These ${\alpha }_i$ functions may be regarded as auxiliary, reference, special “displacements” in two forced harmonic oscillators. Let us underline that the actual motion for a specific transport process is described by the $q_i$ rather than by the ${\alpha }_i$ (which represent just a particular case of all possible $q_i$). Note by the way that the $q_i$ satisfy the same Newton equations (with the same forces) as the ${\alpha }_i$. However, we shall impose to ${\alpha }_i$ boundary conditions that will guarantee zero final excitations whereas the initial conditions for the $q_i$ are arbitrary.

4.2. Boundary Conditions (BC) for $x\left(t\right)$ and ${\alpha }_i\left(t\right)$

We shall assume a transport from $x\left(0\right)=0$ to $x\left(t_f\right)=d$ with additional smooth boundary conditions for the trolley velocity, $\dot{x}\left(t_b\right)=0$ for $t_b=0,t_f$. We shall further assume that the auxiliary functions ${\alpha }_i$, as well as their first and second time derivatives vanish at boundary times $t_b=0,t_f$. We therefore have in principle a total of sixteen boundary conditions (BC), namely

$$\begin{array}{ccc}\hfill {\alpha }_i\left(t_b\right)& =& {\dot{\alpha }}_i\left(t_b\right)={\ddot{\alpha }}_i\left(t_b\right)=0,\hfill \\ \hfill x\left(0\right)& =& 0;x\left(t_f\right)=d,\hfill \\ \hfill \dot{x}\left(0\right)& =& 0;\dot{x}\left(t_f\right)=0.\hfill \end{array}$$

(30)

These boundary conditions guarantee that each invariant $I_i$ coincides with the corresponding Hamiltonian $H_i$ at initial and final times, see (27),

$$\begin{array}{c}\hfill H_q\left(t_b\right)=H_1\left(t_b\right)+H_2\left(t_b\right)=I_1\left(t_b\right)+I_2\left(t_b\right)=I\left(t_b\right).\end{array}$$

At these boundary times, and due to the $\dot{x}\left(t_b\right)=0$ boundary condition, the Hamiltonian represents the total mechanical energy of the system, i. e. $E\left(t_b\right)=H_q\left(t_b\right)$. If a fast finite-time process is designed so that the auxiliary functions ${\alpha }_i$ satisfy the imposed boundary conditions, the energy at final and initial times -regardless of the initial conditions of the hook and load, that is, regardless of the initial conditions set for $q_i\left(0\right)$ and its derivatives- will coincide since

$$E\left(0\right)=H_q\left(0\right)=I\left(0\right)=I\left(t_f\right)=H_q\left(t_f\right)=E\left(t_f\right).$$

Note that in principle the only conditions needed to guarantee $I\left(t_b\right)=H\left(t_b\right)$ are the ones for $\alpha \left(t_b\right)$ and $\dot{\alpha }\left(t_b\right)$. The others have a physical motivation as the desired boundaries for the trolley motion (on $x\left(t_b\right)$ and $\dot{x}\left(t_b\right)$) or are a consequence of the former ones (the ones on $\ddot{\alpha }\left(t_b\right)$ because of (28) and (29)).

In the following subsection we will show how to construct the trolley trajectory $x\left(t\right)$ so that the desired conditions in (30) are satisfied.

4.3. Inverse Engineering

We start by proposing the following ansatz for the trolley velocity $\dot{x}\left(t\right)$, symmetric with respect to $t_f/2$,

$$\dot{x}\left(t\right)=\sum _{j=1}^3{a}_{j}sin\frac{(2j-1)\pi t}{t_f},$$

(31)

with three free parameters $a_1$, $a_2$, and $a_3$ that will be determined from the following three conditions (the second line involves two conditions, one for each frequency, as justified in the Appendix):

$$\begin{array}{ccc}\hfill {\int }_0^{t_f}\dot{x}\left(\tau \right)d\tau & =& d,\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill {\int }_0^{t_f}\dot{x}\left(\tau \right)cos[{\omega }_j(\tau -\frac{t_f}{2})]d\tau & =& 0,\hfill \end{array}$$

(33)

for $j=1,2$.

Different functional forms are possible, but this ansatz is chosen for simplicity and because of very useful properties discussed in the Appendix A (it avoids resonance and boundary effects). It is also remarkable that an ansatz with only three free parameters satisfies the full set of sixteen boundary conditions in (30), see further details in the Appendix A.

The three free parameters can be therefore written in terms of the system physical parameters as

$$\begin{array}{ccc}\hfill a_1& =& \frac{75\pi d\left({\omega }_1^2{t}_f^2-{\pi }^2\right)\left({\omega }_2^2{t}_f^2-{\pi }^2\right)}{128t_f^5{\omega }_1^2{\omega }_2^2},\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill a_2& =& -\frac{75\pi d\left({\omega }_1^2{t}_f^2-9{\pi }^2\right)\left({\omega }_2^2{t}_f^2-9{\pi }^2\right)}{256t_f^5{\omega }_1^2{\omega }_2^2},\hfill \end{array}$$

(35)

$$\begin{array}{ccc}\hfill a_3& =& \frac{15\pi d\left({\omega }_1^2{t}_f^2-25{\pi }^2\right)\left({\omega }_2^2{t}_f^2-25{\pi }^2\right)}{256t_f^5{\omega }_1^2{\omega }_2^2}.\hfill \end{array}$$

(36)

These parameters determine completely the velocity of the trolley by (31), and its trajectory is simply the integral

$$\begin{array}{c}\hfill x\left(t\right)={\int }_0^t\dot{x}\left(\tau \right)d\tau ,\end{array}$$

(37)

which gives an explicit but lengthy expression. See some trolley trajectories and velocities in Figure 2. For long transport times (${\omega }_j{t}_f\gg \pi$) the trolley trajectory becomes independent of the masses and lengths of the pendulum and tends to

$$\begin{array}{ccc}\hfill x_{\infty }\left(t\right)& =& d\left[\frac{1}{2}-\frac{75}{128}cos\left(\frac{\pi t}{t_f}\right)+\frac{25}{256}cos\left(\frac{3\pi t}{t_f}\right)-\frac{3}{256}cos\left(\frac{5\pi t}{t_f}\right)\right].\hfill \end{array}$$

(38)

This trajectory implies a maximal velocity $v_{max}=(15\pi /16)(d/t_f)$ at $t=t_f/2$. In this asymptotic scenario there is only one acceleration time segment up to $t_f/2$ and a subsequent braking segment.

For short times compared to eigenperiods there are several segments of acceleration and braking. In any case this regime is less interesting in practice since the system deviates from the harmonic regime.

5. Numerical Results

5.1. Time Evolution of Suspension Angles

Once the trolley trajectory is designed, the dynamical evolution of the system can be found by numerically integrating the Euler-Lagrange equations of motion using either the exact Lagrangian (3) or the approximate Lagrangian in the harmonic (small oscillations) approximation (5). In Figure 3 some examples of the time evolution of the suspension angles ${\theta }_1$ and ${\theta }_2$ during transport are shown. The initial and final angles are not equal (unless the system is initially at equilibrium), but this is not a requirement for ending with the initial energy. The calculation has been done using the exact Lagrangian, but the results are undistinguishable in the scale of the figure when using the approximate Lagrangian since the involved angles are small throughout the whole transport process. For larger transport distance d or smaller process time $t_f$ these differences will increase and will lead to some errors due to the anharmonicity of the exact model as will be discussed in the following section.

5.2. Anharmonic Effects

For rapid transport operations, the involved angles are larger and the harmonic approximation breaks down, see Figure 4. Therefore, some deviations from the ideal results (i.e., equal final and initial energies) should be expected.

To quantify the excitation at final time in a way that is easy to understand and visualize, we measure the final energy $\Delta E$ in terms of a fictitious angle ${\theta }_f$. This angle is defined as follows: (i) the load and hook are initially in equilibrium (at rest in the vertical position); and (ii) the final energy is artificially interpreted as pure potential energy for a configuration where load and hook are at rest along a line with ${\theta }_f={\theta }_1={\theta }_2$. In other words: ${\theta }_f$ is the final angle when the final energy is considered to be purely potential and the two suspension angles coincide. Using (2) we may write

$$\begin{array}{c}\hfill \Delta E=-2E_0{sin}^2\frac{{\theta }_f}{2},\end{array}$$

(39)

with $E_0=-m_{1}g{l}_1-m_{2}g(l_1+l_2)$ being the energy for the equilibrium configuration. In Figure 4, this fictitious angle is plotted as a function of the process duration time $t_f$.

5.3. Stability

The stability of the proposed transport protocol can be studied by allowing some initial deviations of the angles ${\theta }_1\left(0\right)$ or ${\theta }_2\left(0\right)$ from the equilibrium positions. In Figure 5a, the final time energy excitation, measured in units of the fictitious angle ${\theta }_f$ (39), is plotted as a function of these deviations.

We will compare the resulting excitation with that for a simple third order polynomial ansatz for the trolley trajectory,

$$\begin{array}{c}\hfill x\left(t\right)=3d{\left(\frac{t}{t_f}\right)}^2-2d{\left(\frac{t}{t_f}\right)}^3,\end{array}$$

(40)

which satisfies the four BCs in (30) for $x\left(t\right)$ but not those for the auxiliar functions ${\alpha }_i$. As shown in Figure 5b (which should be compared with Figure 5a), the excitation at final time using this simple trajectory is much larger that the one using the inverse engineered trajectory. Our inverse engineering method leads to much more robust results.

5.4. Example Limiting the Maximal Trolley Speed

The engine power and safety considerations imply limits to the trolley speed. In this example we test the effect of such a limit. We set a load $m_2=1000$ kg transported a distance $d=40$ m. We also set a hook mass $m_1=150$ kg, $l_2=5$ m, and $l_1=40$ m. A maximum velocity of 2 m/s is assumed.

With this data, two transport protocols are compared in Figure 6: (i) inverse engineered trolley trajectory (37) and (ii) directly postulated cubic trajectory (40). For initial conditions at equilibrium and the same final process time $t_f$, our inverse engineering protocol involves higher maximum velocities but the crane ends with much lower energy, almost ending in equilibrium. In the dotted part of the curves the limit of 2 m/s is surpassed.

6. Conclusions

We have applied an invariant based inverse engineering STA method to design fast trolley trajectories of a double pendulum overhead crane. In the small oscillations regime these trajectories guarantee that the transport does not induce any energy excitation, regardless of the initial condition of the double pendulum. We have first found the normal modes and from them the dynamical invariants. Using these invariants, it is possible to inverse engineer STA trolley trajectories. We have performed the numerical simulations with the exact dynamics to see the parameter intervals where the protocol is accurate. Comparisons are also made with less sophisticated trolley trajectories that demonstrate the advantage of the STA approach. We have worked out a particularly simple design for the trolley speed with three sine terms, (31). It should be clear that we have not really optimized the trolley trajectory. One of the interesting facts about STA methods is that the solutions to the inverse problem are not unique. That means that there is much room for finding specific trajectories that optimize variables of interest, or are robust with respect to specific perturbations or parametric uncertainties [35]. STA combine well in particular with optimal control theory [3]. Thus STA provide a useful avenue to minimize the sensitivity to parameter uncertainties, one of the weak points of open-loop approaches. Other possible extension of this work may be to tackle combined or sequential operations with transport and hoisting [25].

Compared to previous work on methods without feedback [26,28,36], this paper exemplifies and introduces the use of shortcuts-to-adiabaticity in mechatronics for multimode systems. We refrain from performing a numerical comparison with “input shaping” methods because virtually any result would be possible given the flexibility of both input-shaping and STA methods to accomodate a vast family of possible designs for the trolley motion, corrections for increased robustness with respect to parameter uncertainties or noise. Nevertheless we would like to underline the simplicity of the basic invariant-based engineering for the moving double pendulum crane, compared to input-shaping approaches [26,28,36]. Even if the choice among methods may be a matter of taste and previous experience, we would like to argue that STA should be in the the toolbox of control methods, if only because STA are well tested and have been intensely developed theoretically along different approaches and applied to many experiments in AMO (atomic, molecular, and optical), and solid state physics [3]. Thus engineering applications may benefit from an important framework of techniques and concepts. By the way, a positive influence in the opposite direction, from mechatronics to AMO physics, is also expected. For example, state manipulation in AMO science has much to learn from a long experience on control with feedback in mechatronics.