Eduardo Queiroga
email: eduardoqueiroga@id.uff.br

Ruslan Sadykov
email: ruslan.sadykov@inria.fr

Eduardo Uchoa
email: uchoa@producao.uff.br

A modern POPMUSIC matheuristic for the capacitated vehicle routing problem

This work proposes a partial optimization metaheuristic under special intensification conditions (POPMUSIC) for the classical capacitated vehicle routing problem (CVRP). The proposed approach uses a branch-cut-and-price algorithm as a powerful heuristic to solve subproblems whose dimensions are typically between 25 and 200 customers. The whole algorithm can be seen as the application of local search over very large neighborhoods, starting from a single initial solution. The main computational experiments were carried out on instances having between 302 and 1000 customers. Using initial solutions generated by some of the best available metaheuristics for the problem, POPMUSIC was able to obtain consistently better solutions for long runs of up to 32 hours. In a final experiment, starting from the best known solutions available in CVRP library (CVRPLIB), POPMUSIC was able to find new best solutions for several instances, including some very large ones.

Introduction

The Capacitated Vehicle Routing Problem (CVRP), introduced by [START_REF] Dantzig | The truck dispatching problem[END_REF], is one of the most widely studied problems in combinatorial optimization and operations research. The CVRP is the prototypical vehicle routing problem. New ideas are often first proposed and tested on CVRP and then generalized to other routing variants. It can be defined as follows. Let G = (V, E) be a complete undirected graph, such that V = {0, 1, . . . , n} is the set of vertices and E is the set of edges, where vertex 0 represents a depot and V + = {1, . . . , n} a set of customers. There is a non-negative cost c ij for each edge {i, j} ∈ E and a demand d i for each customer i ∈ V + . The vehicle capacity is denoted by Q. A route is a path that begins and ends at the depot. A solution consists of a set of routes that respect the following constraints: (i) each customer must be visited exactly once by one of the routes; (ii) the sum of the customer demands in a route can not exceed the vehicle capacity. The objective is to find a set of routes with the minimum total cost.

Given that CVRP is NP-hard, most of the algorithms proposed for this problem are heuristics and metaheuristics [START_REF] Laporte | Chapter 4: Heuristics for the Vehicle Routing Problem[END_REF]. The best performing published algorithms are: the iterated local search with set partitioning (ILS-SP) [START_REF] Subramanian | A hybrid algorithm for a class of vehicle routing problems[END_REF], knowledge-guided local search (KGLS) (Arnold and Sörensen, 2019), hybrid genetic search (HGS) [START_REF] Vidal | A hybrid genetic algorithm for multidepot and periodic vehicle routing problems[END_REF], slack induction by string removals (SISR) [START_REF] Christiaens | Slack induction by string removals for vehicle routing problems[END_REF], and fast ILS localized optimization (FILO) [START_REF] Luca | A fast and scalable heuristic for the solution of large-scale capacitated vehicle routing problems[END_REF]. ILS-SP combines the wellknown ILS [START_REF] Ramalhinho | Iterated local search: Framework and applications[END_REF] with a set partitioning (SP) model. The SP model attempts to build unexplored solutions from the set of routes associated with the local minima found by previous runs of the local search. KGLS presents an efficient GLS with three complementary operators using ideas from sequential search and pruning, as well as a problem-specific knowledge to penalize "bad" edges. HGS is a population-based evolutionary search that also makes use of local search (in a step called education) and a sophisticated mechanism for controlling population diversity. Among the key components of HGS, we can mention the management of a subpopulation with infeasible solutions, as well as the individual evaluation (a.k.a. fitness) driven by the solution cost and its contribution to population diversity. SIRS is a ruin & recreate local search guided by simulated annealing (SA) [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]. The ruin procedure removes strings (sequence of consecutive customers) from routes (inducing a capacity slack ), whereas the recreate procedure reinserts the removed customers in the ruined solution in a greedy manner. Finally, FILO is a scalable metaheuristic that employs novel and existing acceleration techniques during the main iterative part based on ILS, whereas it uses an SA-based acceptance criterion to get a continuous diversification.

On the other hand, the exact methods for CVRP have advanced considerably in recent years [START_REF] Poggi | Chapter 3: New exact algorithms for the capacitated vehicle routing problem[END_REF][START_REF] Costa | Exact branch-price-and-cut algorithms for vehicle routing[END_REF]. The state-of-the-art results are achieved by branchcut-and-price algorithms [START_REF] Pecin | Improved branch-cut-and-price for capacitated vehicle routing[END_REF](Pecin et al., , 2017a;;[START_REF] Pessoa | A generic exact solver for vehicle routing and related problems[END_REF], which combine column and cut generation with several additional mechanisms. According to the experiments carried out in [START_REF] Uchoa | New benchmark instances for the capacitated vehicle routing problem[END_REF], this type of algorithm is able to produce optimal solutions for almost all instances with up to 250 customers, and in some cases, it can solve even larger instances (the largest one already solved has 654 customers). An important observation on the behavior of modern branch-cut-and-price algorithms for CVRP, explored in this work, is the following: while instances with more than 200 customers usually take hours or even days to be solved, many instances with up to 150 customers can be solved in few minutes, and many instances with up to 100 customers can be solved in seconds.

The algorithms that hybridize metaheuristics with mathematical programming approaches [START_REF] Jourdan | Hybridizing exact methods and metaheuristics: A taxonomy[END_REF] are often known as matheuristics. Such methods have already been proposed for several optimization problems, including vehicle routing [START_REF] Archetti | A survey on matheuristics for routing problems[END_REF][START_REF] Leggieri | A matheuristic for the asymmetric capacitated vehicle routing problem[END_REF]. According to [START_REF] Archetti | A survey on matheuristics for routing problems[END_REF], one of the types of matheuristics is based on the decomposition of the original problem into smaller subproblems that can be solved (optimally or sub-optimally) through mathematical programming models. This work proposes a simple Partial OPtimization Metaheuristic Under Special Intensification Conditions (POPMUSIC) [START_REF] Éric | Popmusic -Partial Optimization Metaheuristic under Special Intensification Conditions[END_REF] for the CVRP that uses a modern branchcut-and-price algorithm to solve subproblems (exactly or heuristically). The general idea of POPMUSIC is to optimize subproblems, defined by parts of a solution until a local minimum is reached. This type of algorithm has been shown to be effective for different problems [START_REF] Éric | Popmusic -Partial Optimization Metaheuristic under Special Intensification Conditions[END_REF], including vehicle routing variants [START_REF] Ostertag | Popmusic for a real-world large-scale vehicle routing problem with time windows[END_REF][START_REF] Lalla | A POPMUSIC approach for the multi-depot cumulative capacitated vehicle routing problem[END_REF] and the famous traveling salesman problem [START_REF] Taillard | POPMUSIC for the travelling salesman problem[END_REF].

The remainder of this paper is organized as follows. In Section 2, the proposed POPMUSIC matheuristic for the CVRP is presented. Section 3 describes the modifications to the published branch-cut-and-price algorithm used for solving the subproblems. Section 4 presents and analyses the results of extensive computational experiments. Finally, in Section 5, the final conclusions are presented, as well as suggestions for future work.

A POPMUSIC matheuristic for the CVRP

Algorithm 1 shows the pseudocode of the proposed POPMUSIC matheuristic for the CVRP, which has four inputs: (i) an initial solution S; (ii) an algorithm A to solve subproblems;

(iii) initial value α for the current target dimension dim sp (upper limit on the dimension of subproblems); (iv) step size δ to increase dim sp . The algorithm's output is a (possibly) improved solution S obtained after solving a sequence of subproblems. A solution S is a set {r 1 , . . . , r m } of m routes, whereas the set of customers visited by a route r is denoted by C(r). A set V sp ⊆ V + represents the CVRP subproblem associated with the subgraph G[{0} ∪ V sp ]. We will refer to solutions for subproblems as subsolutions. In addition, in the description of the algorithm, we will consider that c ji = c ij , ∀{i, j} ∈ E, and c ii = 0, ∀i ∈ V + .

Algorithm 1: A POPMUSIC matheuristic for the CVRP The algorithm keeps the current target dimension dim sp , which is the upper limit for |V sp |. dim sp is initialized to α at line 4. The set Π, initialized at line 5, keeps all subproblems already explored during the search together with their subsolutions. Formally, Π is a set of all pairs (V , S ), such that subproblem with set V ⊂ V + of vertices is already solved, and S is its subsolution. At first, a random permutation of the customers in V + produces the array L (line 7). For a given value of dim sp , each customer i ∈ V + is used as a seed to construct a subproblem V sp at lines 11-19. A subproblem V sp with at most dim sp customers is constructed iteratively by including the routes in the current solution S that are closest to vertex i. The distance from i to each route r ∈ S is determined by the smallest cost of edges connecting i and the vertices in r (line 14). The routes already included in V sp are stored in R to avoid repetitions (line 17). The construction of subproblem V sp is finished when the next selected route r cannot be included in subproblem due to the upper limit dim sp on the subprobem dimension (line 19). Figure 1a illustrates the construction of the subproblem for an instance having 109 customers and the current target dimension dim sp = 30. First, the route containing the seed (in black) is included in the subproblem. The second selected route is the red one, while the third is the blue one, and the fourth is the purple one. Adding a fifth route would exceed dim sp , so the obtained subproblem has 24 customers.

1 Data: V, E, c, d, Q 2 Input parameters: initial solution S, algorithm A, α, δ 3 Output: (Possibly) improved solution S 4 dimsp ← α 5 Π ← ∅ 6 while
(a) Initial solution and a constructed subproblem. Seed customer is marked in black.

(b) Improved solution after finding a better subsolution Figure 1: Constructing and solving a subproblem. Depot is the yellow square, and customers are circles with diameter proportional to its demand. For the sake of visualization, the edges adjacent to the depot are not depicted.

The algorithm solves generated subproblem V sp only if it is neither equal nor contained in any subproblem V already solved before (line 21). Indeed, the Π-based condition avoids wasting time on subproblems, i.e., current subsolutions of which are unlikely to be improved because the same or a larger subproblem has been solved already. The solved subproblems together with their solutions are added to set Π at line 23. At line 24, the algorithm A tries to improve the subsolution S sp of S for current subproblem V sp . As algorithm A, we use a branch-cutand-price based heuristic described in Section 3. It is important here to use the cost of the known solution S sp for subproblem V sp to improve the performance of the branch-cut-and-price algorithm. Finally, if the solution S sp found by A is better than S sp , then S is updated, and the search is restarted for the same target dimension dim sp : all customers will be used again as seeds without increasing dim sp . Figure 1b depicts an example of such an improved solution. If all seeds fail to produce an improving subsolution, then the target dimension dim sp is increased by δ, so that larger subproblems can be explored (line 31). The algorithm is interrupted when the time limit is reached or when the target dimension exceeds the number of customers (line 6). From now on, we refer to Algorithm 1 as POP.

A branch-cut-and-price heuristic to solve subproblems

The algorithm A in POP, used for solving the subproblems, is an adaptation of the generic Branch-Cut-and-Price (BCP) algorithm proposed by [START_REF] Pessoa | A generic exact solver for vehicle routing and related problems[END_REF], which is a state-ofthe-art exact algorithm for many VRP variants, including the CVRP. BCP is a well-known technique that incorporates column and cut generation in a branch-and-bound procedure. In particular, the BCP by [START_REF] Pessoa | A generic exact solver for vehicle routing and related problems[END_REF] includes advanced elements, such as: (i) ng-path relaxation [START_REF] Baldacci | New route relaxation and pricing strategies for the vehicle routing problem[END_REF]; (ii) rank-1 cuts with limited memory [START_REF] Jepsen | Subset-row inequalities applied to the vehicle-routing problem with time windows[END_REF][START_REF] Pecin | Improved branch-cut-and-price for capacitated vehicle routing[END_REF]Pecin et al., , 2017b;;[START_REF] Bulhões | On the complete set packing and set partitioning polytopes: Properties and rank 1 facets[END_REF]; (iii) path enumeration [START_REF] Baldacci | An exact algorithm for the vehicle routing problem based on the set partitioning formulation with additional cuts[END_REF][START_REF] Contardo | A new exact algorithm for the multi-depot vehicle routing problem under capacity and route length constraints[END_REF]; (iv) rounded capacity cuts [START_REF] Laporte | A branch and bound algorithm for the capacitated vehicle routing problem[END_REF]; (v) bucket graph based bi-directional labeling algorithm [START_REF] Sadykov | A bucket graph based labeling algorithm with application to vehicle routing[END_REF]; (vi) edge elimination based on reduced costs [START_REF] Irnich | Path-reduced costs for eliminating arcs in routing and scheduling[END_REF][START_REF] Sadykov | A bucket graph based labeling algorithm with application to vehicle routing[END_REF]. The reader is referred to [START_REF] Pessoa | A generic exact solver for vehicle routing and related problems[END_REF] for more details about the BCP algorithm.

Since the methodology proposed in this work is a matheuristic one, optimality does not need to be preserved by the BCP. Thus, we turn the BCP algorithm into a heuristic (named BCP H ) by:

• Imposing a branch-and-bound node limit of 10 and time limit of 3,600 seconds;

• Using the false gap mechanism, described next;

• Using a restricted master heuristic, described below.

As mentioned above, the BCP algorithm uses an elimination procedure that removes edges from graph G by exploiting reduced cost arguments. In particular, if the minimum reduced cost of a path passing by an edge e ∈ E is not smaller than the gap between the current upper bound and the lower bound obtained by the column generation procedure, then edge e can safely be removed from the graph G, as no improving solution contains this edge. Removing edges makes subsequent calls to the labeling algorithm used for solving the pricing problem faster.

In addition to the edge elimination, path enumeration is also dependent on the gap. This procedure tries to enumerate all possible paths with reduced cost smaller than the current gap between upper and lower bounds. If path enumeration is successful (i.e., the number of enumerated paths is less than, say, one million; enough to store them in a table), the pricing problem from now on is solved by inspection. The inspection of enumerated paths is usually much faster to perform than to call the labeling algorithm. If the number of enumerated routes is sufficiently small (less than 10,000), the current node in the search tree can be finished by adding all enumerated routes to the restricted master and solving it as an IP (using a general solver like CPLEX).

The previous two paragraphs show the importance of having very good upper bounds (and, therefore, smaller gaps) for reducing the running time of the BCP algorithm. In fact, that is why POP solves smaller subproblems first (easy even with not so good upper bounds), so the solution of larger subproblems can benefit from already improved upper bounds. To further reduce the running time, the false gap mechanism artificially decreases the gaps when performing edge elimination and route enumeration. The false gap is defined as F G = (U B -LB)/F GF , where the false gap factor F GF > 1 is a parameter. Application of the false gap mechanism can result in removing edges or paths which participate in an improving solution. However, experiments indicate that such an outcome occurs rarely when one uses a moderate value for F GF (we tested F GF = 3).

Another difference from the default BCP algorithm by [START_REF] Pessoa | A generic exact solver for vehicle routing and related problems[END_REF] consists in using an additional heuristic (similar to the one proposed in [START_REF] Pessoa | A robust branch-cut-and-price algorithm for the heterogeneous fleet vehicle routing problem[END_REF]). It is called after the convergence of column and cut generation at every node of the search tree. The idea is to further decrease the false gap (dividing it by two in each iteration) until it is possible to complete the path enumeration. Then, the 10,000 routes with smaller reduced costs are used to create an IP that is solved by CPLEX.

Computational experiments

The proposed algorithm POP was coded in Julia language version 1.4.2. The algorithm BCP H to solve subproblems was obtained by parameterizing the CVRP demo application of VRPSolver [START_REF] Bulhões | VRPSolver: a branch-cut-and-price based exact solver for vehicle routing and some related problems[END_REF]. The parameters are described in Section 4.3. VRPSolver, freely available for academic use, implements the generic BCP algorithm proposed by [START_REF] Pessoa | A generic exact solver for vehicle routing and related problems[END_REF]. It makes use of the BaPCod C++ library [START_REF] Vanderbeck | BaPCod -a generic branch-and-price code[END_REF] as a BCP framework combined with the C++ implementations by [START_REF] Sadykov | A bucket graph based labeling algorithm with application to vehicle routing[END_REF] for solving pricing problems, route enumeration, and separation of rank-1 cuts. It also uses CVRPSEP package [START_REF] Lysgaard | Cvrpsep: A package of separation routines for the capacitated vehicle routing problem[END_REF] for separating rounded capacity cuts. Finally, VRPSolver uses CPLEX 12.9 to solve the LP relaxations and the MIPs over the enumerated paths.

All experiments with POP were performed on a 2 Deca-core Haswell Intel Xeon E5-2680 v3 server with 2.50 GHz and 128 GB of RAM, where each algorithm was executed on a single thread for each instance. Parallel runs for several different instances were performed on the same machine to speed up the experiments, effectively reducing the amount of memory allocated to each process.

Benchmark instances

The tests were performed on the 57 largest instances of the benchmark set X [START_REF] Uchoa | New benchmark instances for the capacitated vehicle routing problem[END_REF], ranging from 303 to 1001 vertices. Indeed, set X is currently the main benchmark used to assess the performance of all recent exact and heuristic algorithms for the CVRP. We skipped the 43 instances with less than 300 customers because most of those instances are now relatively easy for modern heuristics and even for modern exact algorithms. In fact, 39 of them have proved optimal solutions.

For a deeper analysis of some experiments, we split the 57 instances into two subsets: the subset X S of 29 instances with n/K min ≤ 10.8 (i.e., instances with short routes), and the subset X L composed by the other 28 instances (i.e., instances with long routes). The K min value is an instance attribute that means the minimum possible number of routes that a solution can have. For example, the instance X-n561-k42 belongs to X L because n/K min = 561/42 = 13.6 > 10.8. Extensive experiments presented in Pecin et al. (2017a) indicate that modern branch-cut-andprice algorithms for CVRP, like the one we use to solve subproblems in POP, perform considerably better on instances with shorter routes. Therefore, route size is a factor that is likely to affect the overall performance of POP.

Moreover, in the preliminary experiments used for calibration, we consider a small representative subset X R having only seven instances. The choice of X R is described in A.

The gap of a solution S is calculated as 100 • ((cost(S) -BKS)/BKS), where BKS is the best known solution in the CVRPLIB1 , only disregarding the solutions found by executions based on the proposed POP approach. Several optimization groups compete for improving the best known solutions for the instances in CVRPLIB. In fact, there were 24 updates in 2020 by seven distinct groups. Updating a BKS in CVRPLIB does not require the publication of an article; one only has to send the improved solution to be checked, even if the improvement is by only one unit. It is not necessary to describe how the solution was obtained. The competing groups may perform long runs of their methods, try several random number seeds, and even resort to special calibration. Thus, those BKSs are likely to be very close to optimum values.

Obtaining an initial solution

The initialization of POP is a critical issue. Preliminary experiments showed that it did not work so well as a stand-alone algorithm. It means that if it is initialized with a low-quality solution S obtained by a simple constructive heuristic, the overall performance of POP is not competitive with the best existing heuristics. In fact, we are proposing POP essentially as an effective way of improving solutions that are already reasonably good, possibly obtained by running some metaheuristic.

We report results obtained by different variants of POP 1 , which uses the HGS metaheuristic by [START_REF] Vidal | A hybrid genetic algorithm for multidepot and periodic vehicle routing problems[END_REF] to obtain the initial solution. However, as shown in B, the HGS is more effective if the entire algorithm is restarted (with a different random number seed) after 50, 000 iterations without any improvement (a method hereafter called HGS r ). Notation POP 1 t defines the variant that starts POP with the solution obtained by HGS r in t hours. We tested 4 values for t: 0.01 (36 seconds), 0.125 (450 seconds), 0.5 (1800 seconds), and 2 (7,200 seconds). Of course, the initialization time is included in the overall time. For example, in variant POP 1 0.5 , which is run for 32 hours, HGS r obtains the initial solution in 0.5 hours, and then POP spends 31.5 hours improving the initial solution.

The results obtained by several variants POP 1 over the time horizon of 32 hours are compared with those by HGS r itself. We also perform some comparisons with a second metaheuristic, the ILS-SP proposed by [START_REF] Subramanian | A hybrid algorithm for a class of vehicle routing problems[END_REF]. As shown in [START_REF] Uchoa | New benchmark instances for the capacitated vehicle routing problem[END_REF], although the HGS is on average substantially better than the ILS-SP, there are some instances (usually those with very short routes) where the ILS-SP is superior.

Parameterization of the subproblem solver

The default parameterization of the VRPSolver CVRP demo is calibrated to find optimal solutions for hard instances having around 200-300 customers. As POP needs to solve many smaller problems, we propose an alternative parameterization that works better inside POP. C presents the default and the proposed parameterizations, whose performances are compared in Figure 2 by running POP 1 0.5 on the X R instances over 8 hours. The convergence curves of the algorithms show the average gap found at different times.

The figure shows the superior performance obtained when the heuristic version BCP H is used to solve the subproblems. BCP H is obtained from the exact BCP using the proposed parameterization and by setting three additional parameters: RCSPfalseGapFactor to 3 (this activates the false gap mechanism described in Section 3), MaxNbOfBBtreeNodeTreated to 10 (maximum number of nodes in the branch-and-bound tree), and GlobalTimeLimit to 3600 seconds (maximum time for solving a subproblem) in the proposed parameterization. All POP results hereafter are obtained using algorithm BCP H as the subproblem solver.

Calibrating parameters α and δ

Table 1 shows the performance of POP 1 0.5 for different values of parameters α and δ. Each setting was applied to the X R instances over the horizon of 8 hours. The setting (α = 50, δ = 40) achieved the best performance for 2, 4, and 8 hours. Therefore, all POP results below are obtained with parameterization (α = 50, δ = 40). The performance of POP 1 0.01 deserves a separate analysis. It illustrates the behavior of POP as an "almost stand-alone" heuristic, starting from a medium quality solution. Such solutions can be rapidly obtained by any modern metaheuristic for the CVRP. The initial solutions provided by running HGS r for 36 seconds have an average gap of about 2% from the BKS.

• The performance of POP 1 0.01 on instances with shorter routes (set X S ) is very good. After 900 seconds, it already provides solutions that are significantly better than those from HGS r . After 4 hours, it is also consistently better than POP 1 0.125 , POP 1 0.5 , and POP 1 2 and reaches the excellent average gap of 0.017% in 32 hours. It is quite interesting to note that the final gap after 32 hours obtained by POP 1 t on instances X S gets worse as t increases. It seems that worse initial solutions used in POP 1 0.01 are still flexible enough to be transformed into good final solutions by the POP local search mechanism. On the other hand, the much better initial solutions used by POP 1 0.5 , and POP 1 2 seem to be biased towards certain local minima that may not be so globally good.

• On the other hand, POP 1 0.01 performs poorly on instances with long routes (set X L ). It takes 4 hours to obtain an average gap of 0.492%, and it reaches the performance of HGS r only after 16 hours. It is also consistently worse than POP 1 t , for t ∈ {0.125, 0.5, 2}.

The variants POP 1 0.125 , POP 1 0.5 , and POP 1 2 have a more robust performance. When the complete instance set X is considered, all of them are consistently better than HGS r alone (i.e., after POP starts, their average gaps are smaller at all times). This is also true when X S and X L instances are considered separately. The only exception is the variant POP 1 0.125 , which requires four hours to overcome HGS r on X L instances.

Table 3 reports the best solutions found by the algorithms ILS-SP, HGS r , and POP 1 in 32 hours. BKSs marked with a * are proven optimal solutions. Solutions marked in bold are improvements over the BKSs. The variant POP 1 2 achieved the best average and median final gaps, with the exception of the average gap for instances X S , where it is worse than the variants POP 1 0.01 , POP 1 0.125 , and POP 1 0.5 .

4.6

Comparison of the algorithms HGS20 and POP2 over 32 hours

When the work described in this article was already advanced, we were told 2 about the existence of a new implementation of HGS. The new version, specialized to CVRP, is faster and includes one additional neighborhood called SWAP*. We will refer to that yet unpublished algorithm as HGS20. In fact, the performance of HGS20 is much superior to HGS r , and thus it can definitely be considered as a state-of-the-art metaheuristic for CVRP. In this section, we test if POP can still improve HGS20 solutions.

On September 17, 2020, Thibaut Vidal kindly sent us the detailed results of ten 20-hour runs of the algorithm HGS20 on each of the X instances. Those runs are performed on Intel Xeon Gold 6148 @2.40GHz processors (PassMark single thread rating 2056) that are roughly equivalent to our processors (PassMark single thread rating 1840). Moreover, we have also received the solutions obtained after 0.125, 0.5, and 2 hours. Thus, we use them as initial solutions in the variant POP 2 t . Figure 4 depicts the performance of the HGS20 and the variants POP 2 over 32 hours (note that HGS stops at 20 hours). We now analyze these results.

• For each running time, HGS20 obtains solutions with about half of the average gap of the solutions obtained by HGS r , which is a remarkable improvement.

• Considering all X instances, the variant POP 2 0.125 is consistently worse than HGS20 alone, producing inferior solutions for all times.

• Considering all X instances, the variant POP 2 0.5 is slightly better than the HGS20 alone. On X L instances, it only starts to be better at 16 hours.

• Finally, the variant POP 2 2 is clearly better than the HSG20 alone, even on X L instances. This indicates that the proposed approach POP is indeed powerful. It is able to improve solutions obtained by a highly performing metaheuristic, at least in long runs (more than 2 hours).

Table 5 reports instance-by-instance statistics on the solutions found by the HGS20 (after 20 hours) and the variants POP 2 (after 20 and 32 hours). For the HGS20, the average cost and the best cost among ten runs are provided. For the variant POP 2 2 , we give the average cost and the best cost for three runs. We did not have computational resources for running each instance ten times. In order to provide a direct comparison between methods, we also computed the average cost and the best cost of the first three runs of the HGS20.

It is obvious that one can always obtain better solutions by performing multiple runs of any randomized method and picking the best one. But it is still interesting to note that the variant POP 2 2 seems to be particularly well suited for performing multiple runs. In fact, for this variant, the average gap for a single 32-hour run is 0.042%, while the average gap for the best of only three runs decreases to 0.018%, a very substantial decrease. For the instances in the set X S , the approach produces a remarkable gap of -0.001%.

Concluding remarks

In this work, we propose POP, a POPMUSIC matheuristic for the classical and highly competitive CVRP. The algorithm is designed to improve a reasonably good initial solution given as an input.

The results show that our approach outperforms one of the best published metaheuristics for the CVRP in medium and long runs. POP matheuristic is also competitive in long runs with a state-ofthe-art unpublished heuristic, which is specialized to the CVRP. The results are especially good for instances with relatively short routes. Moreover, several best known solutions were improved for literature instances with up to 20,000 customers. This shows a very good scalability of the approach.

POP matheuristic exploits a characteristic of the modern exact algorithms for vehicle routing problems, and in particular, those for the CVRP. If a tight upper bound on the optimum value is provided, those exact algorithms are usually capable of solving to optimality medium-size instances with up to 100-150 customers in a few minutes. Instances with less than 100 customers are usually solved in seconds. Thus, exact (or almost exact, the case of BCP H ) approaches become competitive with the best heuristics for solving such instances. An important advantage of exact approaches is that they "know" when to stop after proving that an improving solution does not exist, whereas traditional heuristics do not possess that information.

POP has interesting features that may be explored in future works:

• It is very different from other existing and well-performing heuristics for the CVRP. This means that their strengths and weaknesses may be complementary. This opens possibilities for many types of hybridization. It could be something simple, like just determining the instance characteristics (besides route size) that make it more or less suited to POP, in order to decide which method should be applied. But it could be something deeper, a full integration where POP and some traditional metaheuristic could take turns on improving parts of a solution and exchange information. This seems to be quite a promising direction for research.

• It is easy to implement, provided that an exact (or nearly exact) code for solving the subproblems is at hand. Given that VRPSolver branch-cut-and-price algorithm is available for academic use and can solve many routing variants other than CVRP, it is natural to try algorithms similar POP on those variants. Of course, there is no guarantee that a straightforward adaptation will obtain good results. Thus, there is room for research on extensions of POP that are more suited for other particular routing problems.

• It is naturally parallelizable. The current sequential version of POP may take a few hours to obtain high-quality solutions and can not be used in practical situations that require faster solutions. That limitation could be much reduced by solving several subproblems in parallel, a natural feature of any POPMUSIC approach. In fact, the larger the instance, the more parallelizable the method becomes.

The last remark is that the underlying implementation of the BCP solver used in POP was not changed by us other than by modifying external parameters. Thus, there is a large potential to improve the efficiency of POP by "going inside the black box". A property that may be exploited is that the subproblems to be solved are often very similar to some already solved subproblem. Keeping information from previous runs, like the generated columns and cuts, may accelerate the algorithm.

C VRPSolver Parameterizations

Table 7 shows the default VRPSolver CVRP parameterization, as well as the changed parameterization in the version used to solve subproblems in POP. The reader is referred to the documentation of [START_REF] Bulhões | VRPSolver: a branch-cut-and-price based exact solver for vehicle routing and some related problems[END_REF] to obtain the description of each parameter. The parameters which have special notation indicated in the second column of Table 7 are also described in [START_REF] Pessoa | A generic exact solver for vehicle routing and related problems[END_REF]. Default values for the last four parameters are not defined because these parameters are not active when the restricted master heuristic is not used.

RCSPmaxNumOfEnumSolsForEndOfNodeMIP is a previously undocumented VRPSolver parameter.

If the number of enumerated routes gets becomes smaller than RCSPmaxNumOfEnumSolutionsForMIP, at any point of a node solution, then the node is immediately solved by MIP. Table 8 shows the additional parameters used for obtaining BCP H . While they reduce running times, the optimal solution of a subproblem may be missed.

Figure 2 :

 2 Figure 2: POP 1 0.5 with three different parameterizations of VRPSolver. The time axis is on a log 2 scale.

Figure 3 :

 3 Figure 3: Convergence curves of POP 1 and HGS r .

Figure 4 :

 4 Figure 4: Convergence curves of POP 2 and HGS20.

  time limit is not exceeded and dimsp ≤ |V+| do 7 L ← a vector of vertices in V+ in random order

	11	Vsp ← ∅
	12	R ← ∅
	13	while |Vsp| < dimsp do
	14	r ← argmin	min j∈C(r) cij
		r∈S,r / ∈R
	15	if |Vsp| + |C(r)| ≤ dimsp then
	16	Vsp ← Vsp ∪ C(r)
	17	R ← R ∪ {r}
	18	else
	19	Go to the line 21
	20	/* If the same or a larger subproblem has not yet been solved, solve Vsp */
	21	if Vsp ⊆ V for all (V , S ) ∈ Π then
	22	Let Ssp be the subsolution for Vsp in S
	23	Π ← Π ∪ (Vsp, Ssp)
	24	Solve Vsp with the algorithm A using cost(Ssp) as the initial upper bound
	25	Let S sp be the subsolution found by the algorithm A, if any
	26	if S sp is found and cost(S sp ) < cost(Ssp) then
	27	Replace (Vsp, Ssp) by (Vsp, S sp ) in Π
	28	Update S by replacing subsolution Ssp by S sp
	29	Go to the line 7
	30	/* Increase the current target dimension */
	31	dimsp ← dimsp + δ

8 for z = 1, 2, . . . , n do 9 i ← L[z] 10 /* Build the subproblem for seed i */

Table 1 :

 1 Avg. gap (%) of POP 1 0.5 on X R instances for different values of α and δ.

Table 2 :

 2 Comparison of the algorithms ILS-SP, HGS r , and POP 1 over 32 hours Figure 3 and Table2show the gap convergence curves for HGS r and POP 1 over the horizon of 32 hours. Average gap (%) of HGS r and POP 1 executions at different times.

		.399	0.380	0.448	0.387	0.390	0.383	0.475	0.460	0.440
	2	0.327	0.303	0.294	0.310	0.308	0.263	0.344	0.328	0.326
	4	0.272	0.245	0.244	0.274	0.242	0.209	0.241	0.236	0.248
	8	0.198	0.170	0.165	0.213	0.193	0.162	0.183	0.214	0.198
	4.5 Instances Time (h)	HGS r	POP 1 0.01	POP 1 0.125	POP 1 0.5	POP 1 2
			0.01	1.996		1.996	-	-	-
			0.125	0.836		1.180	0.836	-	-
			0.25	0.626		0.931	0.564	-	-
			0.5		0.484		0.708	0.457	0.484	-
		All	1 2		0.396 0.330		0.579 0.420	0.355 0.271	0.317 0.240	-0.330
			4		0.283		0.293	0.206	0.182	0.171
			8		0.236		0.201	0.164	0.138	0.126
			16		0.210		0.137	0.117	0.099	0.090
			32		0.184		0.091	0.084	0.076	0.064
			0.01	1.790		1.790	-	-	-
			0.125	0.704		0.754	0.704	-	-
			0.25	0.547		0.524	0.398	-	-
			0.5		0.425		0.350	0.303	0.425	-
		X S	1 2		0.345 0.298		0.268 0.196	0.228 0.163	0.232 0.178	-0.298
			4		0.262		0.101	0.114	0.131	0.108
			8		0.196		0.058	0.069	0.084	0.081
			16		0.176		0.039	0.048	0.056	0.055
			32		0.162		0.017	0.022	0.036	0.037
			0.01	2.209		2.209	-	-	-
			0.125	0.973		1.621	0.973	-	-
			0.25	0.708		1.354	0.735	-	-
			0.5		0.546		1.080	0.617	0.546	-
		X L	1 2		0.448 0.363		0.902 0.652	0.487 0.383	0.404 0.305	-0.363
			4		0.306		0.492	0.301	0.235	0.237
			8		0.278		0.349	0.263	0.193	0.173
			16		0.244		0.238	0.190	0.144	0.127
			32		0.206		0.167	0.149	0.117	0.092

Table 3 :

 3 Best solutions found by ILS-SP, HGS r , and POP 1 after 32 hours.

	Instance	BKS	ILS-SP	HGS r	POP 1 0.01	POP 1 0.125	POP 1 0.5	POP 1 2
	X-n303-k21	21736	21840	21739	21863	21837	21750	21751
	X-n308-k13	25859	25881	25861	25876	25876	25876	25862
	X-n313-k71	94044	94105	94046	94053	94053	94046	94046
	X-n317-k53	78355 *	78355	78355	78355	78355	78355	78355
	X-n322-k28	29834 *	29872	29848	29887	29887	29887	29880
	X-n327-k20	27532	27743	27555	27573	27576	27573	27576
	X-n331-k15	31102 *	31108	31103	31103	31103	31103	31103
	X-n336-k84	139135	139253 139210	139164 139111	139175 139125
	X-n344-k43	42056	42096	42069	42055	42050	42050	42056
	X-n351-k40	25919	26131	25935	25896	25896	25896	25896
	X-n359-k29	51505	51997	51521	51583	51583	51583	51505
	X-n367-k17	22814	22912	22814	22814	22821	22821	22814
	X-n376-k94	147713 *	147713 147713	147713	147713	147713	147713
	X-n384-k52	65941	66382	66048	65941	65999	65956	65947
	X-n393-k38	38260 *	38273	38260	38260	38260	38260	38260
	X-n401-k29	66163	66614	66222	66181	66220	66257	66156
	X-n411-k19	19718	19811	19717	19712	19718	19712	19712
	X-n420-k130	107798 *	107798 107813	107798	107798	107798	107798
	X-n429-k61	65449	65759	65489	65527	65467	65467	65455
	X-n439-k37	36391 *	36402	36395	36395	36395	36395	36395
	X-n449-k29	55233	56131	55336	55332	55236	55259	55258
	X-n459-k26	24139	24421	24184	24193	24208	24209	24160
	X-n469-k138	221824 *	221940 222203	221824	221824	221824	221824
	X-n480-k70	89458	89821	89542	89449	89449	89449	89449
	X-n491-k59	66510	67128	66633	66555	66539	66572	66514
	X-n502-k39	69230	69315	69254	69232	69232	69232	69232
	X-n513-k21	24201	24275	24201	24248	24249	24201	24201
	X-n524-k153	154593 *	154698 154774	154593	154593	154593	154593
	X-n536-k96	94921	95697	95059	94948	94915	95205	95205
	X-n548-k50	86700 *	86710	86737	86701	86701	86701	86701
	X-n561-k42	42717	43016	42744	42758	42773	42758	42758
	X-n573-k30	50673	51074	50782	50807	50882	50742	50735
	X-n586-k159	190423	190767 190581 190365 190340 190375 190379
	X-n599-k92	108489	109147 108781	108498	108558	108517 108462
	X-n613-k62	59535	60318	59671	59561	59544	59606	59656
	X-n627-k43	62164	62762	62369	62182	62213	62245	62266
	X-n641-k35	63694	64449	64019	63773	63989	63919	63863
	X-n655-k131	106780 *	106780 106810	106780	106780	106780	106780
	X-n670-k130	146332	147286 147144	146346	146340	146411	146461
	X-n685-k75	68205	68682	68436	68260	68315	68318	68354
	X-n701-k44	81934	82907	82310	82085	81984	81970	82021
	X-n716-k35	43412	44091	43572	43443	43491	43498	43489
	X-n733-k159	136250	136900 136365 136245 136237	136278 136223
	X-n749-k98	77365	78177	77706	77380	77399	77360	77342
	X-n766-k71	114454	115413 114701	114573	114707	114640	114682
	X-n783-k48	72445	73627	72809	72696	72592	72605	72704
	X-n801-k40	73305	73939	73548	73446	73445	73368	73362
	X-n819-k171	158247	159249 158696 158128 158191 158222 158211
	X-n837-k142	193810	194901 194264	193820 193793 193800	193822
	X-n856-k95	88965	89143	89062	89030	89030	89030	89030
	X-n876-k59	99299	100357	99748	99583	99437	99479	99428
	X-n895-k37	53860	54777	54266	54112	54080	54125	54045
	X-n916-k207	329247	330773 329902	329305 329213	329305	329289
	X-n936-k151	132725	134564 133440	132859	132882	132863	132942
	X-n957-k87	85465	85887	85633	85485	85468	85492	85473
	X-n979-k58	118987	120015 119339	119430	119059	119073	119040
	X-n1001-k43	72359	73810	72766	72714	72506	72460	72486
	Avg. gap (%)		0.629	0.184	0.091	0.084	0.076	0.064
	Median gap (%)		0.607	0.117	0.027	0.031	0.032	0.009
	Avg. gap (%) in X S		0.468	0.162	0.017	0.022	0.036	0.037
	Median gap (%) in X S		0.474	0.117	0.005	0.000	0.002	0.000
	Avg. gap (%) in X L		0.796	0.206	0.167	0.149	0.117	0.092
	Median gap (%) in X L		0.779	0.138	0.150	0.135	0.091	0.073

Table 4 :

 4 Average gap (%) of HGS20 and POP 2 executions at different times.

	Instances Time (h) HGS20	POP 2 0.125	POP 2 0.5	POP 2 2
		0.125	0.407	0.407	-	-
		0.25	0.298	0.320	-	-
		0.5	0.222	0.273	0.222	-
		1	0.177	0.231	0.173	-
	All	2 4	0.136 0.115	0.191 0.154	0.145 0.127	0.136 0.095
		8	0.095	0.128	0.096	0.075
		16	0.083	0.094	0.071	0.056
		20	0.077	0.084	0.064	0.052
		32	-	0.072	0.050	0.039
		0.125	0.352	0.352	-	-
		0.25	0.258	0.242	-	-
		0.5	0.206	0.201	0.206	-
		1	0.166	0.167	0.145	-
	X S	2 4	0.127 0.110	0.139 0.105	0.124 0.104	0.127 0.070
		8	0.079	0.088	0.077	0.048
		16	0.074	0.063	0.055	0.035
		20	0.066	0.047	0.051	0.029
		32	-	0.038	0.033	0.018
		0.125	0.465	0.465	-	-
		0.25	0.339	0.400	-	-
		0.5	0.239	0.346	0.239	-
		1	0.188	0.298	0.201	-
	X L	2 4	0.145 0.120	0.245 0.204	0.167 0.149	0.145 0.121
		8	0.111	0.169	0.116	0.104
		16	0.093	0.126	0.087	0.079
		20	0.090	0.122	0.078	0.077
		32	-	0.107	0.068	0.060

Table 5 :

 5 Detailed statistics for HGS20 and POP 2 2 . Best gaps for 20 hours are underlined.

Table 5 -

 5 continued from previous page

	Instance	BKS
		0.079	0.043	0.080	0.060	0.052	0.027	0.042	0.018
	Median gap (%)	0.047	0.008	0.045	0.014	0.016	0.000	0.011	0.000
	Avg. gap (%) in X S	0.068	0.031	0.068	0.048	0.029	0.007	0.021	-0.001
	Median gap (%) in X S	0.042	0.008	0.040	0.010	0.009	0.000	0.002	0.000
	Avg. gap (%) in X L	0.092	0.055	0.093	0.073	0.076	0.048	0.064	0.037
	Median gap (%) in X L	0.049	0.007	0.051	0.026	0.038	0.018	0.034	0.010

Table 7 :

 7 Default and used parameters of the VRPSolver CVRP application used in

	Parameter	Notation Default value Used value
	RCSPhardTimeThresholdInPricing	τ hard	25 secs	8 secs
	RCSPstopCutGenTimeThresholdInPricing	τ soft	10 secs	3 secs
	RCSPnumberOfBucketsPerVertex	τ hard	25	50
	RCSPmaxNumOfLabelsInEnumeration	ω labels	5 • 10 6	3 • 10 5
	RCSPmaxNumOfEnumeratedSolutions	ω routes	5 • 10 6	10 6
	RCSPmaxNumOfEnumSolutionsForMIP	ω MIP	10 4	5 • 10 3
	RCSPmaxNumOfEnumSolsForEndOfNodeMIP	-	10 4	10 4
	RCSPuseBidirectionalSearch	φ bidir	2	1
	RCSPrankOneCutsMemoryType	θ mem	0	0
	CutTailingOffThreshold	δ gap	0.015	0.03
	StrongBranchingPhaseOneCandidatesNumber StrongBranchingPhaseOneTreeSizeEstimRatio StrongBranchingPhaseTwoCandidatesNumber StrongBranchingPhaseTwoTreeSizeEstimRatio MaxTimeForRestrictedMasterIpHeur	ζ num 1 ζ estim 1 ζ num 2 ζ estim 2 χ rm	100 0.2 5 0.02 -1 (off)	50 0.2 3 0.02 40
	CallFrequencyOfRestrictedMasterIpHeur	-	-	1
	MIPemphasisInRestrictedMasterIpHeur	-	-	1
	RCSPmaxNumOfLabelsInHeurEnumeration	-	-	10 5
	MaxNumEnumSolsInRestrictedMasterIpHeur	-	-	10 4

Table 8 :

 8 Additional parameters for obtaining BCP H

	Parameter	Notation Default value Used value
	GlobalTimeLimit	-	∞	3600 secs
	MaxNbOfBBtreeNodeTreated	-	∞	10
	RCSPfalseGapFactor	-	0 (off)	3

BKSs available in the CVRP Library (http://vrp.atd-lab.inf.puc-rio.br/) on October 31,

Personal communication from Thibaut Vidal.

Coefficient of variation

Acknowledgements

This study was financed in part by CAPES -Finance Code 001, by Capes PrInt UFF no 88881, by CNPq grant 313601/2018-6, and by FAPERJ grant E-26/202.887/2017. This research was also supported by Programa Institucional de Internacionalização (PrInt) from CAPES as part of the project REMATCH (process number 88887.310261/2018-00).

The experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr/).

We would like to thank Anand Subramanian for kindly sending us his code.

We are grateful to Thibaut Vidal for kindly providing us the results of his new heuristic and for insightful discussions.

A A representative small subset of X A minimum subset of the instances X which covers all the characteristics considered in [START_REF] Uchoa | New benchmark instances for the capacitated vehicle routing problem[END_REF] is described below:

• Route size (interval for n/K min ):

• Depot positioning:

• Customers demands:

-Unitary: X-n801-k40 -Small values, large CV 3 : X-n561-k42 -Small values, small CV: X-n393-k38 -Large values, large CV: X-n716-k35 -Large values, small CV: X-n469-k138 -Depending on quadrant: X-n979-k58 -Many small values, few large values: X-n670-k130 Thus, the subset X R is composed by: X-n393-k38, X-n469-k138, X-n561-k42, X-n670-k130, X-n716-k35, X-n801-k40, X-n979-k58. The reader is referred to [START_REF] Uchoa | New benchmark instances for the capacitated vehicle routing problem[END_REF] for a detailed description of all characteristics.

B Comparison of HGS and HGS r

Figure 5 compares the performance of a single execution of HGS and HGS r for all X instances over 8 hours. The convergence curves of both algorithms report the average gap (considering the gap obtained for each instance) found at different times. The final average gaps obtained by HGS and HGS r were 0.281% and 0.236%, respectively.


A Ruiz-Valencia

D Benmeziane

N Pen

E Petit

V Bonniol

M P Belleville

D Paolucci

José Sanchez-Marcano

L Soussan

J Sanchez- Marcano

CO2 valorization by a new microbiological process

Keywords: reduction, formate, bio-catalytic process, microbial process

Introduction

Carbon dioxide (CO2) accounts for more than 75% of global greenhouse gas (GHG) emissions, with annual anthropogenic emissions of about 25 to 35 Gt [START_REF] Giec | Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF]. These emissions, which today represent around 65% of GHG emissions in CO2 equivalents [START_REF] Giec | Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF], derive mainly from heavy and energy producing industries (cement plants, aluminum and steel production sites, coal and oil-fired power plants), as well as transport [START_REF] Giec | Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF][START_REF] Moret | [END_REF]. Use of CO2 as a raw material at the industrial level is currently being studied as a possible solution to reduce levels of this gas in the atmosphere. However, today, industrially, CO2 is mainly used as a fluid rather than a reactant and it originates mainly from natural sources or fermentation plants, its industrial use representing less than 1% (i.e. 100 Mt/yr) of annual emissions [START_REF] Moret | [END_REF]. Indeed, industrial use of the gas as a chemical reactant is currently very limited and based on the affinity of some nucleophilic molecules for CO2 to form urea, carbamates and cyclic organic carbonates [3].

In this context, over the last fifty years direct conversion of anthropogenic CO2 into alternative fuels (such as methane, C1-C5 alcohols or formic acid) has become a major research topic because it could help solve two major future challenges, namely climate change and fossil fuel shortages. Nevertheless, activation of CO2 molecules is very difficult because dissociation of the C=O bond requires a large amount of energy (i.e. 795 kJ.mol -1 ). Activation of CO2 therefore requires efficient and selective catalysts but also, quite often, relatively high temperatures and pressures.

Both homogeneous and heterogeneous catalysts have been reported to catalyze hydrogenation reactions [START_REF] Moret | [END_REF]4,5], artificial photosynthesis [6][7][8][9][10] and electrochemical reduction [3,11,12], giving rise to fuels and molecule platforms from CO2. However, despite recent progress, these chemical conversion routes still exhibit low durability and/or selectivity [6,[13][14][15].

The biological processes for CO2 activation are interesting because they are usually selective and require mild conditions. These biological processes are considered sustainable as they are low-energy-consuming and reagent use and/or by-product release that could have a negative environmental impact is generally limited.

Currently, only processes implementing photosynthetic microorganisms such as microalgae or cyanobacteria have been developed at the industrial scale. The main products recovered from CO2 valorization by these processes are lipids which can be then converted into biodiesel [16,17]. Open raceway technology has allowed the scaling-up of microalgae processes but their development is still limited due to (i) nutrients and mass transfer limitations, (ii) the need for good sunlight capture with special bioreactors requiring very large areas and low depths, (iii) the need to maintain a basic pH (8)(9)(10), and (iv) the need to harvest the microorganisms in order to extract the compounds of interest. All these constraints limit productivity and generate relatively high costs [18]. At the laboratory scale, cyanobacteria have been genetically modified to produce alcohols directly from CO2 [18][19][20][21]. Very recently, the carbon assimilation pathway of the cyanobacterium Synechococcus elongatus PCC 7942 has been directed towards the production of 2,3-butanediol from CO2 and glucose in the dark [22].

Such an advance overcomes the limitations related to the need for sunlight, but the need to add glucose still makes the bioprocess unattractive from an economic point of view. In addition, the biological stability of the modified cyanobacteria has yet to be improved for industrialization purposes [18][19][20][21][22]. electrolyzer tests were carried out with the bacterium to try and generate formate in situ at the cathode by reduction of CO2 [34], but electrochemical reduction of CO2 still remains difficult.

Finally, methanotrophic bacteria (Methylosinus trichosporium IMV 3011 and Methylosinus sporium) have been reported to selectively convert CO2 into methanol with their enzymes (including FDH first that can reduce CO2 into formate) and their intracellular PHB stock being used as a source of electrons and protons [35][36][37]. Contrary to the other biocatalysts, with these bacteria, there is no need to add costly chemicals such as cofactors or H2 to the reaction medium, or to provide the reaction system with light, or to genetically modify the metabolic pathway of the bacteria to recover products directly from CO2. Moreover, methanotrophic bacteria use methane (CH4) as a carbon and energy source to grow, meaning the global carbon balance of the CO2 valorization bioprocess is positive as our city and agricultural wastes are a renewable source of methane. Methanotrophic bacteria are thus interesting biocatalysts for CO2 reduction. However, the methanol production obtained was very low (3.6 µmol methanol•g dry cell -1 ).

In this work, an environmental bacterial consortium rich in methanotrophic bacteria was shown to directly and selectively catalyze the reduction of CO2 into formate, a chemical which could be considered as a future energetic vector since its acidic form, i.e. formic acid, can be used as fuel in Direct Formic Acid Fuel Cells (DFAFCs) or as an H2-storage molecule. Formate could also be employed as an electron donor for oxidation of CH4 into methanol by methanotrophic bacteria [38,39]. Significant amounts of formate were able to be produced (1.06 mmol formate•g dry cell -1 ) without the addition of costly chemicals to the reaction medium or in situ generating compounds for CO2 reduction. In addition, the native consortium was grown on methane, a renewable source of carbon. This new CO2 valorization route appears very promising and a patent application was therefore submitted [START_REF] Soussan | BRED-SATT, Procédé de valorisation du CO2 par réduction biologique -CO2 valorization process based on biological reduction[END_REF], justifying the composition of the bacterial consortium remaining confidential until publication of said patent.

Materials and methods

Bacteria and culture conditions

A bacterial consortium [START_REF] Soussan | BRED-SATT, Procédé de valorisation du CO2 par réduction biologique -CO2 valorization process based on biological reduction[END_REF] containing methanotrophic bacteria, the composition of which remains confidential, was grown in a modified nitrate mineral salts (NMS) medium enriched with copper and iron. The bacteria were grown in sealed vials incubated at 30°C on a rotary shaker (Unimax 1010, Heidolph) operated at 160 rpm. The headspace of the vials was filled with a mixture of CH4 and air (1:1 v/v) that was used as the source of carbon and energy; the gas volume was 3-fold that of the liquid. The NMS medium was composed of 1.06 g/L KH2PO4, 4.34 g/L Na2HPO4·12H2O, 1.7 g/L NaNO3, 0.34 g/L K2SO4 and 0.074 g/L MgSO4·7H2O; the pH of NMS medium was adjusted to 7.0 ± 0.1 with NaOH, 0.1 M or HCl, 0.1 M and the medium was then autoclaved.

Mineral, copper and iron solutions were prepared independently and sterilized by filtration on a 0.2-μm acetate cellulose membrane (Sartorius) before being added to the NMS medium. The final concentrations of the mineral solutions were 0.57 mg/L ZnSO4·7H2O, 0.446 mg/L MnSO4·H2O, 0.124 mg/L H3BO3, 0.096 mg/L Na2MoO4·2H2O, 0.096 mg/L KI, 7.00 mg/L CaCl2·2H2O, that of the copper solution 0.798 mg/L CuSO4, and that of the iron solution 11.20 mg/L FeSO4·7H2O.

Preparation of the bacterial suspension

At the beginning of the stationary phase, the cultured cells were harvested by centrifugation at 4000 g for 20 min at 4°C and then washed once with phosphate buffer 20 mM at pH 7.0. A reaction medium containing 20 mM of phosphate buffer at pH 7.0 and 5 mM of MgCl2 was used to re-suspend the cells to reach a final optical density at 600 nm (OD600) of about 10. Part of the resulting suspension was stored at -20°C with glycerol (20 % w/w) for further biochemical analysis.

CO2-reduction tests in batch mode

The batch reactors were 60 mL-glass vials sealed with aluminium caps. The bacterial suspension was distributed into the batch reactors, with 6 mL being added per reactor. To assess the ability of the consortium to reduce CO2, several experiments were conducted with different headspace atmospheres (Table 1) that were first sterilized by filtration on 0.2-μm PTFE filters (Sartorius). For experiments A, B and C, a large set of vials was filled with the same bacterial suspension and incubated at 30°C under constant rotation (160 rpm). Each experiment was monitored for 20 days and at each sampling point, a vial was taken for analytical measurements.

Different measurements were performed: optical density at 600 nm (OD600), pH, salts analysis by ionic chromatography (IC) and quantification of the volatile organic compounds (VOC) contained in the liquid samples by headspace GC-MS.

Each experiment was conducted in triplicate. CO2-reduction tests were also carried out with 13 CO2 to confirm the origin of the products arising from CO2 reduction by carbon NMR analyses.

1 v/v)
Air was taken from the atmosphere while CO2 (99.995 %) and N2 (99. 995 %) were obtained from Air Liquide and 13 CO2 (99.9 %) was obtained from Sigma Aldrich.

Sample processing

At various reaction times, the bacterial suspensions in the reaction vials (i.e. 6 mL) were recovered and divided into different fractions. A first fraction was used to measure pH and OD600. A second fraction was centrifuged at 10 000 𝑔 at 4°C for 10 min; the supernatant was then recovered to be immediately analyzed (by GC-MS or ionic chromatography). In case of NMR characterization, a portion of the freshly sampled suspension (0.5 mL) was transferred into NMR tubes for direct analysis.

Analytical methods

Optical density (OD600) was measured at 600 nm with a UV-2400 spectrophotometer (Shimadzu, Japan) and pH was measured with a C831 pH-meter (Consort, Belgium).

Volatile organic compounds, including methanol, ethanol, formaldehyde and acetaldehyde, were quantified by gas chromatography (Clarus 580, Perkin Elmer), coupled to mass spectrometry (Clarus SQ-8-MS, Perkin Elmer). A headspace system (Turbomatrix 16S, Perkin Elmer) was used for sample injection. Prior to analysis, acetonitrile was used as an internal standard at a final concentration of 39.3 ± 0.2 mg•L -1 . The method included a 15-min sample heating step to 75°C followed by injection into the GC system. Separation was carried out using a capillary column (Rt-Q-Bond Plot, Restek) under isothermal conditions (150°C for 25 min). Calibration curves were obtained independently for each VOC.

Salt quantification of supernatants was carried out by ionic chromatography (IC). Sample portions of 25 µL

were injected into a Dionex ICS-1000 system (Thermo Scientific™) with an IonPAc AS19 (0.4 x 250 mm) capillary column and elution was carried out under the following conditions: 10 mM KOH (0 to 10 min), then 10-45 mM KOH (10 to 30 min).

For NMR analysis, 450 µL of the sample (complete suspension) was poured into a 5 mm NMR tube with 50 µL of D2O. The analyses were performed using a NMR BRUKER Avance III -500MHz spectrometer -CryoProbe Helium system for four hours. Different standards prepared in the reaction medium were analyzed independently by NMR to obtain the carbon spectra of these molecules.

Results and Discussion

As detailed in Table 1, different experimental conditions were used to check the ability of the consortium to experiments A and B (Table 1); for each experiment, three independent runs were performed in parallel and each point represents the mean values obtained. Regarding OD600, a nearly 25 % decrease in bacterial concentration was observed during the first few days which then remained constant (from day 5), with no differences being observed between headspace atmospheres A and B (Fig. 1). As the initial OD600 was high (7.6 ± 0.6), the biomass concentration still remained significant during the experiments despite the decrease. The observations regarding OD600 were similar for experiment C (data not shown). The fact that methanotrophic bacteria have been reported to be very resistant and able to form spores in stressing conditions [START_REF] Semrau | [END_REF] could also explain the preservation of OD600 over 20 days despite the fact that the bacteria were exposed to unusual conditions (i.e. absence of their substrate CH4 and high CO2 content).

In the reactors where CO2 was added (experiment A), pH decreased from 7.0 to 6.4 ± 0.2 within 5 days and then remained at that value which corresponds to the first acid dissociation constant (pKa) of CO2 in water. A similar change in pH was observed in the reactors of experiment D as they also contained high levels of CO2 (data not shown). In the vials containing only air (experiment B), pH decreased slightly to 6.8 ± 0.1 due to the presence of traces of CO2 (i.e. 400 ppmv). When the headspace contained N2 (experiment C), pH remained constant at its initial value, i.e. 7.0 ± 0.1 (data not shown).

The control tests carried out with the consortium and N2 (experiment C) showed no evidence of microbial byproduct release or formation by GC-MS or ionic chromatography throughout the duration of the experiment.

No products were detected either in the control tests performed without bacteria but with an equimolar mixture of CO2 and air (experiment D), which confirmed that spontaneous CO2 reduction did not occur.

In experiments A and B (consortium exposed to an equimolar mixture of CO2 and air or to air), no VOC were detected by GC-MS. Nevertheless, ionic chromatography analyses demonstrated the appearance of formate after 8 days in both of the experiments (Fig. 1). The consortium was able to produce up to 750 mg•L -1 of formate in the CO2:air mixture (1:1) at day 15 of the reaction, corresponding to a formate production of 280 ± 10 mg•gdrycell -1 (taking into account the bacterial mass initially introduced into the reactors).

When the consortium was exposed to air (experiment B), which contains a small fraction of CO2 (400 ppmv), formate production was lower than that of experiment A (Fig. 1). In fact, 330 mg•L -1 of formate were produced corresponding to a production of 145 ± 15 mg•gdry cell -1 . These results show that formate production is CO2-and bacteria-dependant and suggest that formate could be formed from CO2 present in the atmosphere.

Nevertheless, the amount of formate obtained (i.e. 330 mg•L -1 x 6. 10 -3 mL = 1.98 mg) in experiment B (only air) was much higher than the expected theoretical quantity taking into account the CO2 content of air (54 mL of air with 400 ppmv of CO2 that can be transformed to a maximal formate mass of 0.04 mg). A possible explanation for the appearance of formate in excess in the reaction medium could be mixed acid fermentation.

This bacterial phenomenon occurs under oxygen limitation and could give rise to products such as acetate and formate by cleavage of intracellular pyruvate by pyruvate formate lyase [42]. In this work, anoxic conditions could have occurred due to a rapid depletion of the oxygen initially present in the reactors since oxygen is the usual final electron acceptor for the consortium; the bacterial consortium is indeed first grown under aerobic conditions before being exposed to CO2 (section 2.1).

In order to confirm the actual ability of the consortium to reduce CO2 into formate, a new set of experiments was carried out using a mixture of 13 CO2:air (1:1). The reaction was then monitored by 13 C NMR. Nuclear magnetic resonance spectroscopy (NMR) is only sensitive to 13 C since it has a spin number of ½. 13 CO2labelling was consequently used to monitor any potential 13 C-containing product. In addition, the technique can be used to detect products irrespective of their localization (i.e. inside or outside the cells) as the NMR magnetic field can penetrate bacterial cells. The NMR spectrum obtained after 8 days of reaction is shown in Figure 2. Resonance peaks for 13 CO2 (δ = 124.6 ppm) and its corresponding ionized form H 13 CO3 -(δ=160.2 ppm) were visible at the beginning. After 8 days of reaction, a new peak appeared (δ = 170.9 ppm), corresponding to the standard fingerprint of 13 C-labelled sodium formate. Since no other labelled compound was detected in the reaction medium, the labelled formate is thus likely to result from direct 13 CO2 reduction. This result demonstrates that the consortium is thus able to reduce CO2 into formate. However, in this experiment, it was not possible to recover semi-quantitative information about the labelled formate concentration produced since the 13 CO2 and H 13 CO3 -signals were saturated. The only available information was that the formate concentration on day 8 was higher than the low NMR threshold for detection of labelled formate (i.e. 6 mg•L - 1 ). If it is considered that bacteria exposed to CO2:air (Fig. 1, run A) were able to produce 35 ± 2 mg•L -1 of formate on day 9, it can be assumed that at least 17 % (= 6 / 35 x 100 %) of the formate produced resulted from CO2 reduction.

Reduction of CO2 by the consortium may be explained by several possible enzymatic cycles inside the bacteria.

Among the bacteria of the consortium, methanotrophic bacteria are supposed to be of CO2 reduction catalysts.

These bacteria use CH4 as an energy and carbon source. As described in Figure 3.A, the first step involves methane oxidation to methanol by methanol monooxygenase (MMO); methanol is then successively converted into formaldehyde and then formate by methanol dehydrogenase (MDH) and formaldehyde dehydrogenase (FADH), respectively. Finally, formate dehydrogenase (FDH) catalyzes the oxidation of formate to CO2. The methane consumed is fixed in the carbon cycle for protein and biomass production, and the CO2 produced from CH4 oxidation is reported to be also mainly fixed in the carbon cycle through carboxylase enzymes incorporating CO2 into organic molecules already present inside the cells [43]. In addition, formate can also enter the carbon cycle for the production of bacterial material, and PolyHydroxyButyrate (PHB) in particular.

PolyHydroxyButyrate is an energy storage polymer accumulated during the culture phase under oxygen and nitrogen limitation [44,45]. A potential route for CO2 assimilation by the bacteria is described in Figure 3.B. The CO2 could indeed be reduced into formate by the direct action of FDH and in that case, the equation balance of the reaction would be:

CO 2 + NADH ⟶ HCOO -+ NAD + (Eq. 1)
where NADH is the reduced form of the Nicotinamide Adenine Dinucleotide (NAD) cofactor present in living cells. However, the stock of endogenous NAD is finite. For the CO2 reduction reaction to occur (Eq. 1), the oxidized form of NAD (i.e. NAD + ) has to also be continuously reduced into NADH (Fig. 3.B). Some intracellular enzymes such as β-hydroxybutyrate dehydrogenase are able to recycle NAD+ into NADH by simultaneously oxidizing an electron and proton donor source (Fig. 3.B).

Compared to the work of Xin (2007) [36], no methanol traces were evidenced in the samples even if the compound could have been potentially present in the reaction medium as a product of successive formate reductions catalyzed by FADH and then MDH (Fig. 3.B). One possible explanation for this phenomenon could be that formate production was so high that it led to inhibition of FADH.

Furthermore, formate was assayed in the supernatant of the samples, meaning that the formate was external to the cells. At this stage, it was assumed that the formate was produced inside the cells by the enzymatic machinery required for both CO2 and NAD+ reductions and then released out of the cells. The probability that these reactions were catalyzed by free enzymes resulting from cell lysis is low since the meeting probability between FDH, CO2, NAD, and the enzyme regenerating NAD is low. Regarding formate excretion from the cells, two phenomena could be assumed: (1) formate release could be due to a bacterial regulation process when intracellular formate concentrations become too high for cell viability as formate does indeed have antimicrobial properties [46] and it's intracellular accumulation could lead to cell death, and (2) population differentiation could occur among the bacteria meaning that some bacteria could play the role of formate factories for the survival of other bacteria consuming the formate. Formate is effectively reported in the literature as a precursor used by methanotrophic bacteria for cell maintenance reactions [47,48].

A 5-day period seems to correspond to the adaptation phase for bacteria and their enzymes to their new carbon substrate, i.e. CO2 instead of CH4, before formate production occurs. A decrease in formate concentration (by 27 ± 1%) was then observed after 15 days (Fig. 1). This could be due to formate consumption. Bacterial regulation or changes in bacterial population interactions or substrate limitation could be at the origin of this decrease. However, the theoretical amount of CO2 available in the initial CO2:air mixture of the headspace was calculated according to the gas volume inside the vial (i.e. 54 mL) and the ideal gas equation. The calculated CO2 available (i.e. 1.09 mmol) was nearly 10 times higher than the maximum theoretical amount of CO2 consumed for formate production (i.e. 0.10 mmol) assessed using Eq. ( 1). This means that there was no CO2 limitation and that the bacteria had an unlimited capacity for exogenous formate production; the result therefore indicates that the limiting factor was rather the source of electrons and protons (i.e. the internal PHB) required for NAD recycling and thus for CO2 reduction into formate (Fig. 3.B).

The initial PHB dry mass content of the bacterial suspension was found to be about 15 ± 5% (data not shown).

A theoretical calculation, considering that PHB was composed of 1 000 mononer units, showed that a PHB mass content of approximately 10% would indeed make it possible to reach the maximum concentration of formate produced (i.e. 750 mg•L -1 ). However, the exact number of units making up intracellular PHB is unknown and PHB could also be used in other pathways such as mixed acid fermentation. This hypothesis therefore needs further investigation. Much more experimental work is also required to clarify the explanations provided above and to specify the metabolic ways of this interesting CO2 activation bioprocess.

Conclusions

A bacterial consortium containing methanotrophs was able to reduce CO2 into formate with a positive carbon balance as bacterial biocatalysts were grown on methane, which is also a greenhouse gas. This bioprocess thus constitutes a new promising way to valorize CO2. The role of each bacterium composing the consortium will be further investigated to both fully understand the mechanism and optimize formate production from CO2. In addition, formate production was shown to be enhanced by another phenomenon that is likely to be mixed acid fermentation. At this stage, it is assumed that 17% of the total formate produced resulted from CO2 reduction; however, this content needs to be precisely quantified. The impact of operating conditions should also be studied in the future to control the contribution of CO2 reduction on mixed acid fermentation. Finally, another way to optimize the process will focus on continuous formate recovery from the reaction medium and bipolar electrodialysis will also be tested.

  reduce CO2 into formate. The conditions of experiment A correspond to the reaction conditions for formate production, i.e. bacterial consortium in presence of a CO2:air mixture (1:1 v/v). Experiments B, C and D, performed without adding CO2 (B, C) or without bacteria (D), were blanks or control tests. Experiments B and C aimed to assess whether there was any possible microbial by-product release when the bacteria were placed in presence of air (B) or N2 (C). Experiment D was performed without bacteria and in presence of CO2 to check that spontaneous CO2 reduction could not occur.

Figure 1

 1 Figure 1 presents the results obtained in the batch reactors where formate appeared in the reaction medium, i.e.

Figure 1 .

 1 Figure 1. Reaction profiles of formate production, pH and OD600/(OD600)i obtained for experiments A and B. A: bacterial consortium with an equimolar CO2:air mixture and B: bacterial consortium with only air; (OD600)i means initial OD600; (OD600)i = 7.6 ± 0.6

Figure 2 .

 2 Figure 2. NMR spectrum at (a) t = 0 days, where CO2 and HCO3 -are presents; (b) t = 8 days showing a new compound with a chemical shift corresponding to formate, which is enlarged (c).

Figure 3 .

 3 Figure 3. Methanotrophic bacteria: (A) pathway for CH4 assimilation and (B) suggested pathway for CO2 reduction into formate.

Table 1 .

 1 Experimental conditions for the CO2 reduction tests

	Experiment	Presence of bacteria	Headspace atmosphere
	A	Yes	CO2:Air (1:1 v/v)
	B	Yes	Only Air (100 %)
	C	Yes	Only N2 (100%)
	D	No	CO2:Air (1:

Acknowledgments

A. Ruiz-Valencia acknowledges the Conacyt from Mexico for his PhD scholarship. Mrs Virginie Latil is thanked for reading and correcting the manuscript. SATT is thanked for partial financial support.


Nicolas Rouger
email: nicolas.rouger@laplace.univ-tlse.fr

Yazan Barazi

Marc Cousineau

Fréderic Richardeau

F Boige

Ing J M Blaquiere

Dr Jc L T Le

Modular Multilevel SOI-CMOS Active Gate Driver Architecture for SiC MOSFETs

Keywords: active gate driver, multi-level, SiC gate drivers, SOI CMOS. I

High Voltage SiC power MOSFETs have specific driving challenges such as a reduced short circuit capability (compared to Si IGBTs) and a weak field oxide layer, a large gate to source driving voltage (compared to GaN FETs), a high electric field under negative gate bias in off-state and a high switching speed. The negative bias in off-state creates a high stress which reduces the reliability of the SiC MOSFET. The high positive gate bias can generate large drain saturation current in case of short circuit events. We propose a modular multilevel architecture, which takes benefits of SOI isolation and low voltage CMOS transistors to drive SiC MOSFETs and to improve their reliability using an active and dynamic multilevel-selection on the switching sequences and on/off states.

INTRODUCTION

Wide bandgap power transistors such as SiC MOSFETs show improvements on the efficiencypower density trade-off for next generation power converters. However, SiC MOSFETs have specific behaviors which impact performances, reliability and their dedicated gate drivers: the drain current is governed by majority carriers and channel characteristics [START_REF] Hosoi | Reliability-aware design of metal/high-k gate stack for high-performance SiC power MOSFET[END_REF][START_REF] Wirths | Vertical 1.2kV SiC Power MOSFETs with High-k/Metal Gate Stack[END_REF][START_REF] Fayyaz | Aging and failure mechanisms of SiC Power MOSFETs under repetitive short-circuit pulses of different duration[END_REF], the gate to source voltage swing is rather large to achieve a low on state resistance (high positive VGS, e.g. +18V or +20V) and to guarantee a sufficient immunity to higher dVDS/dt and Miller parasitic coupling (negative VGS, e.g. -4V). To improve the overall performances and reliability of SiC MOSFETs and to offer new gate driving degrees of freedom, we propose an original multilevel gate driver architecture (building-block approach) on a single CMOS IC. The principle is to use the highly efficient low voltage 5V transistors and capacitors for the gate driver, and to connect several stages in series, each stage being isolated via Deep Trench Isolation (DTI) and integrated on the same chip. As a consequence, the applied gate-to-source voltage of the SiC power MOSFET can be easily shaped, with new functionalities.

II. CLASSICAL TOPOLOGY

Typically, the output buffer stage of gate drivers for SiC MOSFETs is based on a split output current amplification push-pull (Fig. 1,2). However, this simple circuit suffers from important limitations: a) the SiC gate driving voltages are linked with the IC technology and gate driver supply levels, hence neither their values nor the switching speed can be dynamically and quickly tuned.

High voltage CMOS MOSFETs (e.g. >40V) are required to be compatible with SiC driving voltages. b) Under hard switching fault (HSF) type I short-circuit, the SiC MOSFET's drain current is only limited by the high saturation current (between 7 and 10 times the rated current), which is related to the large applied gate-to-source voltage during turn-on and the low VGSth compared to Si MOSFETs (Fig. 3). c) to increase the safety margin of parasitic Miller coupling during dVDS/dt, the negative VGS during off-state creates a high stress on the SiC gate stack and also reduces the performances in III rd Author version of the accepted paper. Published in the proceedings of IEEE ISPSD 2020https://doi.org/10.1109/ISPSD46842.2020.9170181

To offer an improved driving scheme for SiC MOSFETs, we propose a Modular Multilevel Active Gate Driver (MMAGD).

Other attempts for multilevel gate drivers have been presented in the literature, mainly with discrete solutions [START_REF] Zhao | Adaptive Multi-Level Active Gate Drivers for SiC Power Devices[END_REF][START_REF] Ning | Feasibility and Limitation of DC/DC Multilevel Converter Power ICs Using Standard CMOS Transistors[END_REF][START_REF] Dong | A gate drive circuit with mid-level voltage for GaN transistors in a 7-MHz isolated resonant converter[END_REF][START_REF] Dymond | Multi-level active gate driver for SiC MOSFETs[END_REF][START_REF] Dearien | Slew Rate Control of High-Voltage SiC MOSFETs using Gate Resistance vs[END_REF][START_REF] Zeng | Comparative Study on Multiple Degrees of Freedom of Gate Drivers for Transient Behavior Regulation of SiC MOSFET[END_REF][START_REF] Idir | Active gate voltage control of turn-on di/dt and turn-off dv/dt in insulated gate transistors[END_REF], while active gate drivers focused on improved switching loss -EMI tradeoff [START_REF] Bau | Sub-Nanosecond Delay CMOS Active Gate Driver for Closed-Loop dV/dt Control of GaN Transistors[END_REF][START_REF] Chen | A Smart IGBT Gate Driver IC with Temperature Compensated Collector Current Sensing[END_REF] or short circuit protection [START_REF] Sadik | Short-Circuit Protection Circuits for Silicon-Carbide Power Transistors[END_REF]. The core idea here is to cascade 5V CMOS push-pull circuits in series, each being supplied by the energy stored in 5V floating capacitors (Fig 4 .). This technique can also be applied with 1.8V or other voltages, which takes benefits of low voltage transistors and integrated capacitors with better figures of merit than high voltage CMOS transistors (e.g. >40V) for 2-level gate drivers. However, it is required to isolate each low voltage transistor, due to the non-zero common mode voltage (e.g. higher than 20V) and floating potentials as a function of the state of the other series-connected transistors. As a consequence, the XFAB XT018 SOI CMOS technology (180nm)

was chosen, thanks to its isolation capability with DTI, efficient 1.8V digital, 5V-40V analog CMOS, high density passive components, high temperature capability and its availability through MPW. Each floating push-pull, pre-amplification circuits and

B. Preliminary results

Fig. 7. presents the characterization of one 5V submodule, which was designed and previously validated [START_REF] Rouger | CMOS SOI gate driver with integrated optical supply and optical driving for fast power transistors[END_REF]. The output 5V

PMOS and NMOS of one submodule have a RDSon of 0.26Ω (25°C), +/-5A source/sink capability and a maximum slew rate of 5V/ns. Other specifications can be found in [START_REF] Rouger | CMOS SOI gate driver with integrated optical supply and optical driving for fast power transistors[END_REF].

C. Design and layout of the SOI-CMOS MMAGD

Each submodule is based on the previously proposed 5V-0.26Ω split output buffer. Hence the complete MMAGD exhibits a maximal total gate driver resistance of 1.3Ω for 5 submodules connected in series. In each submodule, the 5V split output buffer

Fig. 2 .

 2 Fig. 2. Simulated waveforms of Fig. 1, with a constant Ciss (1nF). Two key challenges are highlighted in red, with a typical +16V/-4V VGS driving.

Fig. 1 .Fig. 3 .

 13 Fig. 1. Classical CMOS push-pull buffer with split output, driving SiC MOSFET with Kelvin Source.

Fig. 4 .

 4 Fig. 4. Proposed multi-level modular SOI-CMOS gate driver (MMAGD), based on five submodules. All functions are integrated on one IC (insert).

Fig 5 .

 5 shows the simulation of the MMAGD, where each floating capacitor is previously charged at 4V and an improved SiC driving is proposed. The PWM0..4 signals are generated based on the input PWM signal and the desired SiC VGS voltage sequence. Several voltage levels can also be obtained by a different combination of PWM0..4 signals (e.g. VGS =0V, +4V..+12V). Benefits of the MMAGD are highlighted in Fig5:a reduced HSF short circuit current thanks to a reduced VGS only during turn-on (albeit higher than the worst case Miller plateau voltage), a high VGS for low SiC conduction losses, a dynamic negative VGS to improve immunity to dVDS/dt (high/low side crosstalk), and a zero VGS in off-state to reduce the SiC field oxide stress[START_REF] Hosoi | Reliability-aware design of metal/high-k gate stack for high-performance SiC power MOSFET[END_REF][START_REF] Wirths | Vertical 1.2kV SiC Power MOSFETs with High-k/Metal Gate Stack[END_REF]. Even if the VGS voltage is reduced to +12V during turn-on, the peak gate current is maintained high to guarantee high switching speeds and low SiC switching losses thanks to a lower gate resistance than the one in classical 2-level operation. One has to note that the source internal resistance distributed on each submodule naturally contributes to this apparent gate resistance. The floating capacitors can be recharged by two ways: each submodule has a bootstrap-like circuit which can recharge its floating capacitor as long as the floating ground GND_HV0..4 is connected to the IC GND, and/or by harvesting a part of the SiC gate charge during switching transitions.

Level 1 :Fig. 5 .

 15 Fig. 5. VGS simulated waveform with the multilevel architecture from Fig.4., highlighting the benefits of the proposed multilevel driving scheme on the SiC MOSFET.

Fig. 7 .

 7 Fig. 7. Left, optical microscope picture of a 5V split output buffer. Right, Experimental waveforms on a 1nF constant capacitor, at 1MHz, VGS = 0V...+4V.

  quadrant operation. d) the stored SiC gate Author version of the accepted paper. Published in the proceedings of IEEE ISPSD 2020https://doi.org/10.1109/ISPSD46842.2020.9170181 charge is dissipated through the RGoff gate resistor during turn-off, whereas recycling such a gate charge would offer the possibility to downsize the isolated gate power supply. As a consequence of these limitations, one would like to actively tune the gate to source driving voltage of SiC MOSFETs, to combine a high switching speed capability, low on-state resistance, high immunity to dVDS/dt, improved short circuit capability, and to modify the gate driver architecture.

III. PROPOSED ARCHITECTECTURE

A. Modular Multilevel Active Gate Driver


Yazan Barazi
email: barazi@laplace.univ-tlse.fr

Nicolas Rouger
email: rouger@laplace.univ-tlse.fr

Fréderic Richardeau
email: frederic.richardeau@laplace.univ-tlse.fr

Laplace

CMOS Gate Driver with fast short circuit protection for SiC MOSFETs

Keywords: SiC MOSFET, Active Gate Driver, CMOS, Short-Circuit, Detection and Protection I

This article presents an alternative solution to the short circuit challenges commonly faced by SiC MOSFETs power transistors. In response to this issue, a dedicated fast CMOS active gate driver AGD is designed to detect the short circuits and to protect SiC MOSFETs, using only low voltage analog functions (5V and 40V transistors). The external high voltage diode used in the desaturation monitoring technique is no longer required and the short circuit detection can be much faster based only on signals observed by the gate driver. The IC prototype is based on XFAB XT018 0.18um CMOS SOI technology. The novel detection circuit is based on a 2D diagnosis, without any time dependency, fully integrated, with a low monitoring voltage. Detection of SC events is between 35ns and 130ns depending on the input capacitance of the power transistor, the internal/external resistors, detection threshold levels and other parameters.

INTRODUCTION

Silicon Carbide MOSFETs have been used in power converters to improve the efficiency and power density. SiC MOSFETs compared to Silicon IGBTs have better figures of merit, where SiC MOSFETs can have similar conduction losses with much improved switching losses. As the SiC MOSFETs gain maturity and their cost reduces with time, SiC MOSFETs are expected to be adapted and used in more and more applications [START_REF] Mukherjee | Silicon Carbide: Materials, Processing and Applications in Electronic Devices[END_REF].

The SiC MOSFETs allows higher voltage operations with lower conduction loss, higher switching speed operations with lower switching losses and higher temperature operation capability. In order to ensure these operations, the gate driver should adapt to the specificities of SiC MOSFETs. Essentially, the driving peak current, small propagation delay, high frequency operation, and small power consumption [START_REF] To | Modeling and characterization of 0.35 μm CMOS coreless transformer for gate drivers[END_REF]. Another performance feature to be taken into account by the gate driver is the reliability of detection and protection against short circuit (SC) events. Moreover, this detection and protection should be fast since SiC MOSFET short-circuit withstand capability is within tSC=2µs [START_REF] Romano | Short-circuit failure mechanism of SiC power MOSFETs[END_REF], [START_REF] Romano | Influence of design parameters on the short-circuit ruggedness of SiC power MOSFETs[END_REF] instead of tSC=10µs for Si IGBTs. (tSC = Short-Circuit Time capability). In addition, the SiC Mosfet has a much lower endurance for repeated short-circuit events than a silicon device. This endurance capability is lower, the longer the short-circuit application time, typically for tFaulted_pulse > Tsc/2. This is therefore clearly an additional reason to detect and turn-off softly the device with the shortest delay while providing robust detection to switching noise (typically tdetec. delay << tSC/2). Therefore, the SiC MOSFETs gate driver compared to the Si IGBT ones have lesser time to detect the SC and protect the power device and to avoid ageing of the chip.

Recent trends for integrated gate drivers are to integrate a variety of complex functions to address the challenges and to further improve the device performance. One of the main addressed functions, is the SC detection [START_REF] Yin | A Reliable Gate Driver with Desaturation and Over-Voltage Protection Circuits for SiC MOSFET[END_REF]. In the literature, several detection Author version of the accepted paper. Published in the proceedings of IEEE ISPSD 2020https://doi.org/10.1109/ISPSD46842.2020.9170164 methods were discussed [START_REF] Mocevic | Comparison between desaturation sensing and Rogowski coil current sensing for shortcircuit protection of 1.2 kV, 300 A SiC MOSFET module[END_REF]. The vds desaturation monitoring method is classically used for SiC MOSFETs and Si IGBTs [START_REF] Bertelshofer | Design Rules to Adapt the Desaturation Detection for SiC MOSFET Modules[END_REF].

However, this method requires a long blanking time, a high voltage diode, cross talk between power circuit and gate driver, and therefore impossible to integrate within a low voltage CMOS technology. Moreover, this monitoring strategy is not well suited, with MOSFET device, under high current overload operations because of its ohmic on-state greatly dependent of the junction temperature.

In this article, an active gate driver (AGD) is designed to detect the SC and to protect SiC MOSFETs. In order to do so, two original detection approaches are integrated: the occurrence (or not) of the Miller plateau, and the slope differences at the gate voltage, between normal operation and SC events, which offer different solutions and better performances than the prior art. These SC detection techniques are based on low voltage signals available in the gate driver.

II. SHORT CIRCUIT BEHAVIOR AND TESTBENCH

A. Studied Short Circuit

First and foremost, it is mandatory to understand the behavior of the power semiconductor device under short circuit, before designing its associated active gate driver. There are two main types of SC, Hard switch fault (HSF) (type I), which is the investigated fault and Fault Under load (type II, not discussed in this article). HSF occurs when the switch turns on with a permanent full bus voltage across itself typically caused by faulted control signals. FUL occurs when the device is already in onstate and an external short-circuit occurs [START_REF] Ji | Short-Circuit Characterization and Protection of 10-kV SiC mosfet[END_REF]. Here, the fault event will be detected at the low side of a switching cell to further protect the low side power MOSFET during its turn ON.

Fig. 1 presents the qualitative switching waveforms of a SiC MOSFET under Normal turn-on (NTO) and Hard Switch Fault (HSF). Briefly, the absence of the Miller plateau, changes on the slopes at the gate-source voltage, changes in the dvds/dt and the over-current at the drain are observed under HSF-SC. As it will be presented hereinafter, this distinctive behavior is used to detect and further protect the SiC MOSFET in case of HSF.

III. THE INTEGRATED CIRCUIT: DESIGN AND LAYOUT

A. Gate driver

Although the focus of this work is to propose new SC detection techniques, it is important to offer simple solutions to study robustness of the detection techniques with variable parameters. Among the parameters which affect the power switching speed and turn-ON of SiC MOSFET, the main parameters are the input capacitance of the power device, the internal/external resistor [START_REF] Zhao | An Intelligent Model-Based Multi-Level Active Gate Driver for Power Semiconductor Devices[END_REF] or the gate current [START_REF] Yu | A smart gate driver IC for GaN power transistors[END_REF]. In this work, the AGD was adapted to provide variable gate currents, and selectable gate resistors, to further assess the robustness of the detection methods. The CMOS IC is designed using X-Fab XT0.18 SOI technology, which combines the benefit of SOI wafers with Deep Trench Isolation (DTI) and six metal layers 0.18 µm process. as well as high voltage transistors up to 200V and high temperature operation. These properties facilitate the design of this IC circuit, to source high current and drive the power transistor with the required gate voltage [START_REF] Vogler | Fully integrated high voltage high current gate driver for MOSFET-inverters[END_REF]. fig. 3 3).

 An active Miller Clamp AMC with RDSon<1 ohm, effectively protecting the device from high dv/dt changes and safely shut down the power transistor in the case of negative power supply is removed.

 Last and not least, two short circuit detection approaches and protection functions. These functions will be developed in the section below.

B. Detection and Protection Circuit

The two SC detection approaches that are included in the IC circuit are based on the derivation method. The first approach is the dip approach, where the presence or non-presence of the Miller plateau is detected. The second one is the slope approach, based on the differences between the slope before and after the Miller plateau under NTO; or the unchanged slope under SC.

Once the SC has been detected, the protection circuit is similar to the literature (Soft Shut Down -SSD), where a larger resistance is used to discharge smoothly the gate of the power device. The Rds(on) of the integrated SSD is around 100 ohms, assuring a soft shut down (small size transistor within the IC). Once the SC has been detected, the segmented buffer goes to a high impedance state HZ (HZ+ signal in fig 3 .), where the gate of the power device becomes floating (PMOS and NMOS in output buffer are in OFF state).

Fig.

  Fig. 1. Characteristics switching waveforms of SiC MOSFET under NTO and SC conditions.

Fig. 2 .

 2 Fig. 2. Experimental waveforms for C3M0120090J device under NTO and SC conditions. Top view of the used evaluation board.

Fig. 3 .

 3 Fig. 3. IC schematic including all the integrated functions

Fig. 4 .

 4 Fig. 4. Microscopic view of the IC.

Fig. 5 .

 5 Fig. 5. Detection schematic based on the dip approach.

Fig. 7 -

 7 Fig.7-8 shows the AGD under simulation. Note that the experimental dvgs/dt is presented in the figures to illustrate the behavior and used as input signal into Cadence. vgs and the flags are simulated waveforms. The protection time is 50ns and 100ns respective to C3M and C2M.

Fig. 7

 7 Fig.7shows a detection time of 20ns, and a delay time behavior for the C3M power device of 30ns. For fast detection and original methods reference[START_REF] Miyazaki | 20-ns Short-circuit detection scheme with high variation-tolerance based on analog delay multiplier circuit for advanced IGBTs[END_REF],[START_REF] Hofstetter | Applying the 2D-Short Circuit Detection Method to SiC MOSFETs including an advanced Soft Turn Off[END_REF] are proposed. Yet for the C2M power device under different condition (RgExt=10Ω, vgs=-5V/20V), shows a detection time longer than the previous figure 75ns plus the delay time behavior of 25ns, fig.8. The difference is

Fig. 6 .

 6 Fig. 6. dvgs/dt characteristics under NTO and SC HSF, dip detection approach. Experimental waveforms further implemented in Cadence.

Fig. 7 .

 7 Fig. 7. Simulated waveforms on Cadence of the device under soft protection using experimental C3M0120090 VGS & VDS waveforms .

Fig. 8 .

 8 Fig. 8. Simulated waveforms on Cadence of the device under soft protection using experimental C2M0080120 VGS & VDS waveforms.

  1. Characteristics switching waveforms of SiC MOSFET under NTO and SC conditions.Author version of the accepted paper. Published in the proceedings of IEEE ISPSD 2020 -

	i Load v Bus i ds , v ds i ds , v ds	SC NTO	i ds = i Load I sat-Channel
	v gsNTO, v gsSC v Drv+ v gsM v th i gNTO, i gSC v Drv-		t t
	i gPeak
	i gM
	0 0		t
	(C gs +C gdmin ) (C gs +C gdmin )	(C gs +C gdmax )

https://doi.org/10.1109/ISPSD46842.2020.9170164 B. Testbench and pre-study Real Gate and Drain waveforms are measured (fig.2.) with Wolfspeed's evaluation-board and are further implemented to design the CMOS analog functions. The CREE evaluation-board includes a half bridge converter, with two C3M 1.2kV SiC MOSFETs, and their respective discrete gate drivers.

  This segmented output buffer allows three source/sink current configurations Author version of the accepted paper. Published in the proceedings of IEEE ISPSD 2020https://doi.org/10.1109/ISPSD46842.2020.9170164 (e.g. +/-1.5 A, +/-3A, +/-4.5A, 25°C), with split PMOS outputs which complementary offers other selectable current values through discrete gate resistors and enable bits (labelled EnP1, EnP2 in fig.

	presents the CMOS IC, which includes several
	functions.
	 A segmented output buffer using 40V-range transistors proposed by the technology, with 10ns and 2ns the respective
	propagation and propagation mismatch delay.

ACKNOWLEDGMENT

The authors thank Europractice-IC MPW services and IMEC for their support in the tape-out and packaging, and B. Lantin at Laplace lab.

Author version of the accepted paper. Published in the proceedings of IEEE ISPSD 2020https://doi.org/10.1109/ISPSD46842.2020.9170164

The detection time presented in fig. 7-8 is the time for vgs to reach vgs-Ref and to detect if dvgs/dt is above or below the reference level.

C3M and C2M have different input capacitances therefore different slope speed. The delay time presented in the figures is the time to process the detection flag, to activate the high impedance state of the output buffer, plus an additional delay before activating the SSD.

CONCLUSION

A fast and fully integrated method was proposed to detect and protect Hard Switch Fault type short circuit events for SiC

MOSFETs. With our techniques, the detection times are reduced to less than 100ns, and external components such as high voltage diodes are no longer required. Experiments are ongoing to assess the robustness of our active gate driver and the integrated protection techniques. High safety and improved endurance short-mode capability could thus be achieved with these two proposed ultra-fast detection concepts.