Hélène Collavizza

Michel Rueher
email: rueher@polytech.unice.fr

Pascal Van Hentenryck

CPBPV: A Constraint-Programming Framework For Bounded Program Verification

This paper studies how to verify the conformity of a program with its specification and proposes a novel constraint-programming framework for bounded program verification (CPBPV). The CPBPV framework uses constraint stores to represent the specification and the program and explores execution paths nondeterministically. The input program is partially correct if each constraint store so produced implies the post-condition. CPBPV does not explore spurious execution paths as it incrementally prunes execution paths early by detecting that the constraint store is not consistent. CPBPV uses the rich language of constraint programming to express the constraint store. Finally, CPBPV is parametrized with a list of solvers which are tried in sequence, starting with the least expensive and less general. Experimental results often produce orders of magnitude improvements over earlier approaches, running times being often independent of the variable domains. Moreover, CPBPV was able to detect subtle errors in some programs while other frameworks based on model checking have failed.

Introduction

This paper is concerned with software correctness, a critical issue in software engineering. It proposes a novel constraint-programming framework for bounded program verification (CPBPV), i.e., when the program inputs (e.g., the array lengths and the variable values) are bounded. The goal is to verify the conformity of a program with its specification, that is to demonstrate that the specification is a consequence of the program. The key idea of CPBPV is to use constraint stores to represent the specification and the program, and to non-deterministically explore execution paths over these constraint stores. This non-deterministic constraint-based symbolic execution incrementally refines the constraint store, which initially consists of the precondition. Non-determinism occurs when executing conditional or iterative instructions and the non-deterministic execution refines the constraint store by adding constraints coming from conditions and from assignments. The input program is partially correct if each constraint store produced by the symbolic execution implies the post-condition. It is important to emphasize that CPBPV considers programs with complete specifications and that verifying the conformity between a program and its specification requires to check (explicitly or implicitly) all executables paths. This is not the case in model-checking tools designed to detect violations of some specific property, e.g., safety or liveness properties.

The CPBPV framework has a number of fundamental benefits. First, contrary to earlier work using constraint programming or SMT [START_REF] Benerecetti | Abstraction Refinement of Linear Programs with Arrays[END_REF][START_REF] Collavizza | Software Verification using Constraint Programming Techniques[END_REF][START_REF] Collavizza | Exploring different constraint-based modelings for program verification Procs of CP[END_REF], CPBPV does not use predicate abstraction or explore spurious execution paths, i.e., paths that do not correspond to actual executions over inputs satisfying the pre-condition. CPBPV incrementally prunes execution paths early by detecting that the constraint store is not consistent. Second, CPBPV uses the rich language of constraint programming to express the constraint store, including arbitrary logical and threshold combination of constraints, the element constraint, and global/combinatorial constraints that express complex relationships on a set of variables. Finally, CPBPV is parametrized with a list of solvers which are tried in sequence, starting with the least expensive and less general.

The CPBPV framework was evaluated experimentally on a series of benchmarks from program verification. Experimental results of our (slow) prototype often produce orders of magnitude improvements over earlier approaches, and indicate that the running times are often independent of the variable domains. Moreover, CPBPV was able to found subtle errors in some programs that some other verification frameworks based on model-checking could not detect.

The rest of the paper is organized as follows. Section 2 illustrates how CPBPV handles constraints store on a motivating example. Section 3 formalizes the CPBPV framework for a small programming language and Section 4 discusses the implementation issues. Section 5 presents experimental results on a number of verification problems, comparing our approach with state of the art modelchecking based verification frameworks. Section 6 discusses related work in test generation, bounded program verification and software model checking. Section 7 summarizes the contributions and presents future research directions.

The Constraint-Programming Framework at Work

This section illustrates the CPBPV verifier on a motivating example, the binary search program. CPBPV uses Java programs and JML specifications for the pre-and post-conditions, appropriately enhanced to support the expressivity of constraint programming. Figure 1 depicts a binary search program to determine if a value v is present in a sorted array t. (Note that \result in JML corresponds to the value returned by the program). To verify this program, our prototype implementation requires a bound on the length of array t, on its elements, and on v. We will verify its correctness for specific lengths and simply assume that the values are signed integers on a number of bits.

The initial constraint store of the CPBPV verifier, assuming an input array of length 8, is the precondition

1 c pre ≡ ∀0 ≤ i < 7 : t 0 [i] ≤ t 0 [i + 1]
where t 0 is an array of constraint variables capturing the input. The constraint variables are annotated with a version number as CPBPV performs a SSA-like renaming /*@ requires (\forall int i; i>=0 && i<t.length-1;t Since l 0 ≤ u 0 , it enters the loop body, adds the constraint m 0 = (l 0 + u 0 )/2, which simplifies to m 0 = 3, and considers the conditional statement on line 6. The execution of the statement is nondeterministic: Indeed, both t 0 [3] = v 0 and t 0 [3] = v 0 are consistent with the constraint store, so that the two alternatives, which give rise to two execution paths, must be explored. Note that these two alternatives correspond to actual execution paths in which t [START_REF] Mantovani | Bounded Model Checking of C Programs using a SMT solver instead of a SAT solver[END_REF] in the input is equal to, or different from, input v. The first alternative adds the constraint t 0 [3] = v 0 to the store and executes line 7 which adds the constraint result = m 0 . CPBPV has thus obtained an execution path p whose final constraint store c p is:

[i]<=t[i+1]) @ ensures @ (\result != -1 ==> t[\result] == v) && @ (\result == -1 ==> \forall int k; 0 <= k < t.length ; t[k] != v) @*/ 1 static int binary_search(int[] t, int v) { 2 int l = 0; 3 int u = t.length-1; 4 while (l <= u) { 5 int m = (l + u) / 2; 6 if (t[m]==v) 7 return m; 8 if (t[m] > v) 9 u = m
cpre ∧ l 0 = 0 ∧ u 0 = 7 ∧ m 0 = (l 0 + u 0 )/2 ∧ t 0 [m 0 ] = v 0 ∧ result = m 0
CPBPV then checks whether this store c p implies the post-condition c post by searching for a solution to c p ∧ ¬c post . This test fails, indicating that the computation path p, which captures the set of actual executions in which t[3] = v, satisfies the specification. CPBPV then explores the other alternatives to the conditional statement in line 6. It adds the constraint t 0 [m 0 ] = v 0 and executes the conditional statement in line 8. Once again, this statement is nondeterministic. Its first alternative assumes that the test holds, generating the constraint t 0 [m 0 ] > v 0 and executing the instruction in line 9. Since u is (re-)assigned, CPBPV creates a new variable u 1 and posts the constraint

u 1 = m 0 -1 = 2.
The execution returns to line 4, where the test now reads l 0 ≤ u 1 , since CPBPV always uses the most recent version for each variable. Since the constraint stores entails l 0 ≤ u 1 , the only extension to the current path consists of executing line 5, adding the constraint m 1 = (l 0 + u 1 )/2, which actually simplifies to m 1 = 1.

Another complete execution path is then obtained by executing lines 6 and 7.

Consider now a version of the program in which line 11 is replaced by u = m-1. To illustrate the CPBPV verifier, we specify partial execution paths by indicating which alternative is selected for each nondeterministic instruction. For instance, T 4 , F 6 , T 8 , T 5 , T 6 denotes the last execution path discussed above in which the true alternative is selected for the first execution of the instruction in line 4, the false alternative for the first execution of instruction 6, the true alternative for the first instruction of instruction 8, the true alternative of the second execution of instruction 5, and the true alternative of the second execution of instruction 6. Consider the partial path T 4 , F 6 , F 8 and let us study how it can be extended. The partial path T 4 , F 6 , F 8 , T 4 , T 6 is not explored, since it produces a constraint store containing

c pre ∧ t 0 [3] = v 0 ∧ t 0 [3] ≤ v 0 ∧ t 0 [1]
= v 0 which is clearly inconsistent. Similarly, the path T 4 , F 6 , F 8 , T 4 , F 6 , T 8 cannot be extended. The output of CPBPV on this incorrect program when executed on an array of length 8 (with integers coded on 8-bits to make it readable) produces, in 0.025 seconds, the counterexample:

v 0 = -126 ∧ t 0 = [-128, -127, -126, -125, -124, -123, -122, -121] ∧ result = -1.
This example highlights a few interesting benefits of CPBPV.

1. The verifier only considers paths that correspond to collections of actual inputs (abstracted by constraint stores). The resulting execution paths must all be explored since our goal is to prove the partial correctness of the program. 2. The performance of the verifier is independent of the integer representation on this application: it only requires a bound on the length of the array. 3. The verifier returns a counter-example for debugging the program.

Note that CBM C and ESC/Java2, two state-of-the-art model checkers fail to verify this example as discussed in Section 5.

Formalization of the Framework

This section formalizes the CPBPV verifier on a small abstract language using a small-step SOS semantics. The semantics primarily specifies the execution paths over constraint stores explored by the verifier. It features assert and enforce constructs which are necessary for modular composition.

Syntax Figure 2 depicts the syntax of the programs and the constraints generated by the verifier. In the following, we use s, possibly subscripted, to denote elements of a syntactic entity S.

Renamings CPBPV creates variables and arrays of variables "on-the-fly" when they are needed. This process resembles an SSA normalization but does not introduce the join nodes, since the results of different execution paths are not merged. Similar renamings are used in model checking. The renaming uses mappings of type V ∪A → N which maps variables and arrays into a natural numbers L : list of instructions; I : instructions ; B : Boolean expressions E : integer expressions; A : arrays; V : variables

L ::= I; L | ǫ I ::= A[E] ← E | V ← E | if B I | while B I | assert(B ) | enforce(B ) | return E | {L} B ::= true | f alse | E > E | E ≥ E | E = E | E = E | E ≤ E | E < E B ::= ¬B | B ∧ B | B ∨ B | B ⇒ B E ::= V | A[E] | E + E | E -E | E × E | E/E | C : constraints E + : solver expressions V + = {v i | v ∈ V & i ∈ N } : solver variables A + = {a i | a ∈ A & i ∈ N } : solver arrays C ::= true | f alse | E + > E + | E + ≥ E + | E + = E + | E + = E + | E + ≤ E + | E + < E + C ::= ¬C | C ∧ C | C ∨ C | C ⇒ C E + ::= V | A[E + ] | E + + E + | E + -E + | E + × E + | E + /E + | Fig. 2

. The Syntax of Programs and Constraints

denoting their current "version numbers". In the semantics, the version number is incremented each time a variable or an array element is assigned. We use σ ⊥ to denote the uniform mapping to zero (i.e., ∀x ∈ V ∪ A : σ ⊥ (x) = 0) and σ[x/i] the mapping σ where x now maps to i, i.e., σ[x/i](y) = if x = y then i else σ(y). stating that configuration γ 1 can be rewritten into γ 2 when the conditions hold.

Conditional Instructions The conditional instruction if b i considers two cases. If the constraint c b associated with b is consistent with the constraint store, then the store is augmented with c b and the body is executed. If the negation ¬c b is consistent with the store, then the constraint store is augmented with ¬c b . Both rules may apply, since the store may represent some memory states satisfying the condition and some violating it.

c ∧ (ρ σ b) is satisfiable if b i ; l, σ, c -→ i ; l, σ, c ∧ (ρ σ b) c ∧ ¬(ρ σ b) is satisfiable if b i ; l, σ, c -→ l, σ, c ∧ ¬(ρ σ b)
Iterative Instructions The while instruction while b i also considers two cases. If the constraint c b associated with b is consistent with the constraint store, then the constraint store is augmented with c b , the body is executed, and the while instruction is reconsidered. If the negation ¬c b is consistent with the constraint store, then the constraint store is augmented with ¬c b .

c ∧ (ρ σ b) is satisfiable while b i ; l, σ, c -→ i; while b i ; l, σ, c ∧ (ρ σ b) c ∧ ¬(ρ σ b) is satisfiable while b i ; l, σ, c -→ l, σ, c ∧ ¬(ρ σ b)
Scalar Assignments Scalar assignments create a new constraint variable for the program variable to be assigned and add a constraint specifying that the variable is equal to the right-hand side. A new renaming mapping is produced.

σ 2 = σ 1 [v/σ 1 (v) + 1] & c 2 ≡ (ρ σ 2 v) = (ρ σ 1 e) v ← e ; l, σ 1 , c 1 -→ l, σ 2 , c 1 ∧ c 2

Assignments of Array Elements

The assignment of an array element creates a new constraint array, add a constraint for the index being indexed and posts constraints specifying that all the new constraint variables in the array are equal to their earlier version, except for the element being indexed. Note that the index is an expression which may contain variables as well, giving rise to the well-known element constraint in constraint programming [START_REF] Vanhentenryck | Constraint Satisfaction in Logic Programming[END_REF].

σ 2 = σ 1 [a/σ 1 (a) + 1] c 2 ≡ (ρ σ 2 a)[ρ σ 1 e 1 ] = (ρ σ 1 e 2 ) c 3 ≡ ∀i ∈ 0..a.length : (ρ σ 1 e 1 ) = i ⇒ (ρ σ 2 a)[i] = (ρ σ 1 a)[i] a[e 1 ] ← e 2 , σ 1 ; l, c 1 -→ l, σ 2 , c 1 ∧ c 2 ∧ c 3
Assert Statements An assert statement checks whether the assertion is implied by the control store in which case it proceeds normally. Otherwise, it terminates the execution with an error.

c ⇒ (ρ σ b) assert b ; l, σ, c -→ l, σ, c c ∧ ¬(ρ σ b) is satisfiable assert b ; l, σ, c -→ ⊥, σ, c
Enforce Statements An enforce statement adds a constraint to the constraint store if it is satisfiable.

c ∧ (ρ σ b) is satisfiable enforce b ; l, σ, c -→ l, σ, c ∧ (ρ σ b)
Block Statements Block statements simply remove the braces.

{l 1 } ; l 2 , σ, c -→ l 1 : l 2 , σ, c
Return Statements A return statement simply constrains the result variable.

c 2 ≡ (ρ σ 1 result) = (ρ σ 1 e) return e ; l, σ 1 , c 1 -→ σ 1 , c 1 ∧ c 2
Termination Termination also occurs when no instruction remains.

Implementation issues

The CPBPV framework is parametrized by a list of solvers (S 1 , . . . , S k ) which are tried in sequence, starting with the least expensive and less general. When checking satisfiability, the verifier never tries solver S i+1 , . . . , S k if solver S i is a decision procedure for the constraint store. If solver S i is not a decision procedure, it uses an abstraction α of the constraint store c satisfying c ⇒ α and can still detect failed execution paths quickly. The last solver in the sequence is a constraint-programming solver (CP solver) over finite domains which iterates pruning and searching to find solutions or prove infeasibility. When the CP solver makes a choice, the earlier solvers in the sequence are called once again to prune the search space or find solutions if they have become decision procedures. Our prototype implementation uses a sequence (M IP, CP ), where MIP is the mixed integer-programming tool ILOG CPLEX 2 and CP is the constraintprogramming tool Ilog JSOLVER. Our Java implementation also performs some trivial simplifications such as constant propagation but is otherwise not optimized in its use of the solvers and in its renaming process whose speed and memory usage could be improved substantially. Practically, simplifications are done on the fly and the MIP solver is called at each node of the executable paths. The CP solver is only called at the end of the executable paths when the complete post condition is considered. Currently, the implementation use a depth-first strategy for the CP solver, but modern CP languages now offer highlevel abstractions to implement other exploration strategies. In practice, when CPBPV is used for model checking as discussed below, it is probably advisable to use a depth-first iterative deepening implementation.

Experimental results

In this section, we report experimental results for a set of traditional benchmarks for program verification. We compare CPBVP with the following frameworks:

-ESC/Java is an Extended Static Checker for Java to find common run-time errors in JML-annotated Java programs by static analysis of the code and its annotations. All experiments were performed on the same machine, an Intel(R) Pentium(R) M processor 1.86GHz with 1.5G of memory, using the version of the verifiers that can be downloaded from their web sites (except for EUREKA for which the execution times given in [START_REF] Benerecetti | Abstraction Refinement of Linear Programs with Arrays[END_REF][START_REF] Mantovani | Bounded Model Checking of C Programs using a SMT solver instead of a SAT solver[END_REF] are reported.) For each benchmark program, we describe the data entries and the verification parameters. In the tables, "UNABLE" means that the corresponding framework is unable to validate the program either because a lack of expressiveness or because of time or memory limitations, "NOT FOUND" that it does not detect an error, and "FALSE ERROR" that it reports an error in a correct program. Complete details of the experiments, including input files and error traces, can be found in [START_REF] Collavizza | Comparison between CPBPV with ESC/Java, CBMC, Blast, EUREKA and Why[END_REF].

Binary search We start with the binary search program presented in figure 1.

ESC/Java is applied on the program described in Figure 1. ESC/Java requires a limit on the number of loop unfoldings, which we set to log(n) CBMC does not support first-order expressions such as JML \f orall statement, we generated a C program for each instance of the problem (i.e., each array length). For example, the postcondition for an array of length 8 is given by

+
(result!=-1 && a[result]==x)|| (result==-1 && (a[0]!=x&&a[1]!=x&&a[2]!=x&&a[3]!=x&&a[4]!=x&&a[5]!=x&&a[6]!=x&&a[7]!=x)
For the Why framework, we used the binary search version given in their distribution. This program uses an assert statement to give a loop invariant.

Note that CPBPV does not require any additional information: no invariant and no limits on loop unfoldings. During execution, it selects a path by nondeterministically applying the semantic rules for conditional and loop expressions.

Table 1 reports the experimental results. Execution times for CPBPV are reported as a function of the array length for integers coded on 31 bits. 3 Our implementation is neither optimized for time or space at this stage and times are only given to demonstrate the feasibility of the CPBPV verifier.

The "Why" framework [START_REF] Fillitre | The Why/Krakatoa/Caduceus Platform for Deductive Program Verification Proc. CAV[END_REF] was unable to verify the correctness without the loop invariant; 60% of the proof obligations remained unknown.

The CBMC framework was not able to do the verification for an instance of length 32 (it was interrupted after 6691,87s).

ESC/Java was unable to verify the correctness of this program unless complete loop invariants are provided 4 .

An Incorrect Binary search Table 2 reports experimental results for an incorrect binary search program (see Figure 1, line 11) for CPBPV, ESC/Java, CBMC, and Why using an invariant. The error trace found with CPBPV has been described in Section 2. The error traces provided by CBMC and ESC/Java only show the decisions taken along the faulty path can be found in [START_REF] Collavizza | Comparison between CPBPV with ESC/Java, CBMC, Blast, EUREKA and Why[END_REF]. In contrast to CPBPV, they do not provide any value for the array nor the searched data. Observe that CPBPV provides orders of magnitude improvements in efficiency over CBMC and also outperforms ESC/Java by almost a factor 8 on the largest instance. The tritype program in Java with its specification in JML can be found in [START_REF] Collavizza | Comparison between CPBPV with ESC/Java, CBMC, Blast, EUREKA and Why[END_REF]. Table 3 depicts the experimental results for CPBPV, ESC/Java, CBMC, BLAST and Why. BLAST was unable to validate this example because the current version does not handle linear arithmetic. Observe the excellent performance of CPBPV and note that our previous approach using constraint programming and Boolean abstraction to abstract the conditions, validated this benchmark in 8.52 seconds when integers were coded on 16 bits [START_REF] Collavizza | Exploring different constraint-based modelings for program verification Procs of CP[END_REF]. It also explored 92 spurious paths.

An Incorrect Tritype Program Consider now an incorrect version of Tritype

program in which the test "if ((trityp==2)&&(i+k>j))" in line 22 (see [START_REF] Collavizza | Comparison between CPBPV with ESC/Java, CBMC, Blast, EUREKA and Why[END_REF]) is replaced by "if ((trityp==1)&&(i+k>j))". Since the local variable trityp is equal to 2 when i==k, the condition (i+k)>j implies that (i,j,k) are the sides of an isosceles triangle (the two other triangular inequalities are trivial because j>0). But, when trityp=1, i==j holds and this incorrect version may answer that the triangle is isosceles while it may not be a triangle at all. For example, it will return 2 when (i,j,k)=(1,1,2). Bubble Sort with initial condition This benchmark (see [START_REF] Collavizza | Comparison between CPBPV with ESC/Java, CBMC, Blast, EUREKA and Why[END_REF]) is taken from [START_REF] Benerecetti | Abstraction Refinement of Linear Programs with Arrays[END_REF] and performs a bubble sort of an array t which contains integers from 0 to t.length given in decreasing order. Table 5 shows the comparative results for this benchmark. CPBPV was limited on this benchmark because its recursive implementation uses up all the JAVA stack space. This problem should be remedied by removing recursion in CPBPV.

Selection Sort

We now present a benchmark to highlight both modular verification and the element constraint of constraint programming to index arrays with arbitrary expressions. The benchmark described in [START_REF] Collavizza | Comparison between CPBPV with ESC/Java, CBMC, Blast, EUREKA and Why[END_REF]. Assume that function findMin has been verified for arbitrary integers. When encountering a call to findMin, CPBPV first checks if its precondition is entailed by the constraint store, which requires a consistency check of the constraint store with respect to the negation of the precondition. Then CPBPV replaces the call by the postcondition where the formal parameters are replaced by the actual variables. In particular, for the first iteration of the loop and an array length of 40, CPBPV generates the conjunction 0

≤ k 0 < 40 ∧ t 0 [k 0 ] ≤ t 0 [0] ∧ . . . ∧ t 0 [k 0 ] ≤ t 0 [39]
which features element constraint [START_REF] Vanhentenryck | Constraint Satisfaction in Logic Programming[END_REF]. Indeed, k 0 is a variable and a constraint like t 0 [k 0 ] ≤ t 0 [0] indexes the array t 0 of variables using k 0 .

The modular verification of the selection sort explores only a single path, is independent of the integer representation, and takes less than 0.01s for arrays of size 40. The bottleneck in verifying selection sort is the validation of function findMin, which requires the exploration of many paths. However the complete validation of selection sort takes less than 4 seconds for an array of length 6. Once again, this should be contrasted with the model-checking approach of Eureka [START_REF] Benerecetti | Abstraction Refinement of Linear Programs with Arrays[END_REF]. On a version of selection sort where all variables are assigned specific values (contrary to our verification which makes no assumptions on the inputs), Eureka takes 104 seconds on a faster machine. Reference [START_REF] Benerecetti | Abstraction Refinement of Linear Programs with Arrays[END_REF] also reports that CBMC some assumptions to mean that there is no overflow into the sums, in order to prove the correct version of tritype with this same option.

takes 432.6 seconds, that BLAST cannot solve this problem, and that SATABS [START_REF] Kroening | SATABS: SAT-Based Predicate Abstraction for ANSI-C[END_REF] only verifies the program for an array with 2 elements. Sum of Squares Our last benchmark is described in [START_REF] Collavizza | Comparison between CPBPV with ESC/Java, CBMC, Blast, EUREKA and Why[END_REF] and computes the sum of the square of the n first integers stored in an array. The precondition states that n is the size of the array and that t must contain any possible permutation of the n first integers. The postcondition states that the result is n × (n + 1) × (2 × n + 1)/6. The benchmark illustrates two functionalities of constraint programming: the ability of specifying combinatorial constraints and of solving nonlinear problems. The alldifferent constraint [START_REF] Régin | A filtering algorithm for constraints of difference in CSPs[END_REF] in the precondition specifies that all the elements of the array are different, while the program constraints and postcondition involves quadratic and cubic constraints. The maximum instance that we were able to solve with CPBPV was an array of size 10 in 66.179s.

CPLEX, the MIP solver, plays a key role in all these benchmarks. For instance, the CP solver is never called in the Tritype benchmark. For the Binary search benchmark, there are length calls to the CP solver but almost 75% of the CPU time is spent in the CP solver. Since there is only path in the Buble sort benchmark, the CP solver is only called once. In the Sum of squares example, 80% of the CPU time is spent in the CP solver.

Discussion and Related Work

We briefly review recent work in constraint programming and model checking for software testing, validation, and verification. We outline the main differences between our CPBPV framework and existing approaches.

Constraint Logic Programming Constraint logic programming (CLP) was used for test generation of programs (e.g., [START_REF] Gotlieb | Automatic Test Data Generation using Constraint Solving Techniques[END_REF][START_REF] Jackson | Finding Bugs with a Constraint Solver[END_REF][START_REF] Sy | Automatic Test Data Generation for Programs with Integer and Float Variables[END_REF][START_REF] Godefroid | Automated Whitebox Fuzz Testing[END_REF]) and provides a nice implementation tool extending symbolic execution techniques [START_REF] Botella | Symbolic execution of floating-point computations[END_REF]. Gotlieb et al. showed how to represent imperative programs as constraint logic programs and used predicate abstraction (from model checking) and conditional constraints within a CLP framework. Flanagan [START_REF] Flanagan | Automatic software model checking via constraint logic[END_REF] formalized the translation of imperative programs into CLP, argued that it could be used for bounded model checking, but did not provide an implementation. The test-generation methodology was generalized and applied to bounded program verification in [START_REF] Collavizza | Software Verification using Constraint Programming Techniques[END_REF][START_REF] Collavizza | Exploring different constraint-based modelings for program verification Procs of CP[END_REF]. The implementation used dedicated predicate abstractions to reduce the exploration of spurious execution paths. However, as shown in the paper, the CPBPV verifier is significantly more efficient and often avoids the generation of spurious execution paths completely.

Model Checking It is also useful to contrast the CPBPV verifier with modelchecking of software systems. SAT-based bounded model checking for software [START_REF] Clarke | Bounded Model Checking using Satisfiability Solving[END_REF] consists in building a propositional formula whose models correspond to execution paths of bounded length violating some properties and in using SAT solvers to check whether the resulting formula is satisfiable. SAT-based modelchecking platforms [START_REF] Clarke | Bounded Model Checking using Satisfiability Solving[END_REF] have been widely popular thanks to significant progress in SAT solvers. A fundamental issue faced by model checkers is the state space explosion of the resulting model. Various techniques have been proposed to address this challenge, including generalized symbolic execution (e.g., [START_REF] Khurshid | Generalized Symbolic Execution for Model Checking and Testing[END_REF]), SMTbased model checking, and abstraction/refinement techniques. SMT-based model checking is the idea of representing and checking quantifier-free formulas in a more general decidable theory (e.g. [START_REF] Ganzinger | DPLL(T): Fast Decision Procedures[END_REF][START_REF] Dutertre | A fast linear-arithmetic solver for DPLL(T)[END_REF][START_REF] Nieuwenhuis | Challenges in Satisfiability Modulo Theories[END_REF]). These SMT solvers integrate dedicated solvers and share some of the motivations of constraint programming. Predicate abstraction is another popular technique to address the state space explosion. The idea consists in abstracting the program to obtain an abstract program on which model checking is performed. The model checker may then generate an abstract counterexample which must be checked to determine if it corresponds to a concrete execution path. If the counterexample is spurious, the abstract program is refined and the process is iterated. A successful predicate abstraction consists of abstracting the concrete program into a Boolean program (e.g., [START_REF] Ball | Rajamani Boolean and Cartesian Abstraction for Model Checking C Programs[END_REF][START_REF] Kroening | A Tool for Checking ANSI-C programs[END_REF][START_REF] Kroening | Predicate abstraction of ANSI-C Programs using SAT[END_REF]). In recent work [START_REF] Mantovani | Bounded Model Checking of C Programs using a SMT solver instead of a SAT solver[END_REF][START_REF] Benerecetti | Abstraction Refinement of Linear Programs with Arrays[END_REF], Armando & al proposed to abstract concrete programs into linear programs and used an abstraction of sets of variables and array indices. They showed that their tool compares favourably and, on some of the programs considered in this paper, outperforms model checkers based on predicate abstraction. Our CPBPV verifier contrasts with SAT-based model checkers, SMT-based model checkers and predicate abstraction based approaches: It does not abstract the program and does not generate spurious execution paths. Instead it uses a constraint-solver and nondeterministic exploration to incrementally construct abstractions of execution paths. The abstraction uses constraint stores to represent sets of concrete stores. On many bounded verification benchmarks, our preliminary experimental results show significant improvements over the stateof-the-art results in [START_REF] Benerecetti | Abstraction Refinement of Linear Programs with Arrays[END_REF]. Model checking is well adapted to check low-level C program and hardware applications with numerous Boolean constraints and bitwise operations: It was successfully used to compare an ANSI C program with a circuit given as design in Verilog [START_REF] Kroening | A Tool for Checking ANSI-C programs[END_REF]. However, it is important to observe that in model checking, one is typically interested in checking some specific properties such as buffer overflows, pointer safety, or user-specified assertions. These properties are typically much less detailed than our post-conditions and abstracting the program may speed up the process significantly. In our CPBPV verifier, it is critical to explore all execution paths and the main issue is how to effectively abstract memory stores by constraints and how to check satisfiability incrementally. It is an intriguing issue to determine whether an hybridization of the two approaches would be beneficial for model checking, an issue briefly discussed in the next section. Observe also that this research provides convincing evidence of the benefits of Nieuwenhuis' challenge [START_REF] Nieuwenhuis | Challenges in Satisfiability Modulo Theories[END_REF] aiming at extending SMT6 with CP techniques.

Perspectives and Future Work

This paper introduced the CPBPV framework for bounded program verification. Its novelty is to use constraints to represent sets of memory stores and to explore execution paths over these constraint stores nondeterministically and incrementally. The CPBPV verifier exploits the fact that, when variables and arrays are bounded, the constraint store can always be checked for feasibility. As a result, it never explores spurious execution path contrary to earlier approaches combining constraint programming and predicate abstraction [START_REF] Collavizza | Software Verification using Constraint Programming Techniques[END_REF][START_REF] Collavizza | Exploring different constraint-based modelings for program verification Procs of CP[END_REF] or integrating SMT solvers and the abstraction/refinement approach from model checking [START_REF] Benerecetti | Abstraction Refinement of Linear Programs with Arrays[END_REF]. We demonstrated the CPBPV verifier on a number of standard benchmarks from model checking and program checking as well as on nonlinear programs and functions using complex array indexings, and showed how to perform modular verification. The experimental results demonstrate the potential of the approach:

The CPBPV verifier provides significant gain in performance and functionalities compared to other tools.

Our current work aims at improving and generalizing the framework and implementation. In particular, we would like to include tailored, light-weight solvers for a variety of constraint classes, the optimization of the array implementation, and the integration of Java objects and references. There are also many research avenues opened by this research, two of which are reviewed now.

Currently, the CPBPV verifier does not check for variable overflows: the constraint store enforces that variables take values inside their domains and execution paths violating these constraints are thus not considered. It is possible to generalize the CPBPV verifier to check overflows as the verification proceeds. The key idea is to check before each assignment if the constraint store entails that the value produced fits in the selected integer representation and generate an error otherwise. (Similar assertions must in fact be checked for each subexpression in the right hand-side in the language evaluation order. Interval techniques on floats [START_REF] Botella | Symbolic execution of floating-point computations[END_REF] may be used to obtain conservative checking of such assertions.

An intriguing direction is to use the CPBPV approach for properties checking. Given an assertion to be verified, one may perform a backward execution from the assertion to the function entry point. The negation of the assertion is now the pre-condition and the pre-condition becomes the post-condition. This requires to specify inverse renaming and executions of conditional and iterative statements but these have already been studied in the context of test generation.

Fig. 1 .

 1 Fig. 1. The Binary Search Program

  These mappings are used by a polymorphic renaming function ρ to transform program expressions into constraints. For example, ρ σ b 1 ⊕ b 2 = (ρ σ b 1 ) ⊕ (ρ σ b 2 )(where ⊕ ∈ {∧, ∨, ⇒}) is the rule used to transform a logical expression. Configurations The CPBCV semantics mostly uses configurations of the type l, σ, c , where l is the list of instructions to execute, σ is a version mapping, and c is the set of constraints generated so far. It also uses configurations of the form ⊤, σ, c to denote final states and configurations of the form ⊥, σ, c to denote the violation of an assertion. The semantics is specified by rules of the form conditions γ1 -→γ2

  ǫ, σ, c -→ ⊤, σ, c The CPBPV Semantics Let P be program b pre l b post in which b pre denotes the precondition, l is a list of instructions, and b post the post-condition. Let * -→ be the transitive closure of -→. The final states are specified by the set SFN (b pre , P) = { f, σ, c | i, σ ⊥ , ρ σ ⊥ b pre * -→ * f, σ, c ∧ f ∈ {⊥, ⊤} } The program violates an assertion if the set SFE (b pre , P, b post ) = { ⊥, σ, c ∈ SFN (b pre , P)} is not empty. It violates its specification if the set SFE (b pre , P, b post ) = {⊤, σ, c ∈ SFN (b pre , P) | c ∧ (ρ σ ¬b post ) satisfiable} is not empty. It is partially correct otherwise.

  See http://kind.ucd.ie/products/opensource/ESCJava2/. -CBMC is a Bounded Model Checker for ANSI-C and C++ programs. It allows for the verification of array bounds (buffer overflows), pointer safety, exceptions, and user-specified assertions. See http://www.cprover.org/cbmc/. -BLAST, the Berkeley Lazy Abstraction Software Verification Tool, is a software model checker for C programs. See http://mtc.epfl.ch/softwaretools/blast/.

-EUREKA is a C bounded model checker which uses an SMT solver instead of an SAT solver. See http://www.ai-lab.it/eureka/.

-Why is a software verification platform which integrates many existing provers (proof assistants such as Coq, PVS, HOL 4,...) and decision procedures such as Simplify, Yices, ...). See http://why.lri.fr/.

Of course, neither the expressiveness nor the objectives of all these systems are the same as the one of CPBPV. For instance, some of them can handle CTL/LTL constraints whereas CPBPV dos not yet support this kind of constraints. Nevertheless, this comparison is useful to illustrate the capabilities of CPBPV.

Table 1 .

 1 Comparison table for binary search

	1 which is the

Table 2 .

 2 Experimental Results for an Incorrect Binary Search

	CPBPV ESC/Java CBMC WHY with invariant BLAST
	length 8 0.027s	1.21 s	1.38s	NOT FOUND	UNABLE
	length 16 0.037s	1.347 s	1.69s	NOT FOUND	UNABLE
	length 32 0.064s	1.792 s	7.62s	NOT FOUND	UNABLE
	length 64 0.115s	1.886 s 27.05s	NOT FOUND	UNABLE
	length 128 0.241s	1.964 s 189.20s	NOT FOUND	UNABLE
	CPBPV ESC/Java CBMC Why BLAST
	time 0.287s	1.828s	0.82s 8.85s UNABLE

Table 3 .

 3 Experimental Results on the Tritype ProgramThe Tritype Program The tritype program is a standard benchmark in test case generation and program verification since it contains numerous non-feasible paths: only 10 paths correspond to actual inputs because of complex conditional statements in the program. The program takes three positive integers as inputs (the triangle sides) and returns 2 if the inputs correspond to an isosceles triangle, 3 if they correspond to an equilateral triangle, 1 if they correspond to some other triangle, and 4 otherwise.

Table 4 .

 4 Table 4 depicts the experimental results. Execution times correspond to the time required to find the first error. The error found with CPBPV corresponds to input values (i, j, k) = (1, 1, 2) mentioned earlier. Once again, observe the excellent behavior of CPBPV compared to the remaining tools. 5 Experimental Results for the Incorrect Tritype Program

	CPBPV ESC/Java	CBMC	WHY
	time 0.056s s 1.853s NOT FOUND NOT FOUND
		CPBPV ESC/Java CBMC EUREKA
	length 8	1.45s	3.778 s 1.11s	91s
	length 16 2.97s UNABLE 2.01s UNABLE
	length 32 UNABLE UNABLE 6.10s UNABLE
	length 64 UNABLE UNABLE 37.65s UNABLE

Table 5 .

 5 Experimental Results for Bubble Sort

We omit the domain constraints on the variables for simplicity.

See http://www.ilog.com/products.

The commercial MIP solver fails with 32-bit domains because of scaling issues.

a version with loop invariants that allows to show the correctness of this program has been written by David Cok, a developper of ESC/Java, after we contacted him.

For CBMC, we have contacted D. Kroening who has recommended to use the option CPROVER assert. If we do so, CBMC is able to find the error, but we must add

See also[START_REF] Ït-Kaci | Satisfiability Modulo Structures as Constraint Satisfaction : An Introduction[END_REF] for a study of the relations between constraint programming and Satisfiability Modulo Theories (SMT)

Acknowledgements Many thanks to Jean-Franois Couchot for many helps on the use of the Why framework.


Francesco Ragone
email: francesco.ragone@uclouvain.be

F Bouchet

Rare event algorithm study of extreme warm summers and heatwaves over Europe

The rare event algorithm increases by several orders of magnitude the number of warm summers and heatwaves sampled by the model.

• Warm summers over either France or Scandinavia are linked to wavenumber 3 hemispheric teleconnection patterns. • Warm summers in Scandinavia show bimodality due to different distribution of subsequent subseasonal heatwaves.

Introduction

Recent decades have seen a number of exceptionally warm summers and record breaking heatwaves at the Northern hemisphere midlatitudes [START_REF] Luterbacher | European seasonal and annual temperature variability, trends, and extremes since 1500[END_REF][START_REF] García-Herrera | A Review of the European Summer Heat Wave of 2003[END_REF][START_REF] Barriopedro | Redrawing the temperature record map of europe[END_REF][START_REF] Otto | Reconciling two approaches to attribution of the 2010 russian heat wave[END_REF][START_REF] Agha Kouchak | Extremes in a changing climate detection, analysis and uncertainty[END_REF][START_REF] Ipcc | Managing the risks of extreme events and disasters to advance climate change adaption: special report of the Intergovernmental Panel on Climate Change[END_REF][START_REF] Russo | Top ten european heatwaves since 1950 and their occurrence in the coming decades[END_REF][START_REF] Coumou | A decade of weather extremes[END_REF]. The Intergovernmental Panel on Climate Change (IPCC) has concluded that hot days and heavy precipitation events have become more frequent since 1950 [START_REF] Ipcc | Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF]. The precise extent to which heat extremes will become more common in the future and under the target scenarios of +1.5 o C and +2 o C global surface temperature with respect to preindustrial levels is an active field of research [START_REF] Mueller | Historically hottest summers projected to be the norm for more than half of the world's population within 20 years[END_REF][START_REF] Dosio | Extreme heat waves under 1.5 c and 2 c global warming[END_REF][START_REF] Suarez-Gutierrez | Internal variability in european summer temperatures at 1.5 c and 2 c of global warming[END_REF].

A phenomenological theory exists on the causes of heatwaves [START_REF] Schubert | Northern Eurasian Heat Waves and Droughts[END_REF][START_REF] Perkins | A review on the scientific understanding of heatwaves-Their measurement, driving mechanisms, and changes at the global scale[END_REF][START_REF] Horton | A review of recent advances in research on extreme heat events[END_REF], including large scale atmospheric circulation patterns, [START_REF] Della-Marta | Summer heat waves over western Europe 1880-2003, their relationship to large-scale forcings and predictability[END_REF][START_REF] Cassou | Tropical atlantic influence on european heat waves[END_REF][START_REF] Jezequel | Role of circulation in european heatwaves using flow analogues[END_REF], lack of precipitation and soil moisture [START_REF] Vautard | Summertime european heat and drought waves induced by wintertime mediterranean rainfall deficit[END_REF][START_REF] Zampieri | Hot European summers and the role of soil moisture in the propagation of mediterranean drought[END_REF]D'Andrea et al., 2016), and sea surface temperature anomalies [START_REF] Della-Marta | Summer heat waves over western Europe 1880-2003, their relationship to large-scale forcings and predictability[END_REF][START_REF] Cassou | Tropical atlantic influence on european heat waves[END_REF], with different levels of importance in different regions [START_REF] Stéfanon | Heatwave classification over europe and the mediterranean region[END_REF]. Long lasting heatwaves and warm summers are among the most relevant events in terms of impacts [START_REF] Camargo | Hottest summers the new normal[END_REF][START_REF] Trenberth | Framing the way to relate climate extremes to climate change[END_REF]. Although they are known to be related to persistent weather regimes [START_REF] Lau | The 2010 pakistan flood and russian heat wave: Teleconnection of hydrometeorological extremes[END_REF][START_REF] Teng | Probability of us heat waves affected by a subseasonal planetary wave pattern[END_REF][START_REF] Hoskins | Persistent extratropical regimes and climate extremes[END_REF][START_REF] Petoukhov | Role of quasi-resonant planetary wave dynamics in recent boreal spring-toautumn extreme events[END_REF][START_REF] Boers | Complex networks reveal global pattern of extreme-rainfall teleconnections[END_REF][START_REF] Kornhuber | Extreme weather events in early summer 2018 connected by a recurrent hemispheric wave-7 pattern[END_REF], their dynamics is not fully understood. Some studies show stronger increase with global warming of the frequency of persistent temperature extremes [START_REF] Pfleiderer | Summer weather becomes more persistent in a 2 c world[END_REF]. However, the response of the global circulation to climate change is extremely complex and not well understood [START_REF] Hoskins | Persistent extratropical regimes and climate extremes[END_REF], in particular for the type of rare dynamics that lead to extreme events.

The impact of climate extremes is often dominated by the rarest events. For instance, the death toll of the Western European heatwave in 2003, about 70.000 casualties, exceeded the sum of the fatalities due to all other heatwaves during the last three decades [START_REF] Barriopedro | Redrawing the temperature record map of europe[END_REF]. Moreover, extreme events with very long return times (also called return periods) always occur. For instance, an extreme event with a return time of 1000 years, has a chance 1/1000 to occur next year. Among the huge number of possible extreme events with a return time of 1000 years, a few of them will certainly occur somewhere sometimes next year. It is thus crucial to study the most extreme and rare events However, the study of such events faces strong scientific limitations: we can not rely on historical data for events with return times of 100 years or more, because often no similar events have ever been observed, and it is extremely difficult to study these events using climate models because of the required computational costs.

In this paper we address the problem of computational cost limitations. In order to face this problem, we use a rare event algorithm that concentrates ensemble simulations on trajectories of importance for the extreme events, optimizing the use of the computational resources. Rare event algorithms have recently been applied to turbulence problems [START_REF] Grafke | The instanton method and its numerical implementation in fluid mechanics[END_REF][START_REF] Laurie | Computation of rare transitions in the barotropic quasigeostrophic equations[END_REF][START_REF] Ebener | Instanton base importance sampling for rare events in stochastic pdes[END_REF][START_REF] Bouchet | Rare event algorithm links transitions in turbulent flows with activated nucleations[END_REF][START_REF] Lestang | Numerical study of extreme mechanical force exerted by a turbulent flow on a bluff body by direct and rare-event sampling techniques[END_REF] and climate applications [START_REF] Ragone | Computation of extreme heat waves in climate models using a large deviation algorithm[END_REF][START_REF] Webber | Practical rare event sampling for extreme mesoscale weather[END_REF][START_REF] Ragone | Computation of extreme values of time averaged observables in climate models with large deviation techniques[END_REF][START_REF] Plotkin | Maximizing simulated tropical cyclone intensity with action minimization[END_REF]. Based on the phenomenology of the dynamics for each family of extreme events, an appropriated type of rare event algorithm should be carefully chosen. Here we use a genealogical algorithm [START_REF] Ragone | Computation of extreme heat waves in climate models using a large deviation algorithm[END_REF][START_REF] Ragone | Computation of extreme values of time averaged observables in climate models with large deviation techniques[END_REF] adapted from [START_REF] Del Moral | Genealogical particle analysis of rare events[END_REF][START_REF] Giardina | Simulating Rare Events in Dynamical Processes[END_REF], that is efficient to study long lasting events. This rare event algorithm has been used to study heatwaves with an intermediate complexity model [START_REF] Ragone | Computation of extreme heat waves in climate models using a large deviation algorithm[END_REF] in perpetual summer conditions. Here we use it to sample rare heatwaves in CESM 1.2.2 [START_REF] Hurrell | The community earth system model: A framework for collaborative research[END_REF], a fully realistic climate model, in presence of daily and seasonal cycles.

The goal of this paper is to demonstrate the applicability of rare event algorithms to state-of-the-art climate models, and to highlight properties of extreme long lasting heatwaves that can not be assessed with traditional sampling strategies. We study extreme warm summers and heatwaves over France and Scandinavia in CESM1.2.2, in present-day climate. We compute return times orders of magnitude larger than what feasible with direct sampling, and statistically significant composite maps of dynamical quantities. Our results show that extreme warm summers are associated to wavenumber 3 teleconnection patterns in the Northern hemisphere, and suggest that the most extreme summers are related to the succession of multiple rare subseasonal heatwaves.

Data and methods

Model. All simulations are performed with the Community Earth System Model (CESM) version 1.2.2 [START_REF] Hurrell | The community earth system model: A framework for collaborative research[END_REF]. We use an atmosphere and land only setup, whose active components are the Community Atmospheric Model version 4 (CAM4) and the Community Land Model version 2 (CLM2). The model is run at statistically stationary state, with sea surface temperature (SST), sea ice cover, and the concentration of atmospheric CO2 and other greenhouse gases fixed at values representative of present day climate (year 2000). Contrary to [START_REF] Ragone | Computation of extreme heat waves in climate models using a large deviation algorithm[END_REF], the model features daily and seasonal cycle, and reproduces a fully realistic climate. See the SI for more details. We study the statistics of a 1000 years long control run, and the statistics of several sets of simulations using the rare event algorithm.

Heatwaves and physical observable definition. We consider the surface temperature T s ( r, t), where r is the space variable and t is time, and define the mean surface temperature E (T s ) ( r, t), which varies in time (because of the seasonal cycle) and in space. In practice, E (T s ) will be approximated by the climatological average of T s computed from a 1000 years long control run. We study the statistics of the spatially and temporally averaged surface temperature anomaly defined by

a(t I ) = 1 T t I +T t I A(t)dt where A(t) = 1 |D| D (T s -E (T s )) ( r, t) d r, (1)
where A(t) is the instantaneous spatial average, D is a spatial domain area, t I is the starting date of the time average, T is the averaging time, for instance several days up to one season.

Heatwaves will be defined as extreme values of the observable a(t I ), as heatwave indices developed for dynamical studies often use anomalies rather than absolute values [START_REF] Perkins | A review on the scientific understanding of heatwaves-Their measurement, driving mechanisms, and changes at the global scale[END_REF]. We are specifically interested in long events, for instance warm summers with T =3 months. For summer anomalies, we consider a JJA defined by ( 1) with t I being June 1st and t I + T August 29. In the following, we consider D as the area over either France or Scandinavia (see the SI for the definition of the areas).

Rare event algorithm experiments and importance sampling. Rare event algorithms allows a numerical model to produce very efficiently rare trajectories of importance for the study of extremes. We simulate an ensemble of N model trajectories, and at constant intervals of a resampling time τ we perform trajectory selection, pruning and cloning trajectories in order to favour those leading to the extreme events. See the SI for a more detailed description and [START_REF] Ragone | Computation of extreme heat waves in climate models using a large deviation algorithm[END_REF]. The algorithm also allows to compute probabilities for the obtained trajectories. Let X(t) be the vector of values of all the model variables at time t (which includes but is not limited to the temperature T s defined above). We consider the trajectory X(t) , with all values of X along all the interval t a ≤ t ≤ t a +T a , where t a is the starting date and T a the trajectory duration in the algorithm. We denote P k X(t) the probability distribution function of the trajectories obtained in the algorithm. Then the rare event algorithm produce trajectories distributed according to

P k X(t) = e k ta +Ta ta A(u) du Z P 0 X(t) , (2)
where P 0 X(t) is the probability distribution function of the trajectories in the model climate, Z is a normalization term computed by the algorithm, A is the selection function (with the dependance on the trajectory implicit for simplicity, A(t) = A X(t) ), and k controls the strength of the algorithm selection. Equation 2 is called an importance ratio formula. Using equation 2 we can compute the actual probabilities of the rare trajectories obtained in the simulations with the algorithm. See the SI and [START_REF] Ragone | Computation of extreme heat waves in climate models using a large deviation algorithm[END_REF] for more details.

Results

Importance sampling of extreme warm summers

The main goal of the algorithm is to perform importance sampling for the distribution of the summer temperature anomalies over the target region. Figures 1a and1b show the probability distribution functions of the summer temperature anomalies a JJA for the control run and the rare event algorithm. The control distribution shows a similar variance of about 1 o K for both France and Scandinavia. The algorithm is very effective in performing importance sampling and populating the upper tails of the distributions. The typical value of the seasonal anomaly a JJA in the rare event algorithm experiments is around 4 o K for both cases, which are values never observed in the control run. Thanks to the algorithm most of the computational power is indeed used to simulate extremely warm summers, rather than trajectories belonging to the bulk of the distribution. In the case of Scandinavia the algorithm statistics in figure 1b shows a bimodality that we will discuss in section 3.3.

In figures 1c,d we compare return times of a JJA in the control run and in the experiments with the rare event algorithm. A description of the computation of return times both for direct sampling and importance sampling using formula (2) is presented in [START_REF] Lestang | Computing return times or return periods with rare event algorithms[END_REF]. The black curves are obtained from the control run using 1000 years of data. In order to estimate uncertainty ranges, we compute also an estimate dividing the control run in K=10 samples of 100 years each, computing 10 estimates of the return times curve. We then take their average (blue curves), and compute the empirical standard deviation (blue shaded areas). The red line and shaded areas are obtained in the same way, but using the estimates from K = 10 independent experiments with the algorithm.

With the rare event algorithm we reach return times up to 10 5 -10 6 years with uncertainty ranges comparable with the ones with the control run for return times of order 10 2 years. Given the large value of k we chose, we have a small range of overlap for estimates of the return times from the control run and the rare event algorithm. When they do overlap, the values are consistent with each other within the uncertainty ranges. After about 10 6 -10 7 years the return time curves reach a plateau. Such plateaux are due to undersampling, as discussed in the SI and in [START_REF] Lestang | Computing return times or return periods with rare event algorithms[END_REF]. A full demonstration of the reliability of the results obtained with the rare event algorithm in similar simulations can be found in [START_REF] Ragone | Computation of extreme heat waves in climate models using a large deviation algorithm[END_REF][START_REF] Ragone | Computation of extreme values of time averaged observables in climate models with large deviation techniques[END_REF].

Teleconnection patterns for extreme warm summers

We study the dynamical properties of extremely warm summers with return time larger than 100 years (called 100-year warm summers or seasonal heatwaves from now on). We compute composite maps of anomalies of the local JJA surface temperature and 500 hPa geopotential height conditional on the occurrence of 100-year warm summers, for the control run and the rare event algorithm experiments. One of the key advantage of the rare event algorithm is that it gives much better results than the control run for composite statistics for large return times. Moreover the rare event algorithm gives access to composite statistics for return times larger than 1000 years, which is impossible with the control run.

Figure 2a shows composite statistics of surface temperature and 500 hPa geopotential height of 100-year warm summers over France in the control run. The map hints at the presence of a wavenumber 3 teleconnection pattern for both observables. However, how much of these patterns, which were obtained by averaging over only 10 warm summers, is statistically significant? Figure 2b shows the t-value map for the geopotential height (see the SI for details). Only the geopotential height anomalies over Europe pass a statistical significance test with |t| > 2, which means then we can not really assess the reality of the teleconnection pattern using a 1000 years control run: we just do not have enough data. Figures 2c,d show the same composite maps and statistical significance analysis for 100-year warm summers over France, but computed using the rare event algorithm. The algorithm results are globally significant, except over the northern part of the Pacific area. Transition areas between cyclonic and anticyclonic anomalies have a |t| value smaller than 2, because they have a low value of the conditional average. They are however rather narrow, so that their location can be still considered rather precise. The rare event algorithm thus allows to properly assess the existence of teleconnections of warm summers in Europe, North America and Central Asia.

The algorithm also gives a much better estimate of the amplitude of the anomalies, since it gives access to large amplitude events unavailable in the control run, because they are too rare. In the same way, the rare event algorithm gives also access to composite statistics that totally unavailable with the control run. For instance figure 3a shows composite statistics of surface temperature and 500 hPa geopotential height for 1000-year warm summers over Scandinavia, which also show a teleconnection pattern with a different spatial pattern but similar broad stroke features.

Based on the rare event algorithm results, we obtain a better overview of the characteristics of the synoptic dynamics occurring during extremely warm summers in the considered regions. A warming pattern centered over the target region is present, encompassing a larger area on a spatial scale of a few thousands km, which coincide with central-Western Europe for France and Northern-Eastern Europe for Scandinavia. Persistent anticyclonic synoptic scale structures centered over the target area are associated to this warming pattern. Locally these structures are consistent with the observed synoptic conditions for the occurrence of European heatwaves. In particular, for the North Atlantic and European area, the local patterns are very close to the observed patterns for the Western-European and the Scandinavian heatwave clusters obtained from reanalysis data in [START_REF] Stéfanon | Heatwave classification over europe and the mediterranean region[END_REF].

The application of the rare event algorithm allows to assess in a statistically robust way that these local dynamics are part of hemispheric structures of approximately wavenumber 3, which induce a temperature teleconnection pattern, with persistent regional temperature anomalies of alternating signs around the hemisphere. The patterns for the two regions are in broad strokes similar, although they differ in the exact location of positive and negative anomalies. Warm summers over France occur systematically with warm summers over Siberia and North-East America. The corresponding tripolar structure of anticyclonic anomalies is accompanied by a localized low over Central-Asia and a general lower pressure over the Arctic, with minimum over Greenland. The structure related to Scandinavian warm summers is similar on the North-Atlantic sector, but over Asia it is quite different, with a negative temperature anomaly (and related low pressure) extending from Southern Europe to the whole central Asia, with the positive anomaly constrained over the far East. The circulation over the Arctic is also different, with a strong low over the Pole.

Bimodality of warm Scandinavian summers and subseasonal heatwaves

Figure 1b a bimodality of the Scandinavia summer temperatures in the rare event algorithm data. This suggests that two types of distinct dynamical events might lead to Scandinavian extreme summers. We note that, because of the nonlinear relation (2), a bimodality in the algorithm distribution does not necessarily imply a bimodality in the tail of the model distribution. Still it is interesting to test the hypothesis of two types of dynamics.

Figure 3a shows composite statistics of surface temperature and 500 hPa geopotential height for 1000-year warm summers over Scandinavia. Figures 3b and3c show the composite maps from the rare event algorithm computed separately for a JJA < 4.2 o K and a JJA > 4.2 o K, where 4.2 o K corresponds to the local minimum between the two peaks of the distribution in red in figure 1b. While the overall structure is the same as the total composite map in figure 3a, the map of the first range (figure 3b) shows a weaker teleconnection pattern, and an anticyclonic anomaly over the North-Atlantic that is not present in the map of the second range (figure 3c). The Scandinavian warm summer dynamics is compatible with a northward shift of the jet stream along the entire hemisphere. In the case of the first range (figure 3b), this shift is less clear over the North-Atlantic, where a different dynamics seems to be in place. It is however likely that, if two distinct types of dynamics are occurring, a selection based only on the range of seasonal regional temperature fluctuations is not enough to differentiate between different dynamics, and therefore simple composite maps show from that point of view mixed information.

The study of relation between subseasonal fluctuations of surface temperature and extreme warm summers can give an explanation of this bimodality. To visualize this relation, we look at the genealogical structure of the trajectories in an experiment with the rare event algorithm for Scandinavia in figure 4a. The black lines correspond to trajectories belonging to the first range (a JJA < 4.2 o K), while the red lines to trajectories belonging to the second range (a JJA > 4.2 o K). In the experiment represented in figure 4a, trajectories belonging to each of the two ranges of the bimodal distribution coexist. In figure 4b, we show the same genealogical tree, but indicate with different colors the values of the 5-days temperature anomaly in each segment. These values are computed using equation 1 with T =5 days and t I corresponding to the date of each resampling event. In particular red segments indicate a 5day anomaly above 4.5 o K. The chosen threshold of 4.5 o K is the median of the distribution of the 5-day temperature anomalies rare event algorithm experiments for Scandinavia (shown in the SI). This value corresponds to the 96.6th percentile of the distribution of the 5day temperature anomalies in the control run. We can thus consider 5-days periods with temperature anomaly larger than this threshold as heatwave periods, and a succession of 2 or more consecutive heatwave periods as a subseasonal heatwave.

With this definitions, during warmer summers belonging to the second range (red in figure 4a), we see a first subseasonal heatwave that lasts about 20 days in late June-early July, and then a second subseasonal heatwave that last again about 15 days in August. Less extreme summers (black in figure 4a) have instead only one subseasonal heatwave in June-July. The bimodality of the algorithm statistics therefore corresponds to two qualitatively different types of warm summer, with either one or two subseasonal heatwaves. We note that the 2003 warm summer in France was characterized by two separate subseasonal heatwaves in June and August [START_REF] García-Herrera | A Review of the European Summer Heat Wave of 2003[END_REF].

Conclusions

We have shown how simulations with the algorithm produce hundreds of times more extremes than a control run for the same computational cost, and allow us to estimate return times of extreme events orders of magnitude larger. This allows to compute precisely statistically significant composite maps of dynamical quantities for warm summers with extremely strong seasonal surface temperature anomalies. In this way we are able to identify rigorously the occurrence of persistent teleconnection patterns of wavenumber 3 during very warm summers. These patterns and the related dynamics can not be properly studied with direct sampling, as the corresponding local anomalies of temperature and geopotential height are statistically significant only over the target region. The improved statistics given by the rare event algorithm instead allows to obtain statistically significant maps and a very clear signal also away from the target region.

Teleconnection patterns similar to what we have obtained here were observed in the first application of the rare event algorithm to a climate model [START_REF] Ragone | Computation of extreme heat waves in climate models using a large deviation algorithm[END_REF]. In that case the model was an intermediate complexity model, run in perpetual summer setup, and the target area was the whole Europe. In this case we use a much more realistic model with seasonal cycle, calibrated to reproduce present day climate, and we target areas consistent with observed heatwave clusters [START_REF] Stéfanon | Heatwave classification over europe and the mediterranean region[END_REF] and recent cases of very intense heatwaves or warm summer, like the heatwave over France of 2003 and the warm summer over Northern Europe of 2018. It is striking that we obtain similar qualitative results. This supports the idea that these teleconnections at seasonal scale and the corresponding wavenumber 3 dynamics are robust features, mainly dynamical and little affected by the details of the physical parameterizations.

Comparison with observations is not straightforward, due to the lack of data and the fact that most studies have focused on events at subseasonal scale. Several authors have recently highlighted the role of Rossby waves in determining teleconnection of extreme events [START_REF] Schubert | Warm season subseasonal variability and climate extremes in the northern hemisphere: the role of stationary rossby waves[END_REF][START_REF] Lau | The 2010 pakistan flood and russian heat wave: Teleconnection of hydrometeorological extremes[END_REF][START_REF] Petoukhov | Quasiresonant amplification of planetary waves and recent northern hemisphere weather extremes[END_REF][START_REF] Petoukhov | Role of quasi-resonant planetary wave dynamics in recent boreal spring-toautumn extreme events[END_REF], with particular emphasis on strong heatwaves, e.g. cases over the central USA [START_REF] Teng | Probability of us heat waves affected by a subseasonal planetary wave pattern[END_REF], Alberta [START_REF] Petoukhov | Alberta wildfire 2016: Apt contribution from anomalous planetary wave dynamics[END_REF] and Western Europe [START_REF] Kornhuber | Extreme weather events in early summer 2018 connected by a recurrent hemispheric wave-7 pattern[END_REF]. These studies typically find patterns with wavenumber 5 to 7 and analyse subseasonal temperature fluctuations. Several different detection methods are used in the literature to identify the dynamics leading to heatwaves, including projections of heatwave states on typical variability patterns obtained with cluster analysis [START_REF] Cassou | Tropical atlantic influence on european heat waves[END_REF] or empirical orthogonal functions [START_REF] Teng | Probability of us heat waves affected by a subseasonal planetary wave pattern[END_REF], spectral analysis [START_REF] Petoukhov | Alberta wildfire 2016: Apt contribution from anomalous planetary wave dynamics[END_REF][START_REF] Kornhuber | Extreme weather events in early summer 2018 connected by a recurrent hemispheric wave-7 pattern[END_REF] and/or indicators based on resonance models [START_REF] Petoukhov | Quasiresonant amplification of planetary waves and recent northern hemisphere weather extremes[END_REF][START_REF] Petoukhov | Role of quasi-resonant planetary wave dynamics in recent boreal spring-toautumn extreme events[END_REF]. Merging these approaches with rare event simulations should lead to very promising future studies, aimed at investigating the relation between the seasonal scale structures we find here and the possible role of higher wavenumber stationary Rossby waves proposed in the literature, as well as the possible multimodality in the way of occurrence of some extreme events.

Figure 1 .Figure 2 .Figure 3 .Figure 4 .

 1234 Figure 1. Distribution of seasonal JJA temperatures anomalies averaged over France (a) and Scandinavia (b) from the control run (black) and the rare event algorithm (red). Return times for the seasonal temperature anomalies for France (c) and Scandinavia (d), from the control run (black) and the rare event algorithm (red). The shaded areas, light blue for the control run and light red for the algorithm, correspond to 1 standard deviation of the sample used to compute the estimate (see SI).

We use as selection function A the same function used to define a heatwave in formula (1). We see from equation (2) that if the parameter k is positive, then trajectories with large values of the time average of the surface temperature anomaly will be much more probable in simulations obtained with the algorithm rather than in simple simulations with the model. The larger k, the stronger the selection, and the more probable trajectories with extreme values of the time averaged surface temperature. We analyse extremely warm summers over France and Scandinavia. For each case, we perform K = 10 ensemble simulations with the algorithm, each with N =100 trajectories, biasing parameter k=30, with T a =90 days from June 1st to August 29th, and resampling time τ =5 days (see the SI). A 1000 years long control run is used to provide initial conditions for the experiments with the algorithm, and as a benchmark for the statistics. The computational cost of the experiments with the algorithm is equivalent to simulating 1000 summers in the control run, but they allow to gather a much richer statistics for the extreme events of interest.

Acknowledgments

This work was granted access to the HPC resources of CINES under the allocation 2019-A0070110575 made by GENCI. The computation of this work were partially performed on the PSMN platform of ENS de Lyon. This work has received funding through the ACADEMICS grant of the IDEXLYON, project of the Université de Lyon, PIA operated by ANR-16-IDEX-0005. The research leading to these results has received funding from the European Research Council under the European Union's seventh Framework Programme (FP7/2007-2013 Grant Agreement No. 616811). The data related to the figures presented in this paper have been uploaded with the submission for review purposes. We are discussing with CINES the long term storage and open access of the full dataset of our simulations according to the Enabling FAIR data Project guidelines.


Aurélien Froger Nicholas D Kullman

Jorge E Mendoza

Justin C Goodson

Aurelien Froger

frvcpy: An Open-Source Solver for the Fixed Route Vehicle Charging Problem

HAL is

Introduction

Governmental regulations as well as a growing population of environmentally conscious consumers have led to increased pressure for firms to act sustainably. This pressure is particularly high in the logistics domain, which accounts for about one third of emissions in the United States (Office of Transportation and Air Quality 2019). Electric vehicles (EVs) offer a means to more sustainable transportation; however, they present technical challenges to which their conventional (i.e., internal-combustion engine) vehicle (CV) counterparts are immune. For example, because the distance EVs can travel on a single charge is often less than the distance an equivalent CV can travel on a tank of gas, EVs may demand more frequent recharging operations. This difficulty is compounded by the sparseness of EV recharging infrastructure relative to the network of refueling stations available to CVs, potentially forcing EVs to perform longer detours to recharge their batteries. Further, despite recent improvements to battery and charging station (CS) technology, recharging an EV still requires orders of magnitude more time than refueling a CV. The time required to recharge is also nonlinear with respect to the EV's state of charge (SoC), the relative amount of energy left in its battery, posing yet another challenge not applicable to CVs [START_REF] Uhrig | E-mobility in car parks-guidelines for charging infrastructure expansion planning and operation based on stochastic simulations[END_REF]. Companies choosing to adopt EVs require optimization tools capable of handling these additional challenges.

The development of such optimization tools for conventional vehicles has commanded significant attention from the operations research (OR) community in the study of vehicle routing problems (VRPs).

The incorporation of additional constraints that address the challenges posed by EVs has marked a new family of problems within VRPs known as electric vehicle routing problems (E-VRPs). One of the primary tasks in solving E-VRPs is making good charging decisions -namely, where to recharge and how long to do so. This is the crux of the fixed route vehicle charging problem (FRVCP), in which charging operations must be inserted into a fixed sequence of customers being visited by an EV so as to minimize the time for the EV to reach the end of the sequence in an energy-feasible manner. The FRVCP naturally arises as a subproblem in many E-VRPs, since its solution is required in order to determine the true duration or cost of a given route. Having a capable solution method for the FRVCP is thus crucial to the advancement of E-VRP research.

While FRVCP solution methods exist, they tend to suffer from one or more of the following issues: inexactness (e.g., heuristic methods that provide suboptimal solutions), inefficiency (e.g., mixed-integer programs (MIPs) whose solutions require significant run time), or a lack of robustness (e.g., exact algorithms that are limited to simplified versions of the FRVCP). With this work, we offer an implementation of a solution method that suffers from none of these issues. It is based on the labeling algorithm proposed in [START_REF] Froger | Improved formulations and algorithmic components for the electric vehicle routing problem with nonlinear charging functions[END_REF], which, though capable, is notoriously difficult to implement.

In an attempt to remove the burden of implementation for future E-VRP researchers, we offer our implementation in an open-source Python package, frvcpy (a portmanteau of FRVCP and Python; Kullman et al. 2020a). frvcpy is designed to be easily embedded in more complex solution schemes for E-VRPs (such as in a (meta)heuristic or Benders decomposition): it requires minimal dependencies and inputs, can be accessed either via the command line or a Python API, and includes a translator to generate the required inputs from a common instance format in the VRP community (VRP-REP [START_REF] Mendoza | VRP-REP: the vehicle routing community repository[END_REF])). Our aim with frvcpy is to make it easier to solve E-VRPs, thereby stimulating additional research in this field which promises to bring about more sustainable practices in logistics.

The remainder of the paper is organized as follows. We first define the FRVCP in §2, then discuss some of the previous work on FRVCPs in §3. In §4 we give an overview of the algorithm implemented in frvcpy, then describe the package itself in §5. We conclude with brief comments in §6.

Defining the FRVCP

We consider an EV with a fixed route Π = (π 1 , . . . , π R ) that begins at some node π 1 (usually the depot), has a sequence of stops at other nodes (π i ) R-1 i=2 (customers to visit), and terminates at some node π R (also usually the depot). The vehicle begins at π 1 with its battery at some initial energy level q 0 , often taken to be equal to its maximum battery capacity Q. Let the set of nodes in the route be I = {π 1 , . . . , π R }.

We also consider a set of charging stations C at which the EV may recharge between stops in Π. Each CS c ∈ C has some charging type (e.g., fast, slow) associated with a strictly increasing piecewise-linear concave charging function Φ c (t) specifying, for an empty battery, the resulting energy after charging for time t at CS c (see Figure 3). We denote its inverse by Φ -1 c (q), equal to the time required to charge from empty to energy level q, and we define the related function ūc (q 1 , q 2 ) = Φ -1 c (q 2 ) -Φ -1 c (q 1 ) to be the time required to charge from q 1 to q 2 at CS c. Let the set of breakpoints defining the charging function of

CS c be B c , where a breakpoint b i ∈ B c is a (time, charge) pair: (b t i , b q i ).
When the EV travels between nodes i, j ∈ I ∪ C, it incurs some known travel time t ij and energy consumption e ij (we assume the triangle inequality holds for both). At stops in Π, the EV may also incur some processing time (e.g., waiting). The objective of the FRVCP is to determine charging decisions -how much to recharge, at which CSs, between which stops in Π -that minimize the total time for the EV to traverse the route,

Related Literature

FRVCPs fall under the category of EV routing problems, which are themselves part of a larger body of research on VRPs. We focus our review here solely on FRVCPs. For an overview of E-VRPs we refer the reader to [START_REF] Pelletier | 50th anniversary invited article-goods distribution with electric vehicles: review and research perspectives[END_REF], and similarly to [START_REF] Braekers | The vehicle routing problem: State of the art classification and review[END_REF] for an overview of VRPs. [START_REF] Montoya | A multi-space sampling heuristic for the green vehicle routing problem[END_REF] encounter an FRVCP in their work on a green vehicle routing problem. The FRVCP they consider assumes that vehicles may only visit one CS between stops, that they always fully restore their energy when recharging (that is, they follow a "full recharging strategy"), and that doing so requires constant time. To solve this FRVCP, they offer an exact algorithm. [START_REF] Roberti | The electric traveling salesman problem with time windows[END_REF] address an FRVCP in their work on the E-VRP with time windows (E-VRP-TW) and also offer an algorithm that solves this FRVCP exactly. Their solution accommodates a partial recharging policy, assuming that the time required to recharge is linear with the amount of energy. However, unlike in [START_REF] Montoya | A multi-space sampling heuristic for the green vehicle routing problem[END_REF], they assume that the network of CSs is homogeneous; that is, that all CSs have the same charging technology. The FRVCP again arises in related works by [START_REF] Hiermann | The electric fleet size and mix vehicle routing problem with time windows and recharging stations[END_REF][START_REF] Schiffer | An adaptive large neighborhood search for the location-routing problem with intra-route facilities[END_REF][START_REF] Hiermann | Routing a mix of conventional, plug-in hybrid, and electric vehicles[END_REF]. These studies offer exact algorithms for the FRVCP under the assumption that at most one CS may be visited between stops, that the CSs are homogeneous, and that recharging requires linear time. [START_REF] Hiermann | The electric fleet size and mix vehicle routing problem with time windows and recharging stations[END_REF] additionally assume a full recharging strategy while [START_REF] Schiffer | An adaptive large neighborhood search for the location-routing problem with intra-route facilities[END_REF] and Hiermann et al. ( 2019) allow partial recharging. [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] then consider the first FRVCP that accommodates realistic (nonlinear) recharging times. In the study, they also demonstrate that the assumption of linear recharging times can lead to infeasible or suboptimal solutions. Their FRVCP allows for partial recharging and heterogeneous CSs but assumes that at most one CS may be inserted between stops. To solve their FRVCP, Montoya et al.

offer both a heuristic and a MIP formulation. [START_REF] Koç | The electric vehicle routing problem with shared charging stations[END_REF] adopt the heuristic and MIP formulations from [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] to solve a similar FRVCP that arises in their work on the E-VRP with shared CSs and nonlinear charging. [START_REF] Baum | Shortest feasible paths with charging stops for battery electric vehicles[END_REF] then offer a labeling algorithm to solve an FRVCP on real road networks that also accommodates realistic recharging times and allows for multiple CS insertions, although it is restricted to the special case where the route length is two (an origin-destination (OD) pair).

Finally, [START_REF] Froger | Improved formulations and algorithmic components for the electric vehicle routing problem with nonlinear charging functions[END_REF] propose an exact labeling algorithm to solve the FRVCP from [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF]. Their algorithm is not restricted to OD pairs, and it additionally allows the EV to visit multiple CSs between stops in the route, making theirs the richest of the aforementioned FRVCP variants.

Over a testbed of nearly 30,000 instances, they compare their labeling algorithm against a heuristic and a commercial solver for a MIP formulation. They find that the labeling algorithm matches the optimality of the MIP with a runtime comparable to the heuristic. The algorithm is thus state of the art for solving FRVCPs. However, the authors note that its performance is not without cost. They state that E-VRP researchers may ultimately prefer to adopt the heuristic solution, despite its inferior performance, given the complexity of implementing the labeling algorithm. Here, we offer an implementation of the algorithm in frvcpy in an attempt to ensure that its complexity does not prevent its adoption.

The algorithm from [START_REF] Froger | Improved formulations and algorithmic components for the electric vehicle routing problem with nonlinear charging functions[END_REF] has also been adapted to accommodate FRVCPs with additional constraints. For example, in work on a stochastic E-VRP with public CSs, [START_REF] Nicholas D Kullman | Electric Vehicle Routing with Public Charging Stations[END_REF] adapt the algorithm to accommodate an FRVCP with discrete charging decisions and time-dependent waiting times at CSs. Similarly, Kullman et al. (2020b) adapt an early version of frvcpy to accommo- date an FRVCP with customer time windows. In both cases, the algorithm's speed and exactness were required as it was called repeatedly to solve a subproblem in a Benders-based branch-and-cut procedure.

Overview of Labeling Algorithm from Froger et al. (2019)

Given the algorithm's complexity, we provide here a cursory overview and refer the reader to Froger et al.

(2019) for additional details (see, in particular, their discussion of Algorithm 3 in §5.3 and Appendix E).

To find the optimal charging decisions for a given route Π, the FRVCP is reformulated as a resourceconstrained shortest path problem. The algorithm then works by setting labels at nodes on a modified graph reflecting the vehicle's possible movements along Π (Figure 1, Inset 1). Labels are defined by SoC functions -piecewise-linear functions comprised of supporting points z = (z t , z q ) that describe a state of departure from a node in terms of time z t and charge (SoC) z q .

During the algorithm's execution, labels are extended along nodes in the graph. When a label is extended to a CS node c, we create new supporting points for each breakpoint in B c to which we could charge (that is, breakpoints with a higher energy than that with which we arrived). Consider Figure 2 Figure 2 shows this explicitly for the new supporting point z 4 , corresponding to the decision to recharge to the breakpoint b 4 (more specifically, to b q 4 ). We continue to extend labels along nodes in the graph until the destination node π R is reached, whereat the algorithm returns the time of the first supporting point in the label's SoC function. Bounds on energy and time are established in pre-processing and are used alongside dominance rules during the algorithm's execution to improve its efficiency. a black square for the initial supporting point (z t 1 , z q 1 ). We create additional supporting points (z 2 , z 3 , z 4 ) for each breakpoint to which we could charge, shown by circles b 2 , b 3 , and b 4 in the CS's charging function (left). Axis labels for new supporting point z 4 (right) detail its creation from the decision to recharge to breakpoint b 4 .

The frvcpy Package

frvcpy is an open-source Python-based implementation of the labeling algorithm from [START_REF] Froger | Improved formulations and algorithmic components for the electric vehicle routing problem with nonlinear charging functions[END_REF] for solving the FRVCP. In this section we give an overview of its structure ( §5.1), demonstrate its usage ( §5.2), and briefly comment on its performance ( §5.3). The code accompanying this section can be found at Kullman et al. (2020a).

Structure

frvcpy is a small package (approximately 1000 lines of code) built in the Python 3 programming language; it is available on the Python Package Index and can be installed via "pip install frvcpy." It is comprised of three primary modules: core.py, solver.py, and algorithm.py. The core.py module consists of class definitions for ancillary objects required in the algorithm's execution such as nodes, labels, and the FRVCP problem instance. Module solver.py defines the user-facing Solver class which is responsible for pre-processing, calling the algorithm, and writing solutions to file. The algorithm itself and its accompanying functions are contained in algorithm.py. Additionally, the package contains the module translator.py, which provides an instance translator that allows for interoperability between frvcpy and VRP-REP instances. VRP-REP is a community-driven repository for vehicle routing problem data files; see [START_REF] Mendoza | VRP-REP: the vehicle routing community repository[END_REF] for more details.

Input/output. Users can interact with frvcpy using a Python API or via the command-line interface (CLI). As input, frvcpy requires a compatible instance, the fixed route for the EV to travel, and the EV's initial energy. Compatible instances are JSON files (or equivalent Python dictionaries) following the schema available in Kullman et al. (2020a). After execution, the algorithm returns the optimal energy-feasible route and its duration. The returned route is a list of tuples indicating stops' node IDs and the amount of energy to be recharged there (the latter given by the keyword 'None' for non-CS nodes; see Listing 1, Line 24).

Testing the installation. frvcpy provides simple testing to determine if its installation was successful. From the command line, users can execute the command frvcpy-test to run a suite of tests that performs an instance translation and solves 134 FRVCPs from the [START_REF] Froger | Improved formulations and algorithmic components for the electric vehicle routing problem with nonlinear charging functions[END_REF]

Example Usage

We provide an example demonstrating the use of frvcpy through the Python API. Consider a user with the VRP-REP-compliant instance "vrprep-instance.xml," depicted in Figure 4. The instance contains "fast," "normal," and "slow" CSs whose charging functions are shown in Figure 3. An EV, which begins at the depot with full battery, has been assigned the fixed route Π = (0, 40,12,33,38,16,0), depicted by the gray arrows in Figure 4. Because the EV does not have sufficient energy to traverse Π without recharging, we solve an FRVCP to determine the optimal insertion of charging operations. We can do this using frvcpy as follows: [0,40,12,33,38,16,0] # route to make energy feasible The solution to the FRVCP instructs the EV to recharge at CS 48 between customers 33 and 38, as depicted by the black arrows in Figure 4 and the printed output on line 24 in Listing 1. This results in a total route duration of about 7.34 hours. We note that detouring to CS 48 actually requires more travel time than detouring to CS 41; however, given that CS 48 offers a faster charging speed, it is ultimately preferred over CS 41 (recharging at CS 41 instead of 48 results in an objective of 7.44 hrs).

The above example would be accomplished with the CLI via

Performance

We test the performance of frvcpy over the nearly 30,000 instances comprising the testbed from [START_REF] Froger | Improved formulations and algorithmic components for the electric vehicle routing problem with nonlinear charging functions[END_REF]. These instances have a median route length of 10 stops and a median of 18 CSs that may be inserted. Over the tests, the algorithm has an average run time of 5.6 ms. In addition, in the same tests we find that our translator can translate instances from VRP-REP format in an average of 0.1 s.

These results suggest that frvcpy requires sufficiently low runtime so as to be included in larger solution schemes for E-VRPs.

Conclusion

We introduced frvcpy, a Python-based open-source implementation of the labeling algorithm from [START_REF] Froger | Improved formulations and algorithmic components for the electric vehicle routing problem with nonlinear charging functions[END_REF] for the fixed route vehicle charging problem. The algorithm and our implementation are flexible, able to accommodate realistic problem features such as non-linear recharging times, partial charging decisions, and heterogeneous charging station technologies. Because FRVCPs are often encountered as subproblems of more general EV routing problems, we designed frvcpy to be easily embedded in larger solution schemes. To that end, the package offers two modes of interaction, has minimal requirements, and is computationally efficient. Our hope is that frvcpy facilitates the solution of E-VRPs, lowering the barrier to entry in this field, and ultimately helping bring about a faster transition to more sustainable transportation practices.

Figure 1 :

 1 Figure 1: Instance excerpt depicting a fixed route in the original problem graph (left) and the modified problem graph for the FRVCP (Inset 1), which includes dummy nodes for CS insertions.

  , which depicts this process when extending a label along the edge from customer 33 to CS 48. When it arrives at CS 48, its SoC function has only one supporting point z 1 (assuming the EV has not yet stopped to recharge) depicted by the black square in the right graph of Figure 2. Then for each breakpoint in the CS's charging function to which the EV could recharge (b 2 , b 3 , b 4 ), we add a supporting point to the label's SoC function (z 2 , z 3 , z 4 ) whose time and charge reflect the decision to charge to that breakpoint.

Figure 2 :

 2 Figure 2: Depicting the creation of new supporting points at CS nodes for the case of CS 48 between customers 33 and 38 in Figure 1. Right shows the SoC function of the label extended to CS 48, with

8Figure 3 :

 3 Figure 3: Piecewise linear charging functions for example instance "vrprep-instance.xml."

Listing 2 :

 2 Example frvcpy usage with CLI # translate existing VRP-REP instance, write it to file frvcpy-translate ./instances/vrprep-instance.xml new-frvcp-instance.json frvcpy --instance=new-frvcp-instance.json --route=0,40,12,33,38,16,0 --qinit=16000 --output=my-solution.xml # Duration: 7.339 # Energy-feasible route: # [(0, None), (40, None), (12, None), (33, None), (48, 6673.379615520617), (38, None), (16, None), (0, None)]

Figure 4 :

 4 Figure 4: Depiction of example instance "vrprep-instance.xml." We consider an EV given the route shown by the gray arrows. The solution to the FRVCP for this route instructs the EV to recharge at CS 48 between customers 33 and 38 (black arrows).

  testbed. The same test suite can also be run in Python via

	import frvcpy.test
	frvcpy.test.runAll()

Acknowledgements

Nicholas Kullman is grateful for support from the Société des Professeurs Français et Francophones d'Amérique. This research was also partly funded by the French Agence Nationale de la Recherche through project e-VRO (ANR-15-CE22-0005-01). Justin Goodson wishes to express appreciation for the support from the Center for Supply Chain Excellence at the Richard A. Chaifetz School of Business.