Dynamical modelling of boom tower crane rigging systems: model selection for construction

Accurate dynamical models are imperative to the development of accurate monitoring and control systems, which are foundational to safety in construction and infrastructure projects. However, the highly coupled non-linear dynamics of crane systems requires the application of many simplifying assumptions to the dynamical crane model. To achieve accurate control, simplifications should yield minimal error in modelled behaviour for maximal reduction in model complexity. However, limited information is available on the situational suitability of different combinations of simplifications to construction tower crane models. This paper informs designers of the optimal dynamical models to represent boom tower cranes, with respect to the crane characteristics and selection criteria. The optimal models are determined though the comparison of ten 2D and 3D dynamical models in representation of three variations of boom tower crane that are commonly deployed on construction sites. The comparison includes analysis of over 100 simulations and experimentation. The value of the presented optimal model selection framework is in facilitating systems designers to develop accurate crane monitoring and control systems.


Introduction
A dynamical model is the representation of a real-world system to its monitoring and control systems. The state estimator predicts the state of the system by evaluating the dynamical model, and the control system is developed to control the dynamical model. If the dynamical model inaccurately represents the dynamics of the real-world system, then the performance of the monitoring and control system is reduced upon application to the real-world system [1].
Boom tower cranes ( Fig. 1) are deployed in construction operations to lift large and heavy payloads. Current tower crane control systems are insufficiently accurate in the task of payload alignment to the target [2]. Hence, human workers must physically manipulate the suspended payload into alignment. For steel beam erection [3] and curtain wall installation [4] this task is near to a fall-from-height hazard. This methodology is not safe; the most common types of construction accident include workers falling from a height [5] or being struck by a crane payload [6]. Safety in construction can be improved with more accurate automated crane monitoring [7] and control systems [4].
The development of monitoring and control systems for boom tower cranes is particularly difficult. Many limitations in current tower crane modelling and control literature were identified in [8]. Boom tower cranes are under-actuated and exhibit highly coupled nonlinear dynamics during luffing and slewing motions [9]. Hence, many approximations and simplifications are applied during modelling [1]. While simplifications to the dynamical model are a necessity, there is flexibility for the designer to choose which simplifications to apply. An optimal set of simplifications should yield minimal error in modelled behaviour for maximal reduction in model complexity.
The most common lumped parameter models approximate the system as a single or double pendulum. Control of 1 3 162 Page 2 of 17 these systems was comprehensively reviewed most recently in [9]. In the double pendulum model, the hook block and payload (Fig. 3) are separate point or lumped masses, while in the single pendulum model they are combined [9]. The truss structures (tower, boom, jib) of the crane are usually assumed to be inelastic [10], thus, the position of the boom head is a function of only the actuator positions. The actuator dynamics and limits are also often ignored [8], thus, the position of the boom head is effectively assumed to be precisely position controllable by an independent servomechanism. The hoisting rope is almost always considered straight, inelastic, and massless [9]. The degrees of freedom (DOF) of this double pendulum model are sway (hook block swinging about the top of the hoist with two DOF) and roll (payload tilting about the hook block with two DOF).
Evaluation of the double pendulum model in 2D found that the single pendulum model is insufficient when the sway and roll eigenfrequencies differ [11]. The payload pointmass assumption neglects the payload orientation. However, for long or flat payloads, the payload orientation is coupled with the sway dynamics. This was demonstrated for the 2D bridge crane [12], 3D bridge crane [13], 3D single pendulum jib tower crane [14], and the 3D double pendulum jib tower crane [15].
The approximation of inelastic truss structures neglects the complex dynamical coupling between the tower, boom, jib and the active load [10]. The approximation of inelastic rope neglects the heave vibrational mode (linear oscillation along the rope axis with one DOF). Consideration of heave is important in scenarios requiring high positional accuracy if the payload is heavy and the hoisting rope is long [16]. Transverse vibration in the rope has also been considered, however this significantly complicates the model, reducing the viability of modelling other features [17]. Likewise, consideration of wind force and frictional damping significantly complicates the model [18], especially where the payload shape [19] or transverse rope vibration [17] is considered.
Moreover, many rigging configurations of boom tower crane use multiple hoisting ropes. The sufficiency of the pendulum model to represent the dynamics of multi-rope systems is dependent on the specific rigging configuration, crane design, and operational requirements. For example, the four-nonparallel-rope rigging systems for gantry container cranes exhibit highly complex dynamics dependent on the hoisting rope angle and elasticity [20]. Conversely, in the modelling of a four-parallel-rope robotic crane, the pendulum model simplification is justified and supported with experimental validation showing high accuracy during feedback control [21]. Consideration of model complexities can allow creative control strategies. In [22], a novel tworope system utilises the complexities of the model to change the sway dynamics based on luffing angle, moving the system away from resonance. In [23] auxiliary ropes are used to control external disturbances. In [24], the multi-rope model complexities are used to estimate the payload mass. In [25], the stick-slip behaviour of rope over pulleys is modelled with consideration to the resultant form of the dynamical equations, however, the work is limited by simplifying the system geometry. Assumptions to simplify multi-rope geometry [26] and pulley geometry [21] can significantly affect the modelled skew behaviour (twisting of the hook block about the hoist rope with one DOF).
Specifically for rotary boom cranes, control research is mostly focused on operations to transport shipping containers (container cranes) [9]. However, construction cranes and container cranes differ in rigging configuration. Rotary container cranes typically use multiple separate hoist ropes, each terminating at the hook block, whereas construction cranes typically have only a single hoist rope that passes through sheaves on the hook block before returning to terminate at the top of the hoist (Fig. 1). This specific rigging configuration has been modelled in [27], where numerous vaguely justified simplifying assumptions effectively reduced the system to a single-pendulum model, and in [11], where similar assumptions were used in the 3D study.
The existing research does not evaluate which simplifications are appropriate to apply in the modelling of construction boom tower cranes. Therefore, this paper evaluates the question "What are the optimal dynamical models to represent the boom tower cranes that are commonly deployed on construction sites?".
Answering this question raises the prerequisite question of how the established simplified models from literature relate to these cranes in a geometric sense. For example, when attempting to apply a pendulum model to a crane with a returning hoist rope, where on the real crane ( Fig. 1) should be designated as the top of the pendulum? Most literature does not clearly define this relation [1]. Possible locations are where the hoist rope leaves the boom head sheave (although this is variable with the luffing angle and angle of sway), where the hoist rope returns to terminate, or somewhere between.
Having a rigorous geometric description dually reveals the whole family of cranes that share the same optimal dynamical model. This answers the auxiliary question "Should a controller that was designed for a particular crane system be deployed on a similar crane with a different rigging configuration?".
The methods of this research represent the variations of boom tower crane with either a returning hoist rope or a non-returning hoist rope rigging configuration (Figs. 1 and 3). The hook block is optionally rotationally joined-to or separate-from the payload. Five 3D and five corresponding 2D dynamical models are derived in representation of these cranes. The models rigorously related to the geometry of the real-world crane and are derived to be comparable to each other through a common ancestor model.
The dynamical equations are programmatically generated and simulated to evaluate differences in hook block and payload trajectory, for each model, with respect to variation in model parameters (e.g. mass values, hoist rope length, boom length). The trials are made independent from the specifics of any monitoring or control system by not implementing any control system. Instead, the inputs are position controlled on pre-defined trajectories. The results are additionally verified through experimentation on a lab-scale model.
Model accuracy is evaluated against complexity. This reveals which dynamical models are optimal with respect to the crane and payload characteristics, rigging configuration, and operational requirements.
The main contributions of this paper are: • A methodology to directly compare a series of crane models to each other. The methodology considers separating the comparison from any control or monitoring system, rigorously establishing commonality between models, and permitting various rigging configurations of the real-world crane. • A programmatic method to derive dynamical equations from geometric descriptions, while allowing for a holonomic algebraic constraint. The system of differential algebraic equations (DAEs) is generated by applying the Euler-Lagrange formulation with Lagrange multipliers and the decomposition of the rotation matrix derivative. The DAEs are then reformed into ordinary differential equations (ODEs) by solving for the Lagrange multiplier. This reduces numerical error in simulation by six orders of magnitude. • Precise description of models that represent three variations of boom tower crane that are commonly deployed on construction sites. The most complex model accounts for rope interaction with the pulley sheaves, while the most simplified model is a single pendulum. Each model is fully defined in how it relates to the geometry of the real crane. • A decision tree, Fig. 10, which guides system designers to choose the optimal dynamical models to represent the boom tower cranes that are commonly deployed on construction sites. The optimal model is dependent on the crane and payload characteristics, rigging configuration, and operational requirements. Understanding these dependencies both enables choice of the optimal model during systems development, and guides application of the developed system to the whole family of cranes which share the same optimal model. Section 2 describes a method to programmatically formulate and simulate dynamical models from given geometric descriptions. Section 3 describes the models being compared. Section 4 discusses considerations of simulating high-complexity dynamical models. Section 5 compares the models through scaled simulation and experimentation for common construction scenarios. A decision tree to choose the optimal model is presented in Fig. 10. Section 6 concludes the paper and identifies opportunities for future research.

Method of dynamical model formulation
Crane systems can be mathematically represented by systems of differential equations, where the solution to the equations describes the motion of the system [9]. The chosen crane systems are of high complexity, thus manual formulation of the equations of motion is infeasible.
The equations were generated programmatically with the MATLAB Symbolic Math Toolbox. The input to the equation generator was the geometric description of the system. This was realised through symbolic homogeneous coordinate transformation matrices, symbolic constraint equations, and specification to distinguish between the symbols that represent generalised coordinates or constants.
The output of the generator was the system equations, formulated and rearranged for input to the ordinary differential equation (ODE) solver. Following limitations of MATLAB to solve the more complex ODEs generated by the 3D models, the generator was made to cross-generate C++ code for solving with the SUNDIALS CVODE solver [28]. This solver was chosen based on the comparison of many ODE solvers in [29]. The SUNDIALS CVODE is a variable-step, variable-order (VSVO) explicit ODE solver of orders from 1 to 5 [28]. The solver was configured for solving stiff ODEs by using the Backward Differentiation Formulas.
To improve performance in both generation and solving, as much as possible of the equation formulation was delayed to occur numerically during solving. This required the generator to output the matrix coefficients of the semi-formulated system equations. The coefficients could then be evaluated at runtime, and solving completed through applying linear algebra on the numeric equations. To perform linear algebra at runtime in C++, the Armadillo C++ library [30] was used.
The following subsections describe the model specification (geometric description) which is input to the automated generator and solver; the method of equation generation; and the verification of the methodology.

Notation
The notation in this paper writes scalars in non-bold, vectors in bold-lowercase, and matrices in bold-uppercase. Right-superscripts are used only for powers or the matrix transpose. Right-subscripts of i or j denote the index to a parent vector, matrix, or set. Unless specified otherwise, other right-subscripts and left-scripts are used to distinguish between variables.
Centred dots above variables denote the single, ȧ , and double, ä , time derivatives. The Euclidean norm is denoted ‖a‖.
Square brackets notate matrix concatenation. Where otherwise ambiguous, parentheses notate function arguments. Hence, a ⋅ b contextually represents either scalar multiplication or dot product. In unambiguous cases, ab is also used to represent scalar multiplication.
Position vectors are notated as A p B to specify the vector from the origin of frame A pointing to the origin of frame B, as measured in the coordinates of frame A. The inertial frame is designated the letter O. Rotation matrices are notated as new frame

Formulation of geometric description
For each model, the geometric description was formed by approximating the system as a kinematic linkage and assigning a set of linearly independent generalised coordinates.
where q free is the vector of n number of unactuated DOF, and q input is the vector of m number of externally position controllable DOF, as given by Eqs. (2) and (3) respectively.
The kinematic linkage was then described with a set of homogeneous coordinate transformation matrices. To maintain linear independence for the closed chain linkages, additional dependent coordinates were temporarily assigned and then solved for in terms of the independent generalised coordinates by means of inverse kinematics. The solution was found by equating the transformation from one node of the closed chain to another in a clockwise direction, to the transformation in the anti-clockwise direction. Equating these transformations, the matrix elements were simultaneously solved to isolate the dependent coordinates.
The transformations from the centre of mass frames to the inertial frame were obtained, then decomposed into position and orientation transformations, O p i (q) and O i R(q). Directly solving all the system constraints proved difficult for the system models which describe the pulley behaviour. In this case, one more generalised coordinate than degrees of freedom was used, accompanied by a holonomic constraint equation of the form C(q) = 0 . The following formulation allows either zero or one constraint equations, as is sufficient to describe the system.

Formulation of system equations
The system equations were obtained through the Euler-Lagrange formulation with Lagrange multipliers where L(q,q) is the Lagrangian, Q i (q,q) is the generalised force, C(q) is the left-hand side of the constraint equation, and is the Lagrange multiplier. The Lagrangian is formed as where the kinetic and potential energies, K and V, are found through where O p i is the position of the centre of mass frame i, with respect to the inertial frame. This is obtained from the geometric description. m i is the corresponding mass, and I i is the corresponding moment of inertia tensor, measured in frame i. g defines the gravity vector in the inertial frame. i is the angular velocity of frame i with respect to the inertial frame, as measured in the coordinates of frame i. It is calculated by decomposing the angular velocity tensor, i , as where i (row, column) is the matrix element (row, column) of i , which is obtained from the well-known rotation matrix derivative [31] where O i R is the orientation of frame i, as measured in the inertial frame. This is obtained from the geometric description.
The generalised forces are obtained by remapping the external forces, F j , with where O p j is the position where F j acts.

Reforming of system equations for solving
Computing the derivatives in (4) and moving all terms to the same side of the equation results in the form Of the MATLAB ODE and DAE solvers, only ode15i is compatible with system equations of the form 0 = f(t, x,ẋ) . However, ode15i is an inefficient and low accuracy solver. Hence, the system equations resulting from the Euler-Lagrange formulation were reformed for compatibility with the higher accuracy solvers.
The common methodology is to reform (11) into matrix equations and isolate q . However, due to the constraint, the resulting equations would still be index-1 DAEs. DAE solvers are low accuracy and inefficient in general. Hence, we developed the general solution to reform (11) into a system of ODEs by algebraically solving the constraint equation.
Starting from (11), for each i, a f i can be decomposed into Hence, (11) can be reformed into (13) and then rearranged into (14) To solve for , the constraint equation is first double differentiated and then reformed into a matrix equation Substituting (14) into (15) and solving for gives Substituting (16) into (14) gives the system equations in the form The differential order of the system must be reduced for solving. The vector of the inputs, u , is defined to comprise the externally controlled generalised coordinates and their derivatives (18). The state vector, x , is defined to comprise the free generalised coordinates and their first order derivatives (19).
Appending the order reducing relation q free =q free to (17) gives rise to Substituting the generalised coordinates, x and u , into (20) results in the nonlinear state equation To enable the ODE solvers to solve with time varying input, (21) was aliased with where u(t, x) is the function that calculates u at the solution time t. This function can implement state feedback control; or open loop control if u(t, x) = u(t) . To implement this function alias, a wrapper function with inputs x and t is defined to first evaluate the input vector, and then evaluate (21).
Finally, this results in the time and state dependent state equation This form of the system equations is compatible with ODE solvers, for example the MATLAB ode45 and SUNDIALS CVODE solvers.

Algorithm performance and verification
Two tests were conducted to verify the performance and correct functioning of the dynamical equation generator and simulator.
All of the applicable matlab ODE and DAE solvers were tested. The ode23 and ode45 are low-to-medium accuracy explicit ODE solvers that respectively implement Bogacki and Shampine's Runge-Kutta (2,3) pair, and Dormand and Prince's Runge-Kutta (4,5) pair [32]. The ode78 and ode89 are high accuracy explicit ODE solvers that respectively implement Verner's "most efficient" Runge-Kutta (7,8) pair, and Verner's "most robust" Runge-Kutta (8,9) pair [33]. The ode113 is a variable-step, variable-order (VSVO) Adams-Bashforth-Moulton explicit ODE solver of orders from 1 to 13 [32]. The ode15s is a VSVO explicit ODE and semiexplicit index-1 DAE solver of orders from 1 to 5 [32]. The ode23t is a variation of ode23 that can solve explicit ODEs and semi-explicit index-1 DAEs [34]. The ode15i is a VSVO fully implicit index-1 DAE solver of orders from 1 to 5 [35]. The ode23s, ode23t, and ode23tb solvers were not tested because they require that the mass matrix is constant, which is not generally true for this problem.
The first test modelled a 4-parallel-bar linkage to be equivalent to the point mass pendulum (Fig. 2). The linkage was modelled as an open kinematic chain A-B-C combined with the algebraic constraint equation, to realise (15), as The resultant DAEs were simulated both as DAEs, and then as reformed into ODEs. The results were compared to the analytic solution of an unforced point mass pendulum [36]. The simulation duration was 20 s, from initial conditions with the pendulum almost inverted ( 0 = 0.9 ), and with integration tolerances of 10 −10 .
The second test modelled a 2D triple pendulum, where the first link was given as a prismatic joint. The constraint equation realise (15) was given as 0 =L 1 , where L 1 was the displacement of the prismatic joint. The same simulation conditions were used. The reference 'true' solution was the simulation using ode89 with the MATLABs maximum allowable tolerances.
The results (Table 1) were analysed by the time taken to form and simulate the equations, and by the root-meansquare error to the reference solution. The results for the first test closely follow the analytic solution, near to the integration tolerances, verifying the correct function of the dynamical equation generator and simulator. The reformulation reduced the error by six orders of magnitude. The results for the second test have greater error, as attributed to the greater complexity of the problem. The reformulation allowed for solving the system, where the the original system of DAEs gradually diverged.

Dynamical models
Five 3D and five corresponding 2D models were developed to represent three rotary boom construction tower crane systems. The models are full complexity (FC), zero-radius The models are comparable to each other through a common ancestor, the 'full complexity' model, from which the remaining models were derived (Fig. 3).
All models assume that the hoist rope is inelastic and massless; the tower and boom are inelastic, transverse vibration in the rope can be ignored; and wind force and frictional damping can be ignored. We justify that these assumptions are required to avoid greatly increasing the complexity of the model. As discussed in the introduction, relaxing any of these assumptions significantly inhibits modelling other features.

Full complexity model and crane
The full complexity model describes the crane rigging configuration using a returning hoist rope (Fig. 1). Additionally, this model includes an active skew rotary hook block (rotator), per the recommendations of [37] to deploy rotators on construction cranes for increased safety and economic efficiency. The kinematic linkage of this model is shown in Fig. 3. Coordinate frames were attached to this linkage (Figs. 4 and 5). The orientation of each frame is set to make solving the inverse kinematics equations simpler.
The location A is the intersection of the luffing and slewing axes; B is where the hoist rope returns to attach to the boom; C is the centre of the boom head sheave; and D is where the hoist rope leaves this sheave. These locations  E and H are where the hoist rope leaves the hook block sheaves. F and G are the centres of the hook block sheaves. J, K, and L are in-line, defining the axis in which the skew actuator rotates. This axis is assumed to be perpendicular to the line joining F and G, with J located halfway between F and G. Thus, the hook block is symmetric.
K is the centre of mass of the hook block. L is the point where the hook joins to the payload, which is free to tilt with 1 DOF. M is the centre of mass of the payload.
The generalised coordinates are as defined by Figs. 4 and 5. Note that this definition is overdetermined by 1 DOF. Therefore, the constraint Eq. (28) is used to ensure that the coordinates are consistent.
1 is the slewing angle; 2 is the luffing angle; and 10 is the skew actuator angle. 5 is related to the skew of the hook block. 11 is the angle of payload relative tilt.
The rope leaves each sheave in a direction tangent to the contacting surface at the point of departure, where the points of departure vary with 3 , 7 , and 8 . The tangent relations are described with ∠CDE = ∠DEF = ∠GHB = 2 . This definition allows the rope to bend out of plane from the sheaves, at the angles designated 4 , 6 , and 9 .
The length from D to E is designated L DE . L DE and L HB vary with hoisting actuation and rolling of the hook block along the rope. The effective length of the rope, L rope , describes the path length of the rope from partial a turn before it leaves the boom head sheave, to the end of the rope at B.
Change in this length is wholly determined by the input rotation of the hoist actuator, hence L rope is used as a control input. Since this length is easily measurable with encoders, it is reasonable to use (27) as the constraint equation. Moving all terms of (27) to the right-hand side of the equation and double differentiating gives the constraint in the form of (15) The dependent coordinates, 8 , 9 , and L HB , were found in terms of the independent coordinates with inverse kinematics methodology. Using the transformation from B to G as described in terms of the independent coordinates Splitting (29) into its components gives the inverse kinematic solution Fig. 4 Coordinate frames, generalised coordinates, and dependent coordinates attached to the full complexity model, Fig. 3. Upper part

Model with zero-radius pulleys
Using the full complexity model as a starting point, the remaining models were derived by applying geometric simplifications.
An ideal simplification reduces model complexity without introducing error in modelled behaviour. From the Euler-Lagrange formulation, it can be seen that the behaviour of the system masses is central to the modelled behaviour. Therefore, simplifications should minimally alter the manifold of achievable configurations of the masses. Given that crane systems should always be kept near to the equilibrium configuration, then this requirement should be most strictly imposed near to equilibrium. This is analogous to linearising a system about an equilibrium point.
However, the equilibrium configuration varies with the luffing angle. Hence, all our simplifications are based on the requirement to not change the equilibrium configuration of the masses at the luffing angle 2 = 4 . This angle was chosen for being near to the middle of the operational range, thereby providing the most benefit across all angles of operation.
Starting from the full complexity model, the pulley sheave radii are reduced to zero by means of the simplifications and C and D are moved to where D would be when the rope leaves the main hoist sheave vertically with the luffing angle This model is not physically viable to implement on a construction crane as cables are rated with a minimum allowable cable bending radius to prevent fatigue failure.

Triple pendulum model and crane
The triple pendulum model is formed by assuming that a virtual rope joins the locations I and J, where I is exactly midway between B and D (Fig. 3). Thus, the kinematic linkage is the sequence joining A-I-J-K-L-M. I is a three DOF spherical joint (sway and hook block skew), and J a two DOF universal joint (roll).
Combined with simplifications to the pulley radius, the virtual rope model is argued to offer great value in simplification for eliminating complexity while minimally impacting the expected position of the mass, by comparison to the large scale of the crane system [21]. However, it is cautioned that this small change in payload height is the dominant term in the restoring force against hook block skew [21].
Starting from the full complexity model, two underlying assumptions are made in this formulation: • The position of I does not change. • The virtual rope is of constant length.
For effective control of the system, these parameters should be set so that the model matches the true system at equilibrium. However, the values change between the equilibria of different luffing angles and hoist rope lengths. For an overall close approximation, the equilibrium parameters with the luffing angle of 2 = 4 were used to specify the length of the rope and position of I. The hoist rope length was held constant during experimentation.
This model is additionally representative of a crane with a single, non-returning, hoist rope, where the size of the hook block is significant. Hence, bending can occur either at the join of the hoist to the hook block, or at the join of the hook block to the payload.

Double pendulum model and crane
In the double pendulum model, it is assumed that the payload sways with hook block. This represents systems where the hook block and payload are rigidly joined together, as is seen on container cranes, where the spreader locks to the container.
Starting from the triple pendulum model, the only additional assumption is that 11 = 0 . This definition differs from the double pendulum model that is common in crane control literature, in that the hook block mass is not concentrated at the pivot point J [9]. We consider this to be a more accurate representation of the system for no increase in complexity, given that moment of inertia is included in the model.
The common definition of the double pendulum model in literature assumes that the size of the hook block is negligible [9]. This allows for use of the double pendulum model when the hook block and payload are not rigidly joined. Mechanically, it is not possible for bending to occur at the centre of mass of a typically shaped hook block. Hence, the choice of the pivot point is somewhat arbitrary. If the size of the hook block is significant, then for this case, the triple pendulum model may be more appropriate.

Single pendulum model
In the single pendulum model, it is assumed that the payload sways with hook block and hoist rope; that I, J, K, L, and M all remain inline. The separation of the hook block and payload masses is equivalent to the typical realisation in which they are lumped together at the total centre of mass.
No construction crane directly relates to the single pendulum model as the join between the hoist rope and hook block does not sit at the centre of mass of the hook-payload pair. Hence, the system dynamics can induce bending at this join.

3D to 2D model simplifications
Each 3D model has a corresponding 2D model. The 3D to 2D simplification enforces These are the minimum necessary simplifications required to reduce the 3D model to 2D. The skew actuator angle is set to allow relative tilt between the hook block and the payload (37). An alternative simplification to (37) could be 10 = 11 = 0 , although, this prevents the relative tilt.

Considerations in simulation
MATLAB succeeded in the generation and manipulation of the equations, even for the high complexity models. However, it was unable to solve the higher complexity models. Hence, the coefficients of (13) and (15) were exported to C++ and solved with SUNDIALS CVODE. To visualise the size of these equations, a single evaluation of the full complexity system's system equations calls the trigonometric function sin() 1.7 million times.
A key technique to reduce numerical error and make generation and solving feasible was to delay the required matrix inversion until runtime, where the inversion could be completed numerically. In MATLAB this was achieved by storing the partially formed equations in an object and passing a method to evaluate the equations to the MATLAB's inbuilt ODE solver. In C++, this was achieved by exporting the matrix coefficients of (13) and (15), which were evaluated for the final formation to be completed numerically with Armadillo. Another specific code optimisation was found to be necessary to enable compilation of these large C++ source files. Without optimisation, the 3D triple pendulum model took over 3 h to compile. Post optimisation, the same model compiled in under 1 min. The required optimisation was to reduce the number of calls to sin() and cos() by substituting repeated calls with variables evaluated before the body of the ODE.
All the 2D and 3D simulations completed in less time than the length of simulation. Hence, real time evaluation of even the highly complex models is feasible.
In numerically solving differential equations, numerical error accumulates with each solve step, eventually causing the solution to diverge. For the crane models, divergence occurred sooner with the increasing ratio of the smallest length in the kinematic chain to the largest link length and increasing model complexity. The closed chain models suffered significantly greater numerical error.
For the closed chain models and the triple pendulum model, the mass matrix is singular when the moment of inertia of the hook block is set to zero. Similarly, if other mass or moment of inertia values are set to zero, the mass matrix can become singular. Hence, if these values are small, the ODEs can become very stiff. To solve this problem, the stiff model can be simplified by setting the related joint in the kinematic chain as fixed, resulting in a non-stiff model. For example, where no payload is attached, the angle of the payload can be set to a constant without reducing the accuracy of the solution.

Simulated and experimental comparison
Two set of comparisons were conducted. First, a simulated and experimental comparison of the model trajectories under slewing and luffing input. Second, over 100 simulated trails are compared under luffing input, measuring the maximum angular disturbance for variation in model parameters.

Model realisations
The model parameters represent a Liebherr 710 HC-L 32/64 crane, described by the manufacturer as a high capacity construction crane ( Table 2). This crane is capable of multiple rigging configurations, where the rope either terminates at the hook block (the triple pendulum model), or passes through sheaves on the hook block and returns to terminate at the boom (the full complexity model). The hook block was replaced with an active skew rotary hook block, as discussed in Sect. 3 to increase safety in crane operations.
Two payloads were tested, representing a 1000 kg glass curtain wall module, and a 9000 kg prefabricated concrete wall module. Hence, this crane and payload system is representative of dangerous construction crane operations which have need for greater monitoring and control accuracy. Optionally, the payloads may be mechanically orientation locked to the hook block. In total, this results in three experimental rigging configurations (full complexity, triple pendulum, double pendulum), each with a choice from two payloads (curtain wall module, prefabricated concrete wall module).
All the experiments and simulations were scaled down by a length factor of 25. Dimensional analysis relates the mass and time scaling factors as mass ∶ length 3 and time ∶ √ length . Hence, time scaled 5 times faster and mass was reduced by a factor of 15625. This was verified by comparing that a full-scale simulation and a scaled simulation produce identical scaled results. The following figures and equations use the experimental scale.
The experimental setup used a Universal Robots UR5 robot to emulate the movement of the crane boom, with the end effector holding the head of the boom from which the hook block was hung (Fig. 6). The accuracy of this setup to represent a full-scale crane was limited by imperfect replication of the flexural and torsional rigidity of the crane's structure [38].
The robot was servo controlled to smoothly follow the input trajectories. The motion of the system was captured with 0.1 mm accuracy by a Vicon Bonita motion capture system. The motion capture software Tracker 3 directly returns the global position and orientation of each rigid body in the system, as determined by directly resolving the location of markers attached to the body. Comparison requires the model parameters, inputs, and configurations to be in some way equivalent between every model. Hence, across all models, we required that the centres of mass have the same initial and equilibrium locations with zero-velocity initial conditions. These requirements were satisfied by initialising every model in static equilibrium. For any other initial configuration, the requirements are non-trivial to satisfy due to the geometric differences of the closed kinematic chain models to the pendulum models.
From these initial conditions, motion was induced by actuating the inputs, (26), within the operational limits of a Liebherr 710 HC-L 32/64 crane. The trajectories of the inputs were pre-defined. This enables comparison between the dynamic responses of each model in a way that is independent of any control system. From preliminary tests, inputs that induced the most variation in trajectory between models were chosen. Hence, the reported results show where the models behave most differently from each other.
The trajectory comparison trials used input to induce oscillation with gradually increasing amplitude, up to the amplitude limit of the experimental workspace. This is realised as For the 2D trials, the slew and skew inputs were instead held at their initial values.  For the trails with varying system parameters, the input was a smooth upward luffing motion, realised with a cubic spline. The crane was then held still to continue to capture the oscillations after the luffing finished. The trial ended after the largest amplitude oscillation was captured.

Trajectory comparison
Comparison of the simulated and experimental results for 2D motion is Fig. 7. The comparison uses the inertial frame A1, formed by the triad x A1 y A1 z A1 in Fig. 4. The general shape of motion agrees. The larger amplitude in simulation is attributed to the simulator omitting frictional forces and air resistance. This provides strong evidence to the correct functioning of the simulator. Comparison of 3D experimental and simulated trials for different models and rigging configurations is Fig. 8. All experimental models are plotted. The omitted simulated models diverged during solving for reasons discussed in Sect. 4. Similarly to in Fig. 7, the simulated results have greater amplitude. All the pendulum models closely follow the same trajectory in both simulation and experiment, while behaviour of the full complexity model differs more significantly from the pendulum models. This same pattern occurred for both payloads and inputs.
The input prevalently excited the first vibrational mode, which explains the strong similarity between the results of the pendulum models. For this case, the single pendulum model achieves a very accurate and low complexity representation of the single-hoist-rope system. However, for less smooth inputs which excite other vibrational modes, the similarity between the configurations is expected to reduce. This is tested as follows. Figure 9 compares the maximum angular disturbance from an upward luffing input for differing system parameters.

Comparison with varying system parameters
Across all combinations of system parameters, all models show similar behaviour. The angle of the payload shows the most variation between the models, particularly for lighter payloads.
In general, the single pendulum model exhibits the smallest hook block and payload angles. This is because these angels are bound to the angle of the hoist. The great length of the hoist results in a very large moment of inertia about the pivot. The other models can rotate the hook block and payload more freely.
Likewise, the double pendulum model exhibits smaller payload angels than the triple pendulum and pulley models. Mechanically locking the payload and hook block together increases the moment of inertia of the combined mass, attenuating any disturbance that is applied to either component. This is most evident when the payload is much lighter than the hook block, where the large payload oscillations exhibited by the higher complexity models are attenuated.
The triple pendulum and zero radius pulley sheave models exhibit almost identical behaviour. This similarity can be likened to the similarity between a 4-parallel-bar-linkage and a pendulum because the ropes of the pulley model are close to parallel, as is typical of the rigging on construction cranes. It is important to note that the system with significantly nonparallel ropes has already been shown to behave significantly differently from the pendulum model in [22]. The full complexity model is most similar to the triple pendulum and zero radius pulley sheave models, as is most evident when the payload is lightweight.

Optimal rigging choice and model choice
In choosing the physical rigging configuration, mechanically locking the hook block and payload together has the Fig. 7 Trajectory of the prefabricated concrete wall module relative to the inertial frame A1. Comparison of the simulated triple pendulum model to the experimental triple pendulum. The positions K and M coincide with the masses of the hook block and payload respectively desirable effect of reducing payload swing, particularly where the payload is lightweight. This knowledge can then be used by designers to make the system easier to control. Additionally, this design is highly compatible with the implementation of a robotic hook block which can assist in payload alignment to the target. In choosing which model should represent a given physical system, we present Fig. 10, a decision tree to choose the optimal model for different requirements. The design of this decision tree is discussed as follows.
For each leaf of the decision tree, the optimisation followed the methodology:  Simulated maximum angular disturbance for various system parameters, compared between each dynamical model. In each plot, only one parameter is varied, while holding the other values constant as per Table 2 1. The results are interpreted to classify each model's ability to describe system information (payload skew or roll orientation), and classify the compatibility of the geometric features of each model to each variation of realworld crane (parallelism of the hoist ropes, joint between the hook block and payload). 2. Using the model classifications, models that do not satisfy the categorical requirements of the decision tree leaf are discarded. For example, the single pendulum model inadequately describes payload roll, hence, it is discarded for the applications where payload roll is considered important. 3. Of the remaining models, the least complex model is optimal. The order of most-to-least complex model is full complexity, zero-radius pulley sheaves, triple pendulum, double pendulum, and single pendulum.
Following, we discuss the model classifications and decision tree requirements. Pendulum models have previously been shown to be inappropriate for application to non-parallel-rope systems [22]. Many variations of non-parallel-rope rigging configurations exist. Hence, choice of the specific closed kinematic chain model is out of scope.
Next, the combined impact of all model simplifications should be considered. Figure 8 demonstrates a case where the magnitude of error in the location of the centre of mass from choosing a pendulum model over a higher complexity model is comparable to the error from omitting frictional forces. The approximation of this system as a pendulum system could hence be considered appropriate in the same circumstances that justify omitting friction from the model. However, for a more accurate representation of the system that includes frictional forces, the impact of representing the system with a pendulum model should be carefully tested for the exact system parameters.
Omitting frictional losses is common practice, for example in closed-loop control systems design, which requires very low model complexity [1]. In this case, the simplification is acceptable. Thus, controllers designed for low-friction parallel-rope systems and pendulum systems should be interchangeable without significant impact on performance.
In contrast, a state estimator may make use of the higher accuracy of a closed kinematic chain model to improve estimation accuracy, or to estimate otherwise unobservable parameters [24]. The choice of model is then dependent on the specific problem requirements.
Furthermore, consider what information is required from the model. If only the location of the centre of mass is important, then large error in the payload angle may be acceptable for the case where the length of the hoist is much greater than the length of the payload. Hence, a double or single pendulum model may be sufficient.
Finally, it should be considered that, for the higher-complexity models, the mass matrix can become singular in special loading cases, inhibiting solving of the system equations. Furthermore, the equations can become very stiff for other combinations of system parameters which may realistically be expected to be encountered during crane operation. For example, when no payload is attached, the model should fix the angle of the payload as a constant.

Conclusions and future work
The current practice for construction crane operations is not safe and can be improved with more accurate payload monitoring and control. Crane control systems development requires the application of many simplifying assumptions to the dynamical crane model. The optimal dynamical model should achieve maximal reduction in model complexity, with minimal error in modelled behaviour, while satisfying the application requirements.
For the first time, this research presents the decision tree Fig. 10 which guides systems designers to choose the  10 Decision tree for choosing the optimal model for different requirements optimal dynamical model, and guides the application of these systems to the family of cranes which share the same optimal dynamical model. The models represent three variations of boom tower crane that are commonly deployed on construction sites. We also present a methodology to rigorously evaluate crane models through their relation to the real-world crane, and a complete programmatic method to reform the constrained dynamical equations to reduce numerical error in simulation by six orders of magnitude.
The modelling methodology ignored frictional damping, assumed inelasticity and weightlessness of the hoist rope, and assumed infinite stiffness of the tower and boom. Simultaneous inclusion of all these complexities in a single model results in the model becoming too complex to feasibly solve. The trajectory comparison shows the impact of these assumptions on the simulated trajectory. As a result, this work considers some modelling cases to be out of scope of Fig. 10. Future work is recommended to systematically review the literature that addresses these cases in detail. Further expansion of the investigation to more types of crane is also recommended.
In the derivation of the pendulum models, a virtual rope is assumed to replace the complex pulley system. This requires choosing a single nominal luffing angle where the locations of the masses match the real locations at equilibrium. At other luffing angles there is a constant offset error, which is very small by comparison to the large scale of the crane system. Hence, the error has minimal effect on sway dynamics, irrespective of the choice of nominal luffing angle. However, for use in systems to estimate the location of the payload, it should be recognised that the resultant length of the hoist rope has a small offset error that is dependent on the luffing angle.
In future work, we will apply these results to the development of a robotic crane end effector for assembling prefabricated modules on high-rise buildings. Understanding the effects of rigging configuration on system dynamics enables choice of the optimal rigging configuration. Understanding the family of cranes with similar dynamics enables targeting the system to the whole family of cranes instead of just one crane. Furthermore, the situations in which operation of the device will be safe can be identified. Finally, the perception, monitoring, and control systems onboard the device can deploy the optimal dynamical model to achieve optimal performance.