Path Generation for High-Performance Motion of ROVs Based on a Reference Model

This paper deals with the generation of sufficiently smooth position, velocity, and acceleration references for guiding the motion of an ROV along purposefully defined curvature-continuous paths in automated missions. The references are meant to be employed in high-performance trajectory tracking and dynamic positioning applications. The path planning problem is not in the scope of this work. A reference model that synthesises references concerning a single Degree-of-Freedom (DoF) motion is initially described. Then, the use of the synthesised references as the parametrisation for other references concerning multiple DoF motion along curvature-continuous paths is exploited. Results from computer simulations and full-scale sea trials, both based on the NTNU’s ROV Minerva, are presented and discussed.


Introduction
This paper deals with the generation of sufficiently smooth position, velocity, and acceleration references for guiding the motion of Remotely Operated Vehicles (ROVs) along purposefully defined curvature-continuous (G 2 ) paths in automated missions.The references are meant to be employed in high-performance Motion Control Systems (MCSs) with trajectory tracking and Dynamic Positioning (DP) capabilities, e.g.Fernandes et al. (2015), Omerdic et al. (2012), Sørensen et al. (2012), Caccia (2006), and Hsu et al. (2000).Although this work focuses on (typically) fully-actuated Unmanned Underwater Vehicles (UUVs) such as ROVs, some of the ideas presented herein can be adapted to be used as a means of also guiding the motion of (typically) underactuated marine crafts, e.g.ships and UUVs such as Autonomous Underwater Vehicles (AUVs), in path following motion control applications.
ROVs are indispensable workhorses used worldwide for industrial, research, and military activities, e.g.inspection, intervention, mapping, and survey.They are teleoperated from support surface vessels through umbilical cables which transmit power, commands, and data.Detailed information can be found in e.g.Christ andWernli (2014, 2007).Granted that accurate motion control is desirable regardless of the type of automated mission that is performed, an MCS with trajectory tracking and DP capabilities has to incorporate a guidance subsystem capable of generating suitable and sufficiently smooth references, given that only the references possessing such attributes can be exactly tracked (Sørensen, 2013;Fossen, 2011;Slotine and Li, 2005).
Guidance is concerned with the transient motion behaviour associated with the achievement of the motion control objectives (Fossen, 2011;Breivik and Fossen, 2009), so that the mission specifications and the vehicle dynamics are all simultaneously observed.Furthermore, collisions with stationary obstacles are avoided whenever a collision-free path is closely tracked or followed.The path planning problem is not in the scope of this work.The reader is referred to Tsourdos et al. (2011), Kavraki and LaValle (2008), Minguez et al. (2008), andLaValle (2006) when it comes to the path planning problem, where the latter provides a thorough coverage of the subject.The former reference is concerned with cooperative path planning of unmanned aerial vehicles, whereas both middle references address the robotics task of planning collision-free motion.Naeem et al. (2003) reviewed several guidance laws relevant to UUVs, with an emphasis on AUVs, and asserted that 'the main problem in bringing autonomy to any vehicle lies in the design of a suitable guidance law'.Among other conclusions, Naeem et al. (2003) stated that, 'in practice, Line-of-Sight (LoS) guidance is the key element of all guidance systems', given that the closed-loop path following scheme suits best the needs when it comes to guiding underactuated vehicles.That work made evident the fact that the research on guidance has historically been focusing mainly on underactuated vehicles such as missiles and aircrafts, whose complex guidance problems have been dealt with since the World War II (Breivik and Fossen, 2009).Breivik and Fossen (2009) carried out another survey of the subject keeping the emphasis on AUVs, but taking planar and spatial motion into account.An example of early MCS for UUVs based on the LoS guidance is the work by Healey and Lienard (1993).Improvements aiming at compensating for heading disturbances caused by the sea current can be found in e.g.Aguiar and Pascoal (1997).Fossen and Pettersen (2014) proved that the equilibrium point of the proportional LoS guidance law by Healey and Lienard (1993) is Uni-formly Semi-Globally Exponentially Stable (USGES).The work by Caharija et al. (2014) aims at merging intuitive and theoretical perspectives of the integral LoS guidance for current compensation problems of underactuated ships.The technical challenges underactuated vehicles impose on the reference generation, due to their non-holonomic kinematic constraints (Aicardi et al., 2000;Egeland et al., 1996), justify the attention they have received.The LoS and LoS-based guidance laws are still often employed to guide ROVs, e.g.Omerdic et al. (2012), Sørensen et al. (2012), and Caccia (2006), due to their simplicity and ease of implementation (Breivik and Fossen, 2009;Naeem et al., 2003).An example of LoS-based conjoint guidance and control scheme that generates reference heading to steer an ROV towards the destination point employing a Lyapunov-based algorithm to ensure smoothness and convergence is given by Caccia et al. (1998).The approach was refined in Caccia and Veruggio (2000).A similar MCS by Guo et al. (2003), conceived to control the motion of AUVs, employs a sliding mode fuzzy algorithm in place of the Lyapunov-based algorithm.Dukan (2014) proposed a spatial LoS guidance strategy dedicated to guide fully actuated ROVs.The interested reader is referred to Caharija (2014) and Lekkas (2014) when it comes to more recent extensions regarding the LoS and LoS-based guidance laws for underactuated marine vehicles.
ROV-based missions neither typically need high autonomy levels, nor present technical challenges with respect to non-holonomic constraints.On the other hand, repetitive missions, and missions which require accurate motion control, claim some degree of automation.Furthermore, ROV pilots may feel exhausted and less attentive during long-lasting missions, and these factors may lead to riskier and more expensive missions.The reader is referred to Vasilijević et al. (2013), Ho et al. (2011), andHsu et al. (2000) for discussions about the challenges entailed in ROV piloting.Experienced pilots have affirmed that every contribution towards automating ROV-based missions is effective not only in increasing the overall motion control accuracy level, but also in improving the global mission performance (Hsu et al., 2000).Fully-actuated ROVs permit high-performance motion control to be exercised through trajectory tracking and DP.A collection of guidance techniques, ranging from open-loop Filter-Based Reference Models (FBRMs) to closedloop optimisation-based reference generators, suitable for guiding the motion of fully-actuated UUVs can be found in Fossen (2011) and in the references therein.The simplest, yet effective, technique consists of an FBRM built upon a 2nd-or 3rd-order Low-Pass Filter (LPF) transfer function.It can be easily implemented and modified, and runs fast, even in a digital computer with modest hardware.Its drawback lies in its linearity, because, for any fixed tuning, the generated velocities and accelerations are linearly proportional to the distances to be covered (Fernandes et al., 2012;Fossen, 2011).Fossen (2011) presents alternatives to partially remedy the problem.Linear FBRMs stem from the Model Reference Adaptive Control (MRAC) technique (Landau, 1974), and are commonly employed in trajectory tracking control systems to improve the closed-loop transient response ( Åström and Hägglund, 2011;Slotine and Li, 2005).
This paper makes two contributions to the literature on guidance: • A Reference Model (RM) that generates references for steady single Degree-of-Freedom (DoF) motion and is easy, i.e. meaningful, to tune; • A path generation (reference generation) scheme, built upon the RM, that generates references for steady multiple DoF motion.
The motivation for developing the RM was to develop another alternative to the useful, yet simple, FBRM found in e.g.Sørensen (2013) and Fossen (2011).The inspiration that underlies the development came from works on motion optimisation for AUVs, e.g.Chyba et al. (2008) and Kumar et al. (2005), although this work does not necessarily seek to provide optimal references, and from works on reference generation for highly accurate Computer Numeric Controlled (CNC) machines, e.g.Huo and Poo (2012) and Matsubara et al. (2011).It is important to emphasise that the references are generated in open-loop through an FBRM in Matsubara et al. (2011), whereas they are directly synthesised by computers in Huo and Poo (2012).The proposed RM synthesises references for guiding a single DoF motion, either linear or angular, in a suboptimal manner with regard to time.The maximum -or minimum, depending on the motion direction -velocity is kept for the longest time span possible, as exemplified in Figure 1.Such velocity regime is intended to: i) induce steadier hydrodynamic effects; and ii) demand steadier thrust forces and moments from the ROV's propulsion system.Both factors together favour the attainment of small reference tracking errors, and also provide more favourable conditions for adaptive controllers and/or observers, if any, to faster and more robustly estimate and cope with the unknown current-generated perturbations.The RM generates references via direct computer synthesis, similarly to Huo and Poo (2012).Thus, repeatability, finite convergence time, and facilitated interaction with the references while they are being synthesised, e.g.starting, pausing, resuming, and aborting, can be achieved.The shapes of the synthesised references keep much resemblance with those considered better in Kumar et al. (2005).The present work is an extension of that reported in Fernandes et al. (2012).The computer code of the RM has been improved.It is currently more efficient and needs less memory space to be stored than before.
The path generation scheme utilises the references synthesised by the RM as the parametrisation of a group of references which are generated for guiding the four DoF -surge, sway, heave, and yaw -motion of an ROV along pre-defined curvature-continuous paths.
The rest of the paper is organised as follows.Section 2 deals with the RM.Section 3 deals with the path generation scheme.Section 4 presents selected results from simulations and full-scale sea trials, both based on the NTNU's ROV Minerva, which is introduced in Appendix B. Section 5 presents concluding remarks.Appendix A presents the optimal curve shapes (single DoF) along a rectilinear path under the condition of constrained acceleration and velocity.

Description of the basic version
The RM synthesises position, velocity, and acceleration references of classes C 2 , C 1 , and C 0 , respectively, for guiding a single DoF motion in a suboptimal manner with regard to time.The differentiability class C n , where n ∈ N, denotes a function whose n-th derivative with respect to time exists and is continuous on the domain of definition of the function.The condition of suboptimality has harmless implications in practice.It is mostly a consequence of the fact that the RM has to reflect the constraints associated with the dynamics of the guided ROV, e.g.limited acceleration and velocity.This point becomes clearer when the shapes of the references synthesised by the RM, which are exemplified in Figure 1, are compared with the optimal curve shapes shown and discussed in Appendix A. The references are parameterised by the parameter time t ∈ R 0 [s].They are synthesised through a combination of functions distributed into four consecutive phases, as seen in Figure 1, where each function is active only during a certain amount of time, such that a(t) := a 1 (t) + a 2 (t) + a 3 (t) + a 4 (t), (1)  1.

First phase: the acceleration phase
This phase is characterised by the modulus of the velocity reference increasing from zero to virtually the maximum value, see Figure 1.The phase is split into two subphases in which the modulus of the acceleration reference firstly increases from zero to virtually the maximum absolute acceleration, and thence decreases to virtually zero again.During this phase, the references are defined as and where H(•) is the unit step function, and is the instant at which this phase ends, when the acceleration reference reaches the threshold a 1 (t − 2 ) := θ 0 a m , beyond which the acceleration reference virtually vanishes, where t − 2 denotes t ր t 2 , i.e. parameter t tending to t 2 from the left, and θ 0 ∈ R >0 | θ 0 ≪ 1 is another tuning parameter (see Table 1), which is discussed within the next subsubsection.The expression of t 2 is furnished in Subsubsection 2.1.5.

Second phase: the constant velocity phase
This phase is characterised by the (constant) cruise velocity v m , see Figure 1.During this phase, the references are defined as and p 2 (t) : where . The expression of P 2 is furnished in Subsubsection 2.1.5.The function h 2 (t), which is the characteristic function of the half-closed finite interval [t 2 , t 3 ), is defined as Remark: Every reference undergoes a step discontinuity of negligibly small magnitude at the transition from the 1st phase to the 2nd phase due to the presence of the exponential function in (8).The lesser θ 0 , the lesser the magnitudes of the step discontinuities.On the other hand, the lesser θ 0 , the longer the 1st phase lasts, hence causing the cruise velocity v m to be reduced.For instance, if θ 0 = exp(−7) and 12 × 10 −4 m/s (< 0.1 %).Likewise, if τ 12 = 1 s, then p(t 2 ) undergoes a step of magnitude (1 − θ a ) θ 0 τ 12 < 9.12 × 10 −4 m.If θ 0 = exp(−11.11)≈ 1.5 × 10 −5 instead, whereas all other values are kept the same as before, the step discontinuities would virtually disappear in the face of quantisation ( Åström and Wittenmark, 1997), if the RM was implemented in a 16-bit digital system, given that 2 −16 ≈ 1.53 × 10 −5 yields the resolution of ≈ 0.0015 %.Therefore, in this work, it is defined that exp(−13) θ 0 exp(−7).

Third phase: the deceleration phase
This phase is characterised by the modulus of the velocity reference decreasing from the maximum value, i.e. |v m |, to virtually zero, see Figure 1.The phase is split into two subphases in which the modulus of the acceleration reference firstly increases from zero to virtually the maximum absolute deceleration, and thence decreases to virtually zero again.During this phase, the references are defined as and is the maximum -or minimum, depending on the motion direction -deceleration, τ 31 , τ 32 ∈ R >0 [s] are time constants respectively associated with the 1st and 2nd subphases, is the instant at which the 2nd subphase begins, when the velocity reference reaches the thresh- is the tuning parameter (see Table 1) that defines the fraction of v m at which the velocity reference starts getting bent as it proceeds towards zero.Lastly, , where t − 3 denotes t ր t 3 , i.e. parameter t tending to t 3 from the left, and , where t − 4 denotes t ր t 4 , i.e. parameter t tending to t 4 from the left.The expressions of d m , τ 31 , τ 32 , t 3 , t 4 , P 3 , and P 4 are given in Subsubsection 2.1.5.The auxiliary functions f 31 (t) and f 32 (t) are defined as

Fourth phase: the constant position phase
This phase is characterised by the constant position reference, see Figure 1.During this phase, the references are defined as and where L ∈ R >0 [m] is the tuning parameter (see Table 1) that defines the (straight line) path length.The function h 4 (t), which is the characteristic function of the closed semi-infinite interval [t 5 , ∞), is defined as Remark: Every reference undergoes another step discontinuity of negligibly small magnitude at the transition from the 3rd phase to the 4th phase due to the presence of the exponential function in (19).

Pre-computation
Let the time ratio r T ∈ R >0 be defined as where is the tuning parameter (see Table 1) that defines the desired minimum time to reach which is the tuning parameter (see Table 1) that defines the desired cruise velocity, and is the tuning parameter (see Table 1) that defines the desired minimum time to stop moving from V d .These parameters translate the desired maximum -or minimum, depending on the motion directionacceleration and deceleration through the direct relations 'V d /T a ' and '−V d /T d ', respectively.Let the auxiliary constants κ a ,κ d ∈ R >0 be defined as and where ξ a , ξ d ∈ R >0 are constants.The greater ξ a and ξ d , the steeper the slopes of the references during the 1st subphases of the 1st and 3rd phases, respectively, since these constants directly influence the time constants defined in ( 35) and ( 37) ahead.Even though ξ a and ξ d are not primarily designed to be tuning parameters, it can be useful to be able to change their values.Hence, it is defined that 10 ξ a , ξ d 15 in this work.
Then, the candidate absolute value v c ∈ R >0 [m/s] for the cruise velocity is defined as where ǫ L ∈ R >0 | ǫ L < 0.1 is the tuning parameter (see Table 1) that defines the minimum fraction of the path length L that is to be traversed at the cruise velocity v m , which is defined as The adjusted minimum time to stop moving from v m , i.e. t d ∈ R >0 [s], and the adjusted minimum time to reach v m , i.e. t a ∈ R >0 [s], are defined as and The maximum (minimum) acceleration and deceleration are defined as and The time constants τ 11 through τ 32 are defined as and The auxiliary time instants T 2 := −τ 12 ln(θ 0 ), ( 40) and such that the function switching time instants t i , i ∈ {1, 2, 3, 4, 5}, are consecutively defined, with respect to the initial instant t 0 , as Finally, the positions P i , i ∈ {1, 2, 3, 4}, are respectively defined as and where the relation T 1 /τ 11 = ξ a , which is used in (45), stems from (39), the relation T 2 /τ 12 = − ln(θ 0 ), which is used in (46), stems from (40), and the relation T 4 /τ 31 = ξ d , which is used in (48), stems from (42).

Results from full-scale sea trials
See Subsection 4.1 for the results.

Description of an advanced version
The advanced version builds upon the basic version described in Subsection 2.1 and, thus, it bears strong resemblance to that version.The synthesised references are likewise structured as in ( 1)-( 3).Each of the four phases is separately described in Subsubsections 2.2.1-2.2.4.The required pre-computations are introduced in Subsubsection 2.2.5.Table 1 collects and summarises all tuning parameters.The main differences in comparison with the basic version are the following: • the initial velocity and the final velocity can be both different from zero; • the first and third phases are skipped depending on the initial and final velocities; the second phase is always executed; and the fourth phase is executed whenever the final velocity is null; • the cruise velocity is numerically determined.

First phase: the acceleration phase
This phase is characterised by the modulus of the velocity reference increasing from the modulus of the initial value, i.e.
s], which is the tuning parameter (see Table 1) that defines the desired initial velocity, to virtually the maximum value, i.e. |v m | = |v(t 2 )|, where v m ∈ R [m/s] is the cruise velocity.This phase is skipped whenever v i = v m .Likewise in the basic version, this phase is split into two subphases.During this phase, the references are defined as and where a m ∈ R [m/s 2 ] is the maximum -or minimum, depending on the motion direction -acceleration, which is computed as indicated in (33), τ 11 , τ 12 ∈ R >0 [s] are time constants respectively associated with the 1st and 2nd subphases, t 0 ∈ R 0 [s] is the instant at which this phase begins, and is the instant at which the 2nd subphase begins, when the velocity reference reaches the threshold v ) is the tuning parameter (see Table 1) that defines the fraction of (v m − v i ) at which the velocity reference starts getting bent as it proceeds towards v m .Lastly, The expressions of v i , v m , τ 11 , τ 12 , t 1 , P 1 , and t a , the latter being useful to compute a m as in (33), are given in Subsubsection 2.2.5.The functions f 11 (t), f 12 (t), h 11 (t), and h 12 (t) are defined similarly to those in ( 7)-(10).

Second phase: the constant velocity phase
This phase is characterised by the cruise velocity v m .During this phase, the references are defined as and where is the instant at which this phase begins.The expressions of P 2 and t 2 are furnished in Subsubsection 2.2.5.The function h 2 (t) is defined similarly to that in (14).
Notice the similitude between this phase and the second phase of the basic version in Subsubsection 2.1.2.

Third phase: the deceleration phase
This phase is characterised by the modulus of the velocity reference decreasing from the maximum value, i.e. |v m | = |v(t 3 )|, to virtually the modulus of the final value, i.e.
which is the tuning parameter (see Table 1) that defines the desired final velocity.This phase is skipped whenever v f = v m .Likewise in the basic version, this phase is split into two subphases.During this phase, the references are defined as and is the maximum -or minimum, depending on the motion direction -deceleration, which is computed as indicated in (34), τ 31 , τ 32 ∈ R >0 [s] are time constants respectively associated with the 1st and 2nd subphases, is the instant at which this phase begins, and is the instant at which the 2nd subphase begins, when the velocity reference reaches the threshold is the tuning parameter (see Table 1) that defines the fraction of (v m − v f ) at which the velocity reference starts getting bent as it proceeds towards v f .Lastly, 3 ) [m], and . The expressions of d m , v f , τ 31 , τ 32 , t 3 , t 4 , P 3 , P 4 , and t d , the latter being useful to compute d m as indicated in (34), are given in Subsubsection 2.2.5.The functions f 31 (t), f 32 (t), h 31 (t), and h 32 (t) are defined similarly to those in ( 18)-( 21).

Fourth phase: the constant position phase
This phase is characterised by the constant position reference.It is executed whenever v f = 0 m/s.During this phase, the references are defined as and where the function h 4 (t) is defined similarly to that in (25), and the expression of t 5 , which appears in (25), is furnished in Subsubsection 2.2.5.Notice the similitude between this phase and the second phase of the basic version in Subsubsection 2.1.4.

Pre-computation
Let the adjusted minimum time to reach v m from v i , which is t a ∈ R >0 [s], and the adjusted minimum time to reach v f from v m , which is t d ∈ R >0 [s], be given by and and also the maximum (minimum) acceleration and deceleration be defined as and The time constants τ 11 through τ 32 are defined as and The auxiliary time instants , where T i := t − i , are defined similarly to those in (39)-( 40) and ( 42)-( 43).Notice that T 3 has purposefully been skipped.It is separately treated ahead.
The auxiliary lengths L 11 , L 12 , L 2 , L 31 , L 32 ∈ R 0 [m] are defined as (notice the modulus function) and Finally, any robust, fast-execution numerical method can be employed to determine v m through the recursive use of ( 61)-( 73), including the computation of T 1 , T 2 , T 4 , and T 5 , such that the three following conditions Then, the positions P i , i ∈ {1, 2, 3, 4}, are respectively defined as and Lastly, T 3 is defined as and the function switching time instants t i , i ∈ {1, 2, 3, 4, 5}, are consecutively defined, with respect to the initial instant t 0 , in a similar fashion to those in (44).

Simulation results
See Subsection 4.2 for the results.

Further features
Both the basic and the advanced versions can be coded together, such that they can be somewhat interchangeable.Then, the use of the tuning parameters will dictate which version is the best in every case.So far, both versions of the RM work in open-loop.However, the use of reference and state -either directly measured or estimated -feedback opens up new possibilities.For instance, the following can be listed: • Hysteretic waiting function: by monitoring how closely an ROV tracks the synthesised position reference, it is possible to 'wait' for the ROV, for how long it is needed, in case the ROV is lagging behind such reference, until the difference is reduced to an acceptable value.Since the references are directly synthesised, it is possible to hold the value of the position reference without affecting the velocity reference.Such a measure has to be applied preferably during the second phase, because the acceleration reference is null during this case; • Emergency stop function: if a collision against a stationary object, or against the seabed, is in the imminence of happening, for instance, it is a good idea to bring the ROV to a halt in the shortest time possible.Towards this end, it is possible to also code a slightly modified version of the third phase of the basic version, which is to be used in exceptional circumstances, such that both the acceleration and velocity references are quickly, although not instantly, zeroed, and the ROV at risk of collision is quickly kept safe in DP mode.

Path generation 3.1. Introduction
Motion control applications concerning UUVs typically require the use of two reference frames simultaneously (Sørensen, 2013;Fossen, 2011;SNAME, 1950).The locally inertial North-East-Down (NED) reference frame, see Figure 2, serves as the reference for the following vector of desired position and attitude (heading angle)  (t)] and ṙd (t) are respectively related to the motion in surge, sway, heave, and yaw.The origin of the BF frame is fixed at a convenient point of the UUV.The transformation matrix J(ψ d (t)) ∈ SO(4) (Special Orthogonal group of order 4), which is defined as and Algorithm 1: Desired heading accumulation The desired heading angle ψ d (t) has to be accumulated along consecutive turns, e.g. through the use of Algorithm 1, to be in (80), for the sake of continuity beyond the natural supremum -either 2 π or π [rad], depending on the considered range -, and the natural infimum -either 0 or −π [rad], depending on the considered range, i.e. to avoid step discontinuities at the aforementioned infimum and supremum.
The accumulated heading ψ d (q) that comes out from Algorithm 1, q ∈ N \ {0}, when it comes to a discretetime implementation, is based on the two most recently computed values of the heading, ψ c (q) and ψ c (q − 1), which have the same bounded range, and on the variable counter ∈ Z, which accumulates the number of turns taken either Clockwise (CW) (counter 0), or Anticlockwise (ACW) (counter < 0).The variable is to be initialised as counter(1) = counter(0) = 0.
There are different ways and purposes of defining paths.This work is not concerned with the path planning problem, as mentioned before.It is therefore assumed herein that purposeful planar and spatial curvature-and torsion-continuous paths composed of rectilinear and curvilinear parts are previously defined, either through the use of proper path planning techniques, or directly based on the target applications.Curvature and torsion continuity imply that two connected parts share common centres of curvature and torsion at the join point.To illustrate the point of this section, paths which are comprised of straight lines connected through planar clothoids are used, without loss of generality.In fact, any planar curved line that is capable of satisfying the curvature continuity condition can be used instead.For instance, the Fermat's spiral (Lekkas, 2014;Lekkas et al., 2013), the single polar polynomial (Lai et al., 2007;Nelson, 1989), the cubic Bézier curve (Farouki, 2012;Seidel, 1993), and some Pythagorean hodographs (Farouki, 2008(Farouki, , 1997;;Farouki and Sakkalis, 1990), among others, are suitable alternatives.

Straight lines
Straight lines are fundamental building blocks of curvature-and torsion-continuous paths formed by rectilinear and curvilinear parts.
The length L l ∈ R >0 [m] of a straight line is given by where the superscripts ' 0 ' and ' 1 ' are not exponents in this case, but help to identify the start point (n 0 l , e 0 l , d 0 l ) and the endpoint (n 1 l , e 1 l , d 1 l ), whose coordinates in the NED frame are n The elements of the points (n l (̟ l (t)), e l (̟ l (t)), d l (̟ l (t))) [m], which describe the straight line in the NED frame, are given by and where , which is a constant angle contained in the NE-plane that is measured from the N-axis, and , which is a constant angle contained in a vertical plane that is measured from the NE-plane, are given by and The corresponding velocities ṅl ( ̟l (t)), ėl ( ̟l (t)), and ḋl ( ̟l (t)) [m/s], and accelerations nl ( ̟l (t)), ël ( ̟l (t)), and dl ( ̟l (t)) [m/s 2 ], are given by ṅl ( ̟l (t)) = ̟l (t) L l cos(α v l ) cos(α h l ), ( 90) and dl ( ̟l (t)) = ̟l (t) L l sin(α v l ). (95)

Mirror-symmetric twin clothoids
The clothoid, also known as the Euler spiral, or the Cornu spiral, is a plane curve whose curvature changes linearly with the arc length (Levien, 2008).Harary and Tal (2012) extended it to three dimensions, such that the torsion also changes linearly with the arc length.These characteristics can be explored for constructing curvature-and torsion-continuous paths, which are formed by straight lines connected through pairs of mirror-symmetric twin clothoids.Every pair of twin clothoids presents symmetry with respect to the line bisecting the (total) angle span.The (total) arc length can be exactly and easily computed based on the angle span.By knowing the arc length, the angle span, the rotation direction, and the parity, it is possible to determine the endpoint of the pair of twin clothoids, through the accurate computation of the midpoint, based on the start point.Only planar clothoids are considered here.
The arc length L c ∈ R >0 [m] of a pair of mirror-symmetric twin clothoids that winds from the start point, at which the curvature is null, to the midpoint, at which the curvature is maximum, and thence unwinds to the endpoint, at which the curvature is equal to zero again, is given by where    99) where the superscripts ' 0 ', ' 1 ', and ' 2 ' are not exponents, is the tangent angle at the point (n 0 c , e 0 c ), and H(•) is the unit step function.The parameter ̟ c (t) ranges from zero to 0.5 along the first clothoid arc of the pair of twin clothoids, as indicated in ( 97)-( 99).The terms n 1 c (̟ c (t)) [m] and e 1 c (̟ c (t)) [m] are given by , (100) where c ′ c (̟ c (t)) and s ′ c (̟ c (t)) are given by the following absolutely convergent power series and denotes the rotation direction developed from the point (n 0 c , e 0 c ), i.e. λ c = 1 assigns the CW direction, whereas λ c = −1 assigns the ACW direction.
The tangent angles α 1 c (̟ c (t)) [rad] are given by The parameter ̟ c (t) ranges from 0.5 to 1 along the second clothoid arc of the pair of twin clothoids, as indicated in ( 97)-( 99).The terms denotes the parity, i.e. ρ c = 1 assigns the even parity, which results in a symmetric pair of twin clothoids, whereas ρ c = −1 assigns the odd parity, which results in an antisymmetric pair of twin clothoids.The latter pair presents a curvature discontinuity at the midpoint, due to an instantaneous change of the centre of curvature at the midpoint, which causes such pair of clothoids to degenerate into just tangentcontinuous (G 1 ).The argument χ c is defined as and c ′′′ c (̟ c (t)) and s ′′′ c (̟ c (t)) are given by and (109) The linear and angular velocities ṅc (̟ c (t), ̟c (t)) [m/s], ėc (̟ c (t), ̟c (t)) [m/s], and αc (̟ c (t), ̟c (t)) [rad/s] along the pair of twin clothoids are given by ṅc and the corresponding linear and angular accelerations nc where the derivatives with respect to time in ( 110)-( 111) and ( 113)-( 114) are obtained with the help of where , and s ′′′ c (̟ c (t)), in ( 101), ( 102), (107), and (108), ζ is a dummy constant that represents either 1/κ max c or λ c /κ max c , and INum(k, ̟ c (t)) and IDen(k) respectively stand for the numerator and the denominator of φ c (̟ c (t)).
The magnitude of the terms in the sums to infinity (101), ( 102), ( 107), (108), and (116) decreases rapidly as k increases.The use of only a few terms, e.g.k ∈ {0, 1, . . ., 12}, yields very close approximations for the sums in practice.Exploring this fact leads to expedited, yet accurate, computations of those summands.

Path generation
There are essentially two possibilities of combining straight lines and pairs of twin clothoids -or other suitable types of curved lines -to form paths.The first possibility consists of using the pairs of twin clothoids as approximating curves with specified maximum curvatures, thereby resulting in paths which do not pass through all the waypoints which sketch out the paths.There are typically as many pairs of twin clothoids as the total number of waypoints minus two. Figure 4 shows an example in which all clothoids have the same maximum curvature.The second possibility consists of using the pairs of twin clothoids as interpolating curves with specified maximum curvatures, thereby resulting in paths which do pass through all the waypoints which sketch out the paths.Figure 5 shows an example in which all clothoids have the same maximum curvature.Such a kind of path, which encompasses two successive lane changes, is useful, among other things, to avoid collisions.The path in Figure 5 is based on antisymmetric pairs of twin clothoids -odd parity (ρ c = −1) -, which render the path only tangent-continuous (G 1 ), due to the inherent properties of the clothoids.Notice that this fact has nothing to do with the curve interpolation.Alternatives to circumvent the curvature discontinuity problem are: i) to replace each antisymmetric pair of twin clothoids by two symmetric pairs of twin clothoids linked by a short straight line; or ii) to replaced them with curved lines which yield curvature-continuous lane change paths, e.g. the curve reported in Nelson (1989); or iii) to move along the path in such a manner that the velocity and the acceleration are both null at the join point of both clothoids of the antisymmetric pair.
The total path length L p ∈ R >0 [m] is given by where L i ∈ R >0 [m] denotes the length of the i-th part of the path, which is formed by N juxtaposed parts in total.In essence, the path length can be mapped, from the start point to the endpoint, as a strictly monotonically increasing one-to-one function of the parameter Likewise, the length of the i-th part of the path can be mapped as a function of the parameter and ̟l (t) = ̟i (t), or ̟ c (t) = ̟ i (t), ̟c (t) = ̟i (t), and ̟c (t) = ̟i (t), i ∈ {1, 2, . . ., N }, depending on the type of the line -rectilinear or curvilinear, respectively -, see Subsections 3.2-3.3.The parameter ̟ p (t), and its derivatives with respect to time ̟p (t) and ̟p (t), are defined as and where a(t), v(t), and p(t) come respectively from ( 1)-(3).They can be generated by either the basic version of the RM that is described in Subsection 2.1, or an advanced version of it, as described in Subsection 2.2.The relation is utilised to determine the derivatives with respect to time in ( 120)-( 121).Consequently, the N parameters ̟ i (t), i ∈ {1, 2, . . ., N }, and their derivatives with respect to time ̟i (t) and ̟i (t), are defined as where L j = L i ⇔ j = i.The characteristic functions h i (̟ p (t)) of the first N − 1 (half-closed) subintervals of ̟ p (t), i.e. for i ∈ {1, 2, . . ., N − 1}, are defined as where H(•) denotes the unit step function, whereas the characteristic function of the N -th (closed) subinterval of ̟ p (t), i.e. for i = N , is defined as

Simulation results
See Subsection 4.3 for the results.Figure 6 depicts a 50 m-long surge motion along a straight line heading northwards from the local origin of the NED reference frame.This implies that ), 0, 0, 0 ] T , and νd = [ a(t), 0, 0, 0 ] T , in ( 80)-( 83).The maximum absolute (spatial) position error is < 0.5 m, with maximum absolute depth error < 0.1 m.The maximum absolute heading error is < 4 • .Figure 7 depicts the corresponding velocities along the path.

Selected results
The measured and the estimated velocities remained close to the references all the time.Figure 8 depicts the commanded rotations of the propellers.The curves of the starboard and port thrusters are nearly flat, somewhat mimicking the shape of the reference surge velocity u d (t).The curve of the lateral thruster oscillates gently in order to maintain the desired heading and position, compensating for the motion disturbances.The curve of both vertical thrusters is practically flat all the time.Fernandes et al. (2012) also report a simple performance comparison experiment concerning a single surge motion being guided by either the proposed RM and an FBRM (Fossen, 2011), where the proposed RM outperforms the FBRM, as expected.

Reference model: advanced version: simulation results
This subsection presents simulation results regarding the advanced version of the RM described in Subsection 2.2.The RM is studied in isolation, for the sake of clarity.References regarding a single DoF motion were synthesised.The simulation was based on recursive use of the RM thrice.Different features of the RM were explored every time it was run.The constant sampling frequency was 6. 6 Hz.The tuning parameters are collected in Table 2. See Table 1 for details about them.
The results are concentrated in Figure 9.The 250 mlong motion starts and ends at rest.It is worthwhile to realise that v m = V d , v i = V i , and v f = V f in all cases.The vertical dashed red lines separate each case, which begin at 15 s, ≈ 78 s, and ≈ 128 s, respectively.The latter ends at ≈ 224 s.In the first case, the fourth phase is skipped, since the desired final velocity is not zero.In the second case, the first, third, and fourth phases are skipped, because the cruise velocity is held throughout.Finally, in the third case, all phases are executed.In particular, the fourth phase is executed due to the fact that the desired final velocity is equal to zero.

Path generation: simulation results
This subsection presents simulation results based on the ROV Minerva, see Appendix B, regarding the path generation scheme described in Section 3. The four DoF MCS used to carry out the simulation was that reported in Fernandes et al. (2015), with the guidance subsystem built upon the just mentioned path generation scheme.The simulation was based on recursive use of the advanced version of the RM that is described in Subsection 2.2.Notice that the basic version of the RM that is described in Subsection 2.1 could have been used instead for the sake of example.The sampling frequency was 6. 6 Hz.The tuning parameters of the RM are all collected in Table 3. See Table 1 for details on them.The additional tuning parameters concerning the path generation were r min c = 10 m and ρ c = 1.They held for both pairs of twin clothoids.Table 4 gives the waypoints which sketch out the path.The tuning parameters of the MCS, regarding the controller, were (52.6, 83.8, 53.3, 69.2), (5,10,5,20), (420.3, 533.8, 432.5, 182.0), (127) where i, j ∈ {1, 2, 3, 4}, and I ∈ R 5×5 is an identity matrix, whereas the tuning parameters regarding the state observer, were  4. Figure 11 depicts the path curvature.Notice that κ(̟ p (t)) is continuous, and that κ(̟ p The parameter ̟ p (t), and its derivatives with respect to time, are shown in Figure 12.The whole path, whose length is ≈ 184.3 m, has five parts, each of which with a characteristic function h i (̟ p (t)), and parameterised by a parameter ̟ i (t), i ∈ {1, 2, 3, 4, 5}, see Figure 13.    Figure 17 depicts the commanded thrust forces and moment from the propulsion system of the ROV.The curves somewhat mimic the shapes of the reference velocities depicted in Figure 16.When it comes to the sway and heave forces, and the yaw moment, the curves gently oscillate around constant values in order to compensate for disturbances.The heave thrust force kept the desired depth by counteracting the positive buoyancy force of ≈ 100 N of Minerva.

Concluding remarks
The paper dealt with the generation of sufficiently smooth references for guiding the motion of ROVs along purposefully pre-defined curvature-continuous paths.An RM that synthesises references concerning a single DoF motion was initially described.After that, the references synthesised by the RM were used to parameterise other references concerning the motion along curvature-continuous paths comprised of rectilinear and curvilinear parts.Both the proposed RM and the proposed path generation scheme proved to be able to synthesise references which yield ROV motion: i) under less induced plant parameter variations; ii) that favours energy saving along the motion; iii) with high overall accuracy; and iv) with finite convergence time.Moreover, both of them were easily tuned through the use of meaningful tuning parameters.Future works must seek to find suitable numerical methods for determining the cruise velocity in the proposed advanced version of the RM.Other advanced versions can also be devised in the future, in order to serve yet different purposes.Also, spatial paths based on 3D spirals, e.g. as suggested by Harary and Tal (2012), must be treated.

A. Optimal curve shapes
Figure 18 shows the shapes of the optimal acceleration, velocity, and position curves along a straight line path which is to be traversed under the condition of constrained acceleration and velocity.The velocity curve is optimal in the sense that it yields the shortest travelling time T min := t 3 −t 0 [s] between both path's extremities, which are separated by the distance d ∈ R >0 [m].
The path is to be traversed under the additional constraints: i) a(t 0 ) = v(t 0 ) = p(t 0 ) = 0; and ii) a(t 3 ) = v(t 3 ) = 0 and p(t 3 ) = d.Symmetry between the maximum acceleration and deceleration is assumed for the sake of simplicity without loss of generality.The claim in this appendix can be alternatively formulated and proved in the framework of the calculus of variations (van Brunt, 2004;Forsyth, 1960).

B. NTNU's ROV Minerva
Minerva is a SUB-fighter 7500 ROV delivered by Sperre AS in 2003.The NTNU's Research Vessel (R/V) Gunnerus (http://www.ntnu.edu/marine/gunnerus) is the support vessel used to operate Minerva.The ROV is powered from, and communicates with, R/V Gunnerus through a 600 m-long umbilical cable.Minerva has five thrusters with fixed pitch propellers.The starboard and port thrusters are oriented 10 • towards the longitudinal axis.The lateral thruster is the only that has one propeller at each end of its shaft, whereas all the other thrusters have a single propeller each.Figure 19 shows the thruster installation configuration.Table 5 keeps on the specifications of Minerva.Details can be found in Sørensen et al. (2012) and Dukan et al. (2011).
A high-precision hydroacoustic positioning system model HiPAP 500 by Kongsberg Maritime AS determines the north n and east e position components of Minerva relative to the R/V Gunnerus with accuracy better than 0.1 m at update rates 1 Hz ROV(s) = Remotely Operated Vehicle(s) UUV(s) = Unmanned Underwater Vehicle(s)

Figure 1 :
Figure 1: Shapes of the references synthesised by the RM concerning a single DoF motion.
74) are satisfied with the highest |v m |.To begin with, v m = V d , v i = V i , and v f = V f are tested.If the conditions in (74) are not satisfied, then |v m | < |V d | is used.Notice that it may also be needed to adjust v i and v f in order to ensure that |v i | |v m | and |v f | |v m |, in case |v m | drops down to less than |V i | or |V f |.
is the angle span, and κ max c = (1/r min c ) ∈ R >0 [m −1 ] is the maximum admissible curvature, where r min c ∈ R >0 [m] is the minimum radius of curvature, which occurs at the midpoint, i.e. at the point where L c /2 and θ c /2.It is worth knowing that the arc length ratio L c /L r , where L r ∈ R >0 [m] is the length of a virtual arc of circumference -shortest arc -whose start point and endpoint coincide with those of the pair of twin clothoids, is closely approximated by the curve depicted in Figure 3. Notice that L c becomes more than 10 % longer than L r only for θ c > 105.4 • .The radius ratio r r /r min c , where r r ∈ R >0 [m] is the radius of the virtual arc of circumference whose length is equal to L r , is closely approximated by the curve depicted in Figure 3. Notice that for θ c < 20 • , κ max c is approximately twice that of the circumference.The elements of the points (n c (̟ c (t)), e c (̟ c (t))) [m], which describe the pair of twin clothoids in the NE-plane, and the tangent angles α c (̟ c (t)) [rad] at the points, are given by

Figure 3 :
Figure 3: Arc length ratio and radius ratio.

Figure 5 :
Figure 4: Example of path using approximating curves.

Figure 9 :
Figure 9: References synthesised by the advanced version of the RM.

Figure 12 :
Figure 12: Parameter ̟ p (t) and its derivatives with respect to time ̟p (t) and ̟p (t).

Figure 13 :Figure 14 :Figure 16 :
Figure 13: Characteristic functions and parametrisation of the five parts of the reference path.
. The depth d (down position component), see Figure 2, is determined based on the measurements provided by the precision, temperature compensated, piezo-resistive underwater pressure sensor model miniIPS 0760001-100 by Valeport Ltd.It has full-scale span 100 bar, accuracy ±10 mbar, and resolution 1 mbar.Its maximum output update rate is 8 Hz.The heading angle ψ and the yaw rate r are determined through the use of a dedicated complementary filter (Mahony et al., 2008) that treats the measurements provided by the Micro-Electro-Mechanical System (MEMS)-based Inertial Measurement Unit (IMU) model MTi-100 by Xsens Technologies B.V. The gyroscopes of the IMU have typical full-scale spans 450 • /s, maximum bias repeatability 0.2 • /s (one year), and typical in-run bias stability 10 • /h.Its accelerometers have typical full-scale spans 50 m/s 2 , maximum bias repeatability 30 mm/s 2 (one year), and typical inrun bias stability 40 µg.The maximum output update rate is 2 kHz with latency < 2 ms.The velocities u, v, and w are simultaneously measured by a 1200 kHz Workhorse Navigator Doppler Velocity Log (DVL) by Teledyne RD Instruments, Inc.It has full-scale spans

Table 2 :
Tuning parameters of the RM