domain: phys chim spi

text:
Adrien Rouviere
email: rouviere.adrien@onera.fr

Lucas Pascal

Fabien Méry

Ehouarn Simon

Serge Gratton

Neural prediction model for transition onset of a boundary-layer in presence of 2D surface defects

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

In a global context where the general trend is to reduce greenhouse gas emissions, it has become necessary to reduce fuel consumption of future aircrafts. One of the solutions lies in the NLF (Natural Laminar Flow) strategy, which consists on ensuring that pressure gradients on a profile are optimized in order to delay the transition. Despite extensive research is carried out about NLF, a laminar boundary layer is hard to obtained due to its high sensitivity to surface imperfections (rivets, junctions, etc.) inherent to the manufacturing processes.

The e 𝑁 method [START_REF] Smith | Transition, Pressure Gradient and Stability Theory[END_REF][START_REF] Van Ingen | A suggested semi-empirical method for the calculation of the boundary layer transition region[END_REF] is one of the most efficient tool to predict correctly the transition location. The 𝑁 factor represents the amplification of the boundary layer instabilities and transition occurs when the 𝑁 factor reaches a critical value 𝑁 tr . A surface irregularity usually has a localized effect on the vicinity of the defect and is responsible of a 𝑁 factor shift. Many studies have looked for empirical correlations linking the type of defect and its geometric characteristics to the Δ𝑁 shift. Wind tunnel experiments were conducted by Crouch et al. [START_REF] Crouch | Modeling gap effects on transition dominated by tollmien-schlichting instability[END_REF] on gaps and by Wang and Gaster [START_REF] Wang | Effect of surface steps on boundary layer transition[END_REF] on BFS. These studies showed that the Δ𝑁 corresponded to 0.12𝑏 * tanh(36 ℎ * 𝑏 * ) for the gaps and 4ℎ * -1.40 for the BFSs, where ℎ * and 𝑏 * denote the non-dimensional step height and step width respectively. These empirical relations have the advantage of being simple to implement and not requiring additional experimental or numerical calculations in addition to those required for a smooth surface. On the other hand, each correlation corresponds only to a particular defect geometry and often applies only in a certain range of parameters.

Traditional numerical methods used to predict the transition, such as Local Stability Theory (LST) [START_REF] Perraud | Laminarturbulent transition prediction in the presence of surface imperfections[END_REF] or Parabolized Stability Equations (PSE) [START_REF] Thomas | Effect of small surface deformations on the stability of tollmien-schlichting disturbances[END_REF], have given satisfactory results for dealing with smooth cases or with surface defects of restricted dimensions. However, the effect of a surface irregularity on the transition is poorly taken into account by these methods because of the assumptions made on the base flow. To overcome these limitations, Worner et al. [START_REF] Worner | Humps/Steps Influence on Stability Characteristics of Two-Dimensional Laminar Boundary Layer[END_REF] and Edelmann and Rist [START_REF] Edelmann | Impact of Forward-Facing Steps on Laminar-Turbulent Transition in Transonic Flows[END_REF] used Direct Numerical Simulations (DNS) to study the impact of humps and Forward-Facing Steps (FFS) respectively on the transition. More recently, Franco Sumariva et al. [START_REF] Franco Sumariva | On the influence of two-dimensional hump roughness on laminar-turbulent transition[END_REF] introduced the Adaptive Harmonic Linearized Navier-Stokes (AHLNS) method which they coupled with PSE upstream and downstream of the region with strong streamwise flow variations to study the effect of humps on the transition. A similar technique is used by Hildebrand et al. [START_REF] Hildebrand | Predicting Boundary-Layer Transition over Backward-Facing Steps via Linear Stability Analysis[END_REF] to study Backward-Facing Steps (BFS), by coupling PSE at the inlet domain with Harmonic Linearized Navier-Stokes (HLNS) in the vicinity of the defect.

Nowadays, the rise of Neural Networks (NN) makes it possible to predict instabilities by taking into account many input parameters in a robust manner. Crouch et al. [START_REF] Crouch | Transition prediction for three-dimensional boundary layers in computational fluid dynamics applications[END_REF] used NN for determining the instability growth rates for calculating the 𝑁 factors to predict transition caused by crossflow and Tollmien-Schlichting (TS) waves instabilities. Lately, Giannopoulos and Aider [START_REF] Giannopoulos | Prediction of the dynamics of a backward-facing step flow using focused time-delay neural networks and particle image velocimetry data-sets[END_REF] predicted the dynamics of a BFS flow using velocity fields as inputs for a NN. More recently, Zafar et al. [START_REF] Zafar | Convolutional neural network for transition modeling based on linear stability theory[END_REF] proposed a transition model based on convolutional neural networks to predict the growth rates of instabilities in two-dimensional incompressible boundary layers. The same authors also developed a transition model based on recurrent neural networks to predict the 𝑁 factor envelope as well as the transition position for different wing profiles [START_REF] Zafar | Recurrent neural network for end-to-end modeling of laminar-turbulent transition[END_REF]. The use of artificial NN could allow, due to their architecture, to take into account more complex relations between the geometrical characteristics of a defect in the evaluation of the Δ𝑁 compared to the previous empirical correlations.

The aim of this paper is to use NN methods taking into account different geometric and aerodynamic parameters of several types of gap-like surface defects to generate new Δ𝑁 models and to make the prediction of the transition to turbulence of a two-dimensional incompressible boundary layer more accurate.

COMPUTATIONAL STRATEGY

Governing equations

The flow considered is a boundary layer developing on a flat plate with surface defect and governed by the incompressible two-dimensional Navier-Stokes equations:

∇• 𝒖 = 0, ( 1a
)
𝜕 𝑡 𝒖 + (𝒖 • ∇) 𝒖 = - 1 𝜌 ∇ 𝑝 + 𝜈Δ𝒖, ( 1b
)
where 𝒖 is the velocity vector, 𝑝 is the pressure, 𝜌 is the fluid density and 𝜈 is the kinematic viscosity. The state vector 𝒒 = (𝒖, 𝑝) is decomposed into a steady base flow 𝑸 = (𝑼, 𝑃) plus an unsteady small perturbations field 𝒒 ′ = (𝒖 ′ , 𝑝 ′ ) in the form:

𝒒(𝒙, 𝑡) = 𝑸(𝒙) + 𝜀𝒒 ′ (𝒙, 𝑡), 𝜀 ≪ 1, (2)
where 𝑥 and 𝑦 are the streamwise and normal components of the position vector 𝒙 respectively, and 𝑡 represents time.

Introducing this decomposition into (1), the steady Navier-Stokes equations governing the base flow are obtained:

∇• 𝑼 = 0, (3a)
(𝑼 • ∇) 𝑼 = - 1 𝜌 ∇ 𝑃 + 𝜈Δ𝑼. (3b)
At the first order, the LNS equations governing the dynamic of the perturbations developing onto the base flow are written in the form:

∇• 𝒖 ′ = 0, (4a)
𝜕 𝑡 𝒖 ′ + (𝑼 • ∇) 𝒖 ′ + 𝒖 ′ • ∇ 𝑼 = - 1 𝜌 ∇ 𝑝 ′ + 𝜈Δ𝒖 ′ . (4b)
The set of equations ( 4) is used to calculate the evolution of a small perturbation 𝒒 ′ in the boundary layer in a linear regime. Perturbations and forcing term are assumed to be harmonic in time as follow:

𝒒 ′ (𝒙, 𝑡) = q(𝒙) e -𝑖 𝜔𝑡 (5)
where 𝜔 = 2𝜋 𝑓 is the real angular frequency of the perturbations and 𝑓 is the disturbance frequency. Introducing this decomposition into (4), the governing equations of the spatial structure of the perturbations ( û, p) are obtained:

∇• û = 0, (6a)
-𝑖𝜔 û + (𝑼 • ∇) û + ( û • ∇) 𝑼 = - 1 𝜌 ∇ p + 𝜈Δ û. (6b)

N factor and 𝚫𝑵 method

Assuming that the TS modes are dependent on both 𝑥 and 𝑦 directions, their amplification can be quantified by an infinity norm 𝐴(𝑥), i.e. the maximum absolute value of the longitudinal velocity along the wall normal coordinate at each position in streamwise direction, as follows:

𝐴(𝑥) = max 𝑦 | û(𝑥, 𝑦)| . ( 7
)
This method is similar to an e 𝑁 method when taking the logarithm of this normalized amplitude:

𝑁 𝑓 (𝑥) = ln 𝐴(𝑥) 𝐴 0 , (8)
where 𝐴 0 is the initial wave amplitude at point 𝑥 0 (the critical point at which the instability starts to amplify). Each 𝑁 𝑓 factor curve is defined for a given non-dimensional reduced frequency F defined as

𝐹 = 2𝜋 𝑓 𝜈 𝑈 2 ∞ × 10 6 . ( 9
)
Since there is no a priori knowledge on which frequency will be responsible for triggering transition, an envelope curve of the maximum 𝑁 𝑓 factors over a large range of frequencies is defined as:

𝑁 (𝑥) = max 𝑓 𝑁 𝑓 (𝑥). (10)
The 𝑁 factor method assumes that the transition occurs at a position 𝑥 tr for which the envelope curve 𝑁 reaches a certain threshold value 𝑁 tr . The Δ𝑁 method enables the extension of the e 𝑁 method to cases including a surface defect: the 𝑁 factor for a smooth case configuration is artificially shifted by an additional amplification caused by the defect with a value of Δ𝑁:

𝑁 = 𝑁 sm + Δ𝑁, (11)
where 𝑁 sm is the 𝑁 factor evaluated for a smooth surface. The transition position is therefore shifted upstream as the threshold value 𝑁 tr is reached earlier.

Flow configuration and numerical method

Equations in Section 2.1 are solved by the Onera's code Pims2d and the numerical methods used are presented here.

Defect configuration

Different geometries of steps and gaps are studied in this work. The generic parameters defining such defects are two heights ℎ 1 and ℎ 2 , the width 𝑏 and the incompressible boundary layer displacement thickness for a flat plate at zero pressure gradient at the defect location 𝛿 1,d . The value of 𝛿 1 at a given position 𝑥 is given by the semianalytical Blasius solution 𝛿 1 (𝑥) = 1.72 √ 𝜈 𝑥 𝑈 ∞ . This set of 4 parameters is shown in Fig. 1 and aims to geometrically represents any type of gap-like surface defect (cavity, BFS, sequence of BFS and FFS, etc.) in order to standardize the existing correlations linking the geometric parameters to the Δ𝑁. In the rest of this study, the geometrical dimensions of the gap are made non-dimensional by 𝛿 1,d and will be denoted hereinafter ℎ * 1 , ℎ * 2 and 𝑏 * , and the aerodynamic parameter defining the position of the defect 𝛿 1,d will be the Reynolds number Re 𝛿 1,d .

𝑈 ∞ ℎ 1 ℎ 2 𝑏 𝛿 1,d
Figure 1: Surface defect parameters.

Base flow

The computational domain extends in the 𝑥-direction from an abscissa corresponding to Re 𝛿 1,in = 350 to an abscissa corresponding to Re 𝛿 1,out = Re 𝛿 1,d + 1000, and in the 𝑦-direction over a height 𝐻 = 30𝛿 1,out . A boundary layer develops and encounters a gap with non-dimensional heights ℎ * 1 and ℎ * 2 and non-dimensional width 𝑏 * at an abscissa corresponding to Re 𝛿 1,d . No-slip boundary condition is imposed at the wall while a free-slip condition with zero normal-stress is prescribed as an outflow condition. A self-similar boundary layer profile with displacement thickness 𝛿 1,in is imposed at the inlet boundary.

A two-dimensional triangulation of the domain is performed with the FreeFem++ finite element library [START_REF] Hecht | New development in freefem++[END_REF] with a Delaunay-Voronoi algorithm. Equations ( 3) are discretized with Taylor-Hood finite elements P 2 for the velocity field and P 1 for the pressure. The non-linear solution of the base flow is solved with a classical Newton method, by gradually decreasing the value of the kinematic viscosity from 𝜈 = 1 m 2 s -1 until the value corresponding to the desired Reynolds number Re 𝛿 1,d is reached.

LNS computation

Once the base flow is calculated, the TS waves are artificially excited by introducing a volume force term in the LNS equations [START_REF] Crouch | Modeling gap effects on transition dominated by tollmien-schlichting instability[END_REF]. The computational domain for the LNS calculations is similar to the base flow domain, but has a reduced height of 𝐻 = 15𝛿 1,out to avoid having a too large number of vertices. The base flow is thus interpolated on the new mesh and equations ( 6) are discretized with Taylor-Hood finite elements P 2 for the velocity field and P 1 for the pressure. LNS calculations are performed both on a smooth surface and on a surface with defect, for 91 non-dimensional reduced frequencies in the range 𝐹 ∈ [25, 160] to generate the 𝑁 factor curve. This frequency zone is chosen because it delimits the locally unstable region of a Blasius boundary layer.

Code validation

To validate Pims2d in presence of surface imperfection, the configuration studied by Hildebrand et al. [START_REF] Hildebrand | Predicting Boundary-Layer Transition over Backward-Facing Steps via Linear Stability Analysis[END_REF] is reproduced. The boundary layer meets a BFS of slope 𝜃 = 75 • and height ℎ * = 0.72 at the abscissa 𝑥 𝑑 = 0.30 m. The flow has the following characteristics: a freestream velocity 𝑈 ∞ = 28 m s -1 , a displacement thickness at the defect location 𝛿 1,d = 6.90 × 10 -4 m and a freestream unit Reynold number Re ∞ = 1.86 × 10 6 m -1 . Figure 2 compares the 𝑁 factor curves for different frequencies obtained by our method to those obtained by Hildebrand et al. [START_REF] Hildebrand | Predicting Boundary-Layer Transition over Backward-Facing Steps via Linear Stability Analysis[END_REF]. The results match perfectly with the HLNS calculations and validate Pims2d in the presence of surface irregularities. ) for a boundary layer in presence of a BFS.

STUDY OF A CRITICAL STEP AND GAP

The study of a critical gap which parameters are detailed in Table 1 is done in this section in order to detail the database generation process for a particular case. Figure 3 shows the pressure distribution at the wall and exhibits a favorable pressure gradient upstream of the gap, followed by a a region with strong adverse pressure gradient and finally a zone with again a favorable pressure gradient approaching towards zero at the infinite downstream of the defect. These pressure variations are responsible for a separation bubble and have an impact on the boundary layer thickness. This one becomes thinner just upstream of the gap when the pressure gradient is negative, and then thickens in the defect. At a certain distance from the gap, the boundary layer recovers the behavior of an unperturbed Blasius boundary layer.

Once the base flow is calculated, equations ( 6) are solved for 91 non-dimensional reduced frequencies in the range 𝐹 ∈ [25, 160] and the TS wave amplification curves at the specified frequencies are obtained. The maximum of these amplifications gives the 𝑁 factor curve shown in Fig. 4. This curve is compared to the 𝑁 sm curve obtained for a smooth case and allows to visualize quantitatively the effect of a groove on the transition.

At the upstream infinity, the amplification of the pertur- bations is equivalent with or without irregularity. Slightly upstream of the defect, the favorable pressure gradient region tends to stabilize the disturbances by thinning the boundary layer. However, just after the gap, the 𝑁 factor increases abruptly to an 𝑁 max value and recovers the behavior of the flat plate configuration at the downstream infinity but has been shifted by a Δ𝑁 far factor. It is usual to consider that the transition to turbulence occurs when the value of the 𝑁 factor reaches the threshold value 𝑁 tr = 9 in flight conditions for TS transition scenario. The case considered here is therefore critical in that it triggers turbulence at the defect due to the value of the Δ𝑁 max , while the smooth case is always below 𝑁 tr and therefore remains laminar. The different Δ𝑁 values for this case are reported in Table 2.

NEURAL NETWORK MODEL

Definition

An artificial neuron is a non-linear function that associates to an input vector x = (x 1 , . . . , x 𝑛 ) an output y, as follow:

y = 𝜎 𝑛 ∑ 𝑖=1 w 𝑖 x 𝑖 + b 𝑖 , (12)
where 𝜎 is an activation function, w = (w 1 , . . . , w 𝑛 ) is the vector of connection weights and b is a bias. The activation function introduces non-linearity, allowing the ). The figure also shows the amplification curves for all calculated frequencies ( ) as well as that of the most amplified frequency (

).

neuron to represent arbitrarily complex functional relations between the input variables and corresponds in this study to a rectified linear unit (ReLU) function defined below:

𝜎(𝑧) = max(0, 𝑧). (13)
A Neural Network (NN) is a structure composed of successive hidden layers between the input layer and the output layer, where the output of a neuron becomes the input of all the units of the next layer. The learning of the NN consists in iteratively adjusting the weights and biases of the network by minimizing a loss function L.

Generation of training data

Inputs of the NN are Re 𝛿 1,d as aerodynamic variable and ℎ * 1 , ℎ * 2 and 𝑏 * as geometric parameters. The range of parameters detailed in Table 3 have been selected to represent both stable and critical cases that can trigger the transition, according to the criteria defined by Beguet et al. [START_REF] Beguet | Modeling of transverse gaps effects on boundarylayer transition[END_REF]. For each sample of input parameters generated, Pims2d provides as output the values Δ𝑁 max and Δ𝑁 far , which will be the outputs to be predicted by the NN.

Database analysis

Figure 5 illustrates the influence of the input geometrical parameters on the Δ𝑁 for each case in the database. Analysis of Fig. 5(a) shows that the Δ𝑁 max is higher the closer the defect configuration is to a BFS, i.e. the higher ℎ * 1 is and the lower ℎ * 2 is. As illustrated in Fig. 5(b), the highest values of Δ𝑁 max correspond to the highest 𝑏 * when Δℎ = |ℎ * 1ℎ * 2 | < 0.5. On the other hand, beyond this threshold, the width does not seem to play a determining role and a correlation seems to exist between the heights difference Δℎ and the value of Δ𝑁 max . Moreover, the Δ𝑁 max seems to reach a limit around 2 for Δℎ < 0.5, i.e. for cavities with same heights, while a geometry closer to BFS (high Δℎ) is more destabilizing. There also appears to be a strong linear relation between the Δ𝑁 max and Δ𝑁 far values, as shown in Fig. 5(c). This could translate in future studies into the need to know only one of the two Δ𝑁 to predict the other.

Neural network predictions and validation 4.4.1 Implemented neural networks

Neural network with different structures regarding the number of hidden layers and the number of neurons in each layer are considered. The structure of these networks is detailed in Table 4. Each training of the network is done on a normalized training dataset representing 80% of the total dataset (i.e. 600 samples) and randomly selected, while the validation is done on the remaining 20% which have never been seen by the network (i.e. 150 samples).

The loss function to be minimized by the network during the training is defined as the Mean Square Error (MSE) between the real values y (𝑖) and the values predicted by the neural network ỹ(𝑖) :

L = 1 𝑛 𝑛 ∑ 𝑖=1 ỹ(𝑖) -y (𝑖) 2 . (14)
In order to evaluate the model, the metric used is the Mean Absolute Error (MAE) 𝜖 MAE described below:

𝜖 MAE = 1 𝑚 𝑚 ∑ 𝑖=1 y (𝑖) -ỹ(𝑖) , (15)
where 𝑚 is the size of the validation dataset. The results obtained for each network are listed in the Table 4.

Re 𝛿 1,d ℎ * 1 ℎ * 2 𝑏 * Δ𝑁 max

Validation of neural network results

Once validated on the validation dataset, the NN model is compared to the experimental results of Crouch et al. [START_REF] Crouch | Modeling gap effects on transition dominated by tollmien-schlichting instability[END_REF]. These authors experimentally studied the effect of gaps with identical depths on the transition, for ranges of non-dimensional heights up to 5 and widths up to 50. According [START_REF] Crouch | Surface step effects on boundary-layer transition dominated by tollmien-schlichting instability[END_REF], the uncertainty on the experimental Δ𝑁 value is assumed to be ±0.15 (±10 mm on the transition position). An empirical Δ𝑁 is proposed from the experiments, which takes into account both the height of the gap and its width. This correlation has the advantage of capturing the limiting behavior between a deep and shallow gap:

Δ𝑁 Cr ≈ 0.12𝑏 * tanh 36 ℎ * 𝑏 * . (16)
The comparison between the experimental results, the results from the correlation ( 16) and the results from the model generated by the neural network B is shown in Fig. 7. Since the authors of the experimental data did not provide all the data necessary to reproduce their study, the Reynolds number Re 𝛿 1,d is chosen as Re 𝛿 1,d = 1700, and the widths chosen are 𝑏 * = 4.3, 10, 14, 23 and 45 depending on the width range considered. A good agreement is observed between the experimental data and the NN prediction as long as we stay in the range of heights and widths studied by the networks. Outside this range, i.e. in the yellow and hatched zones of the graph, corresponding to ℎ * ≥ 3 and 𝑏 * ≥ 20, the results are more balanced. For widths studied by the network, the predictions for heights greater than those in the dataset match the experimental results. On the other hand, when the width of the surface irregularity moves away from the range of values known by the network, this one does not predict correctly the Δ𝑁 max anymore. However, the experimental points corresponding to Δ𝑁 max > 6 probably do not correspond to the transition due to TS waves as explained by the authors, but rather to a bypass transition phenomenon. These points correspond to defect geometries located in the critical zone defined by Beguet et al. [START_REF] Beguet | Modeling of transverse gaps effects on boundarylayer transition[END_REF]. Moreover, as explained by Crouch and Kosorygin [START_REF] Crouch | Surface step effects on boundary-layer transition dominated by tollmien-schlichting instability[END_REF], a difficulty lies in comparing the values of an experimental Δ𝑁 to the values of a numerical Δ𝑁. Indeed, the latter is determined by the global change of the amplification factor in the vicinity of the defect without the knowledge of the actual transition location, while the experimental Δ𝑁 is related to the specific modes responsible for the transition. Δ𝑁 max results from the NN model seems to reproduce the experimental Δ𝑁 whereas one should expect an better agreement with the Δ𝑁 far . Unfortunately, only one streamwise position of the gap on the flat plate is available in these experiments.

A second position would have give insights on which Δ𝑁 to be considered. Recently, Methel et al. [START_REF] Methel | Experimental investigation on the effect of forwardfacing steps and gaps combined with wall suction on boundary layer transition[END_REF] have investigated the effect of gaps and forward facing steps on transition on a flat plate with wall suction at a position corresponding to Re 𝛿 1,d = 1640. Some cases without suction have been extracted to be compared to the NN predictions. The agreement between the Δ𝑁 Exp obtained experimentally and the Δ𝑁 far predicted from the NN is good considering the 𝜖 MAE value of each neural network. This comparison can be found in Table 5. The last two gaps are out of the training range.

In order to verify the network's ability to correctly predict the Δ𝑁 outside of its training area, Wang and Gaster [START_REF] Wang | Effect of surface steps on boundary layer transition[END_REF] experimental results on BFS are replicated. The parameter ℎ * 2 is set to zero, while the parameter 𝑏 * is set to zero in a first step and set to 20 in a second step, to ap- 16) and by the Δ𝑁 max predicted by the neural network B ( ). The yellow area ( ) and shaded areas correspond to range of heights and widths that has not been trained by the network. A case with 𝑏 * = 0 is geometrically close to a BFS, while the 𝑏 * = 20 case allows us to know if a large distance between a BFS and an FFS amounts to simulating a BFS. Both cases are situated outside the learning area of the networks. The results are plotted in Fig. 8 and show reasonable agreement between the experiment and the neural predictions, keeping in mind that the network has not been trained with a pure BFS configuration. It is interesting to note that the neural results with 𝑏 * = 20 seem slightly closer to the experimental results when ℎ * 1 < 1.2, while beyond this limit the neural results with 𝑏 * = 0 seem better. In both cases, the observed differences are small and suggest that a "pure" BFS and a succession of BFS and FFS separated by a distance 𝑏 * = 20 have a similar impact on the boundary layer stability. The NN results corresponding to Re 𝛿 1,d > 2000, located outside the training region for the Reynolds number, are represented by filled symbols to differentiate them and have a larger error compared to the experimental results. It is interesting to note that our results are located between the curves corresponding to the Wang and Gaster [START_REF] Wang | Effect of surface steps on boundary layer transition[END_REF] and Hildebrand et al. [START_REF] Hildebrand | Predicting Boundary-Layer Transition over Backward-Facing Steps via Linear Stability Analysis[END_REF] correlations on BFS, recalled below:

Δ𝑁 Wang = 4ℎ * 1 -1.4, (17a)
Δ𝑁 Hild = 3ℎ * 1 -0.55. (17b)

CONCLUSION

A new method for predicting the transition to turbulence of incompressible two-dimensional boundary layers in the presence of gap-like surface defects has been presented. This method is based on the e 𝑁 method and a NN is used to determine the value of the Δ𝑁 generated by the surface irregularity. The proposed model relies on four distinct parameters defining the geometry of the defect and its aerodynamic environment which are its two heights, its width and the Reynolds number based on the displacement thickness of a Blasius boundary layer at the defect position.

The database was generated with the Pims2d code whose advantage is to be quite flexible in the creation of the studied groove mesh, allowing to easily implement boundary layer stability calculations, while the currently existing methods for the prediction of the transition in presence of defects are rather complex to implement. This code has been validated in the case of a boundary layer around a BFS, showing a good agreement with results from AHLNS calculations.

The average error committed by the model is around 8 × 10 -2 on the Δ𝑁 , which allows a fairly robust prediction of the criticality of the surface irregularity. Contrary to the empirical correlations already existing in the literature, the model generated by the NN generalizes to any type of gap insofar as the characteristics of the defect are within the range of parameters studied by the network. The parametrization of the surface defects enables to extend the model to other defects, for instance for BFS.

Encouraging results have been presented in this sense. This work can be followed up by a comparison between the results obtained by the NN model and experimental results on step and gap defects carried out in Onera's wind tunnels. The database will be progressively extended to other types of surface defects and will be expanded to take into account the effects of compressibility and pressure gradients.

10 (Figure 2 :

 102 Figure 2: 𝑁 factor curves for the frequencies 𝑓 = 300 Hz ( ), 𝑓 = 500 Hz ( ), 𝑓 = 700 Hz ( ) and envelope curve () for a boundary layer in presence of a BFS.

Figure 3 :

 3 Figure 3: Pressure distribution at the wall in the boundary layer.

Figure 4 :

 4 Figure 4: 𝑁 factor envelope curve with defect () and for the same aerodynamic configuration but without defect (). The figure also shows the amplification curves for all calculated frequencies ( ) as well as that of the most amplified frequency ().

Figure 5 :

 5 Figure 5: (a)-(b) Evolution of Δ𝑁 max as a function of geometric parameters and (c) evolution of Δ𝑁 far as a function of Δ𝑁 max .

Figure 6 :

 6 Figure 6: Neural Model methodology.

5 ]Figure 7 :

 57 Figure 7: Comparison of the Δ𝑁 obtained experimentally by Crouch et al. [4] ( ), by the correlation (16) and by the Δ𝑁 max predicted by the neural network B ( ). The yellow area () and shaded areas correspond to range of heights and widths that has not been trained by the network.

Figure 8 :

 8 Figure 8: Comparison between predictions of network B, wind tunnel experiments of Wang and Gaster [16] and empirical correlations for BFS. The filled symbols () correspond to ranges of Re 𝛿 1,d for which the network has not been trained.

Table 1 :

 1 Aerodynamic and geometric parameters.

	Re 𝛿 1,d	ℎ * 1	ℎ * 2	𝑏 *
	1795.52 1.72 0.47 7.54

Table 2 :

 2 𝑁 max and Δ𝑁 values for a critical case.

	𝑁 max Δ𝑁 max Δ𝑁 far
	10.10	5.59	3.63

Table 3 :

 3 Statistical description of the database with 750 observations.

	Parameter	Min	Max
	Re 𝛿 1,d ℎ * 1 ℎ * 2 𝑏 * Δ𝑁 max Δ𝑁 far	901.53 1999.41 0.10 3 0 2.90 0.52 14.98 0.03 10.71 0.01 6.51

Table 4 :

 4 Details of networks architectures and results. Architecture of the network corresponds to the number of neurons in each layer.

	Network	Architecture	Parameters 𝜖 MAE
	A B C	[4 -80 -2] [4 -20 -20 -2] [4 -15 -15 -15 -2]	562 562 587	0.07 0.09 0.08

Table 5 :

 5 Comparison between Methel et al.[START_REF] Methel | Experimental investigation on the effect of forwardfacing steps and gaps combined with wall suction on boundary layer transition[END_REF] experimental results without wall suction on a gap and the Δ𝑁 far predicted by the neural networks ℎ

	1.90 3.80 1.90 12.70 1.90 22.20 1.90 28.60	0.15 0.34 0.96 2.20	-0.06 0.46 1.24 1.74	-0.06 0.47 1.10 1.12	0.17 0.58 1.50 2.28
	proximate as closely as possible the geometry of a BFS.

* 𝑏 * Δ𝑁 Exp Δ𝑁 far A Δ𝑁 far B Δ𝑁 far C

ACKNOWLEDGMENTS

This work has been partly funded within the frame of the Joint Technology Initiative JTI Clean Sky 2, AIRFRAME Integrated Technology Demonstrator platform "AIRFRAME ITD" (contract CSJU-CS2-GAM-AIR-2020-945521) being part of the Horizon 2020 research and Innovation framework programme of the European Commission."

SUPPLEMENTARY MATERIAL

The three networks developed in this paper are available at the following address: https://doi.org/10.5281/zenodo.6074956 [START_REF] Rouviere | Neural Prediction Model for Transition Onset of aBoundary-Layer in Presence of 2D Surface Defects[END_REF]. This repository also contains a python script Main.py and a jupyter script Main.ipynb containing explanations and examples on how to use these NN.


Xiaoyue Zhang
email: zhangxiaoyue912@126.com

Yuanhui Wang
email: wangyuanhui@hrbeu.edu.cn

Ahmed Chemori
email: chemori@lirmm.fr

Structural Reliability Based Energy-Efficient Arctic Position Mooring Control of Moored Offshore Structures Under Ice Loads ⋆

Keywords:

This work develops an energy-efficient robust control for arctic position mooring systems subject to discontinuous ice loads and system uncertainties. To provide a safe and energy-efficient reference position for moored offshore structures, an improved structural reliability-based setpoint chasing algorithm is proposed. An ice-load disturbance observer is constructed to handle the arctic environmental disturbances, relaxes the continuity, differentiability, and boundedness requirements of disturbances. The system uncertainties are addressed by an adaptive fuzzy estimator. With the reference position and disturbance estimation, an improved dynamic surface control law is derived by introducing a novel filter to improve the convergence rate of the filter error, it enhances also the control performance of the arctic position mooring with fast response. Theoretical analysis proves the stability of the resulting closed-loop system and the boundedness of all its signals with the invariant set mechanism and Lyapunov theory. Comparative numerical simulations are carried out on a turret-moored floating production storage and offloading (FPSO) vessel to demonstrate the performance and effectiveness of the proposed control scheme.

Introduction

The arctic area is rich in oil and gas resources, accounting for 13% of the world's unproved reserves, and 30% of the world's total reserves, respectively. Due to the increasing energy demand, it has become an attractive resource development area [START_REF] Gautier | Assessment of undiscovered oil and gas in the arctic[END_REF]. However, many arctic exploration areas are located in deep waters with harsh environmental conditions and heavy ice loads. The offshore structures, relying only on a mooring system or a dynamic positioning system, can no longer meet the needs of safe and economic development of arctic resources [START_REF] Li | Arctic fpso: Technical feasibilities and challenges[END_REF]. The position mooring system exploits the capacity of mooring cables to reduce vessel energy consumption and adopts the dynamic positioning system to improve cable safety [START_REF] Skjetne | A survey on modeling and control of thrusterassisted position mooring systems[END_REF]. It provides an effective solution for arctic exploration from both safety and energy efficiency perspectives [START_REF] Sangsoo | Coupled dynamic analysis of thrusterassisted turret-moored fpso[END_REF].

Arctic position mooring control is often exposed to harsh arctic environments, including wind, waves, currents, and especially sea ice [START_REF] Rahimi | Reliability prediction of offshore oil and gas equipment for use in an arctic environment[END_REF], which is the main load affecting ship maneuvering and safety [START_REF] Wu | Statistical analysis of ice loads on ship hull measured during arctic navigations[END_REF]. Ice load rejection thereby is a key objective for the arctic position mooring. In [START_REF] Zhou | Numerical simulation of moored ship in level ice[END_REF], a numerical model is extended based on [START_REF] Su | A numerical method for the prediction of ship performance in level ice[END_REF] to evaluate the continuous ice load acting on a moored ship, ignoring wind, wave and current disturbances. [START_REF] Kjerstad | Modeling and control for dynamic positioned marine vessels in drifting managed sea ice[END_REF] adopted a configurable high-fidelity numerical model to develop a dynamic positioning framework while considering complex environmental disturbances in the managed ice . As an effective solution for handling disturbances, the observer is applied to handle ice loads. Jayasiri, Nandan, Imtiaz, Spencer, Islam and Ahmed (2017) developed an unscented Kalman filter-based observer to estimate both the vessel state and the ice load for dynamic positioning vessels in ice. [START_REF] Zhou | Heading control for turretmoored vessel in level ice based on kalman filter with thrust allocation[END_REF] proposed an ice simulator interconnecting ship and ice dynamics. However, all the above disturbance rejectiom methods, as well as commonly used observers, including nonlinear passive observers [START_REF] Fossen | Passive nonlinear observer design for ships using lyapunov methods: full-scale experiments with a supply vessel[END_REF], high gain observers [START_REF] Ngongi | A high-gain observer-based pd controller design for dynamic positioning of ships[END_REF], finite time observers [START_REF] Xia | Finite time output feedback control for ship dynamic positioning assisted mooring positioning system with disturbances[END_REF], etc., are all constructed based on the restrictive assumptions of continuity, differentiability and boundedness of disturbances. But ice loads may have a rapid variation and discontinuous interaction phenomenon, such as intermittent ice crushing, especially in low-speed oil exploitation operations [START_REF] Mååttånen | Iso 19906 ice crushing load design extension for narrow structures[END_REF]. Accordingly, these assumptions are not appropriate for arctic applications, as such external disturbances may lead to an unstable closedloop system. Therefore, disturbance rejection considering the discontinuous and unsmooth behavior of ice loads is a key issue to be addressed in arctic position mooring control.

On account of the long-term offshore operations of ships in the harsh arctic environment, safety and energy efficiency are two issues to be solved for arctic position mooring control. Some studies have been proposed to solve the above issues. [START_REF] Kjerstad | A robust dynamic positioning tracking control law mitigating integral windup[END_REF] proposed an improved pseudo-derivative feedback dynamic positioning control law to provide highly reactive system responses, especially for ship-ice interaction applications. To directly compensate environmental disturbances and system uncertainties, an acceleration-based feedforward controller is proposed for arctic dynamic positioning (Kjerstad, Skjetne and Jenssen, 2011). In [START_REF] Jayasiri | An optimum control-based approach for dynamic positioning of vessels[END_REF], a nonlinear model predictive controller is proposed to generate online optimal commands for arctic dynamic positioning considering both the actuator and state constraints. Considering energy consumption, Kang, Lee, Lim, Lee, Jang, Jung and Lee (2021) evaluated the improvement of the positioning performance in ice conditions by heading control. Besides, Lee, Choi, [START_REF] Lee | Tension based heading control strategy of the arctic fpso with dp assisted mooring system[END_REF] proposed to introduce tension into optimal heading control to quickly determine the target heading and thus reduce the ice loads acting on the vessels. Indeed, the aforementioned position mooring systems excessively rely on thrusters, but conservatively exploit the capacity of the mooring system to maneuver the vessel position and heading. Consequently, it does not completely address the energy efficiency issue.

To achieve energy-efficient position mooring control, [START_REF] Wang | Reliabilitybased robust dynamic positioning for a turret-moored floating production storage and offloading vessel with unknown time-varying disturbances and input saturation[END_REF] and [START_REF] Tuo | Robust adaptive dynamic surface control based on structural reliability for a turret-moored floating production storage and offloading vessel[END_REF] introduced structural reliability as a quantitative security index to make full use of the mooring system with a guaranteed safety. Nevertheless, it is insufficient to reduce the control effort to take into account only one critical cable with minimal structural reliability. Considering all mooring cables, a setpoint chasing algorithm is proposed to minimize thruster usage and maximize the utilization of the mooring system [START_REF] Skjetne | A survey on modeling and control of thrusterassisted position mooring systems[END_REF]. Nguyen and Sørensen (2009a,b) proposed the setpoint chasing algorithm in the form of a lowfrequency filter to compensate for the mean environmental load. The setpoint is calculated offline by switching with the predetermined sea condition level. [START_REF] Fang | Optimal set-point chasing of position moored vessel[END_REF] enhanced the algorithm with the mooring tensionbased objective function, considering all mooring cables, to calculate the optimal setpoints online without environment classification. To better capture the dynamic effects and the safety of the mooring system, the structural reliability index is further introduced into the objective function as a weighting factor or in place of the tension term [START_REF] Fang | Optimal set-point chasing of position moored vessel[END_REF][START_REF] Fang | Position mooring control based on a structural reliability criterion[END_REF]. Recent research on setpointchasing algorithms further validates the effectiveness of the aforementioned approach [START_REF] Fang | Mooring system diagnosis and structural reliability control for position moored vessels[END_REF][START_REF] Wang | Numerical and experimental study on the influence of the set point on the operation of a thruster-assisted position mooring system[END_REF][START_REF] Bjørnø | Modeling, parameter identification and thruster-assisted position mooring of c/s inocean cat i drillship[END_REF][START_REF] Xia | Finite time output feedback control for ship dynamic positioning assisted mooring positioning system with disturbances[END_REF]. Nevertheless, all the above forms of objective functions neglect tension and reliability constraints, which can cause them to exceed their admission limits, and even lead to cable failure and safety issues in the worst cases. In addition, intelligent algorithms like the deep deterministic policy gradient approach [START_REF] Yu | Optimal setpoint learning of a thruster-assisted position mooring system using a deep deterministic policy gradient approach[END_REF] and Q-learning techniques [START_REF] He | Energy-efficient control of a thruster-assisted position mooring system using neural q-learning techniques[END_REF] are also applied to setpoint chasing algorithms, but unlike the above methods, this is a model-free, off-policy approach that trains reinforcement learning agents with rewards obtained by interacting with the position mooring system to obtain the optimal setpoints. As a result, it is computationally complex, time-consuming, and has poor real-time performance. In conclusion, in order to provide energy-efficient and safe reference setpoints for position mooring systems in real time, we concentrate on solving the problems of objective function based setpoint chasing algorithms to improve the safety of mooring systems.

Inspired by the above studies, we propose to focus on safe and energy-efficient arctic position mooring control for moored offshore structures in the presence of ice loads, unknown time-varying environmental disturbances, and system uncertainties. The main contributions of the present work include:

1. An ice-load disturbance observer is proposed for arctic disturbances. Compared with existing disturbance rejection methods for position mooring systems, the proposed disturbance observer relaxes the constraints on continuity, boundedness, and differentiability of disturbances. It is the first time to address the harsh environmental disturbance with discontinuous ice loads in arctic position mooring control. 2. An improved setpoint chasing algorithm is proposed to obtain optimal reference positions for the arctic position mooring system. The improved algorithm introduces structural reliability constraints into the objective function to meet safety and energy consumption requirements. Compared with the conventional algorithm [START_REF] Fang | Position mooring control based on a structural reliability criterion[END_REF], the proposed one provides higher safety for all cables of the mooring system. 3. An adaptive fuzzy dynamic surface controller with disturbance compensation and optimal reference position is proposed. A novel filter is introduced into the dynamic surface controller to provide fast convergence, and an adaptive fuzzy estimator is employed to handle the system uncertainties. The proposed controller addresses both crucial issues of safety and energy efficiency, and achieves robust arctic position mooring control with a fast response and arbitrarily small errors.

The remainder of this paper is organized as follows. The problem formulation is described in Section 2. The control scheme design as well as the theoretical analysis of the system stability are discussed in Section 3. Comparative numerical simulations are conducted to demonstrate the effectiveness of the proposed control scheme in Section 4. The paper ends with some concluding remarks in Section 5.

Problem formulation

A turret-moored fully actuated surface vessel exposed to wind, waves, currents, and ice is considered. The vessel is moored through the mooring system, with symmetrically distributed cables connected to the vessel's turret and fixed on the seabed by anchors (eg. illustration of Fig. 1). Since wave-frequency motion should not enter the feedback loop, only low-frequency motion should be considered in the controller design [START_REF] Fossen | Guidance and Control of Ocean Vehicles[END_REF], in which a three-degree-offreedom (surge, sway, yaw) low-frequency model is adopted to design the position mooring control scheme. The proposed scheme should achieve energy efficient and safety, avoiding excessive control efforts for wave-frequency motions.

Kinematics and kinetics

A three-degree-of-freedom low-frequency motion kinematics and kinetics model of a moored vessel can be described as follows [START_REF] Fossen | Marine control systems: guidance, navigation and control of ships, rigs and underwater vehicles[END_REF]. ] T denotes the control efforts produced by thrusters, 𝜏 𝑒𝑛𝑣 represents the unknown environmental disturbances related to wind, waves, currents and ice, 𝜏 𝑚 is the mooring force. 𝐽 (𝜂) is the rotation matrix satisfying ‖𝐽 (𝜂)‖ = 1 and 𝐽 (𝜂) T 𝐽 (𝜂) = 𝐼 3×3 , defined as

η = 𝐽 (𝜂)𝜈 (
𝐽 (𝜂) = ⎡ ⎢ ⎢ ⎣ cos (𝜓) -sin (𝜓) 0 sin (𝜓) cos (𝜓) 0 0 0 1 ⎤ ⎥ ⎥ ⎦ (3)
We consider a mooring system with 𝑘 cables, 𝑇 𝑖 = [ 𝑇 𝑖 𝑋 , 𝑇 𝑖 𝑌 , 0 ] T is the contact force vector of the 𝑖th cable, where

𝑇 𝑖 = 𝑇 𝑖 ( ‖ ‖ 𝑝 (𝑡) -𝑝 𝑖 ‖ ‖ ) ≜ ‖ ‖ 𝑇 𝑖 (𝜂) ‖ ‖ , 𝑝 = [𝑥, 𝑦] T and 𝑝 𝑖 = [ 𝑥 𝑖 , 𝑦 𝑖 ]
T are the positions of the vessel and the 𝑖th anchor. The mooring force exerted on the vessel can then be expressed by low-frequency tension because the system uncertainty has taken into account the cable dynamics. Assumption 1. Considering the nonlinear model ( 2) with the unknown system matrix 𝑀, 𝐶(𝜈), and 𝐷(𝜈), the model uncertainties can be described as Δ = Δ𝑀 + Δ𝐶(𝜈) + Δ𝐷(𝜈). And the system uncertainties, include the model uncertainties and unknown system nonlinear terms, is modeled as 𝑓 = -𝐶(𝜈) -𝐷(𝜈) -Δ.

𝜏 𝑚 = [ 𝑘 ∑ 𝑖=1 𝑇 𝑖 𝑋 , 𝑘 ∑ 𝑖=1 𝑇 𝑖 𝑌 , 0 ] T (4)
Remark 1. All terms, including 𝐶(𝜈), 𝐷(𝜈), and environmental disturbances 𝜏 𝑒𝑛𝑣 are time-varying implicitly obtained unknown terms, and are affected by the environment and ship motion. However, the noise will inevitably exist in position and velocity measurements of ship. Hence, it makes sense to assume that the model uncertainties Δ is related to the system matrix. Assumption 1 is reasonable.

Arctic ice loads

Due to the high energy and time-varying nonlinear characteristics of sea ice behavior, unmanaged drifting ice can generate long-term extreme ice loads that threaten vessel safety [START_REF] Sinsabvarodom | Short-term extreme mooring loads prediction and fatigue damage evaluation for station-keeping trials in ice[END_REF]. For arctic offshore production operations, ice loads can only be reduced to a tolerable level by breaking large pieces of drifting ice into smaller and manageable fractions by icebreakers. Therefore, the arctic ice loads should be considered as managed ice, which can be described as a sequence of impulses exerted on the vessel, mixed with relatively long periods of ice-free load, as illustrated by the curve of Fig. 2 [START_REF] Daley | Gpu modeling of ship operations in pack ice[END_REF]. It has a rapid decay rate during the unloading phase [START_REF] Mååttånen | Iso 19906 ice crushing load design extension for narrow structures[END_REF]. The ice loads have a discontinuous and nonsmooth behavior, with a rapid change rate.

Assumption 2. The arctic environmental disturbances induced by wind, waves, currents, and ice can be described as

𝜏 𝑒𝑛𝑣 = 𝑑 + 𝜎 (5)
where 𝑑 is a smooth function with ‖ ḋ‖ < 𝐷 (𝑡) and ‖𝜎‖ ≤ 𝜎 * . 𝐷 (𝑡) and 𝜎 * are a continuous function and an unknown positive constant, respectively.

Remark 2. Arctic environmental disturbances 𝜏 𝑒𝑛𝑣 include ice loads 𝜏 𝑖𝑐𝑒 and wind, wave and current disturbances 𝜏 𝑑 .

Setpoint chasing

Dynamic surface controller

Ice-load disturbance observer Moored offshore structures

K d K Adaptive fuzzy estimator i G 1 s 2 s d f ˆT W [ , K Q m W Figure 3: Block diagram of the proposed arctic position mooring.
It is discontinuous and non-differentiable due to the nonsmooth ice load behavior. Such disturbance can be decomposed as ( 5), with the continuous part denoted by the smooth function 𝑑, and 𝜎 represents the remaining discontinuous part.

Moreover, the non-smooth characteristics of the system uncertainty can also be included in the above disturbances. Disturbance as in Assumption 2 such that its derivative ḋ exists satisfying ‖ ḋ‖ < 𝐷 (𝑡), releases the requirement that the disturbance should be discontinuous and differentiable. However, prior knowledge that the disturbance and its derivative are bounded requires further observer design to resolve.

Structural reliability index

To balance the demand for energy saving and mooring system security, [START_REF] Madsen | Methods of structural safety[END_REF] introduced the structural reliability index as a quantitative safety criterion to describe the possibility of cable breakage.

𝛿 𝑖 = 𝑇 𝑏,𝑖 -𝑇 𝑒𝑥𝑖 𝜎 𝑏,𝑖 = 𝑇 𝑏,𝑖 -𝑘 𝜎 𝜎 𝑖 -𝑇 𝑖 𝜎 𝑏,𝑖 𝑖 = 1, ⋯ , 𝑘 (6)
where 𝑇 𝑏,𝑖 is the average breaking strength of the 𝑖th cable, 𝑘 𝜎 is a scaling factor, 𝜎 𝑖 = 7.5%𝑇 𝑖 is the standard deviation of the time-varying mooring force, 𝑇 𝑖 is the low frequency mooring force, and 𝜎 𝑏,𝑖 = 7.5%𝑇 𝑏,𝑖 is the standard deviation of the average breaking strength. 𝜎 𝑖 includes high-frequency mooring forces, a slowly varying measure describing the sea state, which does not induce the wave-frequency motion of the vessel. So structure reliability is introduced considering the dynamic effects of mooring tension fluctuations. It will decrease with increasing mooring tension, while the cable breakage probability will increase.

Safe and energy efficient position mooring control scheme design

The main control objective is to design a safe and energyefficient control scheme for arctic position mooring that maneuvers the vessel toward a setpoint where the mooring system is fully and safely used, while addressing discontinuous environmental disturbances and system uncertainties. The block diagram of the proposed arctic position mooring system is depicted in Fig. 3.

Improved setpoint chasing algorithm

We improve the setpoint chasing algorithm by introducing structural reliability constraints into the objective function as follows

𝑓 ( 𝛿 𝑖 𝑒 ) = 𝑘 ∑ 𝑖=1 𝑎 𝑖 𝛿 𝑖 𝑒 2 = 𝑘 ∑ 𝑖=1 𝑎 𝑖 ( 𝛿 𝑖 -𝛿 𝑖 𝑑 ) 2 ( 𝛿 𝑖 > 𝛿 𝑖 𝑑 ) (7)
where 𝑎 𝑖 is the weighting factor, and 𝛿 𝑖 𝑑 is the desired structural reliability.

This objective function enables the structural reliability index of each cable as close as feasible to the desired value, so that the full capability of the mooring system can be utilized. And the reliability constraints 𝛿 𝑖 > 𝛿 𝑖 𝑑 are introduced to ensure that the structural reliability is greater than the desired value to avoid cable breakage and guarantee the safety of the mooring system. Note that the structural reliability accounts for the dynamic effect of the tension to offer a more conservative solution than the tension-based objective function.

It is evident from ( 4) and ( 6) that the tension and the vessel's position are related to structural reliability. By minimizing the reliability-based objective function ( 7) while satisfying the constraints, the optimal operating position can be determined. To get the optimal position increment, the mooring tension is expressed in terms of the vessel's horizontal motion [START_REF] Strand | Design of automatic thruster assisted mooring systems for ships[END_REF].

𝑇 𝑖 = 𝑇 𝑜𝑖 + 𝑐 𝑖 Δ𝑑 = 𝑇 𝑜𝑖 -𝑐 𝑖 Δ𝑟 sin ( 𝛽 + 𝛽 𝑖 𝑜 ) (8)
where 𝑇 𝑖 and 𝑇 𝑜𝑖 are the mooring forces at the current and initial positions, respectively. Δ𝑑 is the change in horizontal distance between the ship and the anchor. 𝑐 𝑖 denotes the incremental stiffness tension [START_REF] Faltinsen | Sea loads on ships and offshore structures[END_REF]. The changes in ship position and heading are indicated by Δ𝑟 and 𝛽, respectively. The structural reliability-based objective function is connected to the vessel's position increment and heading by substituting ( 6) and ( 8) into (7). The constrained multivariable function solver is then used to determine the best solution for the minimization objective function, as well as the optimal position increment and heading for the vessel. The reference setpoints can be derived from the initial position.

𝜂 𝑠 = 𝜂 𝑜 + Δ𝑟 [ cos 𝛽 sin 𝛽 0 ] T (9)
To generate a smooth path passing through the setpoints, a reference model is designed as follows [START_REF] Sørensen | Marine cybernetics: Modelling and control[END_REF]

η𝑟𝑒𝑓 = -Λ𝜂 𝑟𝑒𝑓 + Λ𝜂 𝑠 η𝑑 = 𝜈 𝑑 ν𝑑 = -𝜆 1 𝜂 𝑑 -𝜆 2 𝜈 𝑑 + 𝜆 1 𝜂 𝑟𝑒𝑓 (10)
where 𝜂 𝑟𝑒𝑓 is the low-pass signal of the setpoints 𝜂 𝑠 , and Λ is a diagonal gain matrix. 𝜂 𝑑 and 𝜈 𝑑 represent the desired position and velocity, respectively. The design parameters 𝜆 1 and 𝜆 2 are the diagonal stiffness matrix and the nonnegative diagonal damping matrix.

Adaptive fuzzy dynamic surface control

Lemma 1. [START_REF] Wang | A course in fuzzy systems[END_REF] For any real continuous function 𝑓 (𝑥) ∈ 𝑈 (𝑈 ⊂ 𝑅 𝑛 ) and any arbitrary 𝜀 * > 0, there exists a fuzzy system f (𝑥) such that:

sup 𝑥∈𝑈 | | | 𝑓 (𝑥) -f (𝑥) | | | ≤ 𝜀 * (11)
Lemma 1 states that a fuzzy system can approximate any real continuous function 𝑓 (𝑥) ∈ 𝑈 with an arbitrary accuracy. Hence, we construct an adaptive fuzzy estimator to estimate the continuous unknown system uncertainties.

𝑓 ( 𝑥 ) = 𝑊 T 𝜉 ( 𝑥 ) + 𝜀 ( 𝑥 ) (12)
where 𝑥 = [𝑥, 𝑦, 𝜓, 𝑢, 𝑣, 𝑟] T is the fuzzy system input, 𝜀 ( 𝑥 )

is the minimum approximation error, satisfying

‖ ‖ ‖ 𝜀 ( 𝑥 ) ‖ ‖ ‖ ≤ 𝜀 * with 𝜀 * > 0,
and 𝑊 is the ideal parameters vector.

𝑊 = arg min { sup | | | 𝑓 ( 𝑥 ) -Ŵ T 𝜉 ( 𝑥 ) | | | } (13)
The vessel dynamics (2) can be expressed as follows by introducing the fuzzy estimator (12) to handle the system uncertainties 𝑓 (𝑥) = -𝐶(𝜈)𝜈 -𝐷(𝜈)𝜈 -Δ:

𝑀 ν = 𝑊 T 𝜉 + 𝜀 + 𝜏 + 𝑑 + 𝜎 + 𝜏 𝑚 (14)
To proceed with the control design, we integrated a novel filter to enhance the control performance of dynamic surface control. It increases the convergence rate to provide a fast response and better adaptation to rapidly changing arctic environments.

Theorem 1. Regard a first-order filter with the following form [START_REF] Wang | Adaptive backstepping sliding mode control of uncertain semi-strict nonlinear systems and application to permanent magnet synchronous motor[END_REF]:

̇φ 1 = 𝜙 1 -φ1 𝑇 𝑓 + 𝜒 1 ( 𝜙 1 -φ1 ) ‖ ‖ 𝜙 1 -φ1 ‖ ‖ 1 2 + 𝜒 2 (15)
If the input's derivative φ1 is bounded, there exist positive constant matrices 𝑇 𝑓 , 𝜒 1 , and 𝜒 2 which make the filter error 𝑒 𝑓 = φ1 -𝜙 1 arbitrarily small in a finite time.

The proof of Theorem 1 is detailed in Appendix A.

Remark 3. Note that the proposed filter (15) can be reduced to a conventional first-order filter when 𝜒 1 = 0. It differs from the conventional filter in that it enables the filtered error signal to converge in finite time with a fast convergence rate. Additionally, the selection range of parameter 𝜒 2 is increased by the denominator's fractional order power.

The filter ( 15) is introduced to complete the dynamic surface control following two steps below.

Step 1: The first dynamic surface 𝑠 1 ∈ 𝑅 3 describing the position error is defined as

𝑠 1 = 𝜂 -𝜂 𝑑 (16) A virtual control law 𝛼 1 ∈ 𝑅 3 is designed to stabilize 𝑠 1 . 𝛼 1 = 𝐽 -1 ( -𝑘 1 𝑠 1 + η𝑑 ) (17)
where 𝑘 1 ∈ 𝑅 3×3 is a positive definite diagonal matrix. Then, 𝛼 1 is fed into the introduced filter. It is practical to replace α1 with ̇ᾱ 1 to avoid multiple differentiation.

̇ᾱ 1 = 𝛼 1 -ᾱ1 𝑇 𝑓 + 𝜒 1 ( 𝛼 1 -ᾱ1 ) ‖ ‖ 𝛼 1 -ᾱ1 ‖ ‖ 1 2 + 𝜒 2 (18)
The filter error is defined as

α = ᾱ1 -𝛼 1 (19)
Step 2: The second dynamic surface 𝑠 2 ∈ 𝑅 3 , describing the velocity error, is defined as

𝑠 2 = 𝜈 -ᾱ1 (20)
It follows from ( 14) that the derivative of 𝑠 2 is

ṡ2 = 𝑀 -1 ( 𝑊 T 𝜉 + 𝜀 + 𝜏 + 𝑑 + 𝜎 + 𝜏 𝑚 ) -̇ᾱ 1 (21)
Finally, we have a robust nonlinear control law for the position mooring system.

𝜏 = -𝑘 2 𝑠 2 + 𝑀 ̇ᾱ 1 -Ŵ T 𝜉 -d -𝜏 𝑚 (22)
where 𝑘 2 ∈ 𝑅 3×3 is a positive definite diagonal matrix. d is the estimation of 𝑑 in (5) by the ice-load disturbance observer designed in the next subsection.

The adaptive law Ŵ is designed to update the parameters of fuzzy estimator (12) as follows

̇Ŵ = Γ ( 𝜉 ( 𝑥 ) 𝑠 2 -𝛾 Ŵ ) ( 23
)
where Γ is a positive definite adaptation gain matrix, 𝛾 is a parameter matrix.

Ice-load disturbance observer design

According to Assumption 2 and Remark 2, an ice-load disturbance observer is constructed to address the discontinuous and non-smoothness ice-loads. It is independent of the prior assumption that the disturbance and its derivative are bounded.

An auxiliary vector of observer is described as

𝜍 = 𝑑 -𝑘 𝜍 𝑀𝑠 2 (24)
where 𝑘 𝜍 ∈ 𝑅 3×3 is a positive definite diagonal matrix.

In the light of ( 21), the derivative of 𝜍 is

̇𝜍 = ḋ -𝑘 𝜍 ( 𝑊 T 𝜉 + 𝜀 + 𝜏 + 𝑑 + 𝜎 + 𝜏 𝑚 -𝑀 ̇ᾱ 1 ) (25)
Then the estimate of 𝜍 is designed as

̇ς = -𝑘 𝜍 ( Ŵ T 𝜉 + 𝜏 + ς + 𝑘 𝜍 𝑀𝑠 2 + 𝜏 𝑚 -𝑀 ̇ᾱ 1 ) (26)
From ( 24) and ( 26), we construct the ice-load disturbance observer.

d = ς + 𝑘 𝜍 𝑀𝑠 2 (27)
Let us define d = 𝑑 -d = 𝜍 -ς = ς as the estimation error of the observer, we have ̇d from ( 21), ( 26) and ( 27)

̇d = ḋ -̇d = ḋ - ( ̇ς + 𝑘 𝜍 𝑀 ṡ2 ) = ḋ -𝑘 𝜍 ( W T 𝜉 + ς + 𝜀 + 𝜎 ) (28)
Remark 4. It is obvious from (28) that a differentiable disturbance is necessary for the ordinary disturbance observer.

In contrast to that, the suggested observer ( 27) effectively handles discontinuous ice-loads owing to disturbance reconstruction (5). Additionally, it requires no prior assumption of the derivative boundedness and further resolves the unknown 𝐷 by the invariant set theory.

Theorem 2. For the arctic moored vessel ( 1) and ( 2) with unknown time-varying environmental disturbances and unknown system uncertainties under Assumption 1 and 2, the proposed arctic position mooring robust control law ( 22) consists of a fuzzy estimator ( 12), a novel first-order filter ( 18), an adaptive law ( 23) and an ice-load disturbance observer ( 27). For all 𝑉 (0) < 𝐵 0 with 𝐵 0 being any positive constant, all signals of the closed-loop system are guaranteed to be uniformly ultimately bounded. And the position error 𝑠 1 , the adaptive gain error W and the observer estimation error d are settled respectively within the following compact sets:

Ω 𝑠 1 = { 𝑠 1 ∈ 𝑅 3 | | | ‖ ‖ 𝑠 1 ‖ ‖ ≤ √ 𝜌 𝜇 } (29) Ω W = { W ∈ 𝑅 3×𝑁 | | | ‖ ‖ W ‖ ‖ ≤ √ 𝜌∕𝜇 𝜆 max ( Γ -1 ) } (30) Ω d = { d ∈ 𝑅 3 | | | ‖ ‖ d‖ ‖ ≤ 2𝐷 * 𝜆 max ( 𝑘 𝜍 ) + 2 ( 𝜀 * + 𝜎 * ) +2 √ 𝜌∕𝜇 𝜆 max ( Γ -1 ) } ( 31
)
Theorem 2 states that the proposed closed-loop control system is stable and that all of its signals are bounded. Appendix B provides a detailed description of the stability analysis.

A pseudo-code description of the proposed position mooring scheme implementation is provided in Table 1. In initialization phase, the simulation time, initial state of ship, desired reliability index, and control parameters are set. Note that the desired structural reliability is determined through simulations under extreme conditions based on the mooring system characteristics to prevent cable breakage. To generate the reference path, we use the improved setpoint chasing algorithm to minimize the objective function (7) while satisfying the constraints. Equation ( 8) is substituted into equation ( 6) to obtain the relationship between position increment, heading, and structural reliability. It follows from the relationship that the objective function with nonlinear constraints can be solved by the interior point method and the optimal solution, position increment and heading, is obtained. The smooth reference path (10) through the ideal setpoints can then be obtained. Subsequently, the fuzzy logic system (12) and the ice-load disturbance observer (27) are designed to handle system uncertainty and external disturbances, respectively. The adaptive dynamic surface control law ( 22) is designed to track the reference path and compensate for the disturbance estimation. In this way, an improved setpoint chasing algorithm provides the reference path that meets the requirements of cable safety and system energy consumption, and the observer-based adaptive fuzzy dynamic surface controller tracks the reference path, so that a safe and energy-efficient mooring positioning operation is achieved.

Remark 5. The proposed safe and energy-efficient arctic position mooring control scheme has the following application scope and limitations:

1. This paper investigates the regional position mooring of the ships. Due to the limitations of the mooring system, the position mooring system's operating water depth typically ranges from 300-400 m to 1500 m (American Petroleum Institute, 1996). The applied water depth in this paper should therefore fall within the aforementioned range. 2. The application scenario is in arctic, and the environmental disturbance contains wind, waves, currents, and sea ice, with ice load accounting for the main part.

We only consider managed ice with discontinuous and non-smooth behavior and a fast change rate. 3. The proposed ice-load disturbance observer is for discontinuous disturbances without requiring bounded derivatives. 4. The interaction of mooring cables places restrictions on the improved setpoint chasing algorithm, making it unable to regulate the structural reliability indices of all cables to the desired value.

Numerical simulation results and discussion

Numerical simulations have been conducted on a turretmoored floating production storage and offloading (FPSO) vessel to evaluate the effectiveness and robustness of the proposed arctic position mooring control system. The comparative simulations with the state of the art demonstrate clearly the advantages of the proposed control scheme.

Table 1 The algorithm implementation.

Algorithm: Adaptive fuzzy dynamic surface control based setpoint chasing position mooring scheme 1:

Initialization Set 𝑡 = 0, 𝜂(𝑡), 𝛿 𝑑 and control parameters; 2: Procedure 3:

While The controller is not stopped do 4:

Calculate the mooring tension 𝑇 𝑖 ; 5:

Solve the optimization objective function ( 7), find the optimal position increment Δ𝑟 and direction 𝛽; 6:

Generate the smooth reference path 𝜂 𝑑 (10); 7:

Generate the estimation using the fuzzy logic system (12):

-

Performance of the designed control law

In this scenario, the effectiveness of the proposed control strategy is illustrated in a time-varying arctic environmental disturbance considering ice loads.

The FPSO vessel considered in simulations has the following characteristics [START_REF] Sangsoo | Coupled dynamic analysis of thrusterassisted turret-moored fpso[END_REF]: vessel length 𝐿 = 310𝑚, vessel width 𝐵 = 47.2𝑚, mass 𝑚 = 240869𝑡𝑜𝑛 and yaw rotational inertia 𝐼 𝑧 = 1.34 × 10 12 𝐾𝑔 ⋅ 𝑚 2 . The dynamic parameter matrices 𝑀, 𝐶(𝜈), and 𝐷(𝜈) are calculated by using the corresponding parameters [START_REF] Fossen | Mathematical models of ships and underwater vehicles[END_REF]. The system uncertainties term is set to Δ = 0.1 sin(0.001𝑡)𝑀 + 0.1 cos(0.001𝑡)𝐶(𝑣) + 0.1 sin(0.001𝑡)𝐷(𝑣). The considered mooring system includes eight cables attached to the turret [START_REF] Wang | Reliabilitybased robust dynamic positioning for a turret-moored floating production storage and offloading vessel with unknown time-varying disturbances and input saturation[END_REF], as illustrated in Fig. 4, symmetrically anchored at (800, 0), (0, 800), (-800, 0), (0, -800), (565.69, 565.69), (-565.69, 565.69), (-565.69, -565.69) and (565.69, -565.69) following the order of cables No. 1 to No. 8.

The time-varying environmental disturbances consist of wind, waves, currents, and ice. The ice loads are imitated with a controllable pack ice model as shown in Fig. 2 [START_REF] Lee | Tension based heading control strategy of the arctic fpso with dp assisted mooring system[END_REF], and the other disturbances are described by relevant empirical formulas [START_REF] Fossen | Environmental forces and moments[END_REF]. We set the sea state as wind velocity 𝑢 𝑤 = 8𝑚∕𝑠, wind direction 𝛼 𝑤 = 45 • , current velocity 𝑢 𝑐 = 0.4𝑚∕𝑠, current direction 𝛼 𝑐 = 60 • , significant wave height ℎ 𝑠 = 2𝑚, wave direction 𝛼 𝑤𝑎𝑣𝑒 = 15 • and wave period 𝑇 = 8.61𝑠.

The control objective is to maintain the vessel within a safe region in an energy-efficient manner. The region radius 𝑟 max = 110.91𝑚 is determined by simulations under extreme environment. The desired heading is 10 • , the critical and the desired structural reliability are set to 𝛿 𝑐 = 3 and 𝛿 𝑑 = 4 while respecting the constraint 𝛿 𝑑 ≥ 𝛿 𝑐 . The reference position of the vessel is generated by the proposed setpoint chasing algorithm. The initial states of FPSO vessel in terms of position and velocity are 𝜂 (0) = [0, 0, 0] T and 𝜈 (0) = [0, 0, 0] T .

The control design parameters are chosen as: [START_REF] Sinsabvarodom | Short-term extreme mooring loads prediction and fatigue damage evaluation for station-keeping trials in ice[END_REF][START_REF] Sinsabvarodom | Short-term extreme mooring loads prediction and fatigue damage evaluation for station-keeping trials in ice[END_REF][START_REF] Sinsabvarodom | Short-term extreme mooring loads prediction and fatigue damage evaluation for station-keeping trials in ice[END_REF], Γ 1 = Γ 2 = 10 3 𝐼 3 6 ×3 6 , Γ 3 = 10 5 𝐼 3 6 ×3 6 , 𝛾 1 = 𝛾 2 = 10 -3 , 𝛾 3 = 10 -5 and 𝑘 𝜍 = 𝑑𝑖𝑎𝑔(5, 5, 5). Remark 6. Parameters selection criteria: (i) The mooring system does not involve special cases, equally treats all cables with 𝑎 𝑖 = 1. And the weighting factor can be adjusted according to the cable importance while considering specific applications such as fault-tolerant control and cable breakage. (ii) As for the filter parameters, 𝜒 1 and 𝜒 2 are chosen large enough and small enough to satisfy Theorem 1, respectively, so that the filter error converges to arbitrarily small and the inequality condition of the filter design is easily achieved. (iii) The control parameters should follow Theorem 2, increasing 𝑘 1 , 𝑘 2 , 𝛾, 𝑘 𝜍 , and decreasing 𝛾. In this way, it is obvious that (B.11) and 𝜇 > 𝜌 bounded and converge to an arbitrarily small neighborhood of zero.

𝑎 𝑖 = 1, 𝑘 1 = 𝑑𝑖𝑎𝑔(2, 2, 2), 𝑇 𝑓 = 𝑑𝑖𝑎𝑔(1, 1, 1), 𝜒 1 = 𝑑𝑖𝑎𝑔(10, 10, 10), 𝜒 2 = [1, 1, 1] T , 𝑘 2 = 𝑑𝑖𝑎𝑔
The obtained simulation results are depicted in Figs.5-10.

Fig. 5 illustrates the position and heading of the FPSO vessel as it follows the reference setpoint. It shows that the vessel is precisely positioned at the reference setpoint, and the dynamic surface control law provides a satisfactory control performance for quickly following the setpoint with small errors under time-varying disturbances. According to the objective function in Fig. 6, it is clear that its value rapidly falls and then maintains the minimum, indicating that the proposed setpoint chasing algorithm effectively exploits the mooring system's capacity. The simulation results of the mooring system depicting the structural reliability index and tension are shown in Fig. 7. Cables 1 and 8 achieve the desired structural reliability index in Fig. 7(a), and provide more mooring forces in Fig. 7(b). It illustrates that the proposed control strategy considers all the mooring system cables in addition to the critical one. However, the cable arrangement in Fig. 4 causes some mooring tension to increase while others decreases as the ship position changes. Similar variations can be found in the structural reliability. It follows that the sepoint chasing algorithm can only regulate part of the cable structural reliability to the desired value, while satisfying the constrained objective function ( 7). Consequently, combined with the position in Fig. 5, only cables 1 and 8 reach the desired structural reliability. Moreover, the safety of the FPSO vessel is ensured since the structural reliability indices are larger than the critical value and the mooring forces are smaller than the mean breaking strength. Fig. 8 illustrates the evolution versus time of the control forces of the proposed controller, maneuvering the positron and heading of the FPSO vessel following the generated setpoint and resisting the environmental disturbances beyond the capacity of the mooring system. It is worth noting from Fig. 7 that most of the disturbance is compensated by mooring cables No.1 and 8 simultaneously. As the tension increases with the number of cables carrying the environmental loads, the required control effort decreases. Fig. 9 and Fig. 10 demonstrate that the proposed fuzzy estimator and disturbance observer are able to estimate the system uncertainties and environmental disturbances. Additionally, Theorem 2's conclusion regarding the boundedness of the adaptive gain and observer parameters is also demonstrated.

To sum up, it can be concluded that the proposed arctic position mooring control scheme achieves the control objective and provides effective, robust, safe and energy efficient control performance in the presence of system uncertainties and environmental disturbances with discontinuous ice loads.

Comparison with other control schemes

The conducted comparative simulations include three case studies: (i) The first simulation compares the suggested control with the typical dynamic surface control to show the benefits of the proposed control law; (ii) The second comparison is conducted with the reliability-based control to demonstrate the energy efficiency of the suggested setpoint chasing based position mooring control scheme; (iii) The last part is dedicated to illustrating the effectiveness and safety of the improved setpoint chasing algorithm. Case Study 1: The dynamic surface control law [START_REF] Du | Robust dynamic positioning of ships with disturbances under input saturation[END_REF] for the comparison is expressed as follows .

𝜏 𝐷𝑆𝐶 = -𝑘 2 𝑠 2 + 𝑀 ̇ᾱ 2 -Ŵ T 𝜉 -d -𝜏 𝑚 ( 32
)
where 𝑇 𝑑 ̇ᾱ 11 + ᾱ11 = 𝛼 1 is the first-order filter of the dynamic surface comparative controller. The control parameters are set the same as our proposed controller. The setpoint chasing error and the filter error are depicted in Figs. 11 and12, respectively. It is evident from Fig. 11 that the reference setting path can be successfully followed by both controllers. However, the proposed controller provides better control performance with a smaller tracking error and a higher convergence rate. Fig. 12 illustrates that the novel filter apparently provides faster convergence speed and smaller overshoot than the state-of-the-art filter. Consequently, the proposed improved dynamic surface control law outperforms the existing one.

Case Study 2: A structural reliability-based control strategy that directly introduces structural reliability as a controlled variable [START_REF] Wang | Reliabilitybased robust dynamic positioning for a turret-moored floating production storage and offloading vessel with unknown time-varying disturbances and input saturation[END_REF] is implemented for this comparison.

𝜏 𝛿 = 𝑀𝑄 -1 ( -𝑘 𝛿2 𝑠 𝛿2 + ̇ᾱ 𝛿1 -Q𝜐 ) + 𝐶 (𝜐) 𝜐 + 𝐷 (𝜐) 𝜐 -𝜏 𝑚 -𝜏 𝑒𝑛𝑣 (33)
The design parameters of this controller are set the same as the proposed controller. The obtained simulation results are depicted in Figs.13-18. It can be clearly seen from Fig. 13 that the reliabilitybased controller enables the critical reliability and heading to reach and maintain their desired values. Although both controllers keep the vessel's trajectory in Fig. 14 within the safe region, the proposed controller's trajectory is shorter and simpler than the state-of-the-art controller's. From the objective function evolution in Fig. 15, the proposed controller provides more stable performance. In that the comparative structural reliability-based controller introduces the structural reliability index rather than position variation, stabilizing the error between the desired value and the minimum reliability index at zero, only the one critical cable with the lowest structural reliability index can be regulated. So that the critical cable of comparative controller in control scheme provides a more effective strategy, as it simultaneously considers all mooring cables. Fig. 17 shows the evolution of the mooring forces, the reliability-based controller regulates only one cable at a time, whereas the proposed controller can regulates multiple cables. Thus, the capacity of the mooring system is used to its full potential, resulting in a reduction in the required thruster force. It can be verified from Fig. 18 that the reliability-based controller requires more control forces than the proposed controller. As a result, the proposed setpoint chasing controller achieves a more energy efficient position mooring control.

Case Study 3: The conventional setpoint chasing algorithm [START_REF] Fang | Position mooring control based on a structural reliability criterion[END_REF] without structural reliability constraints is considered for this simulation scenario. The obtained results are shown in Figs. 19-22. Only the relevant notable cable is depicted, since the others are far from the critical value.

From Fig. 19, it is obvious that the vessel trajectory is constrained within the safe region by the improved algorithm. Contrarily, the conventional algorithm regulates the trajectory exceed the safe margin, and has its final setpoint outside the safe region as seen in Fig. 20. In addition, although the improved algorithm's objective function value in Fig. 21 is higher than that of the conventional algorithm, the lack of reliability constraint in the conventional algorithm leads to tension exceeding the average breaking strength and structural reliability below the available minimal value in Fig. 22. This can result in cable breakage and compromises the security of the operation. Therefore, in order to guarantee the safety of the mooring system, the structural reliability constraints are crucial for the setpoint chasing algorithm.

In conclusion, the improved setpoint chasing algorithm introduces the structural reliability constraints, ensuring that reliability stays above its critical value, maintaining the probability of cable breakage within tolerance and consequently enhancing the safety of the mooring system.

Conclusion and future works

This paper investigated an arctic position mooring control scheme in the presence of system uncertainties and ice load disturbances. The safety and energy-efficiency coordination issues of the position mooring system are handled by the improved setpoint chasing algorithm introducing structural reliability constraints. The restrictive assumption that the disturbance is smooth, continuous, and derivatively bounded is removed by the ice-load disturbance observer. The proposed adaptive fuzzy dynamic surface control law improves the convergence speed and enhances the control performance. Theoretical analysis verified the system stability and the uniform ultimate boundedness of all the closedloop system signals. Comparative simulation studies clearly demonstrate that our proposed control scheme provides an effective, robust, safe, and energy efficient position mooring control. The proposed position mooring control scheme can be safely applied to arctic marine resource exploration. As for future work, input saturation of thrusters can be considered to further improve the safety and practicality.

Appendix

A. Proof on Theorem 1

PROOF. Let us define the following Lyapunov function

𝑉 𝑓 = 1 2 𝑒 T 𝑓 𝑒 𝑓 (A.1)
It follows from (A.1) that the time derivative of 𝑉 𝑓 yields

V𝑓 = 𝑒 T 𝑓 ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ - 𝑒 𝑓 𝑇 𝑓 - 𝜒 1 𝑒 𝑓 ( ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 1 2 + 𝜒 2 ) -φ1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ≤ - 𝑒 T 𝑓 𝑒 𝑓 𝑇 𝑓 - 𝜒 1 ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 2 ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 1 2 + 𝜒 2 + ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ ‖ ‖ φ1 ‖ ‖ ≤ - ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 1 2 ⎛ ⎜ ⎜ ⎜ ⎝ 𝜒 1 ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 3 2 ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 1 2 + 𝜒 2 - ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 1 2 ‖ ‖ φ1 ‖ ‖ ⎞ ⎟ ⎟ ⎟ ⎠ - 𝑒 T 𝑓 𝑒 𝑓 𝑇 𝑓 (A.2) If the bounded input ‖ ‖ φ1 ‖ ‖ ≤ 𝜙 𝑚 , there exists 𝜅 > 1 such that 𝜒 1 = 𝜅𝜙 𝑚 > ‖ ‖ φ1 ‖ ‖ , then (A.2) leads to V𝑓 ≤ - ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 1 2 ⎛ ⎜ ⎜ ⎜ ⎝ 𝜅𝜙 𝑚 ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 3 2 ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 1 2 + 𝜒 2 - ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 1 2 𝜙 𝑚 ⎞ ⎟ ⎟ ⎟ ⎠ - 𝑒 T 𝑓 𝑒 𝑓 𝑇 𝑓 ≤ - ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 1 2 𝜙 𝑚 ⎛ ⎜ ⎜ ⎜ ⎝ 𝜅 ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 3 2 ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 1 2 + 𝜒 2 - ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 1 2 ⎞ ⎟ ⎟ ⎟ ⎠ - 𝑒 T 𝑓 𝑒 𝑓 𝑇 𝑓 (A.3) If inequality ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ ≥ √ (𝜒 2 +1) 2 4(𝜅-1) 2 + 𝜒 2 𝜅-1 holds, then V𝑓 ≤ - 𝑒 𝑇 𝑓 𝑒 𝑓 𝑇 𝑓 - ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ 3 2 𝜙 𝑚 = - 2 𝑇 𝑓 𝑉 𝑓 -2 3 4 𝜙 𝑚 𝑉 3 4 𝑓 (A.4)
Consequently, the filter error can converge to the set

Ω 𝑒 𝑓 = { 𝑒 𝑓 | | | | | ‖ ‖ ‖ 𝑒 𝑓 ‖ ‖ ‖ ≤ √ (𝜒 2 +1) 2 4(𝜅-1) 2 + 𝜒 2 𝜅-1 } in a finite time,
which can be made arbitrary small by appropriately adjusting the parameters 𝜅 large and 𝜒 2 small. In this way, ‖𝑒 𝑓 ‖ can easily fulfill the inequality condition.

B. Stability analysis of the arctic position mooring system

Choose the following Lyapunov function for the whole system.

𝑉

= 1 2 𝑠 T 1 𝑠 1 + 1 2 𝑠 T 2 𝑀𝑠 2 + 1 2 W T Γ -1 W + 1 2 αT α+ 1 2 dT d (B.5)
The time derivative of 𝑉 leads to the following result.

V = 𝑠 T 1 ṡ1 + 𝑠 T 2 𝑀 ṡ2 + W T Γ -1 ̇W + αT ̇α + dT ̇d (B.
6) It follows from ( 16), ( 17), ( 19), (20), ‖𝐽 (𝜂)‖ = 1 and Young's inequality that

𝑠 T 1 ṡ1 = 𝑠 T 1 ( 𝐽 ( 𝑠 2 + ᾱ1 ) -η𝑑 ) = 𝑠 T 1 ( 𝐽 𝑠 2 + 𝐽 ( α + 𝐽 -1 ( -𝑘 1 𝑠 1 + η𝑑 )) -η𝑑 ) = 𝑠 T 1 ( 𝐽 𝑠 2 + 𝐽 α -𝑘 1 𝑠 1 ) ≤ -𝑠 T 1 𝑘 1 𝑠 1 + 𝑠 T 1 𝑠 1 + 1 2 𝑠 T 2 𝑠 2 + 1 2 αT α (B.7)
In the light of ( 21) and ( 22), we have According to (12), ( 23) and Young's inequality, we deduce For harsh marine environments, especially for ice-loads with infinite energy, it should be noted that ‖ ḋ‖ < 𝐷 with bounded 𝐷 cannot be known a priori. Hence, the upper bound of the disturbance derivative here is unknown. A compact set Ω is constructed to address this issue [START_REF] Slotine | Applied nonlinear control[END_REF]. Consequently, the proof of Theorem 2 is concluded.

W T Γ -1 ̇W = W T Γ -1 ( Ẇ -̇Ŵ ) = -W T 𝜉𝑠 2 + 𝛾 W T Ŵ ≤ -W T 𝜉𝑠 2 - 𝛾 2 W T W
Ω = {( 𝑠 1 ,

Figure 1 :

 1 Figure 1: Illustration of a turret-moored surface vessel.

Figure 2 :

 2 Figure 2: Simulated managed drifting ice loads.

Figure 4 :

 4 Figure 4: Mooring cable arrangement

Figure 5 :Figure 6 :

 56 Figure 5: North, east position and heading of the FPSO vessel.

Figure 7 :Figure 8 :

 78 Figure 7: Simulation results of the mooring system.

Figure 9 :Figure 10 :

 910 Figure 9: The parameter estimate vectors of the fuzzy logic system.

Figure 11 :Figure 12 :

 1112 Figure 11: Evolution versus time of the setpoint chasing error of both controllers.

Figure 13 :

 13 Figure 13: Evolution versus time of the control performance of the structural reliability-based controller.

Figure 14 :Figure 15 :

 1415 Figure 14: Trajectory of the FPSO vessel for both control strategies.

Figure 16 :Figure 17 :Figure 18 :

 161718 Figure 16: Evolution versus time of the structural reliability indices of the mooring cables.

Figure 19 :Figure 20 :

 1920 Figure 19: The trajectories of the FPSO vessel under the two setpoint chasing algorithms.

Figure 21 :Figure 22 :

 2122 Figure 21: Evolution of the cost function for the two setpoint chasing algorithms.

(

  𝑊 T 𝜉 + 𝜀 + 𝜏 + 𝑑 + 𝜎 + 𝜏 𝑚 -𝑀 ̇ᾱ 1 )

  𝜌 arbitrarily small, so that V < 0 holds on 𝑉 = 𝐵 0 . Accordingly, Ω = {( 𝑠 1 , 𝑠 2 , W , α, d) ∶ 𝑉 ≤ 𝐵 0 , ∀𝐵 0 > 0 } is an invariant set. This is to say, all the closed-loop system signals are always on Ω according to invariant set theory. Additionally, it is confirmed that the upper bound of disturbance derivative exists on Ω all the time. Now, it is clear from (B.5) that 𝑉 (𝑡) and 𝑠 1 , 𝑠 2 , W , α and d are uniformly ultimately bounded. Then, it follows from (16), (19), (20), (23), and (28) that 𝜂, 𝛼 1 , ᾱ1 , 𝜈, Ŵ and d are uniformly ultimately bounded. Hence, all the signals of the resulting closed-loop system are proved to be uniformly ultimately bounded.According to (B.16), we deduce that 𝑉 (𝑡) ≤ 𝑉 (0) + 𝜌 2𝜇 . In view of 1 2 𝑠 𝑇 1 𝑠 1 ≤ 𝑉 (𝑡) and

								𝑠 2 , W , α,	d)	∶ 𝑉 ≤ 𝐵 0 , ∀𝐵 0 > 0	}	(B.13)	From (B.15), it can be deduced that
	Consequently, there exists a maximum 𝐷 * of the contin-uous function 𝐷 on Ω, such that ‖ ‖ d‖ ‖ ≤ 2𝐷 𝑘 𝜍 + 2 ‖ ‖ W ‖ ‖ +	0 ≤ 𝑉 (𝑡) ≤	𝜌 2𝜇	+	[ 𝑉 (0) -	𝜌 2𝜇	]	e -2𝜇𝑡	(B.16)
	2 (𝜀 * + 𝜎 * ) always holds. In view of 2𝑠 T 2 𝛾‖ W ‖ 2 4 , the derivative V is derived by Young's inequality. ‖ ‖ W ‖ 4𝑠 T 2 𝑠 2 𝛾 + ‖ ≤
	V ≤ -𝑠 T 1 𝑘 1 𝑠 1 + 𝑠 T 1 𝑠 1 +	1 2	𝑠 T 2 𝑠 2 -𝑠 T 2 𝑘 2 𝑠 2
		-+	𝛾 2 1 2	W T W + αT α -1 𝑇 𝑓	1 2 αT α -𝑊 T 𝑊	1 2	dT 𝑘 𝜍	d	𝜌 2𝜇 , and lim 1 2 W T Γ -1 W ≤ 𝑉 (𝑡), we denote 2	[	𝑉 (0) -𝜌 2𝜇	]	e -2𝜇𝑡 by 𝑉 𝑚
	+ 𝑠 T 2 ≤ -𝑠 T 1 𝑘 1 𝑠 1 + 𝑠 T ( 2𝐷 * 𝑘 𝜍 + 2 ‖ ‖ 1 𝑠 1 -𝑠 T W ‖ ‖ + 3 2 𝑘 2 𝑠 2 + 𝑠 T ( 𝜀 * + 𝜎 * ) 2 𝑠 2 + 4 ) 𝛾 𝑠 T 2 𝑠 2 -𝛾 4 ‖ ‖ W ‖ ‖ 2 + 1 2 ‖ ‖ ‖ 2𝐷 * ∕𝑘 𝜍 + 3 ( 𝜀 * + 𝜎 * ) ‖ ‖ ‖ 2 + 1 2 𝑊 T 𝑊 + 1 2 αT α -1 𝑇 𝑓 αT α -1 2 dT 𝑘 𝜍 d (B.14)	and obtain ‖ ‖ 𝑠 1 ‖ ‖ ≤ ‖ ‖ W ‖ ‖ ≤ Recalling (B.10), (B.11) and (B.18), yields √ 𝜌 𝜇 + 𝑉 𝑚 ≤ √ 𝜌 𝜇 √ 𝜌∕𝜇 + 𝑉 𝑚 𝜆 max ( Γ -1 ) ≤ √ 𝜌∕𝜇 𝜆 max ( Γ -1 )	(B.17) (B.18)
	Furthermore, we have V ≤ -[ 𝜆 min ( 𝑘 1 ) -1 ] 𝑠 T 1 𝑠 1	‖ ‖	d‖ ‖ ≤	2𝐷 * 𝜆 min ( 𝑘 𝜍	) +2	( 𝜀 * + 𝜎 * )	+2	√	𝜌∕𝜇 𝜆 max ( Γ -1	) (B.19)
	-	[ 𝜆 min	( 𝑘 2	)	-1 -	4 𝛾	]	𝑠 T 2 𝑠 2
	--	[ 𝜆 min 𝛾 4 W T W + ( 1 𝑇 𝑓	) 2 1	-𝑊 T 𝑊 1 2 ]	αT α -𝜆 min	( 𝑘 𝜍 2	)	dT d	(B.15)
	+	1 2	‖ ‖ ‖ ‖ ‖	𝑘 𝜍 2𝐷 *	+ 3	(	𝜀 * + 𝜎 * ) ‖ ‖ ‖ ‖ ‖	2
	≤ -2𝜇𝑉 + 𝜌
	where 𝜌 = 1 2 design parameter 𝜇 = min ‖ ‖ ‖ ‖ 2𝐷 * 𝑘 𝜍 + 3 (𝜀 * + 𝜎 * ) { 𝜆 min ( ‖ ‖ 2 ‖ ‖ 𝑘 1 𝜆 min ( 1 𝑇 𝑓 ) -1 2 , 𝜆 min ( 𝑘 𝜍 2 ) } , satisfying 𝜇 > 𝜌 + 1 2 𝑊 T 𝑊 , positive ) -1, 𝜆 min (𝑘 2 )-1-4∕𝛾 , 𝜆 max (𝑀) 2𝐵 0

. 𝜆 min and 𝜆 max are the minimum and maximum eigenvalues of a matrix.

Remark 7. It is clear from (B.15) that there exists a compact set Ω such that V ≤ -2𝜇𝑉 + 𝜌 holds. Appropriately adjust the parameters with increasing 𝑘1, 𝑘2, Γ and 𝑘 𝜍 to make 2𝜇 𝑡→∞ 𝑉 (𝑡) =

⋆ This work is the result of the research project funded by the National Natural Science Foundation of China under Grant No.51879049, and the Natural Science Foundation of Heilongjiang Province under Grant No.LH2019E039.