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

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.


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 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(t) 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.

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].

Lagrangian
In terms of the angles θ 1 and θ 2 , see Figure 1, the Cartesian coordinates of each mass in a rest frame are given by These relations can be inverted to have the generalized velocitiesθ i in terms of the generalized momenta p θ i . The Hamiltonian is found from the Lagrangian as where µ = 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 where Whereas the potential matrix K is diagonal, the kinetic matrix T is not, i. e., p θ 1 and p θ 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(t). 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ẋ(t).

Diagonalization of H θ
Let us first define a new set of coordinates u 1 and u 2 by the linear transformation where the A matrix is yet to be determined. The corresponding momenta transform according to Reference [30] where A −T = (A −1 ) T stands for the transpose of the inverse matrix. The Hamiltonian in these variables reads We now look for a transformation matrix A that diagonalizes simultaneously both the ATA T and A −T KA −1 matrices in the expression above. To do so it is useful to define the following matrix which is symmetric and positive definite and therefore can be diagonalized by an orthogonal matrix O. Without loss of generality, this orthogonal matrix O can be parametrized as and choosing the parameter θ (not to be confused with the angles θ i ) by we have that The eigenvalues ω 2 i are positive since T is a positive definite matrix, and have dimensions of (angular) frequency square. The explicit expressions are in agreement with the eigenfrequencies given in Reference [31]. Now, by writing the transformation matrix as both quadratic terms in the transformed Hamiltonian (11) are diagonal since Finally, the Hamiltonian (11) takes the uncoupled form

Lewis-Leach Family of Hamiltonians and Second Canonical Transformation
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, By using this canonical transformation, the new Hamiltonian is a sum of two independent forced harmonic oscillators that belong to the LL family, Mg l 1 cos θ,

Explicit Expression of Normal Mode Coordinates
Taking into account the two canonical transformations, the explicit expression of the normal mode coordinates and momenta where I 2 is the 2 × 2 identity matrix and, using the explicit expression of A in (16), we have

Designing the STA Protocol
We are now ready to define the invariants and design the driving function x(t).

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 The invariants for (23) have the explicit form [32] provided the functions α i satisfy the following Newton equations, These α 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 α 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 α i . However, we shall impose to α i boundary conditions that will guarantee zero final excitations whereas the initial conditions for the q i are arbitrary.

Boundary Conditions (BC) for x(t) and α i (t)
We shall assume a transport from x(0) = 0 to x(t f ) = d with additional smooth boundary conditions for the trolley velocity,ẋ(t b ) = 0 for t b = 0, t f . We shall further assume that the auxiliary functions α 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 These boundary conditions guarantee that each invariant I i coincides with the corresponding Hamiltonian H i at initial and final times, see (27), At these boundary times, and due to theẋ(t b ) = 0 boundary condition, the Hamiltonian represents the total mechanical energy of the system, i. e. E(t b ) = H q (t b ). If a fast finite-time process is designed so that the auxiliary functions α 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 (0) and its derivatives-will coincide since Note that in principle the only conditions needed to guarantee I(t b ) = H(t b ) are the ones for α(t b ) andα(t b ). The others have a physical motivation as the desired boundaries for the trolley motion (on x(t b ) andẋ(t b )) or are a consequence of the former ones (the ones onα(t b ) because of (28) and (29)).
In the following subsection we will show how to construct the trolley trajectory x(t) so that the desired conditions in (30) are satisfied.

Inverse Engineering
We start by proposing the following ansatz for the trolley velocityẋ(t), symmetric with respect to t f /2,ẋ 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): 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 These parameters determine completely the velocity of the trolley by (31), and its trajectory is simply the integral which gives an explicit but lengthy expression. See some trolley trajectories and velocities in Figure 2. For long transport times (ω j t f π) the trolley trajectory becomes independent of the masses and lengths of the pendulum and tends to This trajectory implies a maximal velocity v max = (15π/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.

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 θ 1 and θ 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.

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 ∆E in terms of a fictitious angle θ 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 θ f = θ 1 = θ 2 . In other words: θ 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 with E 0 = −m 1 gl 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 .

Stability
The stability of the proposed transport protocol can be studied by allowing some initial deviations of the angles θ 1 (0) or θ 2 (0) from the equilibrium positions. In Figure 5a, the final time energy excitation, measured in units of the fictitious angle θ 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, which satisfies the four BCs in (30) for x(t) but not those for the auxiliar functions α 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.

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.

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.