A probabilistic approach for Dynamic Positioning capability and operability predictions

both to a reference supply vessel and a pipe-lay vessel shows the flexibility of the proposed approach for site-specific Dynamic Positioning capability predictions.


Introduction
The performance evaluation of any Dynamic Positioning (DP) system is a mandatory step for the design process of an offshore unit (Balchen et al., 1976), as the DP system is one of the fundamental installations equipping ships working in an offshore environment (Kumar, 2020).Such an analysis allows determining the maximum environmental forces that the DP system can counteract, using all the thruster devices mounted on-board and accounting for system failures (Aalberts et al., 1995).The conventional approach for the study of DP system capabilities aims to determine the maximum sustainable wind speed the unit can face in prescribed environmental conditions (ABS, 2014;DNV, 2021), obtaining the so-called DP capability plots (IMCA, 2000) or regulatory indices assessing DP performances (DNV, 2011;BV, 2021;LR, 2021).This target and its representation can be achieved with different prediction methodologies, implying the adoption of quasi-static calculations (Aalberts et al., 1995;Wang et al., 2018) or time-domain simulations (Smogeli et al., 2013;Mauro and Gaudiano, 2018;Martelli et al., 2022).The approximated evaluation of wind speed   in the ''scatter diagram approach'' is a limitation to its full applicability for dynamic positioning performance predictions in a specific sea environment.Furthermore, the ''scatter diagram approach'' is entirely decoupled from the standard capability calculations, forcing the end users to perform conventional analyses besides the quasi-probabilistic approach to visualise the traditional DP capability.The present work fills this gap and presents a new ''fully probabilistic'' methodology to assess site specific DP capability and operability predictions by improving the wind-wave correlation modelling of the ''scatter diagram approach''.
The development and testing of a novel and original prediction process based on a Monte Carlo (MC) integration allow evaluating the DP operability sampling the environmental conditions from site-specific joint distributions for wind and wave characteristics (Johannessen et al., 2001).However, the DP analysis with the base MC process on two reference vessels in two different sea areas, where environment statistics is available (Li et al., 2013), highlights the necessity to improve the sampling process by implementing an enhanced Quasi-Monte Carlo (QMC) method (Niederreiter, 1992) to reduce the variability of the solution and decrease the total amount of calculations necessary to determine the vessel's DP operability with a given accuracy.
Furthermore, besides operability, the enhanced ''fully probabilistic'' method allows deriving a site-specific capability plot, filling the second weakness of the ''scatter diagram approach'', visualising the maximum sustainable wind speed at each vessel heading  and inheriting the probabilistic nature of the novel DP performance assessment.This last aspect increases the feasibility of the new method, providing an additional output which is a state-of-the-art diagram easily interpretable for offshore operators.
The novel contributions provided by the ''fully probabilistic'' DP analysis cover the following main focal points: • Probabilistic modelling of site-specific environmental conditions.
• Operability calculation for a specific operative site through an MC integration process.• Optimisation of the MC process employing an enhanced QMC algorithm.• Evaluation of DP capability directly from operability calculations.
The present work also allows to clearly distinguish between the available processes to evaluate the capability and operability of a DP system, providing a comprehensive overview and comparison of the methods.The results of the reference cases highlight the importance of enhanced modelling of the combined wind-wave characteristics for a more realistic operability prediction of the DP system.The comparison with standard deterministic and standard site-specific DP analyses is self-evident.Table 1 summarises the main differences between inputs and outputs of the existing processes, highlighting the completeness of the newly proposed ''fully probabilistic'' approach compared to other procedures.
The paper has the following structure: Section 2 presents the reference vessels and the considered sea areas that will be employed through all the calculation presented in this work.The standard quasistatic approach is reviewed in Section 3 to explain the background of the adopted thrust allocation algorithm for equating the environmental loads and the applied external forces.
In Section 4 the hypotheses of the methods #1 to #3 in Table are reviewed together with the corresponding environmental conditions used for DP assessment.The necessity of an enhanced modelling of long-term environmental conditions is discussed in Section 5 by presenting both bivariate and trivariate joint distributions of environmental parameters.
Section 6 deals with the fundamentals of the newly developed DP prediction method #4 (see Table 1) based on Montecarlo and Quasi-Montecarlo samplings.This novel method is physically well grounded and can be suitably used to supply both the DP operability evaluation on the scatter diagram statistics.The way to recover the capability plot from the samplings is detailed in Section 7, where an original vision of the standard DP capability plot in terms of probabilistic safety is introduced for the first time.Final remarks and conclusions are given in Sections 8 and 9, respectively.

Reference vessels and sea areas
The novel procedure for the probabilistic site-specific DP analysis is here applied on two reference vessels: one OSV unit and a PLCV one.These two vessels have been selected due to the availability of information on the thruster's size and their layout but also because the environmental loads are available from model experiments.As the developed methodology concerns site-specific calculations, use is made of two reference areas where detailed environmental statistics are available in literature.Hereafter, the vessels and the selected geographic areas are further described.

Reference OSV and PLCV
For sake of brevity, the two reference vessels will be referred through the text as Vessel-1 and Vessel-2, see Table 2 for their main characteristics.Vessel-1 is an example of OSV vessel, which is equipped with two azimuth thrusters and three tunnel ones, two on the bow and one on the stern.The configuration of the thrusters is symmetric and reflects a standard layout for OSVs.On the other hand, Vessel-2 is a PLCV vessel equipped with six azimuth thrusters.The vessel presents a peculiarity compared to the previous layout, as the foremost thruster T1 is not located on the ships centreline, resulting in an asymmetrical thruster configuration.Moreover, the vessel is equipped with a cantilevered tower for J-lay operations located on the vessels side.Table reports the thruster sizes and their locations for both vessels using the reference system of Fig. 1, together with the magnitude of the pipe loads applied to Vessel-2.Environmental loads coefficients for the two reference vessels are available from model experiments.Table 3 reports also the interaction areas between the azimuth thrusters, showing the main interaction angle   and the amplitude of the interaction area .Due to the clustering of thrusters on the bow, Vessel#2 presents two interaction areas for the bow devices, stern devices have only one as T6 is located at a different depth than T4 and T5.

Reference sea areas
This study aims at the determination of a site-specific DP capability prediction.Therefore, to test the procedure, use has been made of two reference geographic sea areas where joint probability distributions representative for   ,   and   is available from the literature (Li et al., 2013).The two sites are representative of one buoy located in the Atlantic Ocean near Cabo Silleiro and the other in the North Sea near Norway's shore.
For brevity, the first sea area is named Site A and the second one Site B. Fig. 2 shows the scatter diagrams of the two reference seaareas (normalised on 1000 wave realisations according to the Global Wave Statistics (Hogben et al., 1986)).A more detailed description of the joint environmental data of the two sea areas will follow in the dedicated sections.

Quasi-static Dynamic Positioning calculations
An exhaustive and detailed analysis of a DP system requires the execution of time-domain simulations, considering both vessel dynamics and the control system.Models for DP simulations generally employ 3 DOF dynamic equations derived from manoeuvring theory under vessel low-speed assumptions.With reference to system in Fig. 1, motion equations for a vessel with mass  and inertia   have the following form: where   ,   are the centre of gravity coordinates.Variables ,  and  are the surge, sway and yaw velocities, respectively, while ,  and  are external forces and moments acting on the vessel, including the thruster forces derived by the control system.
DP time-domain simulations require considerable computational effort and need the availability and knowledge of all the data necessary to set up both the hydrodynamic and controller models.An alternative and less computational demanding approach is the implementation of a quasi-static DP calculation.Even though this alternative neglects the dynamic effects or only accounts for them with an empirical correction on external loads, the quasi-static method allows the investigation of multiple environmental conditions in a reasonably short amount of time.Therefore, such a method is more suitable to study alternative allocation methods and vessel operability prediction without excluding future extension to the time-domain approach.
The quasi-static DP problem, neglecting vessel dynamics, requires the resolution of the forces and moment equilibrium on the horizontal plane.Therefore, adopting the same reference system of Fig. 1, the following equations system needs to be solved: where   and   are the thrust components in  and  directions delivered by the   on-board actuators, located at   and   coordinates.
The right-hand terms of system (2) represent the environmental and external loads.The latter allows the simulation of additional forces and moments acting on the vessel, like mooring lines or pipeline tensioning.
The resolution of the quasi-static problem requires solving system (2), assuming the thrust components as unknowns.

Thrust allocation algorithm
The total number of unknowns in system (2) is    + 2  , with    being the number of fixed tunnel thrusters and   being the number of steerable ones.Thus, except for the unrealistic case of one steerable and one fixed thrusters, the problem admits infinite feasible solutions, having the equilibrium system only three equations.Therefore, a suitable resolution algorithm has to be used to properly evaluate the distribution and orientation of the thrust among the actuators.This study uses a non-linear optimisation method (Mauro and Nabergoj, 2016), employing a non-linear objective function representative of the sum of absorbed power by each thruster: The algorithm can handle as unknowns both the thrusters thrust components or the thrust intensity and its orientation (Mauro et al., 2021).Here, the selection of thrust components as unknowns allows keeping linear the main optimisation equality constraints, i.e. the equilibrium equations of system (2).Besides, additional inequality constraints have to ensure that the actuators deliver thrust up to their maximum limits: Such constraints-set is non-linear.Thus, non-linear programming or constraints linearisation is suitable to solve the problem.Having the two reference ships' interaction zones between azimuth thrusters (see Table 3), the additional constraints are not sufficient to solve the optimisation problem.More advanced thrust allocation algorithms automatically evaluate the thruster-thruster interaction effect (Prpić-Oršić and Valčić, 2020;Mauro and Nabergoj, 2016;Arditti et al., 2019), autonomously giving preference to thruster orientations with lower power demand (Mauro and Nabergoj, 2015a;Arditti et al., 2015).However, such algorithms rarely apply to standard quasi-static predictions, and a conventional allocation method is more suitable to present the novel probabilistic procedure.With the adopted thrust allocation strategy, additional implementation of forbidden zones can model the thruster-thruster interaction between azimuth devices, forcing the thrust outside the interaction area by automatically adding additional equality constraints when needed (Mauro, 2022).

Environmental and external loads
A relevant aspect of a quasi-static DP prediction is evaluating environmental loads acting on the unit.Besides, particular vessel types as PLV may be subjected to additional external loads like pipeline tensioning.A general decomposition of the right-hand loads of equations system (2) is as follows: In the case of pipeline tensioning, the external load is assumed constant with the vessel heading.On the contrary, the estimate of environmental loads is necessary for all encounter angles the unit will face during stationing.The use of non-dimensional coefficients for the environmental loads characterisation grants more flexibility for the reproduction of multiple environmental conditions.This approach is straightforward for wind loads, where the total forces and moment can be expressed by means of non-dimensional load coefficients with the following expressions: where   and   are the frontal and lateral areas of the unit superstructures,   is the overall length,   the mean wind speed and   the air density and  the heading angle.Current loads can be described with expression analogue to (6), where superstructure areas are substituted by wetted surface .Reference is made to current speed   and water density for salt water   : Wave loads are described by means of mean drift forces, representative of an irregular and usually long-crested sea for specific couples of   and   .Mean drift forces can be measured through dedicated model scale experiments or derived from the quadratic transfer functions (QTF) calculated by diffraction theory.Modelling the irregular sea environment with a spectrum it results: where ∇ is the vessel volume of displacement,  the acceleration of gravity and   is the wave amplitude spectrum expressed as a function of circular wave frequency .
The simplified quasi-static DP calculations lead to the overestimation of the final station keeping ability of a vessel.It is possible to decrease this overestimation by adding a dynamic allowance to environmental loads, approximatively accounting for dynamic effects.The execution of time domain simulations on the same vessel allows a possible estimation of the dynamic allowance.However, without performing time-domain calculations, use can be made of guidelines given by regulations.In the present study, a dynamic allowance coefficient   =1.25 is used to represent such allowances (DNV, 2021).According to classification societies rules, it is possible to combine the environmental loads in different ways, e.g. by changing their directions or their relative contributions.Different ways of combining the simultaneous action of environmental loads and definition of reference environmental conditions leads to different kinds of quasi-static DP analyses.This problem will be investigated in the following section.

Environmental conditions in Dynamic Positioning analyses
A quasi-static approach in predicting DP capabilities allows performing several analyses, depending on the environmental conditions modelling.Conventional methods employ fixed combinations of wind speed and wave parameters through a deterministic relationship (windwave correlation).However, long-term environmental data improve the prediction reliability for site-specific calculations.Any more detailed probabilistic definition of the combined wind-wave environment requires a substantial change in the calculation process suitable to evaluate the effectiveness of a DP system.
The present section describes the commonly used environmental models for combined wind-wave action, defining the associated standard and more advanced DP analyses.

Wind-wave correlation
The most common and simple way to model the combined action of wind and waves is the use of a wind-wave correlation, providing a unique deterministic relationship between   ,   and   .Even tough the relationship between wave parameters and wind speed changes with the operating location of the vessel, classification societies (DNV, 2021) and operators associations (IMCA, 2000) suggest the adoption of reference values reported in Table 4.
Another possibility is to derive the wind-wave correlation from an irregular wave spectrum formulation as, for example, the Pierson-Moskowitz (P-M) one, solving the following system: where  and  for the standard PM spectrum are two constants equal to 0.0081 and 0.74, respectively.It is also possible to determine wind-wave correlations for sitespecific environmental conditions (DNV, 2021).In such a case, windwave correlations are expressed through simplified models obtained from a direct comparison of cumulative density functions (CDF) of marginal distributions of   ,   and   .Such an approximation generates triplets   −   −   corresponding to the same occurrence in the given CDF analysis.That represents a simplified characterisation of a sea state, applicable when joint statistics of environmental variables are not available.Therefore, for DP applications, a site-specific windwave correlation offers environmental modelling with the same logical structure as the standard correlations already described, just adding specificity for the sea area of interest.Then, the adoption of site-specific wind-wave correlation does not change the evaluation process needed to assess DP capability with standard correlations.
Fig. 3 shows the comparison between conventional wind-wave correlations and those specific for Site A and Site B. All the correlations identify alternative triplets   −   −   for the modelled environment.
Besides wind and wave loads, also current loads have to be taken into account.Regulations consider the current speed constant (except for DNV suggestions for low wind speeds) and all the environmental loads simultaneously acting on the vessel with the same encounter angle.The adoption of such kind of environment leads to the determination of the DP system capability through consecutive quasi-static calculations.
The representation of DP capability predictions can follow several standards, but all of them are different versions of the so-called Capability plot.This diagram represents the maximum sustainable wind speed    the DP system can face, using the on-board thrust devices, at each encounter angle .Fig. 4 shows the capability plots obtained on the two reference vessels using the DNV, the IMCA and the Pierson-Moskowitz wind-wave correlations together with site-specific ones, with   =0.75 m/s for all cases and considering all the on-board actuators running.The capability plots indicate that the adoption of alternative wind-wave correlations influences the resulting DP capability prediction.The differences depend on the vessel type and actuators power.As site-specific correlations are defined up to the maximum wind detected in the sea area, the associated capability plots are saturated at that specific limit, resulting in an apparent lower capability for head and stern directions.
The site-specific correlations give the highest    between all the tested options, both for Vessel-1 and Vessel-2.Even tough the correlations for Site A and Site B are different (see Fig. 3), the maximum capability limits do not differ too much, especially for Vessel-2.Within standard correlations, for Vessel-1 the resulting capability plots with DNV and Pierson-Moskowitz wind-wave correlations show a higher capability for beam directions (approx.4 kn spreading).The capability of Vessel-2 is less influenced by the wind-wave correlation, but differences still remain for beam weather directions (approx. 2 kn spreading).For head and stern directions the DNV correlation gives lower capabilities compared to the other two also because of the lower maximum   value provided by the regulation.
These differences are due to the specific environmental conditions (i.e. the unique combinations of   ,   and   ), which differ for all the considered wind-wave correlations.

Scatter diagram approach
An alternative method to a deterministic wind-wave correlation is the adoption of the quasi-probabilistic scatter diagram approach (Mauro and Prpić-Oršić, 2020).This approach allows having a more comprehensive and complete overview of the DP system performances in a specific sea area.Instead of performing deterministic DP capability plot predictions, the scatter diagram approach allows the execution of DP calculations for each combination of   and   of the operational area.Calculations can be carried out for each heading , evaluating whether the DP system is able or not to keep the vessel in position with the resulting sea environment.
A wave scatter diagram gives statistical data for   and   only.Therefore, the corresponding wind speed remains unknown and should be derived using specific assumptions.The methodology described by Mauro and Prpić-Oršić (2020) uses as main assumption the adoption of a Pierson-Moskowitz spectrum for each considered sea state combination   ,   ; therefore, the wind speed is derived by the formulations given in system (9).As the two equations of the system give different values for   except for the specific conditions identified by the windwave correlation, the method assumes that for each combination of   and   , the maximum wind speed between the two equations is used as reference   for the calculation.Such an approach allows an unique determination of   for each   −   couple, estimating whether the DP system keeps position or not.
In this context, the scatter diagram approach produces two output types.The first is the quantitative evaluation of the DP operability in a sea area, obtained considering the wave occurrences and vessel headings: where  ℎ is the probability associated to each of the  ℎ vessel headings and   is the probability associated with the specific wave conditions given by the scatter diagram.The function   is equal to 1 when the DP system holds the position, 0 otherwise.Secondly, increasing the granularity of the calculations for   and   , the procedure is capable to identify critical curves for the DP system at the encounter angles  of interest.The critical curve represents the operative limit for the DP system at each encounter angle, means the maximum sustainable   −   for each   .Thus, DP critical curves clearly divide the   −   diagram in two zones, one where   is equal to 0 and the other where the function equals 1.Furthermore, the intersections between the critical curves and the P-M wind-wave correlation represent the couples    −  of the conventional DP capability plot.Fig. 5 reports the correspondence between DP critical curves and the capability plot for the two reference vessels.
However, the main goal of the scatter diagram approach is the evaluation of DP operability in a given sea-area through the   index.
To give an overview of the process, Table 5 reports the operability of the two reference vessels for Site A and Site B. The results refer to intact (all thrusters running) and single thruster failure conditions.As a general consideration, the scatter diagram approach highlights that the operability in Site B is higher than Site A for both vessels and all the investigated cases.Despite Site B presents higher   compared to Site A (Fig. 2), where the density of unfavourable   −   conditions is higher for Site A. Such a conclusion cannot be achieved by analysing a standard deterministic capability plot only.

Enhanced modelling of long-term environmental conditions
Having detailed environmental data at disposal from local measurements or forecast models allows investigating alternative approaches for enhanced modelling of environmental conditions for site-specific DP predictions.Hindcast or forecast raw data include the principal parameters necessary for fitting long-term distributions, namely   ,   and   .
Separate distributions can be derived for wind   and wave   −   characteristics, giving the possibility to extend the capabilities of the   scatter diagram approach, while keeping the same calculation assumptions.However a more useful option is to jointly combine wind and wave statistics to overcome the deterministic wind-wave correlation issue of the scatter diagram approach.Therefore, the possible usage of a trivariate wind-wave distribution is here described, underlining the need for a new dedicated calculation process to predict probabilistic DP operability.

Marginal 𝑉 𝑤 distribution
The starting point for the modelling of a joint long-term environmental condition is the identification of a marginal distribution for   .It has been widely observed that the probability density function    for the wind speed can be modelled with a two parameters Weibull distribution (DNV, 2014).In such a case, considering the random variable   in (0, +∞) for the wind velocity, the following relationship is valid: )   (11) where   and   are the shape and scale parameter of the distribution, respectively.The distribution parameters vary with reference to the sea area, applying a proper fitting procedure to the measurements.

Bivariate joint 𝐻 𝑠 − 𝑇 𝑝 distribution
Having at disposal only wave measurements data allows to establish a joint   −   distribution    ,  .Such a possibility provides an alternative to the scatter diagram approach, resulting in a continuous joint distribution instead of a discrete one.The joint distribution for   and   is composed by a marginal distribution    and a conditional distribution    |  .Making reference to the two random variables ℎ  and   in (0, +∞) for wave height and period the following relationship is valid: There are different probability distributions that can be used to fit    .An option is the adoption of a 2 parameter Weibull distribution (DNV, 2014) for   as in Eq. ( 11), or a hybrid fit between a Weibull and a log-normal distribution (Haver, 1980;Li et al., 2013).Concerning the joint distribution    |  , a log-normal model is widely agreed in the literature (Johannessen et al., 2001;Li et al., 2013;DNV, 2014): where standard deviation    and mean value    of ln   are expressed as smooth functions of ℎ  (Li et al., 2013).Adopting equation ( 13) as joint distribution together with a marginal distribution for   allows to have a continuous definition of the traditionally discrete scatter diagram.
However, it must be underlined that the above modelling of joint   −   distribution remains still fully decoupled from the   statistics, leading to adopt same assumptions of the scatter diagram approach to determine the DP operability.A possible solution to couple the two loads, having detailed wind and wave data at disposal, is now described.

Trivariate joint 𝑉 𝑤 − 𝐻 𝑠 − 𝑇 𝑝 distribution
Once both wind and wave statistics are available for the sites of interest it is possible to describe the long-term environmental statistics by means of a trivariate distribution, i.e. a joint   −  −  distribution    ,  ,  .Such a distribution is composed by a marginal distribution    for the wind speed and two conditional distributions,    |  for the wave height and    |  ,  for the periods.With this enhanced modelling,   is conditional to both   and   and not only to   as in Eq. ( 13).The resulting formulation for the trivariate distribution is: where   , ℎ  and   are the random variables as previously defined.
Once again, the marginal distribution for   is the one described by Eq. ( 11).For the conditional distribution    |  an equivalent two parameters Weibull model can be used: Shape and scale parameter of the joint distribution are modelled as power functions of   , leading to the dependence of Eq. ( 15) to both   and ℎ  : where the regression parameters derive from non-linear fitting of raw data.For the conditional distribution    |  ,  the log-normal model of Eq. ( 13) can be used.However, the distribution parameters are now a function of ℎ  and   and not only of ℎ  .Therefore, to provide    and    the following regression models can be used: )  ] The regression parameters for the trivariate joint distribution described by Eqs.(15) (16) (17) are reported in Table 6 for Site A and Site B according to the analysis performed by Li et al. (2013).

Implications for operability calculations
The adoption of joint-distributions for the environmental modelling influences the way operability should be assessed through equation (10).Considering the case of the bivariate   −   distribution (equation ( 12)), the determination of the DP system operability can be derived as a direct integration of Eq. ( 10), by using a Monte Carlo process (as described in the following sections) for each angle .In such a way, the   index is no more subject to the discretisation given by the scatter diagram.
For the trivariate joint distribution of Eq. ( 14), there is no more an unique correspondence between   ,   and   , which is required to apply the scatter diagram approach.Therefore, the evaluation of the operability according to Eq. ( 10) has to be modified, as already highlighted for the case of the joint   −   distribution.
The following section presents the newly proposed DP operability calculation based on an MC process using joint-distributions for the environmental conditions.

DP operability calculation as a Monte Carlo process
Adopting the enhanced environmental modelling described in Section 5, the determination of the DP operability may follow a nondeterministic multidimensional Monte Carlo (MC) integration process.The general approximation of a MC integral has the following form: where  ⊂ R  is an -dimensional probability space,  ∈  is a set of  independent random variables and   is the number of samples.It is convenient to define  as an unit hypercube (0, 1)  leading to ∫  d = 1, and to use uniform random variables  ∼ U (0, 1) to define .Then, Eq. ( 18) assumes the following form: Eq. ( 10), defining   in the scatter diagram approach, can be converted in the form of a sum of MC integrals, using the joint   −   distribution previously described for the  ℎ heading angles of interest.Then the   formulation becomes: The same process may suit also for the trivariate joint   −   −   distribution described in Eqs. ( 14) to (18).Therefore, the corresponding MC-like formulation for the operability index assumes the following form:  From the above Eqs.( 20) and ( 21), it is evident that the environmental joint distributions are direct functions of the uniform random variables .Therefore, the sampling method adopted to generate the random distributions influences the final integration.The common practice for MC integration is performing a direct sampling of  with pseudo-random numbers, adopting the so-called crude MC method (Hammersley and Handscomb, 1964).As the final joint distributions are non-uniform, the associated random variables should be   derived from the sampled .The most commonly adopted method is given by the inversion of the cumulative density function  ().Therefore, being  uniform in [0, 1],  −1 () is distributed according to  , and, for a generic variable , the cumulative  () is consequently uniform in [0, 1].This property is valid also for multidimensional variables and applies to the joint-distributions described above.Fig. 6 shows for the reference angle  of 60 degrees an example of the application of a crude MC process for   calculation of the two reference vessels in Site A and Site B. The example refers to the bivariate case (12) and described by Eq. ( 20) for  ℎ =1 and  ℎ 1 =1, considering   = 10 4 samples.The plot shows that the cases with   = 1 and the cases where   = 0 lay in two distinct regions of the   −   space.These two zones are divided by the DP critical curve (displayed in red) obtained with the scatter diagram approach, as described in Section 4.2.The wind-wave correlation derives from the same assumptions used to calculate   from the distinct couples   −   .Fig. 7 represents the same scenarios of Fig. 6 but sampling the environmental conditions from the trivariate   −   −   joint distribution given in Eq. ( 14) and   given by ( 21).For such conditions it is not possible to clearly distinguish two regions on the   −   space having   = 1 and   = 0. Consequently, there is no correlation between the DP critical curve given by the scatter diagram approach and the MC process with a trivariate distribution.The main reason for this effect is the presence of the marginal distribution for   instead of the simplified deterministic modelling of the scatter diagram approach.It is also clearly visible that the wind distribution affects the final value of the   .
The cases reported in Figs. 6 and 7 refer to a single run performed with a crude MC integration with 10 4 samples.Since the random nature of the sampling process leads to different couples or triplets of environmental parameters at each run, the final integral value and the convergence history are different at each calculation.As the approximated integral converges to an exact value as   increases without upper bounds, then the process is subject to uncertainties.Therefore, the use of a crude MC method requires to adopt a sufficiently large number of samples that ensures the matching of a required confidence level for the solution.This can be achieved by calculating a Confidence Interval () across multiple repetitions   with different sample size   .
Considering the same conditions of the examples reported in Figs. 6 and 7 a set of   = 20 repetitions with different sample size   ranging from 10 3 to 5 ⋅ 10 4 has been performed.Being   not large enough to ensure that  can be described by a normal distribution, the following formulation has been adopted: where  is the mean value of the   repetitions,  is the repetitions variance and  is the inverse cumulative density function of the Student -distribution with confidence level  and   − 1 degrees of freedom.Fig. 8 shows  with the associated 2 interval obtained for the reference vessels and conditions using the bivariate and the trivariate joint distributions.Table 7 reports  and the 95%  for the considered cases and the operability in percentage.The results show that the mean value for the operability is only slightly affected by   , as the differences are in the range of 0.25%   for both bivariate and trivariate joint-distribution cases.On the contrary,   strongly influences the , highlighting the necessity to perform a large amount of simulation to reach stability on multiple repetitions.Fig. 8 clearly shows the spread in the results obtained by multiple repetitions with same   , and how the spread reduces when increasing the sample size.Furthermore, it is also clear that , and consequently , reduces as  is closer to unit.Also considering the cases where the  is lower, uncertainties always affect the first or second decimal digit of the  value, while considering   as a percentage.Therefore, with the adoption of a crude MC process, the final value of   depends on the convergence and variance of the integration process, thus on the selection of both   and   .
Even if the obtained   values for the bivariate case differ from the trivariate case, they are comparable with the scatter diagram approach.This can be expected for the wave distributions in Fig. 2 and is confirmed by the partial   values derived from scatter diagram calculations (using the granularity of Fig. 2).For Site A the values   for  = 60 deg are 35.70% and 54.90% for Vessel-1 and Vessel-2, respectively; for Site B: 40.51% and 61.71%.Therefore, crude MC with a joint   −   gives no additional information compared to the quasi-probabilistic scatter diagram approach.
Therefore, for practical applications, it is convenient to implement a process that is still capable to sample joint-distribution while reducing the variability of a crude MC integration.

Quasi-Monte Carlo sampling for the joint 𝑉 𝑤 − 𝐻 𝑠 − 𝑇 𝑝 distribution
An alternative sampling strategy, aimed to reduce the variance of MC integration, is the adoption of Quasi-random methods (Niederreiter, 1987), distributing samples in  as uniformly as possible with a deterministic approach.Thanks to the adoption of low-discrepancy sequences it is possible to achieve lower errors than crude MC on practical integration problems (Niederreiter, 1992).Between the different deterministic low-discrepancy sequences, the Sobol chain presents an attractive option, giving a good reproduction of the uniform distribution even with a low sample size and without high computational effort (Sobol et al., 2011).Furthermore, the process needed to generate Sobol sequences can be easily implemented in software programs (Levitan et al., 1988;Bradley and Fox, 1988).When a Quasi-random method is used to sample uniform or generic distributions, it is common to adopt the nomenclature of Quasi-Monte Carlo (QMC) sampling.
Fig. 9 shows the differences between crude MC and QMC with Sobol sequences for the sampling of a bivariate uniform distribution with 10 3 samples.This example highlights the differences in the coverage of a sampling space given by the two procedures.QMC method grants a more uniform coverage with a lower number of samples, avoiding agglomeration of points and void spaces typical of a crude MC approach.The adoption of a QMC sampling to non-uniform distributions follows the same concept of the MC sampling and can be addressed with the inversion of the cumulative distribution.In the specific case of the joint   −   −   distribution, the QMC sampling process to obtain   ,   and   can be summarised in the following steps: from the probability density function (PDF) given in Eq. ( 11).
5. Compute   ,   and   as follows: Applying the above steps, it is possible to generate   samples of the joint environmental characteristics for a specific area.Fig. 10 compares the sampling of 10 3 triplets   −   −   for Site A (using the parameters in Table 6 for the joint distribution) resulting from crude MC and QMC processes.From the reported example it is evident that the QMC process improves the coverage of the sample space also for joint-distributions, avoiding agglomeration of samples and unfilled spaces, as already stated for the uniform distributions.
Adopting a deterministic Sobol sequence of numbers, the process is no more stochastic, thus it allows to determine an unique value across multiple repetitions.Fig. 11 shows, for the same cases reported in Fig. 7, the convergence history of the QMC integration process compared to 5 of the 20 repetitions of MC runs performed with 5 ⋅ 10 4 samples each.The graphs show also  and the interval 2 of the MC runs with different   values.It can be observed that the results of QMC integration start to oscillate around the final value when   is above 10 3 for cases where the partial   is above 0.90 for Vessel-2 in Site A, and above 5 ⋅ 10 3 for the other cases.Moreover, the value of the single QMC run is close to the mean of the 20 MC repetitions, except for   = 10 3 , where  value for MC integrations is significant, leading to high spread of the solutions.
Even though the QMC method does not require multiple repetitions of the same run, it is still necessary to identify the number of samples    that ensures a sufficient level of convergence to the partial   value at a given angle  of the environmental loads.The convergence can be evaluated checking the relative differences between partial   at consecutive samples: Attention should be given to the convergence threshold to be considered for the integral.Being the partial   a quantity defined between 0 and 1 that indicates the fraction of year a vessel could operate in a sea area at a given encounter angle, the convergence should be related to the time unit used to quantify the operability.In case it is quantified in days, considering 1 day as convergence threshold, the convergence can be reached when   approaches 2.74⋅10 −3 .Considering a threshold of 1 h, then the   of reference is 1.14 ⋅ 10 −4 .Fig. 12 shows the   variations with   for the same test cases previously used in this section.For each case it is possible to distinguish Having the integrating function   only two possible discrete values, the convergence curve is not continuous.For such a reason, the level of convergence for the   should be checked on the sequence of points referring to   = 1, the higher of the two sequences shown in the graphs.Fig. 12 reports also the two thresholds representing the accuracy of 1 day and 1 h convergence values.It can be observed that, for all the tested cases, the convergence at 1 day level is reached after 300-400 samples, the 1 h level at about 10 4 samples.When the partial   is high, then the convergence is reached for a higher number of samples.This is due to the higher occurrence of cases with   = compared to conditions where partial   value is below 0.9.
The reported cases refer to  = 60 deg.To check the behaviour of the convergence with   at different headings, an additional set of simulations has been performed, reporting here the case for   − in Site B. Fig. 13 shows the convergence of the partial   for the headings of 30, 60, 90, 120 and 150 degrees.The behaviour of the   at the different encounter angles reflects the considerations given for the analyses provided for  = 60 degrees.Therefore, it can be concluded that a   = 10 4 can be sufficient to reach a level of convergence around 1 h for the partial   index.

DP operability evaluation
The QMC process described in the previous section is oriented to obtain the partial   indices specific for a heading angle .According to Eq. ( 21) the partial indices should be summed and weighted for the  ℎ associated to each specific .In case the operative profile of the offshore vessel requires to operate at preferential  angles,  ℎ factors may be not homogeneous and should be specified by the operator.However, for a global and general assessment of the DP operability in a given sea area, the  ℎ can be assumed homogeneous, thus considering that  angles have the same occurrence 1∕ ℎ .Therefore, Eq. ( 21) can be rewritten in the following form: where   |  is the partial   at a given angle .The direct application of Eq. ( 25) for   predictions is not advisable and should be done with caution, as the usage of a not sufficiently high discrete number of headings may lead to inaccuracy in the final results.In case the partial   are not of specific interest for the end user, the QMC process can also include the heading by sampling  as an additional independent random variable with uniform distribution in [0, 2].
As it is convenient to identify the headings with lower partial operability, it is advisable to maintain the calculation process within a discrete number of  angles, adopting then an approximate integration procedure for the final   .To this end, the adoption of a Simpson integration method could help keeping a discrete number  ℎ of  angles.In such a case, Eq. ( 25) becomes: where  () is the abbreviation for   |  .Eq. ( 26) is valid for  in radiant, but substituting 2 with 360 it becomes valid also for  in degrees.Fig. 14 shows the dependence of the final   from the  value, applying equation ( 26) to the two reference vessels in Site A and Site B with no thruster failures (intact condition) and all possible single thruster failures.Using six different , ranging from 5 to 45 degrees, the variations are negligible between =5, 10 and 15 degrees.Increasing the step, variation in the final   raises up to the previously mentioned accuracy threshold of 1 day for 1 year of operation.This is true especially for failure conditions, where the value of   is lower than for intact case.In general, the lower is the   the higher is the importance to choose a low .
Table 8 reports the   values, in percentage form, for the reference cases obtained with the lower of the tested steps  of 5 degrees, being the one granting the better visualisation of the partial   results in polar form, as depicted in Fig. 15.Data in tabular form refer to the intact condition and all the possible single thruster failures.The graphical plot shows the intact condition and the minimum envelope of the single thruster failures, means the minimum partial   at each angle between the single failure options.
Observing the values provided in Table 8, it is possible to quantify the differences between operability in the two sea sites for the reference vessel.Site A is the most favourable area for both the analysed offshore units, leading to about 8 unworkable days in one year for Vessel-1 and 1.5 days for Vessel-2, both with all thrusters running.The unworkable days in Site B raise to about 29 days for Vessel-1 and 9.5 for Vessel-2.Concerning the single failure cases, the most restrictive failure for both vessels is the loss of the foremost actuator T1.For such a failure, Vessel-1 cannot operate for about 33 days in Site A and 71 days in Site B; Vessel-2 cannot operate for 12 days in Site A and 32.5 days in Site B.
The polar plots in Fig. 15 highlight the differences not only in the global operability, but also in the partial indices.It is possible to identify the  angles that are most critical for the vessel operation, both in intact condition and with thruster failures.For Vessel-1, being the DP system symmetric with respect to vessels' diametral plane, the   plots are also symmetric between starboard and port side headings.This is true also for the single failure envelope, as also failures are intrinsically symmetric.For this vessel, the most critical headings are around 60 and the symmetric 300 degrees.Vessel-2 presents an asymmetric DP system layout with additional asymmetry induced by the external pipe load.This reflects on the resulting operability plots, where 60 degrees is also here the most critical heading for operability.However, for port headings the worst case is around 240 degrees both for intact and failure minimum envelope conditions.
Comparing the results of Table 8 with the ones from the scatter diagram approach in Table 5, a difference between 20 and 40% can be observed in the final operability index   for all the tested cases, which is essentially due to the different assessment of wind speed   .This underlines even more the importance of a proper modelling of environmental conditions for DP predictions.
Concerning the computational time, the evaluation of partial   for a single heading  with   = 10 4 takes in mean 1 min on a regular laptop without parallelisation.Therefore, the cases reported in this example with 73 headings take about 1 h and 12 min each, but time can be significantly reduced varying the  and running calculation in parallel.The calculation time is higher than in the scatter diagram approach, as QMC process evaluates more environmental conditions, but the global calculation time makes this approach still applicable for design and analysis purposes.

Recovering the capability plot
The developed site-specific probabilistic operability prediction method based on the sampling of a multivariate distribution presents the DP capability in terms of operability in a given area.This is the same vision of the scatter diagram approach.However, offshore operators are more familiar with conventional capability plots for DP system evaluation, thus by associating DP performance limits with a maximum sustainable wind speed varying with .
The scatter diagram approach, as highlighted in Mauro and Prpić-Oršić (2020), is not capable to reproduce an equivalent capability plot representing the wind limit, being the wind directly evaluated by   and   according to Eqs. (9).The most similar output to a capability plot provided by the scatter diagram approach is the determination of the   associated to the minimum   at which the DP system cannot keep position at the different  angles (see Fig. 5).
With the adoption of the joint   −   −   distribution for the environmental parameters, there is no more a single   for each couple   −   , since   follows the marginal distribution given by Eq. ( 11).Analysing the results of the QMC process for the partial   at different  angles it is possible to distinguish between favourable (  = 1) and unfavourable (  = 0) cases as a function of   .This kind of analysis generates two sub-distributions, with associated probability density  as a function of   , one for   = 1 and another for   = 0.Besides establishing the probability of the two possible values of   , it is convenient to determine the relative frequency  * of favourable and unfavourable cases for specific intervals   : where  0 is the number of cases with   = 0 inside the th interval   ,  1 is the number of cases with   = 1 and   is the total number of samples in the th interval   .The adoption of relative frequency is different from considering an empirical probability density function.The relative frequency  * directly indicates the ratio between favourable and unfavourable cases in specific   intervals.Therefore, it is possible to identify which are the   intervals that are leading to more unfavourable than favourable cases.Being the wind loads function of  2  (as shown in Eq. ( 6)), it is reasonable suppose that increasing   ,  * 0 (   ) increases too, and also the associated   will be higher in average (according to the conditional distribution of Eq. ( 15)) as well as   .This implies that, using a sufficient number of samples,  * 0 (   ) monotonically increases with   .Then, evaluating  * for all the consecutive   , it is possible to identify the   leading to a desired value of  * with a simple interpolation.
Repeating the process for all the  angles it is possible to extract all the   that correspond to a selected relative frequency value.As an example, selecting a  * = 0.5 the limiting wind speed where it becomes less frequent the ability to keep the desired position at the considered  angle can be identify.Reporting the wind speed in a polar plot consequently leads to the determination of a capability plot specific for the considered sea area.
To clarify the procedure, Fig. 16 shows the determination of the two sub-populations of favourable and unfavourable cases and the associated frequencies  * for both the two reference vessels in Site B. The graphic example shows the  angles of 30, 60, 90 and 120 degrees for the intact condition.It can be observed that, being the partial   relatively high for all the headings (see Fig. 15), the population with   = 1 is predominant and follows the behaviour of the Weibull marginal distribution defining   .The population with   = 0 covers the higher wind speeds and has a density function more similar to a Normal distribution.In the same figure it is possible to observe the relative frequency  * that confirms the monotonic behaviour.Comparing the results for the two vessels it is possible to identify that, being   higher for Vessel-2, the population of unfavourable cases presents a narrower peak for Vessel-2 and, therefore, corresponds to an identification of higher   associated to  * =0.5.
Fig. 17 shows the obtained probabilistic-based capability plots for the two vessels in both Site A and Site B, comparing them with the conventional capability plot obtained with IMCA wind-wave correlation and the site-specific wind-wave correlations.The Figure shows the curve referring to  * =0.5 together with the band 0< * <1, thus the whole range of   where DP system may not counteract the external loads.As a general remark, it should be underlined that in case the partial   is equal to 1 for a given , the associated limiting   derived for the site-specific capability plot is equal to the maximum   sampled from the marginal wind distribution.Therefore, for headings close to 0 and 180 degrees the site-specific capability curves are saturated to the relative maximum wind speeds of Site A and Site B. This behaviour is in line with the capability plots obtained with site-specific wind-wave correlations, which are also defined up to a maximum wind speed.
Considering the curve for  * =0.5, the reported example shows an agreement with the shapes of the limiting curve obtained with conventional method.For Vessel-1, the behaviour of the probabilistic sitespecific capability curves is more coherent with the Pierson-Moskowitz and DNV wind-wave correlations shown in Fig. 4 rather then with the IMCA one.The site-specific wind-wave correlation maximum wind limit is higher than the  * =0.5 probabilistic curve for head and stern directions, but the limits are close within each other for beam headings (where the conventional site-specific limit is inside the band 0< * <1).For Vessel-2 the shape is comparable, and the differences between the limiting wind speed values are at most in the order of 3 m/s considering deterministic correlations.The conventional site-specific wind limits are higher than the probabilistic ones for both sites, also considering the band 0< * <1.
It is noteworthy that the predicted ''probabilistic capability plots'' for the two reference vessels are different between Site A and Site B, in line with the   calculation presented in the previous sections.This is an improvement compared to the conventional capability plot provided by site-specific wind-wave correlations (Fig. 4), where the limiting    was almost equal between the two see sites, especially for Vessel-2.
For the reference vessels used in this study, the   vector has been discretised in intervals of 1 m/s, starting from zero wind speed up to the maximum values of 30 m/s.The maximum of 30 m/s has been chosen to include the highest values obtained by sampling   from the marginal wind distribution for Site B, which is more severe then that for Site A as highlighted by   calculations.
It should be noted that the obtained probabilistic site-specific capability plots are fully comparable with the conventional one, as the assumption of collinearity of wind, wave and current loads is the same and also the current modelling is equal.However, even though the results in terms of capability plots are comparable, the new approach based on the probabilistic joint-distribution for wind and waves gives additional information on the effective operability of the DP system in a given area.Thus, the probabilistic method adds to a conventional site-specific prediction all the benefits already provided by the quasiprobabilistic scatter diagram approach and an enhanced definition of the limiting environmental conditions.

Final remarks
The modelling of environmental conditions influences the prediction of DP performances.The adoption of wind-wave correlations (deterministic or site-specific) implies the evaluation of DP performances by means of capability plots and additional quantities associated with the maximum sustainable wind at each encounter angle .The employment of   −   joint distributions (or a scatter diagram) allows for the evaluation of operability of the DP system in a given sea area, which is not directly correlated to DP capability.The newly developed probabilistic method based on trivariate joint   −   −   distributions overcome the disassociated evaluation of capability and operability, allowing for a joint prediction of the two aspects within the same calculation process.
To summarise, the following options can be identified for the quasistatic assessment of DP performances on the base of different modelling for environmental conditions:  1. Deterministic DP capability: standard predictions using prescribed deterministic wind-wave correlations.The calculations determine the DP capability plot for the vessel according to different standards, as IMCA (IMCA, 2000) or DNV Level-1 and Level-2 (DNV, 2021) methods.2. Site-specific DP capability: alternative method to determine DP capability employing site-specific wind-wave correlations instead of deterministic ones.An example of standard for these predictions is the DNV Level-2-site (DNV, 2021).3. Scatter diagram approach: quasi-probabilistic method employing a joint distribution or the scatter diagram for the modelling of wave environment coupled with a PM correlation to determine wind speed.The method does not evaluate a capability plot but predicts the operability of the system for a given site.4. Probabilistic approach: fully probabilistic method for the determination of site-specific wind and wave characteristics with trivariate joint distributions.The method is capable to evaluate DP operability with higher accuracy than the scatter diagram approach.Furthermore, the method provides also a capability plot to comply with site-specific conventional predictions.
The probabilistic approach allows for determining both capability and operability of the DP system for a given sea area, by providing a more detailed and comprehensive overview of DP system behaviour compared to the conventional site-specific prediction and scatter diagram approach.
As a last remark, as the probabilistic methods is an extension of the scatter diagram approach, the new method can be suitable to perform combined station keeping and seakeeping predictions (Mauro and Nabergoj, 2015b).Studies performed on the scatter diagram approach (Mauro and Nabergoj, 2020;Mauro et al., 2021), highlighted that imposing appropriate criteria to ship motions/accelerations in critical locations onboard it is possible to evaluate not only the operability of the system but also of the whole operation.This, in theory, applies also to the probabilistic method, inheriting all the improvements highlighted above.Such consideration further stress the capabilities and future developments of the newly proposed probabilistic process.

Conclusions
The adoption of probabilistic joint-distributions necessitates the implementation of a sampling process to evaluate the wind and waves parameters, adopting a numerical integration based on a Monte Carlo process to determine the DP operability.The present paper compares the performances of a conventional Monte Carlo process with a more convenient Quasi-Monte Carlo sampling based on Sobol sequences.Detailed analysis on two reference vessels (one Offshore Supply vessel and a Pipe-Lay vessel) in two different sea areas highlights the benefit

F
.Mauro and R. Nabergoj

F
.Mauro and R. Nabergoj

Fig. 5 .
Fig. 5. Correspondence between site-specific DP critical curves and deterministic DP capability plots using the Pierson-Moskowitz wind-wave correlation for Vessel-1 and Vessel-2.

F
.Mauro and R. Nabergoj

Fig. 8 .
Fig. 8.  and 2 interval for DP operability at =60 deg using a bivariate (top) and trivariate (bottom) joint-distribution for the environmental conditions.

Fig. 11 .
Fig. 11.Solution history for the partial   at  = 60 deg with MC and QMC process.

Fig. 13 .
Fig. 13.QMC integration convergence for the partial   at multiple  for Vessel-1 in Site-B.

Fig. 15 .
Fig. 15.Partial   polar plots for the reference vessels in Site A and Site B considering all thrusters running (intact ) and minimum envelope of single failures (min.sf.).

Fig. 16 .
Fig. 16.Probabilities and relative frequencies of cases with   = 1 and   = 0 for Vessel-1 and Vessel-2 in Site B.

F
.Mauro and R. Nabergoj

Table 1
Methods for DP analyses according to different environmental modelling.

Table 2
Main characteristics of reference vessels.

Table 3
Thrusters characteristics, thrusters interaction zones and external loads of reference vessels.

Fig. 2. Wave
scatter diagrams for Site A and Site B.

Table 5
DP operability   for Vessel-1 and Vessel-2 according to scatter diagram approach.

Table 7
DP operability at =60 deg using crude MC method with 20 repetitions and increasing sample sizes.

Table 8
DP operability   for Vessel-1 and Vessel-2 according to QMC process.