Michel Vasquez
email: michel.vasquez@mines-ales.fr

Mirsad Buljubasic
email: mirsad.buljubasic7@gmail.com

Saïd Hanafi
email: said.hanafi@uphf.fr

B Michel Vasquez

An efficient scenario penalization matheuristic for a stochastic scheduling problem

Keywords: Mixed integer programming, Local search, Stochastic scheduling, Mean-risk, Quantiles

We propose a new scenario penalization matheuristic for a stochastic scheduling prob-lem based on both mathematical programming models and local search methods. The application considered is an NP-hard problem expressed as a risk minimization model involving quantiles related to value at risk which is formulated as a non-linear binary optimization problem with linear constraints. The proposed matheuritic involves a parameterization of the objective function that is progressively modified to generate feasible solutions which are improved by local search procedure. This matheuristic is related to the ghost image process approach by Glover (Comput Oper Res 21(8):801-822, 1994) which is a highly general framework for heuristic search optimization. This approach won the first prize in the senior category of the EURO/ROADEF 2020 chal-lenge. Experimental results are presented which demonstrate the effectiveness of our approach on large instances provided by the French electricity transmission network RTE.

Introduction

The practice of decision making under uncertainty frequently resorts to two criteria mean-risk models (cf. [START_REF] Markowitz | Mean-variance Analysis in Portfolio Choice and Capital Markets[END_REF] where the mean represents the expected outcome, and the risk corresponds to a scalar measure of the variability of outcomes, in which the mean is maximized and the risk is minimized. The standard mean-variance model, called the Markowitz model, uses the variance as the risk measure in the mean-risk analysis. The mean-variance model under a set of linear constraints can be formulated as a quadratic programming problem [START_REF] Markowitz | Portfolio selection[END_REF]. [START_REF] Sharpe | A linear programming approximation for the general portfolio analysis problem[END_REF] proposed a linear programming (LP) approximation to the mean-variance model for the general portfolio analysis problem. Several authors (cf. [START_REF] Ogryczak | Dual stochastic dominance and related mean-risk models[END_REF][START_REF] Ogryczak | Dual stochastic dominance and quantile risk measures[END_REF] have pointed out that in general the mean-variance model is not consistent with stochastic dominance rules [START_REF] Whitmore | Stochastic Dominance: An Approach to Decision-making Under Risk[END_REF]. To overcome this flaw of the mean-variance model, we consider a mean-risk model involving q-quantiles in risk measure. Note that the measure of value at risk (VAR) is defined as the maximum loss at a level q, consequently, VAR is a widely used quantile risk measure (cf. Jorion 1997, and references therein).

We consider the stochastic scheduling problem (SSP) proposed by Réseau de Transport d'Electricité, usually known as RTE. RTEis the electricity transmission system operator of France, responsible for the operation, maintenance and development of the French high-voltage transmission system, which with approximately 100,000 kms, is the largest of Europe. Maintenance policies have a major economic impact in many areas of the industry such as the electricity sector [START_REF] Froger | Maintenance scheduling in the electricity industry: a literature review[END_REF], the manufacturing industry and civil engineering. The main goal is to build an optimal maintenance schedule to ensure the delivery and supply of electricity. This problem corresponds to the competition of ROADEF/EURO challenge 2020 (ROADEF 2020).

The stochastic scheduling problem consists in determining the start time of maintenance activities (called interventions) in a high-voltage transmission network over a given time horizon. Each intervention needs a certain number of time units to be achieved without interruption that depends on its starting time. All interventions must be planned and finished before the end of the time horizon. A feasible plan must be consistent with all activity related restrictions such as resource constraints (e.g. each intervention consumes resources and the total amount of resources used at each time step is bounded from below and above) and exclusions between interventions (e.g. some interventions cannot take place at the same time). Given a set of scenarios, where the risk value of each intervention is known at each time step and for each scenario, the goal is to minimize a convex combination of the expectation and the quantile of the risk.

The stochastic scheduling problem can be formulated as a non-linear Mixed Integer Programming (MIP) model with a non-linear (generally non convex) objective function and linear constraints. Note that SSP optimization problem is similar to mean-VARportfolio problem. [START_REF] Benati | A mixed integer linear programming formulation of the optimal mean/value-at-risk portfolio problem[END_REF] showed that mean-VARportfolio problem is NP-hard and proposed MIP formulation using Cplex to solve medium size instances. It is well known that in general the exact methods require exponential growth of computational effort for large scale instances. These practical limitations have occasioned a considerable research effort focusing on approximation approaches such as heuristics, metaheuristics and matheuristics (cf. [START_REF] Hashimoto | A grasp-based approach for technicians and interventions scheduling for telecommunications[END_REF][START_REF] Gavranović | An efficient local search with noising strategy for google machine reassignment problem[END_REF][START_REF] Buljubašić | Two-phase heuristic for SNCF rolling stock problem[END_REF][START_REF] Hanafi | Mathematical programming based heuristics for the 0-1 MIP: a survey[END_REF]. [START_REF] Gouvine | Mixed-integer programming for the ROADEF/EURO 2020 challenge[END_REF] proposed a hybrid approach combining a branch and cut algorithm [START_REF] Padberg | A branch-and-cut algorithm for the resolution of large-scale symmetric traveling salesman problems[END_REF] with a constraint generation method based on Benders decomposition [START_REF] Rahmaniani | The benders decomposition algorithm: a literature review[END_REF]. [START_REF] Cattaruzza | A slope scaling/Lagrangean perturbation heuristic with long-term memory for multicommodity capacitated fixed-charge network design[END_REF] presented a finite convergent adaptive scenario clustering algorithm that guarantees an optimal solution but is only useful for problems of small size. They also developed an overlapping alternating direction method [START_REF] Glowinski | Sur l'approximation, par éléments finis d'ordre un, et la résolution, par pénalisation-dualité d'une classe de problèmes de dirichlet non linéaires[END_REF][START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite element approximation[END_REF]) that serves as a primal heuristic for quickly computing feasible solutions of good quality for problems in the size range examined. [START_REF] Zholobova | An industry maintenance planning optimization problem using CMA-VNS and its variations[END_REF] developed a hybrid approach combining the Covariance Matrix Adaptation Evolution Strategy [START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF] and Variable Neighborhood Search [START_REF] Hansen | Variable neighborhood search: basics and variants[END_REF]) metaheuristics.

In this paper, we introduce a new matheuristic for the SSP optimization problem, called Scenario Penalization Matheuristic (SPM). The proposed matheuritic involves a parameterization of the objective function that is progressively modified to generate feasible solutions which are improved by a local search procedure and which is applicable to applications of practical sizes.

The SPM is based on the ghost image process (GIP) approach by Glover (1994) which is a highly general framework for heuristic search optimization. Note that GIP is a generalization of many relaxation methods (e.g. Lagrangian, surrogate and composite relaxations) and self-organizing neural networks of [START_REF] Kohonen | Self Organization and Associative Memory[END_REF] as applied to optimization problems. GIP has been applied to some interesting optimization problems. For example, [START_REF] Woodruff | Ghost image processing for minimum covariance determinants[END_REF] developed a GIP application to the problem of computing the minimum covariance determinant estimators. With regard to the parameterized objective function, the GIP approach is related to the concept of slope scaling applied within the context of fixed-charge networks (see e.g. [START_REF] Kim | A solution approach to the fixed charge network flow problem using a dynamic slope scaling procedure[END_REF][START_REF] Gendron | A Tabu search with slope scaling for the multicommodity capacitated location problem with balancing requirements[END_REF][START_REF] Cattaruzza | A slope scaling/Lagrangean perturbation heuristic with long-term memory for multicommodity capacitated fixed-charge network design[END_REF][START_REF] Glover | Parametric ghost image processes for fixed-charge problems: a study of transportation networks[END_REF][START_REF] Gendron | Matheuristics based on iterative linear programming and slope scaling for multicommodity capacitated fixed charge network design[END_REF]. This approach SPM won the first prize in the senior category of the EURO/ROADEF 2020 challenge. Experimental results are presented which demonstrate the effectiveness of our approach on large instances provided by RTE.

The remainder of this paper is organized as follows. A description of the stochastic scheduling problem is provided in Sect. 2 and the scenario penalization matheuristic development for the stochastic scheduling problem is presented in Sect. 3. Computational results obtained from the available set of instances, provided by RTE, are provided in Sect. 4. Finally, Sect. 5 summarizes our conclusions.

Stochastic scheduling problem

In this section, first we will describe the input parameters of the stochastic scheduling problem (SSP). We refer the reader to the challenge subject in ROADEF (2020) for a more thorough description of the maintenance of electricity transmission lines and the implications of the problem from business and environmental perspectives. Additional details of the global RTE strategy and risk management can be found at [START_REF] Crognier | Grid operation-based outage maintenance planning[END_REF]. Next we propose a mathematical programming formulation of this problem.

Input parameters of the problem

We consider a finite set of maintenance activities I (called interventions in the challenge application) to be scheduled over a discrete time horizon T . Interventions are not equal in terms of duration or resource requirements. Because of days off (weekends...), the duration of a given intervention is not fixed in time and depends on when it starts. Therefore, δ i,t ∈ N denotes the actual duration of intervention i ∈ I if it starts at time t ∈ T .

To carry out the different interventions, a workforce is necessary, split into teams (or resources) of various sizes and specific skills. Each team has different specific skills and can potentially be required for any task. The available resources are always limited and vary over the time horizon. Let J denote the set of resources. The amount of resource j ∈ J used at time t ∈ T by intervention i ∈ I starting at time t is given by a j,t i,t . So for every resource j ∈ J and every time step t ∈ T , the amount of resource consumed must be between the lower bound l j t and the upper bound u j t . When an intervention is being performed, the power lines involved must be disconnected, causing the electricity network to be weakened at this time. This implies a certain risk for RTE, which is highly linked to the grid operation: if another nearby site were to break down (due to extreme weather for example), the network may not be able to handle the electricity demand correctly. Even if such events have an extremely low probability of occurring, they must be taken into account in the schedule. In order to financially quantify these risks, RTE previously conducted simulations for various scenarios with different time steps. Let S t be the set of scenarios at time t ∈ T and S = ∪ t∈T S t be the set of all scenarios. The risk value depends on the intervention concerned and on the time period, as it is often much less risky to perform interventions in summer (when there is less demand on the electricity network) rather than in winter. So the risk value (expressed in Euros in the challenge application) is denoted by risk s,t i,t ∈ R for time period t ∈ T , scenario s ∈ S and intervention i ∈ I when i starts at time step t ∈ T . Further, some maintenance activities cannot be performed at the same time. This exclusion between interventions is given by a set E of triplets (i, i , t) such that interventions i and i cannot be both in process at time t.

Intervention preemption is not allowed, and each intervention must be terminated at time |T |. To more formally derive the optimization models and their associated solutions, we require some additional notation. Let T (i) ={t ∈ T : t +δ i,t ≤|T |+1} denote the set of feasible starting times of intervention i ∈ I , and let + (i, t) ={t ∈ T : t ≤ t < t + δ i,t } denote the set of times for which the intervention

T i is in process if it starts at time t. We further let T -(i, t) ={t ∈ T : t ≤ t < t + δ i,t }={t ∈ T : t ∈ T + (i, t )}
denote the set of starting times of intervention i for which the intervention is in process at time t, and let I(t) ⊆ I denote the set of interventions in process at time t ∈ T . For any real value v ∈ R,let v =min{v ∈ N : v ≥ v}.

Mathematical programming formulation

A schedule (or solution) consists of a list of starting times of interventions. Hence a solution will be represented by an integer vector σ ∈ N |I| where each component σ i represents the starting time of intervention i ∈ I. An alternative is to define a solution by a binary matrix x ∈ {0, 1} |I|×|T | such that x i,t = 1 if intervention i ∈ I starts at time step t ∈ T 0 otherwise for i ∈ I and t ∈ T . The compact integer representation σ of a feasible schedule will be used in heuristic procedures and the binary representation x will be exploited in the mathematical programming formulation.

Non linear objective function

The evaluation score of a feasible schedule depends only on the distribution of risks. Two criteria are taken into account: the average and the excess of the risk values. The excess is defined from the quantile value of the risk distribution. More formally, given a feasible schedule represented by the integer vector σ or the binary matrix x, the mean risk Mean(x) = Mean(σ ) and the expected excess risk Excess(x) = Excess(σ ) are evaluated as follows.

Mean cost The cumulative planning risk at t ∈ T for a scenario s ∈ S t , denoted by risk s,t (x) or risk s,t (σ ), is the sum of risks in scenario s over the in-process interventions at t:

risk s,t (x) = i∈I(t) t ∈T -(i,t) risk s,t i,t × x i,t
Then the mean cost overall planning risk is

Mean(x) = 1 |T | t∈T 1 |S t | s∈S t risk s,t (x)
Excess cost The planning quality also takes into account the cost variability. Computing the mean risk over all scenarios induces a loss of information, and critical scenarios inducing extremely high costs may not be adequately captured by the mean. To prevent this kind of outcome from happening, a metric exists to quantify the variability of the scenarios. The expected excess indicator relies on the τ quantile values where τ ∈]0, 1]. For every time period t, we define the quantile value Q t τ as follows:

Q t τ (x) = min{q ∈ R : ∃E ⊆ {1, . . . , |S t |} : |E|≥ τ × |S t | ; ∀s ∈ E, risk s,t (x) ≤ q}
Note that if we sort elements of the set S t in increasing order, then Q t τ (x) will be equal to the element at position τ × |S t | in this sorted set. The expected excess of a planning is: Two risk vectors: r 1 , r 2 and a linear combinaison of them:

Excess(x) = 1 |T | t∈T max ⎛ ⎝ 0, Q t τ (x) - 1 |S t | s∈S t risk s,t (x) ⎞ ⎠ Table 1 Example of non-convexity: 4 scenarios S = {s 1 , s 2 , s 3 , s 4 }, α = 0.5 and τ = 0.7 scenarios r 1 r 2 r 3 = 0.5 × (r 1 + r 2 ) s 1 1.000 1.
r 3 = 0.5 × (r 1 + r 2 ). Bold values correspond to the Q τ values. Obj(r 3 ) = 0.5 × (Obj(r 1 ) + Obj(r 2 ))
The two metrics Mean(x) and Excess(x) described above cannot necessarily be compared directly, as they depend on risk aversion (or risk policies). That is why a scaling factor α ∈ [0, 1] is needed. Then the final score of a planning is:

Objective(x) = α × Mean(x) + (1 -α) × Excess(x).
In general, the quantile function Q t τ (x) is non-convex as shown by [START_REF] Gouvine | Mixed-integer programming for the ROADEF/EURO 2020 challenge[END_REF]. Hence, the objective function Objective(x) is non-convex. Moreover, the objective function Obj(r ) of the vector of decision variables risk r associated with a solution x (i.e. r = (risk s,t (x)) t∈T ,s∈S t ) is a non-convex function as can be observed from Table 1.

The goal is to find a feasible schedule which minimizes the objective function Objective(x). A schedule x is feasible if satisfies all the linear constraints presented below.

Linear constraints

Schedule Interventions have to start at the beginning of a period. Moreover, as interventions require shutting down some lines of the electricity network, once an intervention starts, it cannot be interrupted. More precisely, if intervention i ∈ I starts at time t ∈ T , then it must end at t + δ i,t . All interventions must be executed and completed no later than the end of the horizon. If intervention i ∈ I starts at time t ∈ T , then t + δ i,t ≤ |T |+1. Hence, the following multiple choice constraints must be satisfied:

t∈T x i,t = 1, ∀i ∈ I (1)
From the definition of this schedule problem, we have the reduction

x i,t = 0, ∀i ∈ I, t ∈ T -T (i)
which can be incorporated in the preprocessing step.

Resource limitation

The resources needed cannot exceed the resource capacity but have to be at least equal to the minimum workload, and hence the resource constraints are:

l j t ≤ i∈I t ∈T a j,t i,t × x i,t ≤ u j t , ∀ j ∈ J , t ∈ T (2)
Exclusion between interventions The exclusion constraints can formally be written as:

i ∈ I(t) ⇒ i / ∈ I(t) ∀(i, i , t) ∈ E
This implication can formulated using the binary variables x i,t as follows:

t ∈T -(i,t) x i,t + t ∈T -(i ,t) x i ,t ≤ 1, ∀(i, i , t) ∈ E (3)

MIP_full formulation

Finally, the stochastic scheduling problem SSP can be expressed as a non-linear binary MIP model, denoted MIP_full, given by

(MIP_full) ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ min obj = α |T | t∈T 1 |S t | s∈S t risk s,t (x) + (1-α) |T | t∈T max(0, Q t τ (x) -1 |S t | s∈S t risk s,t (x)) s.t.: (1), (2) and (3), x ∈ {0, 1} |I|×|T |
The quantile value Q t τ (x) for each scenario s in time step t can be expressed by introducing a binary variable y s,t which indicates whether Q t τ (x) is greater than or equal to the total risk for scenario s at time step t, i.e.

y s,t = 1 if Q t τ (x) ≥ risk s,t (x) 0 otherwise , ∀t ∈ T , s ∈ S t (4)
with additional linear constraints requiring that the number of y s,t binary variables that take the value 1 must be greater than or equal to τ × |S t | , i.e.

s∈S t y s,t ≥ τ × |S t | , ∀t ∈ T (5)
The implicit definition (4) of auxilary binary variables y for quantile can be formulated as quadratic constraints

Q t τ (x) ≥ y s,t × risk s,t (x), ∀t ∈ T , s ∈ S
Linear MIP formulations are proposed by [START_REF] Gouvine | Mixed-integer programming for the ROADEF/EURO 2020 challenge[END_REF] and [START_REF] Cattaruzza | A slope scaling/Lagrangean perturbation heuristic with long-term memory for multicommodity capacitated fixed-charge network design[END_REF]. To strengthen the linear relaxation of the quantile model, they generate a set of valid inequalities derived from the expression of the quantile function and linear programming duality theory. The existence of binary variables x and y with max(0, •) function and the non-convexity of quantile function Q t τ (x) makes the MIP_full model very difficult to solve optimally with state-of-the-art MIP solvers like Cplex and Gurobi except for small size instances.

Scenario penalization matheuristic

Our new scenario penalization matheuristic (SPM) for solving the SSP optimization problem is based on four main steps. Each solution obtained throughout Step 1 and

Step 2 of the SPM (Algorithm 1) is evaluated as a candidate for the best solution x * currently found.

Algorithm 1 SPM()

Step 0: Create an initial approximation MIP_mean(c 0 ) of MIP_full model by ignoring the excess cost and set c = c 0 and Objective(x * ) = ∞. while termination criterion is not met do

Step 1: Solve MIP_mean(c) yielding a solution x.

Step 2: Starting from x, apply a Local Search to obtain a better solution x .

Step 3: Update the coefficient matrix c to take into account the excess from x and update the best solution x * = argmin{Objective(y) : y ∈ {x, x , x * }}. end while return x *

In the following, we give details of these steps as adapted to the present context.

Step 0: initial MIP_mean(c 0 ) approximation generating initial solution x 0

As noted above, the MIP_full model is too hard to solve for real world instances. Hence, the initial approximation we consider corresponds to the mean model ignoring the excess. Formally, let c 0 i,t denote the initial resulting mean risk if intervention i ∈ I starts at time step t ∈ T , defined as follows

c 0 i,t = 1 |T | t ∈T + (i,t) 1 |S t | s∈S t risk s,t i,t (6)
Then the initial linear approximation of the MIP_full model without excess cost is given below. 2) and (3)

(MIP_mean(c 0 )) ⎧ ⎨ ⎩ min Mean(x) = i∈I t∈T c 0 i,t × x i,t s.t.: (1), (
x i,t ∈ {0, 1} ∀ i ∈ I, t ∈ T
Required computational time to solve MIP_mean(c 0 ) varies across the data sets: it is usually quite fast but can take up to a few minutes for largest instances. Consequently, solving MIP_mean(c 0 ) optimally or approximatively with an MIP solver will generate an initial feasible solution x 0 for the MIP_full model (i.e. x 0 corresponds to the first solution generated by SPM() Algorithm.

Step 3: scenario penalization heuristic updating the current MIP_mean(c) approximation

Our goal in this step is to modify the current objective coefficient matrix c of the MIP_mean(c) model, stating with c = c 0 to make the resulting model closer to the original MIP_full model. To achieve this, let x be an optimal or best solution of the current MIP_mean(c) model and define the excess risk cost at time step t ∈ T by:

Excess t (x) = max(0, Q t τ (x) - 1 |S t | s∈S t risk s,t (x)) (7) Hence, Excess(x) = 1 |T |
t∈T Excess t (x). We first provide observations that describe the behavior of the function Excess t (x). For any step time t ∈ T , we define the set of scenarios

S + t (x) = {s ∈ S t : risk s,t (x) > Q t τ (x)}.
It is obvious that this modification increases the mean value 1

|S t | s∈S t risk s,t (x).
Moreover, this modification involves only binary variables y s,t = 0 since S + t (x) = {s ∈ S t : y s,t = 0}, and the constraints (5) imposed on the quantile remain satisfied. Consequently, the quantile value Q t τ (x) is unchanged and the value Excess t (x) will decrease after this change. Hence for any time step t ∈ T , increasing the risk value risk s,t (x) for any scenario s ∈ S + t (x) decreases the value Excess t (x). Let Risk t mean (x) (resp. Risk t max (x)) denote the mean (resp. max) risk at time step t ∈ T :

Risk t mean (x) = 1 |S t | s∈S t risk s,t (x) (8) Risk t max (x) = max{risk s,t (x) : s ∈ S t } (9)
Then we can state

t (x) = Q t τ (x) - 1 |S t | s∈S t risk s,t (x)
Moreover, from the definition of the τ quantile Q t τ (x) in Sect. 2.2.1, in this case we have

Q t τ (x) = min{q ∈ R : ∀s ∈ S t , risk s,t (x) ≤ q}
Or equivalently

Q t τ (x) = Risk t max (x) = max{risk s,t (x) : s ∈ S t }
This completes the validation of Eq. ( 10).

Next, observe that the objective coefficients c i,t , can be written:

c i,t = 1 |T | t ∈T + (i,t) s∈S t risk s,t i,t × 1 |S t | (11)
We replace 1 |S t | values in formula ( 11) by non-constant values μ s,t , to yield the following formula for calculating c i,t values:

c i,t = 1 |T | t ∈T + (i,t) s∈S t risk s,t i,t × μ s,t (12)
The μ s,t values, for s ∈ S and t ∈ T , are set to 1 |S t | initially and penalizing scenario s at time step t will be done by increasing μ s,t value. Note that the sum of these values over all scenarios in the time step equals 1 initially ( s∈S t 1 |S t | = 1) and we will keep the sum fixed at 1 when doing the modifications. This condition is imposed because we want the sum of values/coefficients to be the same in each time step in order to not prioritize some time steps over others-only the importance of scenarios inside the time step will change.

Formally, we calculate μ s,t values in the following way. Let Penalty s,t be a current penalty for scenario s at time step t and let Freq s,t denote how many times scenario s at time step t is penalized. Penalty and frequency values are initialized to 0 and they will be updated in each iteration i.e. after each call of MIP_mean(c) as follows. We increase the penalties (in formula ( 13)) such that no excess appear in the solution obtained by MIP_Mean(c) (i.e. Q t τ (x) = Risk t mean (x)). Note that if |S + t (x)|= 0, no penalties are given to the scenario coefficients of this time step.

Step 2: local search (LS)

Step 2 constitutes an improvement method that transforms a feasible solution into one or more best solutions. The efficacy of Local Search (LS) (also called neighborhood search) for solving a wide variety of optimization problems (cf. [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF][START_REF] Glover | Tabu Search[END_REF][START_REF] Hoos | Stochastic Local Search: Foundations and Applications[END_REF][START_REF] Hansen | Variable neighborhood search: basics and variants[END_REF]) depends strongly on the characteristics of one or more neighborhood structures involved and their combination. We are interested in two well-known moves, namely Shift and Swap moves. The shift neighborhood is the set of solutions that can be obtained by changing the start time of one intervention, while the swap neighborhood is the set of solutions that can be obtained by interchanging the start times of two interventions. Formally, given a feasible schedule σ , for any intervention i ∈ I and any period t ∈ T , the shift neighbor solution σ = Shift(σ, i, t) is defined by

σ i = t and σ i = σ i for i ∈ I -{i} and for two interventions i 1 , i 2 ∈ I the swap neighbor solution σ = Swap(σ, i 1 , i 2 ) is defined by σ i 1 = σ i 2 , σ i 2 = σ i 1 and σ i = σ i for i ∈ I -{i 1 , i 2 }
Consequently, the Shift and Swap neighborhood sets are defined by

Shift(σ ) = {σ = Shift(σ, i, t) : σ is feasible, i ∈ I, t ∈ T } Swap(σ ) = {σ = Swap(σ, i 1 , i 2 ) : σ is feasible, i 1 , i 2 ∈ I}
where the statement that σ is feasible means that the solution σ satisfies the linear constraints (1), ( 2

) and (3).

There are several ways to combine two or more neighborhoods and computational results show that certain combinations are superior to others (cf. Di [START_REF] Gaspero | Neighborhood portfolio approach for local search applied to timetabling problems[END_REF][START_REF] Lü | Neighborhood analysis: a case study on curriculum-based course timetabling[END_REF][START_REF] Mjirda | Sequential variable neighborhood descent variants: an empirical study on the traveling salesman problem[END_REF][START_REF] Hansen | Variable neighborhood search: basics and variants[END_REF]). The three basic combinations of two neighborhood structures are strong neighborhood union, selective neighborhood union and token-ring search. For strong neighborhood union, the LS algorithm picks each move (according to the algorithm's selection criteria) from all the Shift and Swap moves. For selective neighborhood union, the LS algorithm selects one of the two neighborhoods to be used at each iteration, choosing the neighborhood Shift with a predefined probability p and choosing Swap with probability 1p. Note that an LS algorithm using only Shift or Swap is a special case of an algorithm using selective neighborhood union where p is set to be 1 and 0 respectively. In token-ring search, the neighborhoods Shift and Swap are alternated, applying the currently selected neighborhood without interruption, starting from the local optimum of the previous neighborhood, until no improvement is possible.

For the SSP challenge problem, the token-ring search implemented in the LS(σ, TL, TL 1 , TL 2 ) procedure is described in Algorithm 2. Specifically, the LS procedure starts from a given feasible solution σ , and uses one neighborhood until a best solution is determined, subject to time limits imposed on the search (TL 1 (resp. TL 2 ) for Shift (resp. Swap) neighborhood exploration. Then the method switches to the other neighborhood, starting from this best solution, and continues the search in the same fashion. The search comes back to the first neighborhood at the end of the second neighborhood exploration, repeating this process until time limit TL is reached. Neighborhoods Shift and Swap are explored by randomly selecting a candidate feasible move which is performed if it improves the cost. However, performing only the moves that improve the objective function can quickly lead to a local optimum. Escaping a local optimum is therefore achieved by occasionally allowing non-improving moves. This is done in a simple way: from time to time (for example once in 1000 iterations) a 'bad' move is accepted.

Algorithm 2 LS(σ, TL, TL 1 , TL 2 ) Input: feasible solution σ , running time limits TL, TL 1 and TL 2 while time limit TL not exceeded do while time limit TL 1 not exceeded do explore Shift neighborhood select random σ ∈ Shift(σ ) if σ is accepted then σ = σ end if end while σ = best solution found in SHIFT while time limit TL 2 not exceeded do explore Swap neighborhood select random σ ∈ Swap(σ ) if σ is accepted then σ = σ end if end while σ = best solution found in SWAP end while return σ

Full SPM algorithm

In this section, we specify the full description of the SPM algorithm. The preliminary initialization (Step 0) constructs the coefficient matrix c = c 0 of the initial MIP_Mean(c) approximation ignoring the excess cost and set the penalties and frequencies to 0 (i.e. Penalty s,t = Freq s,t = 0, ∀t ∈ T , s ∈ S t ). A each iteration, the current MIP_Mean(c) approximation is solved to generate a current feasible solution x (i.e. Step 1). The LS heuristic procedure tries to generate an improved solution x (i.e. Step 2). Step 3 updates the coefficient matrix c of the current MIP_Mean(c) approximation. The process is repeated until the time limit allocated is exceeded. The full algorithm pseudo-code is given in Algorithm 3.

Algorithm 3 SPM(TL) for i ∈ I do

Step 0: Initialization

for t ∈ T do c i,t = c 0 i,t
end for end for for t ∈ T do for s ∈ S t do Penalty s,t = Freq s,t = 0 end for end for while time limit TL not exceeded do x = MIP_Mean(c)

Step 1: Solve the current approximation MIP_Mean(c) x = LS(σ (x))

Step 2: Local Search for t ∈ T do

Step To clarify the procedure, we consider the small example provided in the next section.

scenario ( s,1 1,1 ) ( s,2 1,2 ) ( s,3 1,3 ) ( s,1 2,1 ) ( s,2 2,2 ) ( s,3 2,3 ) ( s,1 3,1 ) ( s,2 3,2 ) ( s,3 3,3 ) ( s,1 4,1 ) ( s,2 4,1 ) ( s,2 4,2 ) ( s,

A small example

We illustrate Algorithm 3 for an instance with four interventions I = {i 1 , i 2 , i 3 , i 4 }, three time steps T = {1, 2, 3}, seven scenarios in each time step S = {s 1 , . . . , s 7 }, τ = 0.7 and α = 0.5. The duration of the first three interventions is 1 (i.e. δ i 1 ,t = δ i 2 ,t = δ i 3 ,t = 1) and the duration of the last intervention is 2 (δ i 4 ,t = 2). All risk values are given in Table 2. For example, starting at time step 1, the risk s 4 ,2 i 4 ,1 value of intervention i 4 for scenario s 4 at time step 2 is equal to 8 (bold cell). The initial objective function coefficients c i,t are computed using formula (6): we obtain the following value for intervention i = 1 starting day t = 1: c 1,1 = 1 3 × 1 7 × (10 + 12 + 8 + 5 + 12 + 7 + 10) = 64 21 = 3.048. The evolution of these coefficients during 8 iterations of SPM procedure is given in Table 3.

Solving MIP_Mean(c) with these c i,t values produces the solution with optimal mean risk. Since we do not have any resource or exclusion constraints, the optimal solution can be obtained by simply selecting the time step with the smallest c i,t value for each intervention i. Selected time steps are bold at iteration 1 given the solution x 1,1 = x 2,3 = x 3,3 = x 4,1 = 1, and in italic at iteration 2 given the solution x 1,2 = x 2,3 = x 3,1 = x 4,2 = 1. We have calculated new scenario coefficients, update penalties and frequencies only for the two scenarios with greatest risk (|S t |-τ ×|S t | = 7-5 = 2)-achieved in Algorithm 3 by adding 1 to the previous frequency values and adding penalty =

(Q t τ (x)-Risk t mean (x))×|S t | | t |-τ ×| t |
to the previous penalty values, giving us the values shown in line 2 of S Table S 3. In the last line of Table 3 bolditalic values correspond to the solution x 1,3 = x 2,3 = x 3,1 = x 4,1 = 1 which is optimal.

Experimental results in the context of ROADEF-EURO challenge

Our SPM matheuristic obtains the best results in the final phase of the ROADEF/EURO 2020 challenge. In this phase, 13 teams (i.e. algorithms), denoted by A = {J 3, J 24, J 43, J 49, J 73, S14, S19, S28, S34, S56, S58, S66, S68}, were qualified as finalists over 74 registered junior and senior teams. The SPM matheuristic corresponds to the algorithm proposed by team S34 composed by the two first authors of this paper.

The computer used to evaluate the programs of the teams is a Linux OS machine with 2 CPU, 16 GB of RAM, and the allowed list of MILP solvers are CPLEX, Gurobi and LocalSolver. Our winning SPM algorithm was implemented in C++ language and used Gurobi as an MIP solver to solve the MIP_Mean(c) approximations of MIP_full.

Dataset To evaluate the 13 finalist algorithms, a set of 30 industrial problem instances was provided by RTE, divided into two data sets (Data set C and Data set X ) with 15 instances for each one. Those data sets are available at the website of the challenge (ROADEF 2020). The characteristics of these 30 instances are given in Table 4.

Time limit Two stopping criteria of each algorithm on each instance are used to differentiate the 13 algorithms by imposing two time limits of 15 and 90 min.

Objective evaluation Given a time limit TL and a dataset K = C or X . For each instance k ∈ K and each algorithm a ∈ A, the evaluation function of the returned solution x by algorithm a is given by obj TL k,a = Objective(x) if x is feasible, and otherwise obj TL k,a = ∞. The final score attribuated to an algorithm a ∈ A on an instance k ∈ K uses the convex weighting obj k,a = 0.8 × obj 15 k,a + 0.2 × obj 90 k,a .

Finally, given a dataset K , for each algorithm a ∈ A, the number of times where algorithm a fails to find a feasible solution or crashes with time limit TL is given by Crash TL K ,a = |{k ∈ K : a is infeasible or crashes with time limit TL}|.

Ranking method Let Better k,a be the number of algorithms with a result strictly better than the result of algorithm a on instance k, i.e.

Better k,a = |{a ∈ A : a = a, obj k,a < obj k,a }|.

For each dataset K , we compute the sum of Better k,a over k ∈ K , i.e. The score of an algorithm a for the dataset K is given by Score * K ,a = |{k ∈ K : Score k,a = Better * }|.

Given a time limit TL and a dataset K , for each algorithm a ∈ A, the number of times where the algorithm a is the best in terms of the relative gap value, is defined as follows

Best TL K ,a = |{k ∈ K : TL k,a = 0}| where TL k,a = 1000 × obj k,a -obj * k obj k,a
and where obj * k is the best value over the algorithms a ∈ A, i.e. obj * k = min{obj k,a : a ∈ A}.

Table 5 Our SPM algorithm (i.e. a * = S34) strictly dominated all the other algorithms a ∈ A -{a * } on the 18 comparison criteria (Score K ,a , Score * K ,a , Better K ,a , Best TL K ,a ) for K ∈ {C, X , C + X } and TL ∈ {15, 90}, (which justifies the word "efficient" in the title of this paper!). However, there is no dominance among the remaining 4 best algorithms, i.e. A * = {a ∈ A : Score * C+X ,a ≥ 150} with or without the 6 criteria Crash TL K ,a for K ∈ {C, X , C + X } and TL ∈ {15, 90}. For example, the second algorithm a 2 = S66 is dominated by the third algorithm a 3 = S56 on the two criteria Best 15 K ,a ) and Crash 15 K ,a . The fourth algorithm a 4 = S19 dominates the second algorithm a 2 and the third algorithm a 3 on the 4 criteria (Score * C,a , Better * C,a , Crash 90 X ,a , Crash 15 C+X ,a ). Basic lower and upper bounds: Since for any feasible solution x the value Excess(x) is non negative and α ∈ [0, 1] we have Objective(x) ≥ Mean(x). Consequently, a lower bound on the optimal value of MIP_full can be obtained by solving the MIP_mean(c 0 ), i.e.

LB = α × v(MIP_mean(c 0 )) ≤ v(MIP_full)
where v(P) is the optimal value of a given optimization problem P.

From Proposition 2, for any t ∈ T such that |S t |= τ ×|S t | , we have Excess t (x) = Risk t max (x) -Risk t mean (x) which is valid for τ = 1, the objective function can be expressed as follows: In this case (i.e. τ = 1), the non-linear optimization MIP_full can be stated as the linear mixed integer programming problem:

Objective(x) = 2 × α -1 |T | t∈T Risk t mean (x) + (1 -α) |T | t∈T Risk t max (x)
(MIP_Q1) ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ min obj = 2×α-1 |T | t∈T Risk t mean (x) + (1-α) |T | t∈T z t s.t.: (1), (2) and (3), z t ≥ Risk s,t (x), ∀t ∈ T , s ∈ S t , x ∈ {0, 1} |I|×|T | , z ∈ R |T |
where z t represents Risk t max (x) for t ∈ T . Since Risk t max (x) ≤ Risk t mean (x), a basic upper bound

UB = v(MIP_Q1) ≥ v(MIP_Full).
Table 6 gives the lower bound (LB) and the upper bound (UB) for C and X datasets. MIP solver has been used to compute the bounds: we have set parameters gap to 0.0005 and time limit to 3600 s. All solutions produced by SPM algorithm have a better objective than UB. But non exact approaches can produce worse solutions than UB.

Behavior of SPM algorithm Results have been obtained by running the SPM algorithm with 10 seeds. Last column gives the best values obtained by the challenge competitors. Table 6 below summarizes the results obtained by our approach.

Conclusion

Our new scenario penalization matheuristic (SPM), for the mean-risk model involving quantiles related to value at risk measure (VAR) is an instance of the highly general framework of ghost image processing (GIP) proposed by [START_REF] Glover | Optimization by ghost image processes in neural networks[END_REF], and is based on both mixed integer programming (MIP) models and local search (LS) methods. In this paper, we considered the application of the mean-risk model arising in the stochastic scheduling problem (SSP) formulated as a MIP model (MIP_full) with binary decision variables where the objective function is a non-linear (generally non-convex) function and the constraints are linear. The SPM matheuritic involves a parameterization of the objective function that is progressively modified to generate feasible solutions which are improved by local search. The initial parameterization provides an approximation MIP_mean(c) that corresponds to the MIP_full model by ignoring the excess cost (the non linear part of the original objective function). At each iteration, from an optimal or best solution of the current approximation MIP_mean(c), the parameter cost matrix c is modified in order to improve the next approximation by penalizing the scenario risks. The current feasible solution is improved by a local search procedure based on the combination of shift and swap neighborhoods. This approach won the first prize in the senior category of the EURO/ROADEF 2020 challenge. Experimental results are presented which demonstrate the effectiveness of our approach on large instances provided by the French electricity transmission network RTE. 11,202.71 33,514.49 33,513.15 33,511.70 33,513.48 33,512.63 33,512.26 33,795.93 C04 11,076.18 37,592.72 37,589.07 37,585.73 37,590.55 37,587.86 37,586.31 12,923.01 43,343.11 43,342.48 43,342.49 43,342.15 43,341.83 43,341.84 ,726.23 12,720.54 12,721.13 12,725.77 12,719.95 12,718.79 14,062.89 C13 12,730.23 42,490.05 42,488.24 42,487.99 42,486.11 42,484.17 42,484.56 42,722.61 C14 10,696.36 26,484.12 26,476.34 26,467.22 26,475.58 26,462.26 26,457.11 31,291.59 C15 13,913.86 39,759.93 39,759.43 39,758.03 39,758.25 39,757.03 39,757.54

  Better K ,a = k∈K Better k,a and the number of instances with Better k,a = 0 over k ∈ K , i.e. Better 0 K ,a = |{k ∈ K : Better k,a = 0}|. The score of an algorithm a for the instance k is defined by Score k,a = max(0; Better * -Better k,a ) if x is feasible 0 if x is unfeasible or crash where Better * is the maximal score that an algorithm can earn from one instance. During both the qualification and final phases, Better * is equal to 10. The global score of an algorithm a, on a given dataset K denoted as score(a) is defined by Score K ,a = k∈K Score k,a .

  For any t ∈ T such that |S t |= τ × |S t | , we have

	Excess t (x) = Risk t max (x) -Risk t mean (x).	(10)

Proposition 1

Proof If |S t |= τ × |S t | then from the τ quantile constraint (5) we deduce that the indicator binary variables y s,t for all s ∈ S t are set to 1. Consequently, from (4) we have Q t τ (x) ≥ risk s,t (x) for all s ∈ S t and the operator max(0, .) can be dropped from the expression to compute Excess t (x), i.e. now we have Excess

  For each time step t ∈ T with excess i.e. for which Q t τ (x) > Risk t mean (x) we update the penalties and frequencies for |S t |-τ × |S t | scenarios with the highest risk s,t (x) values (those are scenarios with risk values at least Q t

		μ s,t =	β s,t s ∈S t β s ,t	, ∀t ∈ T , s ∈ S t	(15)
	where
	β s,t =	1 1 + Penalty s,t /Freq s,t otherwise i f Freq s,t = 0	, ∀t ∈ T , s ∈ S t	(16)

τ (x)):

Penalty s,t = Penalty s,t + (Q t τ (x) -Risk t mean (x)) × |S t | |S t |-τ × |S t | (13)

Freq s,t = Freq s,t + 1 ( 1 4 ) Then the μ s,t values are calculated by setting:

Table 2

 2

	risk s,t i,t for 7 scenarios
	Intervention i1	Intervention i2	Intervention i3	Intervention i4

Table 3

 3 Objective coefficients (c i,t ) evolution in 8 iterations of SPM procedure without local search: column Full obj. corresponds to α × Mean + (1 -α) × Excess

	Iter.	Intervention i1	Intervention i2	Intervention i3	Intervention i4	Mean Full obj.
	1	3.048	3.095	3.048	3.048	3.143	2.476	3.095	3.190	3.00	7.143	7.381	7.833	8.333
	2	3.396	3.125	3.528	3.115	3.146	2.796	2.760	3.250	3.00	7.604	7.361	8.021	9.000
	3	3.302	3.229	3.428	3.246	3.156	2.644	3.016	3.458	3.17	8.291	8.702	8.590	8.500
	4	3.333	3.491	3.389	3.250	3.351	2.524	2.981	3.134	3.26	8.028	8.558	8.433	8.333
	5	3.280	3.321	3.238	3.043	3.267	2.475	3.135	3.125	3.17	7.492	8.017	8.165	8.238
	6	3.273	3.321	3.323	3.050	3.267	2.593	3.140	3.125	3.10	7.521	7.989	8.245	8.333
	7	3.358	3.323	3.366	3.079	3.268	2.633	3.021	3.108	3.11	7.546	7.832	8.261	8.500
	8	3.364	3.447	3.361	3.077	3.364	2.628	3.018	3.072	3.12	7.690	7.920	8.349	8.167

Table 4

 4 Characteristics of the 30 instances with |J |= 9 and α = 0.5

	Inst	Data set C				Data set X
		|I|	| T |	|E|	| S t |	τ	|I|	|T |	|E|	| S t |	τ
	1	120	53	54	169-207	0.95	120	53	48	169-207	0.80
	2	120	53	43	169-207	0.80	706	53	1234	56-69	0.85
	3	706	53	1223	56-69	0.85	280	53	162	169-207	0.80
	4	706	53	1194	56-69	0.90	426	25	490	175-203	0.80
	5	706	53	1377	56-69	0.95	467	220	604	84-103	0.85
	6	280	53	183	169-207	0.80	528	300	703	45-55	0.95
	7	120	42	38	113-138	0.95	209	300	80	56-69	0.90
	8	426	25	340	175-207	0.80	209	300	57	56-69	0.90
	9	110	53	38	169-207	0.90	548	30	820	141-173	0.80
	10	522	102	705	56-69	0.95	460	35	527	146-173	0.95
	11	89	102	35	171-207	0.90	521	131	725	56-69	0.95
	12	298	191	195	84-103	0.80	522	131	723	56-69	0.95
	13	505	230	533	56-69	0.95	336	212	248	84-103	0.90
	14	465	220	620	84-103	0.85	613	180	951	56-69	0.95
	15	528	300	624	45-55	0.95	613	180	917	56-69	0.95

  presents the values of the 6 parameters (Score K ,a , Score * K ,a , Better K ,a , Best TL K ,a , Crash TL K ,a ) used to evaluate an algorithm a ∈ A over a dataset K ∈ {C, X , C + X }. The values provided in this table are obtained from disggregated information described in the 4 matrices (Score k,a , Better k,a , 15 k,a , 90 k,a ) presented in the Appendix. The final ranking of the ROADEF-EURO challenge is based on the Score * C+X ,a values for the submitted algorithms a ∈ A. The winning algorithm a * is the one that receives the highest score, i.e.

	Score * C+X ,a

* = max{Score * C+X ,a : a ∈ A} .

Table 5

 5 Final results of the ROADEF/EURO 2020 challenge on datasets C, X

	+ X
	and C

Table 6

 6 Behavior of SPM algorithm

	Instance LB	15 min			90 min			UB
			Average Best	Challenge Average Best	Challenge
	C01	2004.38	8515.90	8515.90	8515.90	8515.90	8515.90	8515.90	8677.32
	C02	2197.87	3541.27	3540.24	3541.65	3539.67	3538.16	3539.80	4168.92
	C03