Simulation-based optimisation for stochastic maintenance routing in an offshore wind farm

Scheduling maintenance routing for an offshore wind farm is a challenging and complex task. The problem is to find the best routes for the Crew Transfer Vessels to maintain the turbines in order to minimise the total cost. This paper primarily proposes an efficient solution method to solve the deterministic maintenance routing problem in an offshore wind farm. The proposed solution method is based on the Large Neighbourhood Search metaheuristic. The efficiency of the proposed metaheuristic is validated against state of the art algorithms. The results obtained from the computational experiments validate the effectiveness of the proposed method. In addition, as the maintenance activities are affected by uncertain conditions, a simulation-based optimisation algorithm is developed to tackle these uncertainties. This algorithm benefits from the fast computational time and solution quality of the proposed metaheuristic, combined with Monte Carlo simulation. The uncertain factors considered include the travel time for a vessel to visit turbines, the required time to maintain a turbine, and the transfer time for technicians and equipment to a turbine. Moreover, the proposed simulation-based optimisation algorithm is devised to tackle unpredictable broken-down turbines. The performance of this algorithm is evaluated using a case study based on a reference wind farm scenario developed in the EU FP7 LEANWIND project.


Introduction
The offshore wind industry is growing at a rapid rate. For example, in Europe the installed capacity has remarkably increased from less than 40 MW in 2000to 12.6 GW in 2016(EWEA, 2017. Although the offshore industry gives more benefits such as the electricity generated and the low visual impact compared to its counterpart onshore (Akbari et al., 2017), the major challenge for the industry is to reduce its high operating costs in order to make this industry profitable. One of the major cost components of operating an offshore wind farm is the cost of the operation and maintenance (O&M). This cost is expected to contribute up to a quarter of the life-cycle cost of an offshore wind farm (Snyder and Kaiser, 2009;Musial and Ram, 2010;Garrad-Hassan, 2013;Wiser et al., 2016). Scheduling maintenance for an offshore wind farm is a challenging and complex task (Shafiee, 2015;Shafiee and Sørensen, 2017). The resources needed to maintain the turbines such as vessels, parts, and technicians are commonly located at the nearest port or O&M base. The offshore turbines are relatively more vulnerable to breakdown than the onshore ones. Moreover, the accessibility to the wind farm sites (using vessels) are restricted and uncertain. In addition, as the number of turbines along with their capacity increases, new wind farms are planned to be located even further offshore. This situation increases the limitation to access the windfarm site.
One way to reduce the maintenance cost is to design the optimal route for each selected vessel to access a set of turbines in the wind farm. This also includes the number of each type of technicians required by each vessel. The principal skill division when considering wind turbine maintenance is between electrical and mechanical expertise. However, there also exists an intersection of tasks that require both electrical and mechanical expertise. Therefore, in this paper, we use three types of technicians: electrical, mechanical and electrical-mechanical. When designing the routes, several factors need to be considered including weather/time window, availability of various resources (e.g. vessels, technicians, and spare parts) and the loss of electricity generation. The weather/time window is defined as the period of time in which the weather is sufficiently clement to allow the maintenance operation to be performed. In our case, the time window is represented by the interval time between the time for a vessel to leave the O&M base and the time for the vessel to return to the O&M base. Two types of maintenance, namely corrective and preventive maintenance, are considered in this paper. Preventive maintenance (PM) is desirable in order to ensure that offshore wind farms operate in an efficient way. PM aims to prevent failures and is comprised of monitoring, detection and restoration actions in order to avoid future breakdowns. Ideally, PM in offshore wind farms will be performed at a time when any associated shutdown will have the least impact on the net energy production capacity (Zhang et al., 2012). On the other hand, corrective maintenance (CM) is conducted following the detection of an unforeseen failure (Halvorsen-Weare et al., 2013).
There are various sources of stochasticity involved in the maintenance of an offshore wind turbine that require to be considered as uncertain rather than deterministic. These principally include turbine breakdown and weather related factors, although these two sources are related on a longer term basis. The uncertainty related to breakdown depends on the past preventative and corrective maintenance history of the turbine as well as factors such as age and location. With a growing size of data and hence evidence regarding offshore wind turbine breakdown, this uncertainty is becoming more quantifiable as time progresses. The second main factor causing uncertainty is the short-term weather conditions at the desired time of maintenance. These have an impact on factors such as travel times and time to complete tasks, as well as ruling out maintenance completely if parameters such as wind speed and wave height surpass pre-defined thresholds. Whilst short-term weather is not predictable in the medium or long term, sufficient meteorological data exists to give quantifiable probabilities for the above weather factors in a given season and location.
In this paper, we first propose an efficient solution method based on large neighbourhood search to solve the deterministic maintenance routing in an offshore wind farm. The solution method is then incorporated in the proposed simulation-based optimisation framework for the stochastic maintenance routing problem where uncertain conditions are taken into account. The uncertain conditions considered are the travel time of the vessel, the transfer time for technicians and equipment to a turbine, and the required time to maintain a turbine. The main contributions of this paper include: (i) A new efficient metaheuristic approach to solve the deterministic maintenance routing problem; (ii) A novel simulation-based optimisation framework for stochastic maintenance routing incorporating Monte Carlo simulation and a metaheuristic.
The remainder of this paper is organised as follows: Section 2 presents a review on maintenance routing in offshore wind farms. Section 3 describes the optimisation model and the metaheuristic technique for the deterministic maintenance routing problem. Section 4 presents a set of computational experiments to analyse the performance of the proposed metaheuristic method in solving the deterministic problem. In Section 5, the simulation-optimisation based model for dealing with uncertainty in the stochastic maintenance routing problem is discussed. Section 6 presents the results of computational experiments to evaluate the performance of the proposed simulation-based optimisation algorithm in solving the stochastic maintenance routing problem. The last section highlights the findings and some avenues for future research.

Literature review
A variety of models ranging from one period to multi-period have been developed.  proposed the maintenance routing for an offshore windfarm for a one-day planning horizon.The model developed is a Mixed Integer Linear Programme (MILP) considering preventive and corrective maintenance activities along with the capacity of vessels transporting technicians.
The model is reformulated as a path-flow model where Dantzig-Wolfe decomposition (Dantzig and Wolfe, 1960) is applied. In their work, a heuristic is proposed to generate a subset of possible routes. Dai et al. (2015) introduced an MILP model that integrates maintenance scheduling and routing problems. The model aims to optimise the schedule and the route for each vessel to maintain a set of wind turbines over a planning period of several days. In other words, the model proposed considers the multiple period problem. In addition, the model takes into account the capacity of vessels in transporting spare parts. This problem was solved using an exact method.
An optimisation model for maintenance scheduling and routing for multiple periods, multiple O&M bases and multiple offshore wind farms was investigated by . The model is an enhancement of the model proposed by Dai et al. (2015) where the aim is to find the optimal schedule for maintaining the turbines and the optimal routes for the vessels to access the turbines.
They designed an efficient algorithm based on the Dantzig-Wolfe decomposition method where an MILP for each subset of turbines is optimally solved to generate all feasible routes and maintenance schedules. Raknes et al. (2017) studied an optimisation model for the maintenance schedule and routing problem taking into account multiple work shifts. Moreover, the model deals with the situation where vessels are able to stay offshore for several shifts. Schrotenboer et al. (2018) investigate the problem proposed by  considering sharing of technicians between wind farms. The authors propose an Adaptive Large Neighborhood Search heuristic for solving the problem. Stock-Williams and Swamy (2019) study the daily maintenance planning problem where metaheuristic optimisation model based on a Genetic Algorithm is combined with a simulator to evaluate solutions.
Most of the papers cited above deal with the deterministic case where all parameters are treated as certain parameters. In this paper, we first propose an efficient metaheuristic approach based on Large Neighborhood Search to solve the maintenance routing problem for a one-day planning horizon. The simulation-based optimisation algorithm combines the metaheuristic with Monte-Carlo simulation to solve the stochastic maintenance routing problem. In the proposed simulation-based optimisation algorithm, Monte Carlo simulation is developed to handle the stochastic parameters including the travel time of each vessel, the required time to maintain a turbine and the transfer time for technicians and equipment to a turbine. The stochastic parameters can be affected by many factors including weather conditions, condition of the turbine, technicians skills, vessel conditions, and the weight of equipment/parts. When the travel time or maintenance time are randomly generated at very high values, the turbines may not be able to be maintained or visited for safety reasons or because of a prohibitively high total maintenance cost.

The optimisation model and solution method for the deterministic maintenance routing problem
This section presents the optimisation model for the deterministic maintenance routing problem (MRP) for an offshore windfarm. It is followed by the proposed metaheuristic method for solving the MRP to determine the best route of each vessel.

The optimisation model
The optimisation model for the deterministic Maintenance Routing Problem (MRP) for an offshore windfarm is developed based on the model proposed by . The model requires information regarding the turbines that need to be maintained and the resources (e.g. vessels, spare parts and technicians) to do so. The information required from the turbines includes: the type of maintenance task (PM/CM), the maintenance/repair time, the number of technicians required, the weight of spare parts needed, downtime cost and turbine penalty cost. The downtime cost is incurred due to maintenance activities as the turbine stops operating. In other words, a revenue loss occurs. The turbine penalty cost exists when the turbine cannot be visited/maintained on a given day due to various factors/constraints including limited resources (e.g. vessels, parts and technicians) and bad weather conditions. For a turbine that needs a preventive maintenance task, the turbine penalty cost occurs when the turbine has to operate at a reduced performance (derated). For a turbine that requires a corrective maintenance task, the turbine penalty cost could be calculated based on the revenue lost for one day as the turbine is not operating (broken-down).
This cost is also determined based on the electricity price on that day.
The required information on vessels includes the travel cost/time for each vessel, the transfer time for technicians and equipment from a vessel to the turbine, the vessel capacity (number of technicians on board and total weight of spare parts/equipment) and the information whether a vessel needs to be present during the maintenance operation on a turbine. The information on the weather window for each vessel includes the earliest time for a vessel to leave from and return to the O&M base. In other words, a vessel is not allowed to stay overnight at sea. In addition, the availability of technicians in the O&M base is also considered. The aim of the models is to minimise the total cost including travel, downtime (PM and CM) and turbine penalty costs. Here, the technician cost is not considered as this cost is relatively small compared to the other costs. The following notations are used to describe the sets and parameters of the deterministic maintenance routing model (MRP).

Sets
V : set of vessels with v as its index I: set of turbines with i as its index and n = |I| I p : set of turbines that need preventive maintenance (I p ⊂ I) I c : set of turbines that need corrective maintenance (I c ⊂ I) I v : set of turbines that can be maintained by vessel v ∈ V I: set of turbines whose require a given vessel during the maintenance operations V i : set of vessels that can maintain turbine i, V i ⊆ V, i ∈ I P : set of technician types (such as electrical, mechanical, and electromechanical) with p as its index Let N = {0, ..., 2n + 1} where node 0 and 2n + 1 represent an O&M base for the beginning and ending of a vessel's route. Each turbine consists of two types of nodes, namely a delivery and a pickup node, where turbine i has delivery node i and pick-up node (n + i). Therefore, a set of delivery nodes is similar to a set of turbines (i ∈ I, i = 1, . . . , n) whereas a set of pickup nodes is denoted by i ∈ I (i = n + 1, . . . , 2n). For modelling purposes, let ρ (n+1)p = −ρ ip for i = 1, . . . , n and p ∈ P with ρ ip = 0 for i = {0, (2n + 1)} and p ∈ P .

Decision Variables
T vi : the time when vessel v ∈ V visits (delivery/pick-up) node i ∈ N Q vip : the number of technicians of type p ∈ P on vessel v ∈ V after leaving node i ∈ N

The objective function
The objective function aims to minimise the total maintenance cost (Z) including travel, downtime (due to preventive and corrective tasks) and turbine penalty costs. The objective function is expressed as follows: where Constraints (6) and (7) ensure that each turbine is visited only once for delivery and only once for pick-up. Constraints (8) and (9) ensure that a vessel leaves and returns to the port only once.
Constraints (10) ensure that a vessel returns to the port from a pick-up node. Flow conservation at each node is expressed by Constraints (11). Here, if a vessel visits a node, this vessel must leave this node as well. Constraints (12) ensure that both delivery and pick up at a turbine are done. Vessel capacity is described by Constraints (13). These ensure that the weight of spare parts transported by a vessel does not exceed its capacity. Constraints (14) -(16) guarantee that the number of technicians onboard does not exceed the vessel's capacity. The number of technicians onboard the vessel when arriving at each node is tracked by Constraints (17) and (18). Constraints (19) and (20) ensure that a vessel must leave after the start of, and return within, the weather window. The travel time compatibility of the vessel is maintained by Constraints (21) where M is an arbitrarily large constant. These constraints also determine the arrival time of a vessel at a node.
Constraints (22) guarantee that the time between the delivery and the pickup is greater than the time required to perform maintenance at the turbine. Constraints (23) ensure that the number of technicians required is less than or equal to the number of technicians available in the O&M base.
Constraints (24) make sure that the vessel travels directly from the delivery node to the pickup node (the same location) if the vessel needs to be present during the maintenance activity. Constraints (25) -(28) also ensure flow conservation at each node. The ability of vessels to maintain a turbine is represented by Constraints (29) and (30). Constraints (31) and (33) indicate the type of the decision variables.

Metaheuristic Approach for the Deterministic Problem
The proposed solution method benefits from the large neighborhood search principle to explore the solution space. The underlying principle of the Large Neighborhood Search (LNS) is to destroy and reconstruct iteratively the current solution in order to improve it (Eskandarpour et al., 2017).
The proposed solution method (LNS) aims to obtain a good solution that minimises the total maintenance cost expressed in Equation (1) within an acceptable computational time. The main steps of the LNS are depicted in Algorithm 1. We first initialise the required parameters (line 1).
Set S * v consists of a set of turbines that will be visited by vessel v and set R stores all the turbines that are not visited by any vessel. We call the latter the request set. At the beginning, all the turbines are assigned to set R = I. The determination of the initial solution (line 2) is detailed in Section 3.2.1. We use simply a maximum number of iterations (IterMax) as the termination criterion (line 6). At each iteration, a so-called removal operator is randomly selected to eliminate a turbine from the current solution and insert it into set R (line 7-9). The resulting solution is called a partial solution S v , after eliminating this turbine. These removal operators are described in Section 3.2.2. Each turbine belonging to the set R will be reinserted to the partial solution S v to construct a new solution (line 11-12). This process is described in Section 3.2.3. Finally, the algorithm returns the best solution along with the turbines visited by each vessel and the ones that are not visited by any vessel.

Initial solution
The initial solution assigns randomly σ turbines to each vessel. σ is set to |I| |V | − 2 . This value is determined based on numerous computational tests and ensures that each vessel is able to visit all the nodes (delivery and pick up) with respect to all the constraints such as those concerning time windows and personnel. Once the sets have been populated, the order of all the nodes (pick up and delivery) visited by the vessel are determined using the reconstruct operator explained in Section 3.2.3. Randomly choose a Removal operator 8:

Algorithm 1 The proposed LNS for the deterministic model
Select a turbine to be removed from S v and move it to set R 9: 10: Reconstruct the solution by inserting turbine R i to solution S v 12:

13:
Denote z and z * the objective values of solutions S v and S * v , respectively 14: if z < z * then 15:

Removal operators
The goal of these operators is to select a turbine to be eliminated from the current solution. For this purpose, five operators are defined to explore the solution space from different perspectives. All the operators work in a hierarchical approach in the way they first determine which vessel should be targeted and then which turbine should be eliminated. To diversify the search, all the removal operators (except Random removal) use the biased roulette wheel selection principle introduced by Prescott-Gagnon et al. (2009). Without loss of generality, let us describe this principle applied to the vessel selection part only. The reasoning is exactly the same for selecting the turbine to remove. The underlying principle is to randomly select a vessel, favoring those providing a greater proximity measure considering the operator of concern. All vessels are sorted in non-increasing order according to the desired proximity measure. The vessel at position N ρ D is chosen, where N is the number of vessels, ρ a number generated randomly between 0 and 1, and D a constant greater than or equal to 1. For these experiments, the value of D was set at 2.
1. Travel cost removal: Select the vessel with the highest traveling cost and then select a turbine that imposes the highest travelling cost.
2. Corrective maintenance time removal: Select the vessel with the highest corrective maintenance cost and then select a turbine that has the most costly corrective maintenance.
3. Preventive maintenance: Select the vessel with the highest preventive maintenance cost and then select a turbine that has the most costly preventive maintenance.

Total cost:
Select the route with the highest total cost and select the turbine with the highest summation of traveling, corrective maintenance and preventive maintenance costs among the ones assigned to this vessel.

Randomly removal:
Randomly select a Route and select a turbine for removal out of it.
Equations 1 are used for the above operators excluding Randomly removal to calculate the sum of traveling, corrective maintenance and preventive maintenance costs associated with each vessel and turbine needed.

Reconstruct operator
The outcome of applying a removal operator is a partial solution, S v , and an updated request set. More precisely, a turbine is selected by the removal operator to be removed from the solution and is inserted into the request set. The goal of the reconstruct operator is to reconstruct the partial solution and find the best permutation of the assigned turbines (delivery and pick up) for each vessel v. Once the order of visiting the turbines is determined x vij , other variables (T vi , Q vip ) can be computed according to Equations 14 -23.
Algorithm 2 details the whole process of reconstructing a partial solution. The algorithm starts with a set of turbines assigned to each vessel and a request set, which contains a set of turbines that are not yet assigned. The idea is to insert each turbine that belongs to the request set into the set of turbines visited by each vessel (S v , v ∈ V ) and explore all the neighboring solutions.
To obtain the neighboring solutions, all the turbines of set S v convert to their equivalent set of delivery nodes and pick up nodes. As far as delivery nodes are concerned, complete permutations are counted. Then, for each permutation of the delivery nodes, all possible permutations of pick up nodes are considered. If the vessel must be present during the maintenance operations, pick up node i + n is placed next to delivery node i. Otherwise, all possible places are iteratively considered as a candidate place for pick up node i + n. Of course, pick up node i + n must be visited after delivery node i. Once the permutation of delivery and pick up nodes are determined, the feasibility check is called to ensure all the constraints are respected.

Computational experiments to assess the proposed LNS
In this section, the computational results for the deterministic maintenance routing problem  Nine instances of each case study are used to assess the performance of the LNS where the number of turbines (n) that are required to be visited is set from 7 to 14, which consist of PM and CM turbines. The number of vessels |V | is set from 2 to 4 with different specifications presented in Table, 1, which is taken from . An example of the test data for |V | = 3 and |I| = 9 is presented in Table, 2, where PM and CM activities need to be performed on 7 and 2 turbines, respectively. In these experiments, it is assumed that the vessels are available on the O&M base and all vessels are able to transport the spare parts required for maintaining the turbines.
The time window for each vessel is set to 12 hours which means that a vessel needs to return to O&M base within 12 hours after the vessel left the O&M base. In addition, the transfer time for technicians (and equipment) from a vessel to a turbine is set to 11 minutes and 45 technicians are available at the O&M base.   T1  PM  7  710  7800  650  4  No  T2  PM  7  628  7800  650  2  No  T3  PM  7  882  7800  650  2  No  T4  PM  7  561  7800  650  2  No  T5  PM  7  486  7800  650  4  No  T6  PM  7  552  7800  650  2  No  T7  PM  7  400  7800  650  3  No  T8  CM  3  500  23400  650  4  No  T9  CM  4  320  23400  650  4  No As the LNS is incorporated within the simulation-based optimisation for solving the stochastic problem (SMRP), it is very important to ensure that the LNS can produce good quality solutions within an acceptable computing time. In this paper, the MRP is solved by six methods: For the West Gabbard case study, Table 3 presents the experimental results for the MRP using various methods (EM, DA, SA, VNS, TS and LNS). Based on the table, the EM generated the optimal solutions for relatively small problems (|V | = 2 with |I| = 6 − 7). However, the EM was not able to guarantee optimality for the relatively large problems ( |I| ≥ 8) within 3 hours. Therefore, for this case the solutions generated by other approaches (DA, SA, VNS, TS and LNS) are not guaranteed to be the optimal solutions. It is worthwhile noting that all methods avoid turbine penalty costs in their solutions produced. In general, the transportation cost accounts for approx.
5% of the total maintenance cost, which is the smallest cost component of the maintenance cost.
On the other hand, the downtime cost for PM is the largest cost component (81%) due to the large number of turbines that need to be preventively maintained compared to the corrective ones (14%).
According to Table 3, in general, the metaheuristic methods perform quite well as they generate   For the Thanet case study, Table 4 shows the results of the solution methods. In line with the previous experiments, the EM was able to obtain the optimal solutions for relatively small problems (|V | = 2 with |I| = 6 − 7). Overall, the metaheuristic methods perform well as relatively good solutions with low %Dev are generated within a relatively short time (181 secs on average). This CPU time is much smaller than the ones of EM (8,457 secs) and DA (22,465 secs). Similar to previous results, among the metaheuristic methods, LNS outperforms the other methods in terms of the quality of solutions produced as it yields a relatively small deviation of 0.38% compared to the best solutions obtained. Therefore, this method can be very useful to be incorporated in the proposed simulation-based optimisation method to address the stochastic maintenance routing problem in an offshore wind farm. This simulation-based optimisation approach requires an optimiser (method) to be executed iteratively. Therefore, a powerful optimiser that runs very fast while generating good solutions is needed. The EM and DA are not practical as they require a long computational time to solve the problem, especially for large problems.

Stochastic Maintenance Routing Problem (SMRP)
This section mainly details the proposed simulation-based optimisation method to address two main streams: (i) uncertain parameters (ii) unexpected broken-down turbines. However, in the first subsection, we present a review on the stochastic routing problem as the MRP in an offshore wind farm can be categorised as a vehicle routing problem with pick-up and delivery (VRPDP) based on the classification proposed by Berbeglia et al. (2007).

Overview of the stochastic routing problem
In this subsection, a review on the routing problem with stochastic travel and service times is presented as in the proposed stochastic routing problem the uncertain travel and maintenance times are also considered. Kenyon and Morton (2003)

Simulation-based optimisation method
In the proposed stochastic optimisation model, we consider several uncertain parameters including the travel time of vessel v ∈ V (τ v ), the required time to maintain turbine i ∈ I (τ i ) and the transfer time for technicians and equipment to turbine i ∈ I (τ i ). In this study, we assume that these parameters follow a normal distribution where the realisations have positive outcomes. Here, we propose a simulation-based optimisation method to solve the stochastic maintenance routing for an offshore windfarm. In this method, the hybridisation of the LNS and Monte Carlo simulation is proposed. The LNS presented in the previous section is used to solve the deterministic model to find the route of each vessel to visit the turbines. When solving this deterministic problem, the stochastic parameters are transformed into their deterministic counterparts. When the route of each vessel has been determined, Monte Carlo simulation is used to obtain an estimate of the expected total maintenance cost by considering the stochastic parameters. Monte Carlo simulation is an iterative process where at each iteration, we generate random numbers to represent the stochastic parameters. For each iteration, we propose a vessel penalty cost for the vessel that returns later thant v .
This can be considered as the recourse cost. Therefore, an additional parameter is added denoted byc to represent the vessel penalty cost per hour for the late vessel. In other words, the stochastic problem has two types of penalty cost. The first penalty cost is the penalty cost for not visiting the turbines while the second one is for the vessels that return outside the time weather/window.

Algorithm 3 presents the main ingredients of the proposed simulation-based optimisation ap-
proach to solve the SMRP. Parameter γ needs to be first defined representing the γ%-quantile total cost data used to determine the expected total maintenance cost. Several scenarios can be performed to obtain a solution including 50%, 70% and 90%-quantile. The 90%-quantile represents the worst conditions where the maintenance duration and travel time are long (above the average).
Therefore, a solution based on the 90%-quantile can be treated as a risk-averse solution as it deals with the worst of the three scenarios.
The simulation-based optimisation method consists of two stages. The first stage is an iterative process where the route of each vessel (S v ) is first determined by the LNS. Here, the stochastic parameters τ v ,τ i andτ i are transformed into their deterministic counterparts. In the first iteration, these parameters are based on their mean values, whereas in the remaining iterations, the parameters are based on the θ%-quantile data generated from Monte Carlo simulation as shown in Algorithm 4. It is worth noting that the value of θ is adjusted systematically, where θ = (α − 1) · 10 ·T /100.
Once the route of each vessel (S v ) has been determined, the Monte Carlo simulation is performed to accommodate the uncertainty of parameters τ v ,τ i andτ i . The expected total maintenance cost (z s ) is obtained by running the simulationT times which can be considered as short simulation (e.g. T = 10, 000). This procedure is repeated Γ times and the route of each vessel (S * v ) that gives the smallest expected total maintenance cost z s * is then selected as an input to the next stage. In the Set θ = (α − 1) · 10 ·T /100 and use the θ%-quantile data generated from Monte Carlo simulation in the previous iteration (iteration [α − 1]) to estimate parameters τ v ,τ i andτ i 9: end if (τ i ) based on their distribution.

10:
Solve the MRP using the LNS with deterministic parameters. The route for each vessel is then obtained. Let S v (v ∈ V ) be a set of turbines that will be visited by vessel v.

11:
Run Monte Carlo simulationT times (short simulation) using S v (the route of each vessel obtained by previous step) by calling Algorithm 4 -Monte Carlo Simulation (S v ,T , γ, z s ). Record the expected total cost value (z s ) based on the γ%-quantile total maintenance cost data.

12:
if z s < z s * then 13: 14: end if 15: end for 16: Stage 2: 17: Run Monte Carlo simulationT times (long simulation) using S * v obtained from the previous stage by calling Algorithm 4 -Monte Carlo Simulation. Record the final expected total maintenance cost value z s * and take S * v as the route for each vessel in the final solution.
second stage, Monte Carlo simulation (e.g.T = 100, 000) is performed based on (S * v ) to obtain the final expected total maintenance cost z s * . As this method requires the solution of the deterministic problem iteratively (several times), a powerful optimiser that can solve such a problem within a relatively short time while producing good quality solutions is needed. Therefore, we propose a metaheuristic that will be embedded in the proposed simulation-based optimisation method.
The procedure of Monte Carlo for the SMRP is presented in Algorithm 4 which requires the route of each vessel (S v ) as input. This route is determined by the LNS. The simulation is executedT times, where in each run parameters τ v ,τ i andτ i are generated randomly based on their distribution.
As the route is fixed, the number of technicians required (Q vip ), the drop/pick-up time (T vi ) and the return time (T v(2n+1) ) for each vessel v ∈ V can be determined based on the generated values of τ v ,τ i andτ i . In addition, the total maintenance cost can be calculated taking into account the vessel penalty cost for the vessels that return to O&M base beyond the recommended time (t v ).
This vessel penalty cost can be determined by introducing an additional parameter calledc which represents the vessel penalty cost per hour for the late vessel. OnceT have been executed, the expected total maintenance cost is taken from the γ%-quantile total maintenance cost (z k ).

Algorithm 4 The Proposed Monte Carlo Simulation
Require: S v ,T , γ, z s 1: for k = 1 toT do 2: Generate randomly the travel time of each vessel (τ v ), the maintenance time for each turbine (τ i ) and the transfer time for technicians and equipment to a turbine (τ i ) based on their distribution.

3:
Using the fixed route of each vessel (S v ), determine the number of technicians required (Q vip ), the drop/pick-up time (T vi ) and the return time (T v(2n+1) )

4:
Compute the total maintenance cost of the k th iteration considering the vessel penalty cost for the late vessels returning to O&M base. 5: end for 6: Use the γ%-quantile total maintenance cost data (z k ) as the expected total maintenance cost value (z s ).

The simulation-based optimisation considering unexpected broken-down turbines
An additional consideration when maintenance route planning is that a turbine could suffer an unexpected break down shortly before or during the execution of the planned route. A breakdown here is defined as an event sufficiently severe as to stop the turbine from producing any power.
This type of breakdown occurs according to a continuous statistical distribution, that may or may not be known in advance and is governed by a failure rate parameter as well as information from the maintenance history of the turbine. With growing data and hence knowledge of wind turbine breakdowns, these failure rates are becoming more predictable.
There are several strategies to deal with turbines which break down or are likely to breakdown close to the planned route of the vessel. Sampling, which is one of the common strategies, incorporates stochastic knowledge by generating scenarios based on realisations drawn from random variable distributions (Pillac et al., 2013). The fundamental idea is to include the sampled turbines which have high chances to breakdown to the existing set of turbines to be visited. Therefore, we may end up with three scenarios: (i) A-priori optimal route which ignores turbines likely to breakdown; (ii) Tour with sampled turbines selected based on their failure rate and their maintenance history; (iii) Optimal scenario without sampled turbines which is sub-optimal regarding a cost evaluation, but leaves room to accommodate new turbines at a lower cost Figure 3 illustrates how scenarios are generated where based on the current turbines (Figure 3a), tour with sampled turbines, i.e. turbine T5, (Figure 3b), and a tour without the sampled turbine T5 (Figure 3c).

Figure 3: Scenario generation in sampling approaches
In this study, we utilise Monte Carlo simulation to predict the turbine breakdowns based on the failure rate and the maintenance history of the offshore turbine. We consider the failure rate per turbine per year (λ) together with the failure rate per turbine components per year (λ c , c ∈ C), where set C denotes the components in an offshore turbine including the gearbox, transformer, blades, tower. The number of technicians, the repair time and the weight of spare parts required to maintain the turbine are determined based on the component that needs to be repaired. Let J be the set of turbines in an offshore wind farm whereas set I denotes the set of turbines that will be maintained (PM and CM) where these turbines are determined at the beginning. Let J = J − I be the turbines that will not be visited. However, there is a probability that these turbines might breakdown, where the unexpected broken down turbines need corrective maintenance straight away in order to reduce cost and maintain power output.
The main procedure of Monte Carlo simulation to predict the turbine breakdowns is presented in Algorithm 5. In the first step, based on the maintenance history, the last maintenance date for each turbine j ∈ J is retrieved where t j , j ∈ J indicates that turbine j was maintained t j days ago. Assuming that the failure rate (λ) remains constant over the life of a turbine, the probability of failure before a certain time can be approximated using Exponential Distribution. This means that the reliability for each turbine (R j ) can be calculated as follows: (R j = e −λt j ). The simulation is executedT times, and in each run a turbine is randomly determined whether it will breakdown or not based on its reliability. If it is the turbine is broken down then the component that needs to be repaired is also randomly generated using the roulette wheel selection principle based on the failure rate of the components. Parametersτ ji , w ji , and ρ jp are determined based on the selected component. Then, we select β turbines with the highest number of failures (η j ) to be included as the sampled turbines (Ĵ).

7:
if υ > R j then 8: Select randomly the type of component that need to be repaired using the roulette wheel selection principle based on the failure rate of the components.

9:
Determine the values ofτ ji , w ji , and ρ jp based on the selected component.

end if
12: end for 13: end for 14: Calculateτ j = T i=1τ ji /η j , w j = T i=1 w ji /η j , and ρ j = T i=1 ρ jp /η j 15: LetĴ be the top β turbines that have the highest value of η j whereĴ ∈ J

Computational experiments for the stochastic maintenance routing problem
Extensive computational experiments have been carried out to evaluate the performance of the proposed stochastic model for the maintenance routing problem under uncertainty (SMRP). The set of data used for the SMRP is similar to the one used for the MRP in Section 4. However, the data has been adapted in order to be suitable for the SMRP. Here, the experiments on the stochastic model are performed on the West Gabbard case study only. The specification of the vessels used in the model is based on Table 1. The travel time per km of each vessel follows a normal distribution with mean µ = 1.71 and standard deviation σ = 0.7. We also assume that the required transfer time for technicians and equipment from a vessel to a turbine follows a normal distribution with mean 11 minutes and standard deviation 4 minutes , N (11, 4). In the Monte Carlo simulation, the vessel penalty cost (c) is set to e650 per hour to penalise a vessel that returns after the weather/time window (12 hours) which is almost the same as the downtime cost per hour. The example of the turbines data is similar to the ones in Table 2 where the maintenance time is also assumed to follow a normal distribution. The mean repair time is set to the values of the third column of Table 2 whereas the standard deviation is set to 2 and 3 hours for PM and CM activities respectively.

Computational Results for SMRP
In this study, three scenarios have been considered by varying the value of θ to 50%, 70% and 90%. In other words, the experiments have been performed using the 50%, 70% and 90% quantiles to obtain the expected total maintenance cost. The experimental results are presented in Table   5, which is divided into three sections where the results with 50%, 70% and 90% quantiles (θ) are provided. The table presents the breakdown cost together with turbine and vessel penalty costs which refer to penalty costs for not maintaining a turbine and for not returning to the O&M base within the weather window, respectively. The table also shows the CPU time required to solve each instance, which increases with the number of turbines to maintain. The table also reveals that the total maintenance cost increases with the number of turbines. The transportation cost is relatively small when compared to the total downtime cost. When a higher value of θ is used, the total penalty cost increases. This is reasonable as the value of θ affects the expected maintenance duration and travel time. When the value of θ is high, the maintenance duration and travel times are long. also provided. It can be noted that diff (%) increases when the value of θ increases. Here, 50%, 70% and 90% quantiles yield approximately a difference of 3.5%, 7.9% and 15.6% respectively. The total maintenance costs produced by the proposed stochastic model are slightly higher than the one produced by the deterministic MRP model. However, the solutions on the stochastic model are more realistic as uncertain conditions are considered.  Table 7 shows the different routes for |V | = 3 and |I| = 9 generated by solving the deterministic MRP and stochastic MRP with θ= 50%, 70% and 90%. In the table, the text in red indicates the pick-up nodes as the the vessels will pick up technicians from the turbines. The table also presents the number of technicians in the vessels after leaving the nodes (turbines). According to the solutions obtained, all three vessels are used to maintain the turbines. The solutions generated also prioritise the broken down turbines (Turbines T8 and T9) where these turbines are visited first in order to minimise the total cost. The corrective maintenance (CM) task is needed for these turbines as they are not operating. Surprisingly, in the solution generated for SMRP with θ= 50%,
In the Monte Carlo Simulation of Algorithm 5, the value ofT is set to 10,000 simulation runs.
The value of β is set to min{(|V | − 1), 2} where β turbines with the highest probability value is treated as the sampled turbines. Figure 4 presents the pareto chart of the probability of turbines that will break down generated for the instance with |V | = 3 and |I| = 9. In this case, Turbines T5 and T6 are considered as the sampled turbines and they are predicted to breakdown on that day. Table 9 reveals the experimental results for the SMRP considering the unexpected turbine breakdowns. The table provides the total maintenance costs for the a-priori route, route with The total maintenance cost is still 5.6% higher than the one of the a priori route although the sampled turbines have been removed from the routes. However, the routes based on the optimised scenario without sampled turbines are more robust as they consider the turbines that have a high chance of breakdown. We also observe quite interesting results when considering the turbine breakdowns. Figure 5 shows the routes generated for the instance with |V | = 2 and |I = 8|. As illustrated on Figure   5a, in the a priori routes, all turbines can be maintained and visited by 2 vessels. However, once turbine T9 is considered as a sampled turbine that has high potential to breakdown, and turbine T6 that needs PM cannot be visited by Vessel 1 due to the lack of resources. This leads to the  Figure 5c are more robust and flexible to accommodate turbine T9 that may breakdown.

Sensitivity Analysis
A sensitivity analysis on the vessel penalty cost for late return vessels is carried out to provide useful information to decision makers. We evaluate the effect of the vessel penalty cost on the change in the total lateness (hours) and the routes generated for each vessel. In the previous experiments, a vessel penalty cost of e650 per hour (c) is used which is based on the downtime cost. We vary this cost to e350 and e50 per hour. The former is determined based on the average personnel cost (for 12 technicians) per hour, whereas the latter indicates a relatively small vessel penalty cost.
The experiments are carried out on three problems. The first one is the problem with |V | = 2 and |I| = 7 using θ = 50%, the second one |V | = 3 and |I| = 9 using θ = 70%, and the last one |V | = 4 and |I| = 14 using θ = 90%. Figure 6 shows the total lateness for a vessel penalty cost (c) for a different problem. According to the figure, when a lower vessel penalty cost is used, a solution with a larger total lateness is generated. Whenc is set smaller than the travel cost of a vessel per hour, a turbine tends to be visited by a vessel that has visited a turbine near to it at expense of a longer return time of this vessel. This may lead to a vessel penalty cost for the lateness for a vessel to return to O&M base.
However, this cost is still smaller than the travel cost for other vessel to visit that turbine. In other words, changing the vessel penalty cost per hour affects the turbine allocation configuration. For example, in the problem with |V | = 3 and |I| = 9 using θ = 70%, Vessels 1 and 2 visit 5 and 2 turbines respectively whenc= 50. If we increase the value ofc to e650, it changes the turbine allocation configuration where Vessels 1 and 2 visit 4 and 3 turbines, respectively. According to the results of these experiments, it is more realistic to set the value ofc to e650 as a balance turbine allocation to the vessels can be obtained while yielding the smallest lateness. We also perform a sensitivity analysis on the maintenance/repair time where different statistical distributions are used to approach this stochastic parameter. In the previous experiments, the normal distribution was used for the maintenance time. Here, we also apply exponential and lognormal distributions as these can also be used for this parameter in an offshore wind farm (Seyr and Muskulus, 2016). The total maintenance cost is assessed by the change in the maintenance/repair time generated based on its distribution. The experiments are conducted on three problems using θ = 50% with |V | = 2 and |I| = 6, |V | = 3 and |I| = 10, and |V | = 4 and |I| = 12. The results of our experiments are presented in Figure 7, where the total maintenance cost based on the different distributions for the maintenance time is presented. The figure shows that for each problem, the use of normal distribution generates the lowest total maintenance cost followed by exponential and lognormal distributions. This indicates that a longer maintenance time will be generated when exponential or lognormal distributions are assumed. There are a number of possible extensions of the models developed in this paper to make it more applicable to offshore wind farms in operation. The models can be integrated with the O&M Strategy model for a more holistic decision support framework considering strategic, tactical and operational time scales. Further research into the choice between the three scenario generation techniques illustrated by Figure 3 under different real world circumstances would also be an interesting line of potential future research. The models can also be enhanced for the maintenance scheduling and routing of Service Operation Vessels (SOV) or other "mother vessels" which can stay offshore for multiple days. One challenge is coordinating the operation of such vessels with the use of daughter vessels, ordinary Crew Transfer Vessels and possibly also helicopters.