Application of dynamic optimisation to the trajectory of a cable-suspended load

The paper discusses a dynamic model of the system consisting of an on-board hoisting winch on a moving vessel, cable, and load. The model accounting for large rope sag and hydrodynamic drag force was used to solve the problem of dynamic optimisation. The essence of which is what angle of rotation should be selected for the hoisting winch to ensure that the load during the defined movement of the vessel shall remain at the set distance from the seabed, which may be undulating. To solve this problem, the nonlinear optimisation methods were used. The presented calculation results show the efficiency of the developed 3D system model and its possible application to solve other similar optimisation problems.


Introduction
In contemporary offshore engineering systems, increasingly complicated devices are used that enable operation at great depths and in difficult conditions. To simulate operation of such devices, complex dynamics models, applied in commercial software packages (Rifflex, Orcaflex, Proteus), are required. The models of cables Ł. Drąg (B) University of Bielsko-Biała, ul. Willowa 2, 43-309 Bielsko-Biała, Poland e-mail: ldrag@ath.bielsko.pl or raisers, which frequently occur in the offshore engineering systems, have to account for their large sags or bowing, variable configuration, and loads caused by water impact (wavy motion, buoyancy force, the presence of sea currents, hydrodynamic drag, volume of associated water, etc.). To divide the flexible structure into discrete elements, the following methods are used: finite element method [1,2], lumped mass method [3], segment method [4][5][6]. Generally speaking, for modelling of flexible structure in the offshore engineering systems the methods used include those developed, inter alia, in the analysis of multibody system dynamics with recognition of flexibility. A broad review of methods used is presented in the paper [7]. For modelling of flexible structures the rigid finite element method (RFEM) is also used [8][9][10]. The modification of that method is presented in this paper. A similar approach was applied in the study [11] where the flexible segment model was presented. The model developed there assumes, however, dependency of coordinates of ith element on the coordinates of preceding elements. This assumption is redundant in the proposed RFEM modification due to the formulation of appropriate equations of geometric bonds. The issue of method selection used for dividing susceptible members such as cables or raisers, into discrete elements, gains special importance in solving the problem of controlling the propulsion drives of marine equipment. The numerically effective model of the system dynamics, reflecting its characteristics with acceptable accuracy, enables the formulation and even solution of appropriate optimisation tasks or their sequence for a specific set of control parameters. The results may then be provided as inputs for artificial neural networks and consequently used for control. The questions of dynamic optimisation of cables and raisers are quite rarely presented in the literature. This problem for cranes operated in marine conditions is discussed in papers [12,13]. In the above papers, as in many others, the cable is modelled as an elastic and damping element without mass. Such a simplification is not used in the aforementioned commercial software packages. That method was also applied to partition the raiser into discrete elements [14]. Also the monograph [15] formulated the problem of dynamic optimisation involving the compensation of wave motion during pipe laying using a drum method with the application of rigid finite element method (RFEM).
This paper presents a 3D model of the elastic element with bending and possibly torsional compliance, enabling the mapping of large cable and raiser sags. The model developed through the application of own modification of RFEM eliminates translation rigidity from further consideration (shearing in both directions and longitudinally one) and optionally torsion. The formulation of the rigid finite element method presented in this paper differs considerably from the previous formulations. The motion of each rigid finite element (rfe) is described by three displacements, two rotation angles reflecting bending, and one rotation angle reflecting torsion (when necessary). Thus, each rfe has either five or six degrees of freedom. The elements are connected by means of three geometrical constraint equations due to the necessary continuity of displacements at the points connecting the elements. Such an approach has already been applied by the author but for planar systems with consideration of bending [16,17] or bending and longitudinal flexibilities. In the formulation presented in [16], each rfe has four degrees of freedom: two translations, bending rotation, and elongation, while the elements with three degrees of freedom (two translations and bending) are used in [17]. Constraint equations are formulated for both the planar and spatial models; however, there are 2*(n+1) constraint equations in planar models (2D) described in [16,17], while for the spatial models (3D) considered in this paper, 3*(n+1) constraint equations have yet to be formulated. In both formulations, n+1 is the number of rigid finite elements into which the flexible link is divided.
Thus, the models of spatial slender systems presented in this paper are generalisations of the planar models presented in [16,17]. Similarly to the approach from [17], the longitudinal flexibility is omitted, which enables the equations of motion to be integrated using relatively large integration step. So, the presented modification of RFEM is characterised by high numerical efficiency, which enables its use for dividing the cable into discrete elements while solving a dynamic optimisation problem. The dynamic optimisation problem considered further in this paper is shown in Fig. 1.
It involves determination of the angle of rotation of a hoisting winch ϕ(t), which shall ensure that the load is at a constant distance from the surface of the sea bottom. Figure 2a shows elastic beam element of constant cross section, which during primary division was partitioned into n equal sections ( Fig. 2b) with the length of:

Susceptible element model
During secondary division, in the middle of these sections spring damping elements were placed (sde) ⊗ as shown in Fig. 2c. Parts of unit located between unit ends and sde, and between sde, are hereinafter regarded as rigid elements (rfe). Their number is n + 1 (numbered from 0 to n). As a result of that division, the unit is replaced with the n + 1 system of rigid solids (rfe) assuming mass properties of the unit and n non-dimensional spring damping elements (sde) without mass assuming unit's elastic (bending) properties.
The operational procedure when the unit is of variable section or concentrated loads are attached is described in detail in papers [9,18].
The movement of each rfe is defined by the vector components: where x i , y i , z i are coordinates of point A i and ψ i , θ i , ϕ i are Euler's angles Z Y X , presented in Fig. 3, that are rotations around system axis {i} combined with rfe i. Such angles are used in aviation mechanics and are called heading, attitude, and bank [19]. While applying homogeneous transformations, the coordinates of the point defined in the local coordinate system {i} combined with rfe i may be expressed in the inertial system according to the following formula: where Elements of B i matrix are defined unequivocally by the components of the vector q i . The kinetic energy of rfe i can be defined with following dependence [9]: where H i = m i r i r T i dm i . It should be noted that elements of H i matrix are constant. Under the assumption that axes y , z of the {i} system are main axes, whereas x i axis is symmetrical axis, the H i matrix referred to in (4) shall assume the following form: where Following the procedure referred to in (17), the Lagrange operators from kinetic energy of rfe i may be presented as follows: where The calculation of A i matrix elements and h i vector using formulae referred to in (6) is very timeconsuming in terms of computer processing time. Once derivatives of B i matrix have been calculated, it would be desirable to calculate products and traces of matrices present in formulae of a i, jk and h i, j . The high number of elements present in (5) is at the same time equal to nil. Therefore, in this paper a similar method was used to that referred to in [9]: products and traces of matrices, referred to in (5), were calculated analytically.
Individual rigid finite elements are combined to form a kinematic chain by means of geometric bonds (Fig. 4a). Forces present in links are shown in Fig. 4b.
To account for reactions in the links of rigid finite elements (rfe), it is necessary to introduce generalised forces. Generalised forces imposed on rfe i, originating from F i and from F i+1 , might be presented as follows: Since Finally, the generalised forces caused by forces affecting rigid finite elements rfe i reactions in A i and A i+1 may be noted as: connected with the rigid finite elements (rfe i) link with i −1, after their double differentiation, shall assume the following form: where Given that movement of point A n+1 , r n+1 = r n+1 (t) is known, the additional equation of bonds may be noted in acceleration-related form as: where −l n big(ψ 2 n cψ n cθ n −2ψ nθn sψ n sθ n +θ 2 n cψ n cθ n −l n ψ 2 n sψ n cθ n −2ψ nθn cψ n sθ n +θ 2 n sψ n cθ n l nθ It should be noted here that given that the hoisting winch drum rotates, the cable can reel/unreel on the drum. As a consequence, the length l n of the rfe n changes and the variable figure of n, being the number of elements into which the cable was divided. The potential energy of terrestrial gravity forces of rigid finite elements (rfe i) may be presented as follows: where g-acceleration of gravity, y C i -y coordinate of mass centre C i . Using the denotations in Fig. 2, the following calculation can be made: When flexural (bending) and torsional compliance are considered for discrete element, the moments caused by elastic strain of spring damping elements (sde) ⊗ must be introduced to element motion equations. Sde 1642 Ł. Drąg Generally, it is assumed [9] that moments M i,ψ , In the case of linear physical relationships, it can be assumed that: where k i,ψ , k i,θ , k i,ϕ are flexural and torsional rigidity factors. Angles ψ i , θ i , and ϕ i can be calculated based on the following formula: If moments M i,ψ , M i,θ , and M i,ϕ are known, then their effect causes the occurrence of appropriate generalised forces in the motion equations. According to Fig. 5a, rigid finite elements (rfe i − 1) shall be affected by the moments M i,ψ , M i,θ and M i,ϕ , and rigid finite elements (rfe i) by -M i,ψ , -M i,θ , -M i,ϕ . So spring damping element (sde) ⊗ of i number shall lead to occurrence of generalised forces: A. In rfe i − 1 motion equations and M i,ϕ may be roughly assumed to be equal to M i,ψ , M i,θ , and M i,ϕ , or can be calculated precisely using the following formula: The motion equations of rfe 0, …, n (so far without taking into account the impact of aquatic environment) can be presented as follows: This is a system of 5·(n+1) differentiation ordinary equations of second order with 5·(n+1) unknowns, which are components of the vectors q 0 ÷ q n , and 3·(n+1) unknowns, which are components of reaction vectors F 1 ÷F n+1 . The above equations have to be supplemented with 3·(n+1) bond equations (11) and (12): Before the presentation of the method how to solve the system of equations (20) and (21) the water impact forces were introduced into the system.

Consideration of aquatic environment impact
The aquatic environment impact forces, for examined flexible elements, of large length and small cross section, can be introduced using Morison equations [20]. If forces caused by flow inside the element are disregarded, then consideration of aquatic environment impact requires taking into account the following: -buoyancy force -viscous resistance forces -force of inertia.
Below a method is presented for introducing the above forces into the element motion equations.

Buoyancy force
If ith rigid finite element (rfe) is fully immersed in water, then the buoyancy force ( Fig. 6) can be calculated as: where ρ w is water specific gravity and C i is crosssectional area for rigid finite element (rfe i).
Meanwhile, the generalised force will be added in the motion equations for ith element: The circular unit section which is dξ in length is affected by the hydrodynamic resistance force [21] ( Fig. 6): where D x ,D yz -longitudinal and lateral resistance factors, ⎦ is relative velocity vector for dξ element in the coordinate system for rigid finite element rfe {i} , The generalised force that occurred due to the effect of df H,i force corresponding to coordinate q i,k is presented by the following dependence: The integrals present in (26) may be calculated using Gauss formulae. The generalised forces generated by water resistance shall occur on the right-hand side of rigid finite elements (rfe) motion equations:

Force of inertia
When applying the Morison equation, it may be assumed that the force of inertia influencing ith rfe ele-ment, which is dξ long, is defined in the inertial system as: where df w I,i = e w I a w dξ, e w I = C M C i ρ w , df a I,i = e a I adξ, e a I = −(C M − 1)C i ρ w , C M -coefficient, a w -water acceleration, a-acceleration of rigid finite element (rfe i) with the length of dξ .
The first component in (28) shows forces caused by water movement (acceleration), and the second one by movement (acceleration) of the rigid finite element dξ rfe i.
After df I,i integration over the length of element l i , it may be stated that generalised forces caused by water acceleration are as follows: The generalised forces caused by effect of forces df a I,i can be calculated as: and after the following transformations: where P i,a = ⎡ ⎢ ⎢ ⎣ The generalised forces, referred to in (31), may be presented as follows: where Taking into consideration the generalised forces originating in the aquatic environment, one can record movement rfe 0 ÷ n in the following manner:

The equations of bonds (21) remain unchanged.
Fourth-order Runge-Kutta method was applied for integration of motion equations (33) with the constant time step of numerical integration. Since bond equations were transformed into acceleration-related form, to stabilise them, the Baumgart method was used [22].

Model validation
Before proceeding to solving the optimisation problem, the cable model was validated from the perspectives of static and dynamic calculations. A comparison was made of the results obtained using the model with accurate solution for catenary line, disregarding the rigidity of flexible cable, and the results generated by the Abaqus commercial software package. The impact of torsion on the results of numerical simulations and computation time was also examined.

Static problem
For static calculations, own results were benchmarked against results obtained analytically. Exact values of x, z coordinates and force T in the cables were calculated using the following formula [3]: where w = ρ Ag-weight of cable unit length. The numerical computations were made for the cable with the length of L = 300 m, fixed cross section A = 0.07 2 /4 m 2 , Young modulus E = 10 11 N/m 2 , density ρ = 6.5 · 10 3 kg/m3, and loading with forces: H = 20000 N, V = 50000 N (Fig. 7a). The comparison of the curves z(x) obtained for the analytical solution and for the proposed model, for the line division n = 15 elements, is presented in Fig. 7b. As can be seen, the proposed method correctly reflects the shape of the rope chain.
The values of relative percentage error ε γ were calculated using the following formula: where f c γ -approximate value using the presented method, f a γ -value calculated using formulae (34a) and (34b), for s i = i j=0 l j and i = 0...n, γ ∈ {x, z}. Maximum relative percentage error ε max γ was calculated using the following formula: The analysis of the values of errors shown in Fig. 8a showed that the maximum differences in displacements in x, z directions did not exceed 0.5 %. The highest values of errors ε x , ε z are present at points where rope chain was characterised by the biggest curvature. The division of the rope chain into a higher number of elements reduces computational error (Fig. 8b).

Dynamics problem
In order to validate the dynamics model for large sags, the computational results were compared with the results obtained using Abaqus commercial software package. The calculations made using the commercial software package considered rope chain division into n=100 elements with the same characteristics as under static problem section. Then, load was applied to the system as follows: end E of rope was subject to effect of forces F = F x , F y , F z T (Fig. 9). The curves of F x and F y forces are presented in Fig. 9b, and the component F z = −3000 N was assumed to be constant. Numerical simulations were started in the position of free sag of rope. Total analytical time was 100 s. Figure 10 shows the trajectory of E point in the x y plane.
Similarly to the static problem, very high conformity was attained between own model results and those obtained from commercial package. Maximum relative error for displacement towards x and y for the rope divided into n = 100 elements is less than 0.5 %. The total simulation time using fourth-order Runge-Kutta method after the adoption of integration step h = 2 · 10 − 4 s was 5242 s. The fact that obtained results are highly compliant with those of Abaqus commercial software package proves that the rope model  percentage errors a ε x , ε z , b ε max x ,ε max z Fig. 9 Analysed system a forces applied at the end of rope, b curves of F x and F y forces proposed in this paper may be used in the dynamic optimisation.
The impact of torsion on the E point trajectory and calculation time was examined. An important advantage of the model presented earlier is that omission of torsion requires only the omission of ϕ i (q 1,6 ) coordinate in the formula (2) and then reduction of matrix and vectors used in Eq. (21) and (33) from 6 × 6 and 6 × 1 to 5 × 5 and 5 × 1. It should be also noted that omission of shearing would enable drastic reduction of the integration step. For flexible elements such as a rope, the torsional rigidity factor is calculated using the following formula: where G-modulus of volume rigidity; I 0 i -polar moment of inertia for cross-sectional area; -length of initial element as shown in Fig. 1. Given small values of mass moments of inertia for rigid elements (rfe) of the rope, high values of c i,ϕ lead to the occurrence of high frequency of torsional vibrations. Meanwhile, this necessitates the reduction of the step of integration. Thus, the elimination of coordinates ϕ i enables increase in step of integration of system displacement equation, and thus reduction of calculation time. Table 1 shows the impact of the length of step of integration on cal-   It may be seen that the E point trajectory error even with the step of integration 200 times greater does not exceed 0.15 %. Therefore, the torsion was omitted in the optimisation problem analysis shown below.

Optimisation problem: stabilisation of rope end at set depth
The system model consists of the drive system drum and the rope. It was assumed in the optimisation that vessel displacement at the place where rope hoist- ing winch is located would be known. Assuming that R radius of the drum is small, it would enable calculation of coordinates Because hoisting winch drum may rotate, the rope's length varies, which enables adjustment of load position (height δ). In the examined problem of dynamic optimisation the course of ϕ(x A (t), y A (t)) = ϕ(t) drum rotation was selected, which shall ensure that despite horizontal movement of the vessel, point A 0 shall stay at the same distance δ from the sea bottom (Fig. 11). The vertical displacements of point A were omitted. As shown by the results published in [15], the vertical movement may be compensated by adoption of: where ϕ(z A ) = z A r D , r D is radius of winch drum. The objective function was defined as follows: where z(x, y)-function describing seabed topographical features, z sb -seabed depth for A 0 at t = 0, The function describing seabed topographical features was defined as follows: where x 0 , y 0 , k, p, q are parameters of the function describing sea bottom profile. Calculation of the value of functional involves integration of motion equations (33) for adopted functions defining x A (t), y A (t). To analyse the problem as nonlinear optimisation, it was assumed that ϕ(t) function is defined as third-degree catenation function. The catenation function ϕ i = ϕ(t i ) values are searched in equidistant points of the 0; T interval: where T -motion simulation time, p-number of subintervals in 0; T . It is assumed that value ϕ 0 = ϕ(0) = 0 is known and that it results from the system initial configuration. The objective of minimisation of the function referred to in (39) is to search for such values ϕ 1 ÷ ϕ p being components of vector f= ϕ 1 , ϕ 2 , ..., ϕ p T which fulfil limitations ϕ i,min ≤ ϕ i ≤ ϕ i,max and minimise function = (f) referred to in (39).
To solve the problem of nonlinear optimisation, the creeping simplex method was used. The motion equations were integrated assuming the step of integration h ≤ 0.1 s securing proper stability and accuracy of solution.

Results of numerical simulations
In order to depict displacements towards x axis of E end of the rope sitting on motionless drum of the onboard hoisting winch on the vessel moving towards x,  numerical simulations for data shown in Table 2 were made. As a result of hydrodynamic resistance having impact on the rope, the end of rope rises (Fig. 12). The selection of the angle of rotation of the hoisting winch drum ϕ(t) in the problem of stabilisation of the rope end at a set height above seabed shall concern the compensation of the rope end lifting caused by vessel movement and elimination of changes caused by set seabed profile. It should be noted that when hoisting winch is immobilised, the E line end is displaced along z axis by approximately 7 m. This is caused by the fact that the rope length is maintained unchanged, despite the effect of hydrodynamic resistance.
Total simulation time was T = 210 s. The values of decision variables in the dynamic optimisation problem were selected at equidistant points, dividing 0; T interval into p = 14 subintervals. It was assumed that motion of point A n+1 (vessel) along x and y axes is described as in Fig. 13a. Thus, the trajectory is the same as in Fig. 13b.
Consecutive drawings 14 ÷ 16 show E point trajectory, before and after optimisation for both examined options of seabed profile under the assumption that the rope end is kept at a set distance δ from the seabed.  The value of objective function after optimisation, for A and B options, was A = 0.076 and B = 0.105, respectively, which means that the average deviation of the rope end from set distance δ = 5 m above the seabed was below 0.11 m. The result of optimisation for examined seabed profiles should be considered as good, given the height at which E point of rope is raised as a result of vessel movement (dotted lines in Fig. 14) (Figs. 15, 16).
The range of angles of the rotation of the hoisting winch drum ϕ(t) selected during optimisation process for option A and option B is shown in Fig. 17. In both cases, the angle of drum rotation was selected in such manner as to eliminate changes caused by the raising of E end of the rope due to hydrodynamic resistance and change of seabed profile.
Shape and position of the rope end at selected points in time depending on seabed profile are presented in Fig. 18. To describe vessel movement in Fig. 13, it was noted that rope end rises maximally approx. 12 m from seabed. Following the application of optimisation process (selection of the hoisting winch rotation angle ϕ(t)), the distance is maintained at almost constant level-approximately δ = 5 m above seabed.
The results of dynamic optimisation presented above indicate that through appropriate selection of the hoisting winch drum rotation angle, it is possible to control rope end position, which will guarantee keeping it at a set level above seabed. The major factor influencing the optimisation process time is the number of decision variables. Therefore, it has been planned to subsequently develop an algorithm ensuring the appropriate

Summary
The 3D model presented in the paper concerned with the dynamics of flexible member with flexural compliance was developed through modification of the rigid finite element method. An important feature of the model is a specific selection of generalised (global) coordinates. Each rigid element (rfe) has 5 or 6 DoF. Unlike the conventional rigid finite element method (RFEM) where rigid elements are combined by means of spring damping elements with 6 rigidity coefficients (3 translational, 2 flexible, and 1 torsional), in this paper elastic and damping elements have only 2 or 3 rigidity factors (2 flexural and sometimes 1 torsional). The identical displacements at rigid finite element links were ensured by defining bond equations. This enabled elimination of translational rigidity coefficients of high values from the system of motion equations. Consequently, computation with a longer step of integration was possible. The motion equations of the system were supplemented with the parts connected with the impact of marine environment. The model was validated by comparison of the statics computations with analytical solution (in catenary case), and dynamics computations were compared with the solution developed using MES commercial software package. In both cases, satisfactory conformity of results was obtained already when flexible member was divided into n ≥ 10 elements. High numerical efficiency of the developed model and computer software enabled their application to solve the problem of dynamic optimisation. It is expected that the developed model and computer software shall be used in simulation of dynamic positioning of risers.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.