Dynamic Interaction of a Heavy Crane and a Ship in Wave Motion

Previous work on the dynamics of vehicle-manipulator systems is extended to offshore ships with heavy cranes. The proposed method is based on a Newton-Euler formulation where the forces of constraint are eliminated using projection matrices based on the method of Kane’s equations of motion. This leads to an efficient method for developing the equations of motion of a ship with a heavy crane so that the motion of the crane influences the motion of the ship and vice versa. The calculation of the projection matrices is made efficient and intuitive by observing that the columns of the projection matrices are the screw axes of the joint twists in Plücker coordinates. Wave excitation of the ship is modeled with force RAOs based on established wave spectra. This gives a model that is well suited for design and testing of crane control systems, and for studying the feasibility of demanding crane operations in different weather conditions.The resulting equations of motion have been validated in simulation experiments for a ship with a 3 DOF heavy crane with a payload, where the ship is excited by a JONSWAP wave spectrum using a simple controller based on feedback linearization. The simulations clearly demonstrated that the ship responded in a physically meaningful way to the motion of the crane.


Introduction
Crane operations are important in the offshore industry, where heavy cranes are mounted on floaters and ships.Offshore cranes are used for heavy loads, including Blow-Out Preventers (BOP), which have a mass in the order of 350 metric tonnes.At the same time there is an interest in using smaller ships to reduce costs.This means that the mass of the crane and the load can be significant compared to the mass of the ship, which means that the motion of the crane will significantly influence the motion of the ship.In terms of modeling, such systems have typically been modeled with separate models for the ship and the crane.This was done in Chu et al. (2017), where a model of a ship and a model of a crane was connected using the Func-tional Mock-Up Interface (FMI).This approach is very efficient, but it will not model the inertia coupling of ship and crane.Alternatively, this can be implemented so that the motion of the ship will influence the crane, while the motion of the crane will not influence the motion of the ship.Then the ship motion can be simulated and the results can be used as inputs to the crane model.This was done in Masoud et al. (2003) where the control of payload pendulations on marine cranes was studied, in Love et al. (2004), where tracking control was studied for a manipulator mounted on a ship using a Lagrangian approach.
The resulting motion of the ship can be modeled using transfer functions in the form of motion RAOs or force RAOs (Response Amplitude Operators) that are calculated from the geometry of the hull using detailed hydrodynamic models.The resulting RAOs are transfer functions given in non-parametric form in terms of the numerical values of the transfer function gain and phase as a function of the wave frequency.The use of force RAOs makes it possible to combine this wave model with forces from other sources, like a crane mounted on the ship.Perez et al. (2004) and Fossen and Perez (2007) review the kinematic models commonly used in seakeeping and manoeuvring theories and provide a complete derivation of the transformations that link these descriptions and combining the models to study ship motion.
An offshore crane will have a kinematic arrangement of the same type as an industrial manipulator.This means that the dynamic model of a crane on a ship can be derived in the same way as the dynamics of a manipulator on a moving base, like a spacecraft, an underwater vehicle, or a ship.An important difference compared to the usual robotic dynamics is the orientation of the base which will be given as a rotation matrix R in SO(3), and there will exist no minimum representation of the rotation in terms of generalized coordinates.A model of a manipulator on a spacecraft was presented in Egeland and Sagli (1993) where Kane's equations of motion were used so that generalized speeds could be used instead of generalized coordinates.From et al. (2010From et al. ( , 2014) ) derived a singularityfree dynamic equations of a robotic manipulator on a non-inertial base, and showed how the equations of motion could be developed for vehicle-manipulator systems based on Lagrange's equations using Lie group techniques.
In this paper we derive the equations of motion for the combined ship and crane dynamics based on a Newton-Euler formulation where the forces of constraints are eliminated using the technique of Kane's dynamic equations Kane and Levinson (1985).The proposed method is based on projecting the equations of motion for each link using the partial linear velocities and the partial angular velocities defined by the generalized speeds, which was done in Egeland and Sagli (1993) for a spacecraft-manipulator system.An important improvement in this work is that the partial linear velocities and the partial angular velocities are defined as screws in the form of lines in Plücker coordinates.It turns out that the relevant lines are the lines of the joint axes, and the description of the lines as screws gives useful geometric insight and allows for efficient transformations using the screw transformations.This simplifies the modeling, and reduces the risk for errors in the derivations.In addition, it is shown how the wave forces due to a wave spectrum can be computed using force RAOs for use in the equations of motion.
The proposed equations of motion can be used for controller design and for evaluating the performance of combinations of ship and crane for different operational scenarios.In this paper a simple controller based on feedback linearization is used for validation of the model.The method is demonstrated for a ship with a heavy crane with rotary joints and a heavy load, where the ship is operating in waves described by the JON-SWAP wave spectrum.
The paper is organized as follows: The theoretical background is presented in Sections 2, 3 and 4. First the equations of motion are presented in Section 2 using Kane's technique.Then screws and twists are presented in Section 3. Then wave modeling and force RAOs are presented in Section 4. The kinematics of a ship with a crane with three joints and a payload is presented in Section 5.The equation of motion is presented in Section 6, and a controller is presented in Section 7. Finally, simulation experiments are presented in Section 8, where the proposed equations of motion were validated for a ship with a heave crane with three joints, which was simulated using a simple controller and a JONSWAP wave spectrum.

Kane's equations of motion for a vehicle with a serial-link mechanism
In this section it is reviewed how the equations of motion for a serial link mechanism on a vehicle with 6 degrees of freedom can be developed using a Newton-Euler formulation based on the method of Kane and Levinson (1985).The presentation is based on Egeland and Sagli (1993).Consider a vehicle with 6 degrees of freedom and a serial-link mechanism mounted on the vehicle.The vehicle is called link 0, while the links of the mechanism are called link i for i = 1, . . ., n l .The links are connected with rotary joints with one degree of freedom so that joint i connects link i − 1 and link i.The North-East-Down frame {n} is supposed to be an inertial frame, and frame {i} is fixed in link i.The velocity of link i, i = 0, . . ., n l relative to frame n is denoted v i i , where the superscript i denotes that the vector is given in the coordinates of frame {i}.The angular velocity of link i relative to frame {n} is denoted ω i i .The equations of motion for link i is given by Egeland and Gravdahl (2002) (1) where m i is the mass of link i, I i i/i is the inertia tensor of link i about the origin of frame {i} in the coordinates of frame {i}, r i ig is the position vector from frame {i} to the center of mass for link i, f i i is the active force acting on link i, f i(c) i are the forces of constraint, m i i is the active torque acting on link i, m i(c) i are the torques of constraint, and ŵ is the skew symmetric form of a vector w.
The vector ν i i is defined by and the vector of joint variables is defined by θ = [θ 1 , . . ., θ n l ] T .Next the vector of generalized speeds is defined by Then the vector ν i i can be expressed in terms of the generalized speeds as and the generalized forces corresponding to the generalized speeds are given by According to the principle of virtual work, the forces and torques of constraints can be eliminated using Then, in view of (5) and the fact that the generalized speeds u i are independent variables, this implies that Then the forces and moments of constraint can be eliminated by combining the equations of motion (1,2) using ( 6), ( 8) and the Jacobi identity, which is written âb c + bĉa + ĉâb = 0 for any three vector a, b, c.This gives the equation of motion for the total system in the form where the link mass matrix and the matrix Finally, the equations νi i = P i u + Ṗi u and (5) are inserted, and the equations of motion becomes where the mass matrix of the system is while the Coriolis matrix is This formulation satisfies the condition that Ṁ − 2C is skew symmetric, which is seen from Sagatun and Fossen (1991 where The equation of motion as given by ( 12) is in closed form and has a simple structure.The only challenge is to find expressions for the P i matrices, which are given by and the associated time derivatives Ṗi .In the following it will be shown that expressions for the P i matrices can be found efficiently using screw theory, and that this leads to efficient computation of the P i matrices.Moreover, a simple geometric interpretation of the P i matrices can be given, which is useful to check the correctness of the solutions.

Definition
A screw s /b is an ordered pair of vectors which satisfies the screw transformation when the reference point of the screw is changed from the origin of {b} to the origin of frame {a}.Here p ab is the position vector from {a} to {b} McCarthy and Soh (2011).
In coordinate form the screw is written where u j and w j are vectors given in the coordinates of frame {j}.The screw transformation from reference {b} and coordinates in {j} to reference {a} and coordinates in {i} is given by is a screw rotation matrix, which transforms the coordinates from {j} to {i}, and is a screw transformation matrix, which changes the point of reference from {b} to {a}.Here p j ab is the position vector from {a} to {b} in the coordinates of {j}.It is noted that the screw transformation can be written Ri where R i j pj ab = pi ab R i j .The resulting screw referenced to {a} is

Lines as screws
A line can be described as a screw in terms of the Plücker coordinates of the line.The line is then given by l = ( a, m), where a is the direction vector of the line, and m = q × a is the moment of the line, where q is the vector from the reference point to an arbitrary point on the line.A line is a special type of screw where the two vectors of the screw are perpendicular to each other.
A line l j/j through the origin of frame {j} with reference point {j} will be given by l j/j = ( a j , 0), or in coordinate form in the coordinates of {j} as This line can be referenced to {i} by a screw transformation which gives l j/i = ( a j , p ij × a j ), or in coordinate form in the coordinates of {i} as

The time derivative of a line
Consider the time derivative of the line L i j/i given by ( 27) where it is assumed that the direction vector is fixed in frame {j}, so that da j j /dt = 0.The time derivative of the line is found by taking the time derivative of each of the vectors of the screw.First the time derivative of the direction vector is found to be

Twists
The linear and angular velocity of a rigid body B with a body-fixed frame {b} relative to a frame {a} can be described with the twist Angeles (1988) where ω ab is the angular velocity of frame {b} relative to frame {a}, and v ab/b = v ab is the velocity of frame {b} relative to frame {a} referenced to {b}, which is the velocity of the origin of {b} relative to frame {a}.
A twist is a screw, and is transformed according to the screw transformation when the reference point is changed.The twist can therefore be referenced to frame {a} as which is referenced to frame b.The twist in matrix form can be transformed according to or in vector form according to the screw transformation where Note that the screw transformation in terms of Ad T a b transforms both reference frame and coordinate frame from {b} to {a}, while a general screw transformation (24) may change reference frame independently from the coordinate frame.The twist t a ab/a is called the spatial velocity, while the twist t b ab/b is called the body velocity in Murray et al. (1994).

The twists of a composite displacement
Consider the composite displacement T a c = T a b T b c , and define the twists with vector forms Then the twist of the composite displacement is given in matrix form by tc ac/c = (T a c ) −1 Ṫ a c , which gives It is seen that the twist of a composite displacement is the sum of the screws of the individual displacements, where all screws are referenced to the origin of the same reference frame.Obviously, this applies also for the vector formulation t c ac/c = t c ab/c + t c bc/c , and will still apply after a change of coordinates, e.g., to frame {a} by t a ac/c = Ra b t a ac/c , which gives t a ac/c = t a ab/c + t a bc/c .

Link twists
In the Denavit-Hartenberg convention, the link transformation matrix from link i − 1 to link i is which is a rotation θ i about the z i−1 axis, followed by a translation d i along the z i−1 axis, and then a rotation α i about x i axis and translation d i along x i .It is assumed that joint i is rotary with joint variable θ i .Then the link twist will be is the line through the origin of frame {i − 1} along the z axis of frame {i − 1}.This means that L is the line along the joint axis of joint i.The twist of link i relative to the frame n referenced to i is found by adding the twists of the composite displacement from {n} to {i}, which gives Insertion of (44) then gives where L i j−1/i is found from the line L j−1 j−1/j−1 of joint j − 1 through a screw transformation from {j − 1} to {i}, which gives 4 Wave Modeling

The Wave Spectrum
Ship motion in waves can be computed from a wave spectrum, which describes the frequency distribution of the wave elevation ζ(t).Commonly used wave spectra include the Pierson-Moskovitch spectrum, the ITTC spectrum, and the JONSWAP spectrum Journee and Massie (2001); Perez (2005).In this section, it is shown how the wave spectrum is related to the power spectral density.This is based on the Fourier transformation Oppenheim and Verghese ( 2016) The wave elevation is an ergodic signal, and the autocorrelation is given by The Fourier transform of the autocorrelation is the power spectral density S ζζ (ω), which gives The average value of the square of the wave elevation is ζ(t) 2 = R ζζ (0), and it follows from (52) that The wave spectrum is by convention defined as It is seen from S ζζ (ω) = S ζζ (−ω) that the wave spectrum S(ω) satisfies The significant wave height H s is defined as The physical interpretation is that H s is approximately the average peak-to-peak value of the largest one third of the waves, which has formerly been used as the definition of H s .

Examples of Wave Spectra
The JONSWAP spectrum Hasselmann et al. (1973) describes the waves of the North Sea, and is given by where Here ω p is the peak frequency, and default values are α = 0.0081 and while γ can be set to γ = 3.3.The JONSWAP spectrum can be related to the significant wave height H s by Journee and Massie (2001); Fossen ( 2011) where T p = 2π/ω p is the period corresponding to the peak frequency.The JONSWAP spectrum is expected to give good results for 3.6 < T p / √ H s < 5 DNV GL (2017)

Long crested irregular sea
The standard technique for simulating a given wave spectrum S(ω) is to approximate the wave spectrum with discrete spectrum S N (ω) where the wave elevation ζ(t) is the sum of N single frequency components ζ i (t) according to where Z i is the amplitude and i is a random phase angle of frequency component i.The frequency ω i is generated as a random number in the interval and the amplitudes Z i are selected as Then the resulting discrete wave spectrum will be where δ(•) is the Dirac delta function.It follows that S N (ω) approximates S(ω) in the sense that

Short crested irregular sea
The wave spectrum is further discretized with M different wave directions using a spreading function where χ 0 is the dominant wave propagation direction and χ j is randomly chosen in the interval Then for each frequency ω i the wave elevation is the sum with direction χ j , phase ij and amplitude

Wave forces from force RAOs
The wave load on a ship is described as the superposition of two effects Faltinsen (1990).The first is the wave excitation due to the incoming waves acting on a nonmoving ship.These wave forces include the Froude-Krylov forces and the diffraction forces.These forces are given in the hydrodynamic frame {h}, which has the xy plane in the undisturbed sea surface, and origin determined by the undisturbed position of the ship.
The second effect is due to the ship moving on a sea with no incoming waves.This includes the radiation forces and the restoring forces.
The Froude-Krylov forces and the diffraction forces due to the incoming waves acting on a nonmoving ship are calculated from the wave components ζ ij (t) using force RAOs (Response Amplitude Operators) Faltinsen (1990).The force RAO F k (ω, χ) in the degree of freedom k is a transfer function that is given in terms of its amplitude |F k (ω, χ)| and phase ∠F k (ω, χ), which are calculated from the geometry of the ship hull.The resulting wave forces are the diffraction forces and the Froude-Krylov forces.The wave force in degree of freedom k is found from where the phase is µ ij = ij + ∠F k (ω i , χ j ).The resulting wave forces τ w,k are given in the {h} frame.

Seakeeping Model
The equation of motion for a ship moving in waves can be given in frame {h} as the seakeeping model Smogeli et al. (2005); Ross et al. (2006).
where ξ is a vector of the generalized coordinates for the ship, comprising the three position coordinates and three roll-pitch-yaw Euler angles.M 0 is the associated mass matrix of the ship in frame {h}.The vector τ w is the generalized wave induced Froude-Krylov and diffraction forces given by the force RAOs in equation ( 71), while τ R is the radiation force vector due to hydrodynamic added mass, damping and restoring forces.τ is the vector of generalized control forces, while τ (c) denotes the forces of constraint.
The radiation force vector τ R is given by where A(ω) is the frequency-dependent added mass, and B(ω) is the frequency-depended damping, and Gξ is the restoring force.
The frequency dependence of the added mass and the damping in ( 73) is due to the common procedure in seakeeping analysis to treat single-frequency motion.In Cummins (1962); Ogilvie (1964) it is explained that the frequency dependence is due to the memory effect due to the wave pattern that is set up by the motion of the ship on an undisturbed surface.Moreover, for general motion, this can be modeled by introducing a convolution term µ in the expression for the radiation force, which gives where is a convolution term with kernel K(t).In Kristiansen et al. (2005) it was shown that this convolution term can be represented by a state-space model of the form

Ship model in vessel-fixed frame
The equations of motion for the ship is formulated in the vessel-fixed frame {0}.Following Ross et al. (2006), it is assumed that the rotation from frame {h} to frame {0} is small, that the hydrodynamic coefficients A(ω) and B(ω) are computed in CG, and that the origin of frame {0} is located at the CG.Then, the hydrodynamic parameters and hydrodynamic forces in frame {h} will be the same as in frame {0}, while ξ = ν 0 0 and ξ = ν0 0 .This gives the model in the vessel-fixed frame {0} as where M 0,A = M 0 + A(∞) and D = B(∞), and the generalized forces are given by τ 0 = τ w + τ thr + τ (c) 0 , where τ thr are the controlled thruster forces of the ship, and τ (c) 0 are the forces of constraint due to the crane.The generalized wave forces τ w due Froude-Krylov and diffraction effects are given by the force RAOs in equation (71).

Marine vessel
The position of the vessel frame {0} with respect to the {n} frame is p n n0 = [x 0 , y 0 , z 0 ] T .The rotation matrix from {n} to {0} is given by the roll-pitch-yaw Euler angles Θ = [φ, θ, ψ] T as where R z (ψ) is the yaw rotation about the z axis, R y (θ) is pitch rotation about the current y axis, and R x (φ) is the roll rotation about the current x axis.The linear velocity of {0} relative to {n} in the coordinates of {0} is v 0 0 = [u, v, w] T , while ω 0 0 = [p, q, r] T is the angular velocity of {0} relative to {n}, in the coordinates of {0}.
The position and orientation of the vessel is represented in the coordinates of {n} by the vector while the linear and angular velocity is given in the coordinates of {0} by The kinematic differential equation of where and

Crane with payload
The crane and payload is modelled with 5 degrees of freedom with three joints in the crane, and two degrees of freedom for a swinging load that is attached to the tip of the crane with a wire of constant length.This is described with 5 links that are connected with rotational joints in a serial arrangement so that link i − 1 is connected with link i with joint i.The ship is considered to be link 0, which is to link 1 with joint 1 of the crane.
The position and orientation of link frame {i} relative to link frame {i − 1} is given by the homogeneous transformation matrix coordinates of the origin of frame {i} relative to the origin of frame {i − 1}.In the Denavit-Hartenberg convention this is described as a composite displacement of a rotation θ i about the z i−1 axis, a translation d i along the same axis, a rotation α i about the x i axis followed by a translation a i along the same axis.Note that θ i is the joint variable, while the parameters d i , α i and a i describe the geometry of the link transformation.This gives a homogeneous transformation matrix described by the Denavit-Hartenberg parameters θ, d, α and θ in Table 1.
The configuration of the crane and payload system, is given by the generalized coordinates This is the vector of crane joint angles q c = [q 1 , q 2 , q 3 ] T and angles q p = [q 4 , q 5 ] T describing the payload.
6 Equations of motion

Generalized speeds and projection matrices
The generalized speeds for the system are given by the 11-dimensional vector and is related to The projection matrices P i for i = 1, . . ., 5 are then found from the expressions of the link twists t i ni/i , using the transformation where It is noted that H = H −1 , and that Then, from (47) it is seen that the P i matrices are given by where the expressions for the lines L i i−1/i are given by (48).

Gravity forces and moments
The active forces and moments of link i = 1, . . ., 5 include the gravitation forces which result in the generalized force where and g = 9.81 is the acceleration of gravity.

Equation of Motion
The equations of motion as given by Kane and Levinson (1985) is then found from equation ( 9) to be This can be written where τ contr denotes the torque applied at the crane joints and

Link parameters
The mass matrix M i i of link i is given by ( 10), and the matrix C i i is given by ( 14).The parameters of the mass matrices are the positions of the center of mass, given by and the inertia tensors I i i/i given by where and the parallel axes theorem has been used.

Controller objective
The control objective is to control tip position p n n,tip of the crane relative to {n}.This is done to avoid excessive oscillations of the load.Given a desired position d p n n,tip of the crane tip corresponding to in frame {0}, the error of the tip position given in coordinates of frame {0} is The crane tip position relative to and expressed in frame {0}, is where the forward kinematics is represented as a mapping f : q c → p.

Controller
The equation of motion for q c in joint space is from equation ( 99) where We can rewrite the equation of motion in (110) in terms of p ∈ R 3 , by using the Jacobian of the mapping in (109), where It follows that We can now substitute these expressions into (110) and pre-multiply by J −T = J −1 T to obtain where ) Using the computed torque technique Murray et al. (1994), the controller becomes where M22 is positive definite, and it follows that the error dynamics are asymptotically stable and given by ë + K d ė + K p e = 0 (117)

Simulation
Simulations were performed with significant wave height H s = 4 m, peak frequency ω p = 1.3 rad/s and the mass of the crane given by m 1 + m 2 + m 3 = 296 metric tons.The Denavit-Hartenberg parameters of Table 1 were d 1 = 4 m and a 2 = a 3 = 2.5 m.The simulations presented in Sections (8.1) and (8.2) were done with a payload mass of m 5 = 150 metric tons, while in the simulations reported in Sections (8.3) and (8.4) the payload mass was m 5 = 250 metric tons.The vessel was dynamically positioned with desired values surge u = 0, sway v = 0 and yaw ψ = 0.
In the simulations, the ship was a supply vessel with 82.5 m between the perpendiculars, 6 m draught, 8 m breadth and a mass of m 0 = 6, 362 metric tons.The hydrodynamic coefficients of the vessel, the force RAOs, the parameters of the radiation force model and the rigid-body mass matrix of the vessel was taken from the Marine Systems Simulator (MSS), Fossen and Perez (2004).The Marine Systems Simulator has a SIMULINK model of the ship, whereas the simulator in the present paper is developed in MATLAB.Our MATLAB simulator of the ship without the crane was validated by extensive simulation where it was compared to the MSS SIMULINK simulator, and it was found that the two simulators gave the same results.

Movement of the crane tip parallel to the z-axis of the inertia frame
A simulation study was done with the crane links (the translations a 2 and a 3 ) along the positive sway direction.The simulation studied the roll response of the vessel due to steps in the z coordinate of the crane tip relative to and expressed in {n}.The roll responses are illustrated in the lower windows in Figures 3 and  4, while the z coordinate of the crane tip (relative to and expressed in {n}) is illustrated in the upper windows in Figures 3 and 4. The red graphs in Figures 3  and 4, illustrate a simulation where the desired position d p n n,tip in ( 107), was constant.The blue graphs in Figure 3, illustrate the result from a simulation where the crane tip was lifted up 0.5 m (along negative z-axis of frame {n}) after 10 seconds and then lowered down to initial position after approximately 15 seconds.This caused the roll φ of the vessel to positive peak (increased roll angle) due to the lift of the crane tip, and negative peak (decreased angle) due to the crane tip lowered to initial position.The difference (between the simulations with stationary crane tip and non-stationary crane tip) in roll angles after 15 seconds, was due to small oscillations of the pendulum caused by the steps of the crane tip.The blue graphs in Figure 4, illustrate the result from a simulation where the crane tip was lowered down 0.5 m (along positive z-axis in frame {n}) after 10 seconds and then lifted up to initial position after approximately 15 seconds.This caused the roll φ of the vessel to negative peak (decreased roll angle) due to the crane tip lowered 0.5 m, and negative peak (increased angle) due to the crane tip lifted up to initial position.
The difference in roll angle after 15 seconds, was due to small oscillations of the pendulum caused by the steps of the crane tip.

Movement of the crane tip parallel to the xy-plane in the inertia frame
A simulation study was done with the crane links (the translations a 2 and a 3 ) initially along the positive sway direction.The simulation studied the vessel motions, roll φ and pitch θ, due to the crane moved from a pose where the crane links were along sway direction, to a pose where the crane links were along surge direction.This movement was achieved by moving the crane tip's xy-coordinates linearly (see blue graphs in Figure 5), relative to and expressed in {n}.The z coordinate of the crane tip relative to and expressed in {n}, was kept constant.The vessel motions are illustrated in Figure 6, the vessel's roll angle decreased to φ ≈ 0 and the pitch θ increased slightly due to the movement of the crane links.The red graphs in Figures 5 and  6, illustrate a simulation where the desired position d p n n,tip in (107), was constant.

Increased payload weight
A simulation study was done with the crane links (the translations a 2 and a 3 ) along the positive sway direction and the crane tip kept stationary relative to and expressed in {n}.The simulation studied the vessel motions, heave z and roll φ , due to an increased payload weight.In Figure 7, the blue graphs represent the vessel motions due to a payload weight m 5 = 250 metric tons, while the red graph shows result with payload weight m 5 = 150 metric tons.As seen in Figure 7, the heave z of the vessel increased with approximately 0.1 m and the roll of the vessel oscillated with higher peakto-peak value, with an increased of payload weight of 100 metric tons.

Oscillation of pendulum
The normalized pendulum vector n n 5g = r n 5g / r n 5g from the origin of frame {5} to the center of the payload-mass m p = m 5 expressed in {n}, is from equation (100) described by the Euler angles φ y and φ x (see Figure 8) which gives the Euler angles as φ y = Atan2 n n 5g (1), n n 5g (2) 2 + n n 5g (3) 2 (120)   A simulation study was done with the crane links (the translations a 2 and a 3 ) along the positive sway direction, the crane tip kept stationary relative to and expressed in {n} and the payload weight was m 5 = 250 metric tons.This simulation studied the response of the vessel and the crane tip due to different initial conditions of φ x and φ y .Figure 9 shows the oscillation of the pendulum, Figure 10 shows the crane tip's position relative to and expressed in {n} and Figure 11  As seen in Figure 10, the swinging pendulum (red graph) caused more disturbance on the crane tip than in the case with non-swinging pendulum (blue graph).As seen in Figure 11, the swinging pendulum (red graph) affected heave z and pitch θ of the vessel less, while it affected the roll φ more than in the case with a non-swinging pendulum (blue graph).The swinging pendulum caused a less smooth roll motion of the vessel.

Conclusion
The equations of motion for a ship with a heavy crane has been modeled using a Newton-Euler formulation where the forces of constraint are eliminated using the method of Kane's equation of motion.The use of generalized coordinates were avoided by relying on generalized speeds.The projection matrices associated with the generalized speeds were derived using screw theory, which lead to simple expressions with clear geometric interpretation for the projection matrices.The equations of motion were tested in extensive simulations using a simple controller based on feedback linearization to stabilize the tip of the crane.The simulations clearly showed that the ship responded in a physically meaningful way to the motion of the crane.
the velocity of a point fixed in {b} that passes through the origin of {a}.The coordinate form of the twist is related to the time derivative of the homogeneous transformation matrixMurray et al. (1994)

Figure 3 :
Figure 3: Upper window: The z coordinate of the crane tip, relative to and expressed in inertia frame {n}.The blue graph represents a simulation where the crane tip is lifted up 0.5 m (along negative z-axis of inertia {n}) after 10 seconds, and then lowered to initial position after 15 seconds.The red graph represents a simulation where the crane tip was kept stationary.Lower window: Blue graph represents the roll φ angle of the vessel due to disturbances (lift and lower operation) of the crane tip illustrated as blue graph in the upper window.Red graph represents the roll φ angle of the vessel due to stationary crane tip.

Figure 4 :
Figure 4: Upper window: The z coordinate of the crane tip, relative to and expressed in inertia frame {n}.The blue graph represents a simulation where the crane tip is lowered 0.5 m (along positive z-axis of inertia {n}) after 10 seconds, and then lifted to initial position after 15 seconds.The red graph represents a simulation where the crane tip was kept stationary.Lower window: Blue graph represents the roll φ angle of the vessel due to disturbances (lift and lower operation) of the crane tip illustrated as blue graph in the upper window.Red graph represents the roll φ angle of the vessel due to stationary crane tip.

Figure 5 :
Figure 5: x, y and z coordinates of the crane tip relative to and expressed in coordinates of inertia frame {n}.The red graph represents a simulation where the crane tip was kept stationary.In this case, the crane links (a 2 and a 3 ) lies in the yz-plane of inertia frame {n}.The blue graph illustrates a simulation of the crane tip, where the crane links lies initially in the yz-plane of {n}, from 10 to 30 seconds the crane tip moves linearly to a position where the crane links lies in the xzplane of inertia {n}.

Figure 6 :
Figure 6: Upper window: Roll φ angle of the vessel due to the crane tip trajectories illustrated in Figure 5. Lower window: Pitch θ angle of the vessel due to the crane tip trajectories illustrated in Figure 5.

Figure 7 :
Figure 7: Upper window: The red graph represents the heave z of the vessel due to payload weight m 5 = 150 metric tons.The blue graph represents the heave z of the vessel due to payload weight m 5 = 250 metric tons.Lower window: The red graph represents the roll φ of the vessel due to payload weight m 5 = 150 metric tons.The blue graph represents the roll φ of the vessel due to payload weight m 5 = 250 metric tons.

Figure 10 :
Figure 10: x, y and z coordinates of the crane tip relative to and expressed in coordinates of inertia frame {n}.The desired crane tip is kept stationary relative to and expressed in {n}.The blue graph is the crane tip position due to orientations [φ x , φ y ] of the pendulum with initial conditions [φ x (0), φ y (0)] ≈ [0 • , 0 • ].The red graph is the crane tip position due to orientations [φ x , φ y ] of the swinging pendulum with initial conditions [φ x (0), φ y (0)] ≈ [−45 • , 45 • ]

Table 1 :
DH-parameters -Crane-payload Link a i