Solving a wind turbine maintenance scheduling problem

Driven by climate change mitigation efforts, the wind energy industry has significantly increased in recent years. In this context, it is essential to make its exploitation cost-effective. Maintenance of wind turbines therefore plays an essential role in reducing breakdowns and ensuring high productivity levels. In this paper, we discuss a challenging maintenance scheduling problem rising in the onshore wind power industry. While the research in the field primarily focuses on condition-based maintenance strategies, we aim to address the problem on a short-term horizon considering the wind speed forecast and a fine-grained resource management. The objective is to find a maintenance plan that maximizes the revenue from the electricity production of the turbines while taking into account multiple task execution modes and task-technician assignment constraints. To solve this problem, we propose a constraint programming-based large neighborhood search (CPLNS) approach. We also propose two integer linear programming formulations that we solve using a commercial solver. We report results on randomly generated instances built with input from wind forecasting and maintenance scheduling software companies. The CPLNS shows an average gap of 1.2% with respect to the optimal solutions if known, or to the best upper bounds otherwise. These computational results demonstrate the overall efficiency of the proposed metaheuristic.

1. The Global Wind Energy Council—Global Wind Report Annual Market Update 2015—http://www.gwec.net/wp-content/uploads/vip/GWEC-Global-Wind-2015-Report_April-2016_22_04.pdf, last accessed: 2016-09-15.

2. The ratio of the duration that a generating unit is available to provide energy to the grid, for the time considered, to the duration of the same time period.

3. A wind farm substation collects the electricity produced by all the turbines of the farm and distributes it through the grid.

4. We also implemented the adaptive layer as proposed in Ropke and Pisinger (2006), but after some preliminary experimentation we concluded that the contribution of this component to the accuracy of the method did not payoff the loss of simplicity and the effort needed to fine tune the additional parameters. We therefore limit the discussion in the paper to the basic LNS version.

5. Average of the mean gap found for each instance over 3 runs.

6. Average of the mean percentage of tasks scheduled in the solution found for each instance over 3 runs.

7. Average of the mean gap found for each instance over 10 runs.

8. Average of the mean percentage of tasks scheduled in the solution found for each instance over 10 runs.

9. Average of the best (minimal) gap found for each instance over 10 runs.

10. Average of the worst (maximal) gap found for each instance over 10 runs.

The authors would like to kindly thank two anonymous reviewers for their comments and suggestions. This work was supported by Angers Loire Métropole through its research grant program; and by the Natural Sciences and Engineering Research Council of Canada (NSERC) through a grant that enabled the collaboration with the Canadian company WPred, which we would like to thank for their expertise.

Authors and Affiliations

1. Université Bretagne Loire, Université Catholique de l’Ouest, LARIS EA 7315, Angers, France

   Aurélien Froger & Eric Pinson

2. CIRRELT Interuniversity Research Centre on Enterprise Networks, Logistics and Transportation, and Département de mathématiques et de génie industriel, Polytechnique Montréal, Montreal, Canada

   Michel Gendreau & Louis-Martin Rousseau

3. CNRS, LI EA 6300, ROOT ERL CNRS 6305, Université François-Rabelais de Tours, Tours, France

   Jorge E. Mendoza

4. Centre de Recherches Mathématiques (UMI 3457 CNRS), Montreal, Canada

   Jorge E. Mendoza

Corresponding author

Correspondence to Aurélien Froger.

Appendix 1: Notations

• Time

  • \(\mathcal {T}\) : time horizon

  • \(\mathcal {D}\): set of days

  • \(\mathcal {T}_d\): set of time periods that belong to day \(d\in \mathcal {D}\)

  • \(t^{rest}_d\): last time period \(t\in \mathcal {T}\) before the rest time period following day \(d\in \mathcal {D}\)

• Locations

  • \(\mathcal {L}\): set of locations (wind farms and technician home depots)

  • \(\delta _{ll'}\): binary parameter equal to 1 if and only if locations l and \(l'\) are compatible

• Technicians

  • \(\mathcal {S}\): set of skills

  • \(\mathcal {R}\): set of technicians

  • \(\lambda _{rs}\): binary parameter equal to 1 if and only if technician \(r\in \mathcal {R}\) masters skill \(s\in \mathcal {S}\)

  • \(\pi ^{t}_{r}\): binary parameter equal to 1 if and only if technician \(r\in \mathcal {R}\) is available during time period \(t\in \mathcal {T}\)

  • \(l^t_r\): location of technician \(r\in \mathcal {R}\) when he or she is not available during time period \(t\in \lbrace t'\in \mathcal {T}, \pi ^{t'}_{r}=1\rbrace \)

• Tasks

  • \(\mathcal {I}\): set of tasks to be performed on the turbines

  • \(ov(\mathcal {I})\): set of subsets of non-overlapping tasks

  • \(l_i\): location where task \(i\in \mathcal {I}\) has to be performed

  • \(\mathcal {M}_i\): set of execution modes for task \(i\in \mathcal {I}\)

  • \(m^0_i\): execution mode related to the postponement of task \(i\in \mathcal {I}\) (\(m^0_i\in \mathcal {M}_i\))

  • \(q_{im}\): number of technicians required during each time period to perform task \(i\in \mathcal {I}\) in mode \(m\in \mathcal {M}_i\) (\(q_{im^0_i}=0\))

  • \(d_{im}\): duration of task \(i\in \mathcal {I}\) if performed in mode \(m\in \mathcal {M}_i\) (\(d_{im^0_i}=0\))

  • \(\gamma ^{t}_{i}\): binary parameter equal to 1 if and only if task \(i\in \mathcal {I}\) can be executed during time period \(t\in \mathcal {T}\)

  • \(s_i\): skill required to perform task \(i\in \mathcal {I}\)

  • \(o_{i}\): penalty if task \(i\in \mathcal {I}\) is postponed.

• Turbines

  • \(\mathcal {W}\): set of turbines

  • \(b_{wi}\): binary parameter equal to 1 if and only if the execution of task \(i\in \mathcal {I}\) shuts down turbine \(w\in \mathcal {W}\) when technicians are effectively working on i

  • \(\widetilde{b}_{wi}\): binary parameter equal to 1 if and only if the execution of task \(i\in \mathcal {I}\) shuts down turbine \(w\in \mathcal {W}\) during the rest time periods it overlaps

  • \(g^{t}_{w}\): revenue if turbine \(w\in \mathcal {W}\) can produce electricity during time period \(t\in \mathcal {T}\)

  • \(\widetilde{g}^{d}_{w}\): revenue if turbine \(w\in \mathcal {W}\) can produce electricity during the rest time period following day \(d\in D\)

• Plans

  • \(\mathcal {P}\): set of plans

  • \(\mathcal {P}_i\): set of plans involving task \(i\in \mathcal {I}\)

  • \(i_p\): task involved in plan \(p\in \mathcal {P}\)

  • \(a^{t}_{p}\): binary parameter equal to 1 if and only if task \(i_p\) is executed during time period \(t\in \mathcal {T}\)

  • \(q_p\): number of required technicians if plan \(p\in \mathcal {P}\) is selected

  • \(\mathcal {R}_p = \mathcal {R}_{i_p}\)

  • \(b_{wp} = b_{wi_p}\)

  • \(\widetilde{b}_{wp} = \widetilde{b}_{wi_p}\)

  • \(o_{p} = o_{i_p}\)

Appendix 2: Detailed computational results

1.1 Testbed G1

See Tables 10 and 11.

1.2 Testbed G2

Table 12 reports the results obtained when solving the ILP and CP models and when running the CPLNS. To assess the quality of the results, it is noteworthy that large values for the gap are essentially associated with Type A instances. Because the technician-to-work ratio is tighter, it is more difficult to schedule all the tasks. Scheduling one additional task would drastically reduce this gap.

From the results, we do not see any significant impact of adding corrective tasks on the difficulty to solve the instances. It is not so surprising as corrective tasks usually need to be scheduled at the beginning of the planning horizon in order not to induce too much lost revenue. For our models, this might implies less symmetry on the scheduling of the tasks. Moreover, we can draw identical conclusions on the efficiency of the approaches. Solving the ILP formulations does not lead to an efficient exact approach as only the small-sized instances can be optimally solved, whereas solving the CP formulation yields good results for Type B instances. The CPLNS gives the overall best results and near-optimal solutions for small and medium-sized instances. Lastly, the difficulty of the instances seems also to be related to the same characteristic as for testbed G1.

For each approach, we also report in Table 13 the gap to the best solution that we were able to compute in our tests. We find similar results as for testbed G1. Globally, the CPLNS outperforms the other approaches.

Appendix 3: Instance generation

An instance of the problem is primarily characterized by:

• a finite time horizon (a finite number of time periods)

• a number of time periods per day (yielding the number of days)

• a set of locations (wind farms + home depot)

• a set of wind turbines distributed over the wind farms

• a set of maintenance tasks to perform at the different locations and that impact the availability of the turbines

• a set of technicians to perform the tasks

• wind speed for each time period and location

• postponing penalties

The generator is based on the following parameters:

• \(n_{\mathcal {T}}\), \(n_{\mathcal {D}}\), \( n_{\mathcal {I}}\), \(n_{\mathcal {S}}\) (length of time horizon, number of days, number of wind farms, number of tasks and number of skills)

• \(Dn^{\mathcal {L}}\): probability distribution of the number of locations

• \(Dl^{xy}\): probability distribution of the coordinates associated with each location

• \(Dn^{\mathcal {L}}_{\mathcal {W}}\): probability distribution of the number of turbines per location

• \(Dn^{\mathcal {W}}_{\mathcal {I}}\): probability distribution of the number of tasks per turbine

• \(\varDelta l^{min}\): minimum distance between two locations

• \(\varDelta r^{max}\): maximum distance between two locations such that they can be visited by the a technician during the same day

• K: set of all types of preventive tasks that we consider

• p(k): probability of generating a task of type \(k\in K\)

• \(Di_{impact}(k)\): probability distribution of the impact of each type of preventive task on the wind turbines

• \(Di_{dur}(k)\): probability distribution of the duration of each type of preventive task

• \(Di_{req}(k)\): probability distribution of the number of technicians that can perform each type of preventive task during any time period

• \(Dr_{\#skills}\): probability distribution of the number of skills mastered by a technician

• \(Dr_{\mathbbm {P}(unv)}\): probability that a technician has some unavailability time periods during the time horizon

• \(Dr_{\#unv}\): probability distribution of the number of time periods during which a technician is unavailable

• \(Dr_{d^{unv}}\): probability distribution of the duration of the unavailability of a technician (in man-hours)

• \(Dw_{power}\): probability distribution of the nominal power (in kW) of each turbine

• \(\widehat{\phi }\): average wind speed on each wind farm

• \(\varUpsilon ^{safety}_{max}\): maximum wind speed allowed to perform a task

• \(\varDelta l^{max}\): maximum distances for the spatial correlation of the wind speed

• \(\delta \): number of values used in the moving average for the timewise dependency between the wind speeds

• \(\alpha \): correlation factor between wind speed

We generate an instance following multiple steps. First of all, the length of time horizon, the number of days, the number of tasks and the number of skills are input values. This yields directly the set \(\mathcal {T}\) of time periods and the set \(\mathcal {D}\) of days.

We then start the generation of an instance by building the set \(\mathcal {L}\) of locations whose cardinality is set by sampling the \(Dn^{\mathcal {L}}\) distribution. According to the distance \(\varDelta l^{min}\), we then generate the coordinates of each location by sampling the \(Dl^{xy}\) distribution. Based on these coordinates and on the distance \(\varDelta r^{max}\), we compute the parameters \(\left( \delta _{ll'}\right) _{(l,l')\in \mathcal {L}^2}\) that enable to define the daily location-based incompatibility constraints.

Afterward, we built the set \(\mathcal {W}\) of wind turbines. To this end, according to the target number of tasks, we start by generating a number of wind turbines per locations by sampling the \(Dn^{\mathcal {L}}_{\mathcal {W}}\) distribution. For each location where there is at least one wind turbine (i.e., this location is a wind farm), we then generate a nominal power by sampling the \(Dw_{power}\) distribution and we set the nominal power \(P_w\) equal to this latter value for each wind turbine \(w\in \mathcal {W}\) of the wind farm.

After that, we call procedure \(\mathbf {genTasks()}\) to create the set \(\mathcal {I}\) of tasks. Notice that for each task \(i\in \mathcal {I}\) we build the set \(\mathcal {M}_i\) of execution modes such that it meets the two following requirements:

• \(\forall m,m'\in \mathcal {M}_i\), \(q_{im}\ne q_{im'},\)

• \(\forall m,m'\in \mathcal {M}_i\), \(q_{im}< q_{im'} \rightarrow d_{im}> d_{im'}\).

Arbitrarily, we build \(ov(\mathcal {I})\) considering that overlapping tasks are forbidden on the same turbine. Notice that, according to some experts in the field, it is reasonably realistic to only consider these subsets. After the generation of the tasks, we generate the set \(\mathcal {R}\) of technicians using procedure \(\mathbf {genTechnicians()}\).

The last part of the generator concerns the parameters related to the objective function. For the sake of convenience, we introduce the set \(\mathcal {T}^{+}\) of all time periods formed by the union of set \(\mathcal {T}\) and the set of rest time periods that occur between each day. More specifically, we include a rest time period after every \(\dfrac{\vert \mathcal {T}\vert }{\vert \mathcal {D}\vert }\) consecutive time periods of \(\mathcal {T}\).

As it concerns the wind speed forecast at hub height, the main purpose is to use realistic values. First, we generate wind speed \(\bar{\phi }^t_l\) for every location \(l\in \mathcal {L}\) and every time period \(t\in \mathcal {T}^{+}\) using a Rayleigh distribution with a scale parameter equal to \(\widehat{\phi }\sqrt{\dfrac{2}{\pi }}\) (so that the expected wind speed is \(\widehat{\phi }\)). Since space correlation can be significant, we compute a corrected wind speed \(\bar{\bar{\phi }}^t_l\) for every location l and every time period t as follows:

Wind speeds were generated independently from a time period to another one. However, this timewise independence assumption is unlikely to be verified in practice. To smooth out the speed values, we use a \(\delta \)-weighted moving average that yields wind speed \(\phi ^t_l\) according to the following formula:

The resulting values are rounded to the nearest tenth. From our perspective, they compare well to realistic data.

Afterward, for each task \(i\in \mathcal {I}\) and every time period \(t\in \mathcal {T}\), we compute the binary parameter \(\widetilde{\vartheta }^{t}_{i}\) equal to 1 if and only if \(\phi ^t_l<\varUpsilon ^{safety}_{max}\) (i.e., the task i can be scheduled during time period t according to safety concerns). Arbitrarily, we set each parameter \(\vartheta ^{t}_{i}\) equal to 1 for every task i and every time period t. We point out here that this choice makes the instances more complicated to solve as there is a wide flexibility to schedule the maintenance operations. This also matches field observations.

The last step consists in computing the revenue value \(g^t_w\) for every wind turbine \(w\in \mathcal {W}\) during each time period \(t\in \mathcal {T}^{+}\). We compute the revenue from the nominal power \(P_{w}\) of the wind turbine and from the wind speed \(\phi ^{t}_{l_w}\). We also use an estimation hours(t) of the number of hours during every time period t. More specifically, we compute the revenue \(g^t_w\) generated by each turbine \(w\in \mathcal {W}\) that is available during time period \(t\in \mathcal {T}\) as follows:

where

• 0.08: is an approximation to the selling price in euros of 1 kWh of wind energy (this selling price is guaranteed for the next 10 years in France).

• hours(t): estimation of the number of hours during time period \(t\in \mathcal {T}^{+}\)

• \(P_w\): nominal power of wind turbine \(w\in \mathcal {W}\)

• \(\phi ^t_{l_w}\): wind speed forecast during time period t at the location of turbine \(w\in \mathcal {W}\)

• \(CF(\phi )\): the ratio of the net electricity generated according to a wind speed equal to \(\phi \) to the energy that could have been generated at full-power operation (this ratio is given by a piecewise linear function estimated from real data)

Finally, we compute a single postponing penalty across all tasks. This penalty is equal to the maximum loss of revenue that can be generated by the scheduling of a task of \(\mathcal {I}\) plus one. With this definition, we almost always (if not always) ensure that postponing a task is non-profitable. With this penalty, we therefore almost ensure to schedule the maximum number of tasks according to the total number of technicians and their availability. This is quite in line with the practice in the field.

Table 14 presents the detail parameter setting used in the generation process.

Reprints and permissions