Valentin Antuori Renault

Emmanuel Hebrard

Marie-José Huguet

Alain Nguyen

Valentin Antuori

Combining Monte Carlo Tree Search and Depth First Search Methods for a Car Manufacturing Workshop Scheduling Problem

Keywords: 2012 ACM Subject Classification Mathematics of computing → Combinatoric problems, Mathematics of computing → Combinatorial optimization, Computing methodologies → Planning and scheduling, Computing methodologies → Discrete space search Monte-Carlo Tree Search, Travelling Salesman Problem, Scheduling. Digital Object Identifier 10.4230/LIPIcs.CP.2021.54

Many state-of-the-art methods for combinatorial games rely on Monte Carlo Tree Search (MCTS) method, coupled with machine learning techniques, and these techniques have also recently been applied to combinatorial optimization. In this paper, we propose an efficient approach to a Travelling Salesman Problem with time windows and capacity constraints from the automotive industry. This approach combines the principles of MCTS to balance exploration and exploitation of the search space and a backtracking method to explore promising branches, and to collect relevant information on visited subtrees. This is done simply by replacing the Monte-Carlo rollouts by budget-limited runs of a DFS method. Moreover, the evaluation of the promise of a node in the Monte-Carlo search tree is key, and is a major difference with the case of games. For that purpose, we propose to evaluate a node using the marginal increase of a lower bound of the objective function, weighted with an exponential decay on the depth, in previous simulations. Finally, since the number of Monte-Carlo rollouts and hence the confidence on the evaluation is higher towards the root of the search tree, we propose to adjust the balance exploration/exploitation to the length of the branch. Our experiments show that this method clearly outperforms the best known approaches for this problem.

Introduction

The assembly floor of our car manufacturer partner contains several machines, each producing a certain type of components and as many machines consuming those components. The process of moving components across the workshop, from the point where they are produced to the point where they are consumed is a major bottleneck for the production rate of the plant. The resulting transportation problem can be seen as a repetitive single vehicle pickup and delivery problem with time windows and capacity constraint. The repetitive aspect comes from the fact that over a weekly schedule, the pickups and deliveries between the same pairs of machines is repeated at a given frequency, and for the same reason, both tasks are constrained in time. Finally, the capacity comes from the specific trolleys used by operators,

54:2

Combining MCTS and DFS for a Car Manufactoring Workshop Scheduling Problem which can be stacked in trains of a bounded length.

The method used in the industrial context is a large scale scheduling model solved using local search solver. A range of approaches relying on reinforcement learning (RL) were recently proposed in [START_REF] Antuori | Leveraging Reinforcement Learning, Constraint Programming and Local Search: A Case Study in Car Manufacturing[END_REF]. A simple stochastic branching policy (a linear model over some problem-specific parameters) is learned via RL, and used either to guide a constraint programming approach with rapid restarts, a constraint approach with limited discrepancy search, or a multistart local search method. All three methods vastly outperform the industrial method both on real and synthetic data sets.

In this paper, we introduce a new approach, combining Monte-Carlo Tree Seach (MCTS) with budget-limited Depth First Search (DFS). MCTS was initially designed for solving AI games [START_REF] Coulom | Efficient Selectivity and Backup Operators in Monte-Carlo Tree Search[END_REF], and, over the last few years, MCTS combined with reinforcement learning and deep learning has enabled a breakthrough in the resolution of many combinatorial games (such as Go, with AlphaGo and AlphaGo Zero [START_REF] Silver | Mastering the game of Go with deep neural networks and tree search[END_REF][START_REF] Silver | Mastering the game of go without human knowledge[END_REF]). Monte-Carlo Tree Search [START_REF] Coulom | Efficient Selectivity and Backup Operators in Monte-Carlo Tree Search[END_REF] offers a good generic strategy to tackle combinatorial problems. The expected outcome of a subtree is evaluated via Monte-Carlo simulation: starting from an open node of the search tree, a complete solution is built using a randomized heuristic policy. The outcome of the rollout is back-propagated to that node and all its ancestors down to the root by computing an average. Then, the next node to expand is selected by traversing the search tree from the root using multi-armed bandits algorithms (e.g., Upper Confidence bounds applied to Trees, UCT [START_REF] Kocsis | Bandit Based Monte-Carlo Planning[END_REF]) until reaching a node that has not been expanded yet. Without requiring built-in domain knowledge, Monte-Carlo rollouts provide good guidance, and the expansion phase gives guarantees on the compromise between exploration and exploitation. We show that in our problem, replacing the Monte-Carlo rollouts by randomized, limited-budget DFS is effective.

Several hybridizations with combinatorial optimization frameworks have been proposed. MCTS has been combined with constraint programming (CP) in [START_REF] Loth | Bandit-Based Search for Constraint Programming[END_REF], where the simulation phase stops at first fail, and the authors do not allow backtracking. Moreover, in order to allow restarts, instead of keeping an evaluation of every open node, this is done on pairs variable/value, in a way inspired by the RAVE (Rapid Action Value Estimation) heuristic used in Go [START_REF] Gelly | Combining MCTS and DFS for a Car Manufactoring Workshop Scheduling Problem 9 Jack Goffinet and Raghuram Ramanujan. Monte-Carlo Tree Search for the Maximum Satisfiability Problem[END_REF]. Finally, took advantage of the fact that Gecode [START_REF]Gecode: Generic constraint development environment[END_REF] uses copying instead of trailing, to open every search node visited during a rollout. In the field of Boolean satisfiability (SAT), MTCS has been has been combined with SAT solver [START_REF] Previti | Monte-Carlo Style UCT Search for Boolean Satisfiability[END_REF], however, in this case without including the defining characteristics (clause learning, VSIDS, ect.) of modern SAT solvers. In [9], the authors propose to hybridize MCTS with local search to solve the MAX-SAT problem. They use a fixed limited-budget stochastic local search in place of the rollouts. Finally, in [START_REF] Sabharwal | Guiding combinatorial optimization with UCT[END_REF], the UCT algorithm has been used in mixed integer linear programming (MILP), although replacing Monte-Carlo rollouts by a lower bound obtains with the Linear Programming (LP) relaxation.

We are not aware of MCTS approaches using DFS rollouts. However, this is closely related to the Hybrid Best First Search (HBFS) algorithm introduced in [START_REF] Allouche | Anytime hybrid best-first search with tree decomposition for weighted CSP[END_REF] where limited DFS is interleaved with BFS, although the choice of leaf to expand is not based on the same principles. Besides using DFS, we propose two adaptations of MTCS method designed to be effective on our problem, but directly applicable in any combinatorial problem.

First, since the goal of a Monte-Carlo rollout is to evaluate a single decision, and since each subsequent heuristic decision reduces the relative impact of that first decision, we argue that the definition of the overall reward should reflect this form of "diminishing returns". Therefore, we propose to define the outcome of a rollout as the sum of the marginal increments of the lower bound at each step, weighted by a coefficient in ]0, 1[ that exponentially decreases V. Antuori, E. Hebrard, M-J. Huguet, S. Essodaigui and A. Nguyen 54:3 with the depth. When the coefficient tends towards 1, the outcome tends towards the overall objective value of the solution, and when it tends towards 0, the short term growth of the lower bound weigh more and more. Observe that this scheme is generic, it only requires a lower bound of the objective function which is monotonically non decreasing at each decision.

Second, in the multi-armed bandit algorithm, the tradeoff between exploration and exploitation is controled by a constant factor c for the exploration term. As the tree becomes deeper, the number of iterations of the multi-armed bandit along a branch grows. Therefore, the probability that it will deviate from the best branch so far grows exponentially with the depth of the branch. To offset this, we apply an exponential decay to the parameter c towards the root, so that the likelihood of deviating at the root decreases rapidly when the depth of the tree grows.

The paper is organized as follows. First, in Section 2 we describe the problem of routing vehicle components in car manufacturing workshops and we give a detailed overview of the standard MCTS algorithm in Section 3. Then, we present the novel aspects of our approach in Section 4. Finally, we give the specific implementation details for the considered problem in a MCTS framework in Section 5, and we report the results of extensive experiments on both industrial and synthetic data in Section 6. These experiments show that our adaptations of the MCTS method significantly outperforms previous methods, including the local search approach currently used in the industry.

Problem Description

The industrial assembly line consists of a set of m components to be moved across a workshop, from the point where they are produced to where they are consumed. Each component is produced and consumed by two unique machines, and it is carried from one to the other using four dedicated trolleys. Initially, there are two trolleys standing at the production point and two trolleys at the consumption point. On each side, one of them is full and the other is empty. However, the empty trolley at the production point is being filled, and the full trolley at the consumption point is being emptied. The full trolley at the production point must be brought to the consumption point before the initially full trolley there has been emptied, and symmetrically, the empty trolley at the consumption point must be brought to the production point before the initially empty trolley there has been filled. A production cycle is the time c i taken to produce (resp. consume) component i, that is, to fill (resp. empty) a trolley. The two pickups and the two deliveries (of empty and full trolleys) described above must then be done whithin this time window. The end of a production cycle marks the start of the next, hence there are n i = H ci cycles over a time horizon H for the component i. The problem is illustrated on a small example in Figure 1. In this example, there are 3 components having their own production and consumption machines, denoted by P i and C i in Figure 1(a). The lines between the machines represent the routes in the assembly line.

The time cycles of each component and the time horizon (H) are given in Figure 1 Precedences: Pickups must precede their deliveries in the same period.

χ(pf k i ) < χ(df k i ) ∧ χ(pe k i ) < χ(de k i ) ∀i ∈ [1, m] ∀k ∈ [1, n i ] (1)
Train length: The operator may assemble trolleys into a train (trolleys can be extracted out of the train in any order), so a pickup need not be directly followed by its delivery.

However, the total length of the train of trolleys must never exceed a length T max .

Notice that there are only two possible orderings for the four operations of a production cycle. Indeed, since the first delivery (which can be either the full or the empty trolley since they happen in parallel) and the second pickup take place at the same location, doing the second pickup before the first delivery is dominated: the train will needlessly contain both a full and an empty trolley for the same component, and this delivery will need to be done eventually and can only incurs further time loss.

This industrial problem is a repetitive single vehicle pickup and delivery problem with time windows and capacity constraint. In this problem, the production-consumption cycles of each component entail a very particular structure: the four operations of each component must take place in the same time windows and all of these operations are repeated for every cycle. In addition, all operations are mandatory and there is no objective function for the industrial application, instead, feasibility is hard. As a result, the efficiency of the Large Neighborhood Search approaches proposed in [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF] for such routing problems, are severely hampered since they rely on the length of the tour as the objective to evaluate the moves and the insertion of relaxed requests is often very constrained by the specific precedence structure.

This problem was previously introduced in [START_REF] Antuori | Leveraging Reinforcement Learning, Constraint Programming and Local Search: A Case Study in Car Manufacturing[END_REF], and both exact and heuristic methods were proposed to solve it. These approaches rely on a fine tuned heuristic, and it was observed for some instances that greedy dives of the solvers were able to find a solution. The main motivation for a MCTS approach comes from this observation as the algorithm strongly rely on greedy dives, and is entirely guided by them. The Monte-Carlo Tree Search Method

In this section, we give some overview of the Monte-Carlo Tree Search method, and we introduce notations that will be used in the following.

MCTS is a tree search heuristic method based on multi-armed bandit principles to guide the tree expansion and to ensure a compromise between exploration and exploitation. This method was widely studied in the context of games but also for solving optimization problems [START_REF] Sabharwal | Guiding combinatorial optimization with UCT[END_REF][START_REF] Philip Runarsson | Pilot, Rollout and Monte Carlo Tree Search Methods for Job Shop Scheduling[END_REF][START_REF] Minh | A monte carlo tree search for traveling salesman problem with drone[END_REF][START_REF] Mandziuk | UCT-Based Approach to Capacitated Vehicle Routing Problem[END_REF][START_REF] Matsumoto | Monte-carlo tree search for a reentrant scheduling problem[END_REF][START_REF] Chaslot | Monte-Carlo Tree Search in Production Management Problems[END_REF]. For a detailed survey on the MCTS method, the reader may refer to [START_REF] Browne | A Survey of Monte Carlo Tree Search Methods[END_REF]. In a nutshell, the MCTS method develops a search tree where a node corresponds to a state of a given problem, with final states being solutions. Each node is associated with a set of feasible actions leading to child nodes in the tree. The aim is to find a path from the root node to a final state maximizing a reward. The MCTS method is based on four principles:

1. a reward can be computed at each final state;

2. a simulation process, also called rollout, is used to produce a path from a given node to a final state (for instance based on random sampling);

3. a backpropagaton method to update node information after each new rollouts;

4. a selection mechanism, usually based on multi-armed bandit [START_REF] Kocsis | Bandit based monte-carlo planning[END_REF], for guiding the tree expansion and insuring a compromise between exploitation (select the most promising node) and exploration (visit different parts of the tree).

Let A be a set of actions. A state σ ∈ A * is a sequence of actions, and |σ| denotes its length. We note σ|a the state reached when applying action a in state σ, A(σ) denote the set of possible actions in state σ, and p(σ) the parent state of σ. The MCTS method stores in memory the tree T it has already explored, and for every state σ, it stores the triplet:

⟨N (σ), P r(σ), V (σ)⟩, where N (σ)
is the number of time (σ) has been visited, P r(σ) is the prior probability or prior preferences to choose the state σ from its parent state p(σ), and V (σ) is the expected value of subtrees rooted at σ, and computed by averaging the outcomes of Monte-Carlo rollouts. Notice that P r(σ) was introduce in the MCTS in [START_REF] Silver | Mastering the game of Go with deep neural networks and tree search[END_REF] but was not in the original form of MCTS.

The algorithm iterates over the four following phases until some stopping criteria are met.

Selection

The selection phase begins at the root node of T , and finishes when we reach a node that has not yet been explored. At each node σ ∈ T , an action is selected according to the statistics stored in σ:

a * = arg max a∈A(σ) Ṽ (σ|a) + c * U (σ|a) (2)
where Ṽ (σ|a) is the exploitation term (based on the value of node V (σ|a)), U (σ|a) is the exploration term, and c is a parameter which represents the balance between the two terms.

This process continues from the state σ|a * until a non-visited node is reached , i.e. a leaf of the subtree T .

In adversarial games, the value V of a node is the expected outcome, e.g., 1 for a win and 0 or -1 for a loss. In the context of combinatorial optimisation, however, several definitions have been used. A first possibility is to simply store the expected objective value, although this technique entails that rollouts must be complete, even when they are suboptimal early on. In [START_REF] Minh | A monte carlo tree search for traveling salesman problem with drone[END_REF] and [START_REF] Mandziuk | UCT-Based Approach to Capacitated Vehicle Routing Problem[END_REF], the authors consider a solution whose objective value is within a factor α of the best known solution as a "win" (the effective value is in [0, 1] depending on the C P 2 0 2 1

54:6

Combining MCTS and DFS for a Car Manufactoring Workshop Scheduling Problem quality of solution) and all other outcomes as loss (0). The parameter α must therefore be carefully chosen, and the likelihood of a positive reward decreases when the best known solution improves. In [START_REF] Loth | Bandit-Based Search for Constraint Programming[END_REF], the MCTS is frequently restarted (and hence the MCTS tree lost), then the authors store the outcomes of the rollouts on variable/value pair instead. In this technique the rollouts are depth first search calls stopped on the first fail, and the expected relative failure depth is stored for each variable/value pair instantiated in the selection phase.

Finally, in [START_REF] Sabharwal | Guiding combinatorial optimization with UCT[END_REF], instead of a rollout, the lower bound of the LP relaxation is backpropagated instead.

Observe that it is important to normalize the value V stored on the node, to make the choice of the balance exploitation/exploration parameter c more robust. A state value σ|a ending on the action a can be normalized in [-1, 1] as follows:

Ṽ (σ|a) = 2 * V + -V (σ|a) V + -V --1 if N (σ|a) > 0 0 otherwise (3)
Where

V + = max{V (σ|a) | a ∈ A(σ), N (σ|a) > 0} and V -= min{V (σ|a) | a ∈ A(σ)
, N (σ|a) > 0} are, respectively, the maximum and minimum values of any explored sibling state.

Finally, the exploration term is [START_REF] Silver | Mastering the game of Go with deep neural networks and tree search[END_REF]:

U (σ) = P r(σ) N (p(σ)) N (σ) + 1 (4)
The rationale is to select the action a that maximizes Ṽ (σ|a) plus a bonus that decreases with each visit in order to promote exploration. The prior probability P r(σ) biases the initial exploration by the knowledge we have on the state. The square root term could be replaced by a logarithmic term which is often used in MCTS, without changing this rationale.

Expansion

Let σ be the node returned by the selection procedure, during the expansion phase, for all a ∈ A(σ), a child σ|a is added to σ and initialized its visit counter N (σ|a) with its expected objective value V (σ|a) set to 0. If no prior probability P r(σ|a) is available for this state, the uniform distribution 1/|A(σ)| can be used instead.

Simulation

In the simulation phase, the state obtained by the selection phase is extended to a final state τ via a Monte-Carlo rollout. In the context of combinatorial optimization the final state is a feasible solution. The rollout is typically done by random sampling of the possible actions A(σ) from state σ following a stochastic policy. For instance, one can use the probability distribution given by P r(σ|a) | a ∈ A(σ). Alternatively, this can be done by any randomized greedy heuristic tailored to the problem at hand [START_REF] Mandziuk | UCT-Based Approach to Capacitated Vehicle Routing Problem[END_REF][START_REF] Minh | A monte carlo tree search for traveling salesman problem with drone[END_REF][START_REF] Philip Runarsson | Pilot, Rollout and Monte Carlo Tree Search Methods for Job Shop Scheduling[END_REF]. As mention before, hybridization with existing technologies can take place in this phase, whether it is a linear relaxation [START_REF] Sabharwal | Guiding combinatorial optimization with UCT[END_REF], a local search [9] or a call to a CP solver [START_REF] Loth | Bandit-Based Search for Constraint Programming[END_REF].

V. Antuori, E. Hebrard, M-J. Huguet, S. Essodaigui and A. Nguyen

54:7 Backpropagation

Finally for each node σ traversed during the selection procedure, we update its statistics regarding the final state τ obtained by the simulation phase:

V (σ) ← V (σ) + z(τ, σ) -V (σ) N (σ) + 1 N (σ) ← N (σ) + 1
with z(τ, σ) the outcome of the rollout τ evaluated from node σ. The first update rule allows to maintain the average outcome of the rollouts for each node traversed during the selection step. It is possible to change the rule to only keep the best outcome find when traversing the node, instead of the average [START_REF] Sabharwal | Guiding combinatorial optimization with UCT[END_REF][START_REF] Philip Runarsson | Pilot, Rollout and Monte Carlo Tree Search Methods for Job Shop Scheduling[END_REF]. The rationale is the same as minimax algorithms for games, the optimistic view is that eventually search will find the best completion of a partial solution, and therefore its expected value is closer to the best rollout than to the average of all rollouts. However, the preferred choice may depend on the standard deviation of the outcomes of rollouts under a given node, and on the ratio of the whole search tree that the algorithm will eventually explore. For this reason, for larger problems, and when the heuristic used during the rollouts is robust, the average may be better.

Tailoring Monte-Carlo Tree Search to Combinatorial Optimization

In this section, we introduce three modifications of standard Monte Carlo Tree Search which we empirically found beneficial in the context of optimization problems. These modifications are generic, in the sense that they hold outside of our industrial application, as long as we have a lower bound computation technique for the objective function and a depth first search procedure for the target problem.

Evaluation based on the objective function

In game playing, the outcome of a Monte-Carlo rollout may only be known when the game ends. Typically, the rollout is given a value of 1 for a win, -1 for a loss and 0 for a draw.

Standard adaptations to combinatorial optimization are to normalize the objective value in a way or another as described in Section 3.

When simulating long branches, however, a "mistake" on a single decision along the branch can make the final outcome irrelevant. In fact, look-ahead methods often exhibit diminishing returns. For instance, it was observed in Chess that the rate of wins in self-plays between an algorithm looking k + 1 plies ahead versus the same algorithm looking k plies ahead declines as k grows [START_REF] Heinz | New Self-Play Results in Computer Chess[END_REF]. In the case of a greedy procedure, it is therefore natural to conjecture that as the length of the branch grows, the correlation between the quality of the initial decision and the overall outcome decreases.

In combinatorial optimization problems, however, we usually have a lower bound on the objective that monotonically grows with every decision. Therefore, the evolution of this value can provide a better insight into the quality of an initial decision. Let LB : A * → R be a lower bound on sequences of actions, with LB(σ) equals to the objective value if σ is a final state. Then, for a given node σ we propose to evaluate a state σ ′ reachable from σ as the sum of the marginal increment of the lower bound LB in the path from σ to σ ′ , weighted by an exponetially decaying coefficient γ. Hence we can define this sum recursively as follows:

z(σ ′ , σ) = LB(σ) -LB(p(σ)) if σ ′ = σ γ |σ ′ |-|σ| (LB(σ ′ ) -LB(p(σ ′ ))) + z(p(σ ′ ), σ) otherwise (5)
C P 2 0 2 1 54:8

Combining MCTS and DFS for a Car Manufactoring Workshop Scheduling Problem

The evaluation of a final state τ obtained by a rollout is then simply z(τ, σ) and represents an upper bound of the optimal solution.

Algorithm 1 implements backpropagation following the reward defined in Equation 5.

This algorithm takes as an input the sequence (R) of the lower bound increments given by the rollout, the node selected in the selection phase, and the decay rate.

Algorithm 1 Backpropagation procedure Data: R : sequence of the lower bound increments, σ : selected node, γ : decay rate 1 // Sum of exponentially decaying marginal increment of the lower bound

2 val ← |R| i=1 γ i-1 R i 3 // Backpropagation until the root node 4 repeat 5 val ← γ * val + LB(σ) -LB(p(σ)) 6 N (σ) ← N (σ) + 1 7 V (σ) ← V (σ) + val-V (σ) N (σ) 8 σ ← p(σ) 9 until σ = N il;
The proposed evaluation method puts more weight on the short-term impact of a decision, wagering on it being more reliable than long term observations. For γ = 1, the score reflects the objective value LB(τ ) of the rollout, whereas greater weight is put on short-term impacts when γ tends towards 0.

Moreover, the lower bound computations can be used during the expansion phase to avoid expending into sequences whose objective value cannot be lower than the current upper bound (best known solution). Thus, a node σ ′ that cannot be expanded further (all potential children nodes are suboptimal) is removed from the search tree. In that case, the information is backpropagated along the branch that leads to this node, that is, each node σ containing the deleted node σ ′ in its subtree are updated:

V (σ) ← 1 N (σ) -N (σ ′ ) (V (σ) * N (σ) -V (σ ′ ) * N (σ ′ ))
and

N (σ) ← N (σ) -N (σ ′ )
Indeed, all information contained in the deleted node is now irrelevant for the rest of the search as it is not in the tree anymore. Then previous iterations which have passed throughout this node should not have an impact on the future search.

Then, for the implementation of the proposed evaluation function, we should store LB(σ)

at each node σ in addition to the triplet {N (σ), P r(σ), V (σ)}.

A potential limit with this evaluation method is that it may skew search towards postponing actions that greatly increase the lower bound, but must eventually be done. For instance, consider a Travelling Salesman Problem with an isolated city far away from all other cities. Rollouts where this city is visited last will be preferred to rollouts where it is visited early. Lower bounds that take into account the future decisions in a reasonable way (e.g., minimum spanning tree for the travelling salesman problem, or the prehemptive relaxation in scheduling) may prevent this phenomenon since the cost of an exceptionally remote city or of an exceptionally large task would contribute to the lower bound anyways.

V. Antuori, E. Hebrard, M-J. Huguet, S. Essodaigui and A. Nguyen

54:9

The lower bound we used in our industrial problem, however, is extremelly basic and yet this did not seem to be an issue in our experiments.

Dynamic Exploitation vs Exploration Balance

Since the tree grows deeper as search progresses, the likelihood to deviate from the best branch increases. Therefore, we propose to dynamically adapt the parameter that control the balance between exploration and exploitation, depending on the depth of the tree, in order to promote exploitation on deeper nodes. Let td(T ) be the depth of the tree T , then at step t of the selection phase, the exploitation/exploration coefficient will be

β td(T )-t * c ( 6
)
with β < 1 a parameter. This mechanism has a similar effect as committing to a move at the root node. At the root t = 0 and thus β td(T )-t * c tends towards 0 when td(T ) grows, so the first decision is very unlikely to deviate from the most promising choice once the search tree has sufficiently grown. Conversely, at a leaf, this term tends towards the original value c and hence less promising -but less frequently visited -nodes will be selected more often. In the context of games, when a move is actually made, it makes sense to forget the siblings and parents of the corresponding state. In optimization, this mechanism has been implemented in several approaches in order to limit the combinatorial explosion [START_REF] Bertsimas | A comparison of Monte Carlo tree search and rolling horizon optimization for large-scale dynamic resource allocation problems[END_REF][START_REF] Mandziuk | UCT-Based Approach to Capacitated Vehicle Routing Problem[END_REF]. Since commits are irreversible, the algorithm is no longer complete, and budget parameters controling such commits need to be carefully chosen. Instead, the mechanism we propose has a similar effect but in a "smooth" way: near the root, it is more likely that the best move will be chosen, however other states can still be reached.

Depth First Search as a rollout

Finally we propose to use a Depth First Search procedure instead of a randomized greedy heuristic in the simulation phase. More precisely, in order to intensify the search around promising areas, a budget is defined after a first greedy "dive" and a budget-limited DFS is performed. For this purpose, the simulation is split into three steps:

The first step is a greedy randomized procedure from the selected node σ until a contradiction is detected. This contradiction can happen because a constraint is violated, or because the lower bound exceeds the upper bound. At this point, we define a budget for the DFS by evaluating the current state σ ′ . This budget will be larger if this is a promising state, and maximal if no contradiction was encountered (and hence a new upper bound was found). On the other hand, if the state σ ′ is not promising, then the budget will be smaller or null.

The second step of the simulation is a DFS, from the state reached by the greedy procedure σ ′ . This search is only performed on the subtree rooted at the node σ selected in the selection phase. The DFS algorithm must be able to store the best branch discovered, that is, the best solution or the best partial sequence according to the evaluation we described previously.

The third step begins when the budget is consumed (or the search is complete for the subtree rooted at the selected node σ in the selection phase). If no solution was found during the previous step, the greedy randomized procedure is used to extend the best branch found by the DFS to a complete solution which can be evaluated before backpropagation.

C P 2 0 2 1

54:10 Combining MCTS and DFS for a Car Manufactoring Workshop Scheduling Problem

The evaluation procedures for the states and for the budget will be detailed in Section 5

as their definitions depends on the considered problem.

Adaptation to the industrial Workshop Scheduling Problem

Tree model

In the search tree of the MCTS method, a state σ represents a partial sequence of operations and the set of actions correspond to the set of operations of the routing problem described in Section 2, i.e., actions are operations A = O. In the search tree, a sequence σ|a is the sequence σ extended by the action (operation) a. The set of possible actions A(σ) from a sequence σ contains every operation a such that the (partial) sequence σ|a is feasible with respect to the constraints.

Objective function

Since our industrial application is a satisfaction problem (the existence of a tour without delay), we need to generalize it to an optimization problem to apply MCTS as described in Sections 3 and 4. Therefore, during the simulation phase, we relax the due date constraints and instead we minimize the maximum tardiness:

L(σ) = max(0, max 1≤j≤|σ| (e σ(j) -d σ(j) ))
Since in this case operations can finish later than their due dates, it is necessary to make the precedence constraints due to production cycles explicit:

max(ρ(df k-1 i ), ρ(de k-1 i )) < min(ρ(pe k i ), ρ(pf k i )) ∀i ∈ [1, m] ∀k ∈ [2, n i ] (7)
Furthermore, during the expansion phase we do not add a child node that would violate a due date constraint, as our primary goal is to find a solution σ without any late job, that is, such that L(σ) = 0.

We use a trivial lower bound, which is at state σ the maximum tardiness L(σ) of the associated partial sequence also taking into account tardiness of all pending operations.

Pending operations are all the operations that belong to a production cycle in which at least one operation is available to extend the current sequence, ignoring the train constraint.

Therefore, Equation ( 5) is the sum of exponentially decaying marginal increments of the maximum tardiness with a small look ahead.

Heuristic

For the simulation phase as well as for the probabilities of the expansion phase, we use the heuristic tuned by reinforcement learning proposed in [START_REF] Antuori | Leveraging Reinforcement Learning, Constraint Programming and Local Search: A Case Study in Car Manufacturing[END_REF]. This heuristic is stochastic and provides a probability distribution over the set of available operations for a given state. More precisely, at a given state σ, each operation a ∈ A(σ) is evaluated using a fitness function f (σ, a) defined as a linear combination of four criteria: f (σ, a) = θ ⊺ λ(σ, a). These criteria λ i correspond to:

1. The emergency of the operation: lst(a, σ) -max(r a , e σ(|σ|) + tt σ(|σ|),a ), with lst(a, σ) the latest starting time of the operation a in order to satisfy the due date constrains with respect to the operations belonging to σ and the precedences constraints; 4. The type of operation, equal to 1 for pickups and 0 for deliveries.

The parameter θ is set to the proposed learned values (0.251, 0.576, 0.148, 0.023). Then, a softmax function is applied to turn the fitness evaluation into a probability distribution for guiding the choice of the next node in the greedy heuristic:

∀o ∈ A(σ) π θ (o | σ) = e (1-f (σ,o))/δ o ′ ∈A(σ) e (1-f (σ,o ′ ))/δ (8)
where the parameter δ controls the "greedyness" of the heuristic, that is, a "low" value for δ encourages to select the best choice with high probability, whereas a more "neutral" value of δ produces more randomized choices. In the experiments, we will set a value of δ = 0.005 in the simulation phase, and a value of δ = 0.1 to initialize the prior probabilities of the new nodes in the expansion phase.

Simulation

The greedy procedures before and after the DFS simply consist in taking at random the next operation following the probability distribution defined by equation 8.

For the DFS we define a backtrack budget between 0 and B, depending on when the first tardiness was detected during the first dive. If the first dive finds an improving solution, then the budget is maximum (B), in order to find other related improving solutions. Otherwise, we rely on the rank ϕ where the lower bound became positive to define the budget. Let ϕ * be the highest rank for any previous solution, the backtrack budget is then:

       B if ϕ ≥ ϕ * B( ϕ * -ϕ ϕ * -α * ϕ * ) 2 if ϕ * > ϕ > α * ϕ * 0 otherwise (9)
with α ≤ 1, a threshold parameter.

During the DFS, we define a probability distribution over the children using the softmax function of the greedy heuristic and we limit the breadth of the tree by keeping only actions with a probability greater than 10 -6 , which typically leaves all but 1 to 3 children approximately. Then, those children are sorted by their probabilities, and in order to randomize the DFS, a random child (again, using the same probability distribution) is swapped with the first one, to be branched on first by the DFS. As instances can be very large, this is sufficient to keep variety in the solutions, while removing many "bad" decisions.

This is also why we rely on the backtrack count instead of the fail count to define the budget, as a lot of nodes may have only one child. We add a geometric restart policy in the DFS step, where the search is reset to the node selected in the selection phase. The growth factor is reset at each MCTS iteration. At the end, the DFS returns the longest (potentially partial) sequence for which the lower bound remains null (i.e., for as the largest number of operations). Then, the greedy procedure is called to extend this sequence to a complete solution.

Experimental Evaluation

We report in this section the results of our experiments. First, we assess the respective impact of using the new evaluation policy, the dynamic exploration/exploitation balance and the DFS in the simulation phase. In a second part, we compare our MCTS adaptations to state-of-the-art methods for this problem. The number of components is highly correlated with hardness, and directly related to the branching factor in the Monte-Carlo search tree. Indeed, each node has at most two children per component (ie., from 30 children for instances of category A to 60 children for instances of category D). In addition, the depth of the search tree grows with the number of operations, that depends both on the time horizon and on the number of components. This depth varies from 450 for the "shift" schedules, up to 14500 for the "weekly" schedules.

We ran every method 10 times for each of the 120 instances with a timeout of 1h. All experiments were run on a cluster composed of Xeon E5-2695 v3 @ 2.30GHz processors.

Our methods were implemented using in C++ and compiled with GCC-8.0. The two methods from [START_REF] Antuori | Leveraging Reinforcement Learning, Constraint Programming and Local Search: A Case Study in Car Manufacturing[END_REF] were implementd using JAVA and were run in the same conditions, and Choco-4.10 [START_REF] Prud'homme | Choco Solver Documentation[END_REF] for CP.

Impact of the MCTS adaptations

In the first part of the experiments, the goal is to assess the respective impact of the proposed adaptations for the MCTS method. We evaluated 6 different versions of the MCTS, adding the adaptations we propose one at a time:

MCTS is the standard MCTS method without any of the proposed adaptation. This baseline method uses the value of the objective function as the result of the rollouts, and backpropagates this value through the tree to the root node.

MCTS+DFS is the same algorithm as MCTS except that it uses the DFS in the simulation phase.

SEDMI is the variant of MCTS that uses the sum of exponentially decaying marginal increments of the lower bound to evaluate the nodes.

SEDMI+DFS adds the DFS to SEDMI for the simulation phase.

SEDMI+DFS+DC extends SEDMI+DFS with the dynamic exploitation/exploration compromise.

SEDMI+SAT-DFS+DC is the variant of SEDMI+DFS+DC in which the upper bound on the objective function is fixed to 1 in the DFS, i.e. the DFS tries to solve the satisfaction version of the problem instead of trying to improve the global upper bound. However, the last part of the simulation still provides a complete solution via a greedy procedure, and hence this method also provides an upper bound.

All parameters for the proposed methods are given in Table 1. We recall that c is the exploitation/exploration tradeoff parameter. The higher value for this parameter, the more the MCTS will explore. Then, β is the decay rate for the adaptation of c, and γ is the decay rate of the evaluation function. Finally, α and B are respectively the threshold parameter, and the maximum backtrack budget for the DFS. All the values for these parameters were chosen by preliminary experiments, and the chosen combination appears to give relatively good overall results.

The results are shown in Table 2 and3, in which we report the number of solved runs, and the average maximum tardiness. For all the methods we consider that an instance is solved if and only if the value of the objective function is null i.e. there is no tardiness. Table 2 shows the performance of the different variants of the MCTS averaged by classes of instances, and by time horizons. In this table, for each method, a line corresponds to 100 runs (10 instances and 10 runs for every time horizon), then the number of solved runs is a sum over these 100 runs. In Table 3 the same results are presented aggregated by time horizons, and the number of solved instances is in percentage (over the 400 runs by line and by method).

In these tables, we can see the benefit of using the DFS in the simulation phase. Using DFS, as expected, allows the MCTS methods to solve more instances on the week horizon.

In fact, those instances are too large to be solved via rollouts only, and the DFS allows to intensify the search on the deepest parts of the tree, that are not explored in the MCTS.

Unfortunately, the effect of the DFS is not visible on the shift horizon. We can also see the benefit of using the sum of exponentially decaying marginal increments as node evaluation on day and week horizons in terms of objective value. However, this adaptation slightly degrades the performance on shorter horizon meaning that this time horizon is too short to take advantage of this mechanism. Overall, the combination of both mechanisms outperforms the two versions with only one of these mechanisms. Finally, the effect of the dynamic compromise can be seen on the day horizon. This time horizon is small, but not enough for the MCTS to advance deep enough in the search tree to find solutions. This mechanism forces the MCTS to explore the tree deeper and faster, and as a results, to improve the number of solved instances.

Comparison with previous methods

For the second part of the experiments, we compare the two best MCTS methods, namely SEDMI+DFS+DC (the method leading to the lowest objective function) and SEDMI+SAT-DFS+DC (the method with the highest number of solved instances) to the two best methods introduced in [START_REF] Antuori | Leveraging Reinforcement Learning, Constraint Programming and Local Search: A Case Study in Car Manufacturing[END_REF], that are both based on the stochastic branching policy described in Section 5:

CP: a constraint programming approach with rapid restarts. This method solves the satisfaction version of the problem. As a result, it is slightly better for finding solutions without tardiness.

GRASP: a multi-start local search procedure. This method considers the optimization problem with relaxed due dates as in the MCTS methods, hence we can compare the overall tardiness. The results, given in Table 4, show that overall, the proposed MCTS adaptations outperform the CP and the local search approaches on both criteria: the number of solved instances, and the maximum tardiness. More precisely, the dominance is clear for horizons shift and day in terms of number of instances solved, but we can see that our method does not outperform the CP model on the week horizon. Finally, between the CP approach and the SEDMI+SAT-DFS+DC variant, there is a difference of 9.5% of instances solved in favor of the latter. There is still half of the instances that are not solved to optimality. However, the instances of the data set were randomly generated without a guarantee of satisfiability, and, we believe that the majority of unsolved instances are not satisfiable (especially for the week horizon).

Conclusion

In this paper, we have presented and applied several variants of the Monte Carlo Tree Search method to solve a repetitive single vehicle pickup and delivery problem with time windows and capacity constraint, issuing from car manufacturing assembly lines. We defined a way of evaluating the rollouts based on the growth of the lower bound of the objective function. We also proposed an adaptation of the balance parameter between exploitation and exploration in order to be able to solve larger instances. Moreover, we proposed an hybridization of Monte Carlo Tree Search with Depth First Search used during the simulation phase. The experimental evaluation demonstrates that these proposals allow us to outperform previous approaches on the considered problem, and show the benefit of our contributions.

These three proposals, although well suited to a dedicated problem, are generic. The next step is then to demonstrate the genericity of these Monte Carlo Tree Search variants by considering their application to other combinatorial optimization problems. We also plan to integrate our MCTS method in existing constraint programming solvers to take advantage of their search tree exploration in the Depth Fisrt Search part, further reinforcing the hybrid nature of the approach. Finally, we would like to explore further the learning aspects of the method. Indeed, in the simulation phase, we are repeatedly dealing with similar subproblems in different part of the tree, and the policy used in a subtree could be adjusted after each iteration in order to have different policies adapted to different parts of the tree search.

Figure 1

 1 Figure 1 Illustrative example

2 . 3 .

 23 The travel/waiting time of the operation: max(tt σ(|σ|),a , (r a -e σ(|σ|) )); The (negated) length of the trolley;

C P 2 0 2 1 54: 12 6 . 1

 11261 Combining MCTS and DFS for a Car Manufactoring Workshop Scheduling Problem Experimental protocolWe use the same data set as in[START_REF] Antuori | Leveraging Reinforcement Learning, Constraint Programming and Local Search: A Case Study in Car Manufacturing[END_REF] composed of 120 synthetic instances. The data set is made of four categories characterized by the number of components (15 in category A, 20 in B, 25 in C and 30 in D). Moreover, all of these categories are associated to three time horizons: a work shift of an operator (7 hours and 15 minutes), a work day (made up of three shifts) and a full week (6 days).

  An operation o occurring at period k for component i is given a release date r o = (k -1)c i and a due date d o = kc i , with r o ≤ s o and e o ≤ d o .

	(b). In
	this example, there are 3 time cycles for the yellow component, 4 time cycles for the red
	component and 2 for the blue one.
	For each component i, for each of its cycles k, there are two pickups and two deliveries:
	the pickup pe k i and delivery de k i of the empty trolley from the consumption area to the
	production one, and the pickup pf k i and delivery df k i of the full trolley from production
	to consumption. The processing time of an operation o is denoted pt o and the travel time
	between operations o and o ′ is denoted tt o,o ′ .
	Let O be a set of all pickup and delivery operations with |O| = n. The problem is
	to compute a sequence ω : {1, . . . , n} → O of the operations O, where ω(j) is the j-th
	operation in the sequence, and χ = ω -1 its inverse. The sequence ω must satisfy the
	following constraints:
	Routing: For every 1 < j ≤ n, operation ω(j) must be given a start time s ω(j) (and
	end time e ω(j) = s ω(j) + pt ω(j) ) taking into account duration and travel time: s ω(j) ≥
	s ω(j-1) + pt ω(j-1) + tt ω(j-1),ω(j) (and s ω(1) = 0).
	Time windows:

Table 1

 1

	Parameters value
	c	1
	β	0.995
	γ	0.9977
	α	0.9
	B	50000
	Restart (base)	100
	Restart (factor)	1.2

Table 2

 2 Comparison of the MCTS adaptations

		H	MCTS	MCTS+DFS	SEDMI	SEDMI+DFS SEDMI+DFS+DC SEDMI+SAT-DFS+DC
		#S Lmax #S Lmax #S Lmax #S Lmax #S Lmax	#S	Lmax
		shift 100	0	100	0	100	0	100	0	100	0	100	0
	A	day 90		135	90	133	90	115	90	77 100	0	98	0
		week 68 1996 78 1850 70 1800 80 1839 77	1840	80	1858
		shift 80		420	79 258 90	372	82	353	81	349	87	439
	B	day 50 2954 54 2959 60	2522	60 2439 63	2121	70	2134
		week 10 21070 29 20771 10 20572 32 20541 31 20355 36	20635
		shift 49 1676 48 1708 40	1901	45 1727 40	1824	40	2012
	C	day 10 9503 11 9248 10	8683	26 8656 36	8747	35	9022
		week 0 64442 8 64713 0	64480	9 64584 9	64474 10	64445
		shift 40 2154 33 2146 30	2338	33 2018 30	2304	30	2621
	D	day 0 13659 0 13664 0	12657	0 12723 13 12225 11	12340
		week 0 101474 0 101444 0 100533 0 100760 0	100954	0	100840
	Average 41 18290 44 18241 42 17998 46 17976 48 17933 50	18029

Table 3

 3 Results aggregated by time horizon

	H	MCTS	MCTS+DFS	SEDMI	SEDMI+DFS SEDMI+DFS+DC SEDMI+SAT-DFS+DC
		#S Lmax #S Lmax #S Lmax #S Lmax #S	Lmax	#S	Lmax
	Shift 0.67 1063 0.65 1028 0.65 1153 0.65 1024 0.63	1119	0.64	1268
	Day 0.38 6562 0.39 6501 0.40 5994 0.44 5974 0.53 5773 0.54	5874
	Week 0.20 47245 0.29 47195 0.20 46846 0.30 46931 0.29 46906 0.32	46944
						C P 2 0 2 1

54:14 Combining MCTS and DFS for a Car Manufactoring Workshop Scheduling Problem

Table 4

 4 Comparison with previous methods

		H	CP		GRASP	SEDMI+DFS+DC SEDMI+SAT-DFS+DC
			#S	#S	Lmax	#S	Lmax	#S	Lmax
		shift	90	90	10	100	0	100	0
	A	day	90	90	193	100	0	98	0
		week	80	70	2433	77	1840	80	1858
		shift	60	60	467	81	349	87	439
	B	day	52	46	3218	63	2121	70	2134
		week	35	10	26915	31	20355	36	20635
		shift	40	40	1941	40	1824	40	2012
	C	day	10	10	9498	36	8747	35	9022
		week	10	0	71104	9	64474	10	64445
		shift	19	16	2677	30	2304	30	2621
	D	day	0	0	13994	13	12225	11	12340
		week	0	0	107186	0	100954	0	100840
	Average	40.5	36	19969	48	17933	50	18029


Mickael D Chekroun

Henk Dijkstra

Taylan S ¸eng Ül

Shouhong Wang

TRANSITIONS OF ZONAL FLOWS IN A TWO-LAYER QUASI-GEOSTROPHIC OCEAN MODEL

We consider a 2-layer quasi-geostrophic ocean model where the upper layer is forced by a steady Kolmogorov wind stress in a periodic channel domain, which allows to mathematically study the nonlinear development of the resulting flow. The model supports a steady parallel shear flow as a response to the wind stress. As the maximal velocity of the shear flow (equivalently the maximal amplitude of the wind forcing) exceeds a critical threshold, the zonal jet destabilizes due to baroclinic instability and we numerically demonstrate that a first transition occurs. We obtain reduced equations of the system using the formalism of dynamic transition theory and establish two scenarios which completely describe this first transition. The generic scenario is that two modes become critical and a Hopf bifurcation occurs as a result. Under an appropriate set of parameters describing midlatitude oceanic flows, we show that this first transition is continuous: a supercritical Hopf bifurcation occurs and a stable time periodic solution bifurcates. We also investigate the case of double Hopf bifurcations which occur when four modes of the linear stability problem simultaneously destabilize the zonal jet. In this case we prove that, in the relevant parameter regime, the flow exhibits a continuous transition accompanied by a bifurcated attractor homeomorphic to S 3 . The topological structure of this attractor is analyzed in detail and is shown to depend on the system parameters. In particular, this attractor contains (stable or unstable) time-periodic solutions and a quasi-periodic solution.

Introduction

Baroclinic instability is among the most important geophysical fluid dynamical instabilities playing a crucial role in the dynamics of atmospheres and oceans. In particular, this instability mechanism is the dominant process in atmospheric dynamics shaping the cyclones and anticyclones that dominate weather in mid-latitudes, as well as the mesoscale ocean eddies that play various roles in oceanic dynamics and the transport of heat of salt [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics[END_REF]. Much is known on the linear stability of zonal jets in a horizontally unbounded ocean in the quasi-geostrophic (QG) flow regime. Classical models, such as the continuously stratified Eady model [START_REF] Eady | Long waves and cyclone waves[END_REF] and the two-layer Phillips model [START_REF] Phillips | A simple three-dimensional model for the study of large-scale extratropical flow patterns[END_REF], have lead to a detailed understanding of the mechanism of baroclinic instability of a zonal jet in the inviscid case. Long waves destabilize the zonal jet with maximum growth rates occurring for perturbations having wavelengths on the order of the Rossby deformation radius, typically 50 -100 km for the mid-latitude ocean [START_REF] Talley | Descriptive Physical Oceanography[END_REF].

In case linear friction is included in the two-layer model, the neutral curve has a minimum at (k c , µ c ) where k c is the critical wavenumber and µ c the critical value of the control parameter (e.g. the maximum speed of the zonal jet). The nonlinear development of these perturbations has been extensively analyzed in the weakly nonlinear case [START_REF] Pedlosky | Finite-Amplitude Baroclinic Waves[END_REF][START_REF] Romea | The effects of friction and β on finite amplitude baroclinic waves[END_REF]29,[START_REF] Lovegrove | Bifurcations and instabilities in rotating, two-layer fluids: II. β-plane[END_REF]. In the regime |k -k c | = O( ) and |µ -µ c | = O( 2 ), [27] showed that on a long time scale T = 2 t and large spatial scale X = (x -c g t), where c g is the group velocity of the waves at criticality, the complex amplitude A of the wave packet destabilizing the jet satisfies a Ginzburg-Landau equation, written as

∂A ∂T = γ 1 A + γ 2 ∂ 2 A ∂X 2 -γ 3 A|A| 2
, where the γ i are complex constants. [27] also showed that the fixed point solution of this equation can become unstable to sideband instabilities. Subsequent analysis has shown [START_REF] Dijkstra | On the physics of upgradient momentum transport in unstable eastward jets[END_REF] that upgradient momentum transport can occur due to the self-interaction of the instabilities leading to rectification of the zonal jet.

In reality, the ocean basins are zonally bounded by continents and the midlatitude zonal jets are part of the gyre system, for example the subpolar gyre and subtropical gyre in the North Atlantic, forced by the surface wind stress through Ekman pumping [START_REF] Pedlosky | Geophysical Fluid Dynamics[END_REF]. The problem of baroclinic instability of such non-parallel flows is much more complicated and has so far only been tackled numerically. When the wind-forced QG equations are discretized, the linear stability problem for the gyre flow results in a large-dimensional generalized eigenvalue problem, typically of dimension 10 4 . There are many results for the one-layer single-and double-gyre flows (for an overview, see chapter 5 of [START_REF] Dijkstra | Nonlinear Physical Oceanography: A Dynamical Systems Approach to the Large Scale Ocean Circulation and El Niño, 2nd Revised and Enlarged edition[END_REF]), but in this case there is no baroclinic instability. There are relatively few results for the two-layer case. In [START_REF] Dijkstra | Temporal variability of the Wind-Driven Quasi-geostrophic Double Gyre Ocean Circulation: Basic Bifurcation Diagrams[END_REF], it was shown that in the two-layer case the double-gyre flow becomes unstable through a sequence of Hopf bifurcations. The perturbation flow patterns at criticality are 'bananashaped', locally resembling those of baroclinic instability in the Philips model. Stable periodic orbits result from these Hopf bifurcations, typically given rise to meandering motion of the gyre boundary.

As an intermediate, more analytically tractable case, we consider here the baroclinic instability of a zonal jet for a two-layer QG model in a zonally periodic channel. In this case, the properties of the bounded geometry are somehow represented, as the patterns of the unstable modes are restricted by the periodicity of the channel, so a sequence of Hopf bifurcations is expected just as in the more realistic gyre case. In addition, parallel flow solutions exist in the zonally periodic channel which simplifies the linear stability problem substantially such that a more detailed nonlinear analysis, akin to that in the horizontally unbounded case, can be performed. The parallel flow can also be connected to the surface wind stress, as in the full gyre case, but at the expense of adding an additional linear friction term to the upper layer vorticity equation; for more details, see Section 2 below.

The case specifically studied in the paper is the circulation set up by a time-independent Kolmogorov windstress field (for k = 1, 2, . . .)

τ x (y) = -τ τ 0 kπ cos kπ y L y ; τ y = 0
where τ 0 is a characteristic mid-latitude wind-stress value. This wind stress forces an ocean enclosed in a rectangular basin [0, 2L x ] × [-L y , L y ] on the β-plane. The case k = 1 and k = 2 are often referred to as the single-and double-gyre forcing. The stratification is modeled in terms of a two-layer system and the wind stress only directly forces the upper layer. As a response to this wind stress, the system supports a basic shear flow ψ s . The amplitude τ that controls the wind-stress curl, or equivalently the maximal velocity of the shear flow ψ s is chosen as the bifurcation parameter. We first perform a numerical linear stability analysis of this basic shear flow; for small values of τ , all associated eigenvalues have negative real parts such that the jet is stable. When the aspect ratio of the channel a = L y /L x is large, the eigenvalues remain in the left complex plane regardless of the value of τ . However, when the aspect ratio gets small, the basic shear flow loses stability at a critical τ in the form of a pair of single or double complex eigenvalues crossing the imaginary axis, giving rise to a Hopf or double Hopf bifurcation. We next use the idea and method of the dynamic transition theory [START_REF] Ma | Bifurcation Theory and Applications[END_REF][START_REF] Ma | Phase Transition Dynamics[END_REF], which is aimed to determine all the local attractors near a transition. The approach comes with a classification of all transitions into three classes known as continuous, catastrophic and random types. In this way, our study extends previous results using this approach on the single-layer barotropic case [START_REF] Simonnet | Hopf bifurcation in quasi-geostrophic channel flow[END_REF][START_REF] Dijkstra | Dynamic Transitions of Quasi-geostrophic Channel Flow[END_REF], the two-layer case for constant zonal jet velocities [START_REF] Cai | Baroclinic Instability and Transitions in a Two-Layer Quasigeostrophic Channel Model[END_REF] and the barotropic Munk western boundary layer current profile case [START_REF] Kieu | On the Hopf (double Hopf) bifurcations and transitions of twolayer western boundary currents[END_REF], to the zonally periodically bounded two-layer case.

Using the center manifold reduction, we obtain reduced (ordinary differential equation, ODE) models describing this transition. From the coefficients of these ODEs, the transition numbers in dynamic transition theory can be calculated. The case of a Hopf bifurcation is generic while the case of a double Hopf bifurcation is degenerate and requires fine tuning of the aspect ratio to critical values where two pairs of eigenmodes with consecutive wavenumbers cross the imaginary axis simultaneously. Under standard set of parameters describing the midlatitude ocean, we perform numerical computations of the transition number for the forcing patterns corresponding to k = 1, 2, 3 and the aspect ratios a ≥ 3. We find that in the parameter regimes we are interested with, the Hopf bifurcation is supercritical and a stable limit cycle bifurcates. For the double Hopf bifurcation, we find that after the corresponding transition takes place, the system exhibits a bifurcated local attractor [START_REF] Ma | Bifurcation Theory and Applications[END_REF] near the basic shear flow which is homeomorphic to the 3D-sphere. The topological structure of this attractor is analyzed and depending on the parameters, it is found to contain a combination of limit cycles and a quasi periodic solution.

The paper is organized as follows. In section 2, the quasi-geostrophic model is presented. This is followed by section 3 where the theory and numerical results for the linear stability problem (section 3.1), the Hopf bifurcation case (section 3.2) and the double-Hopf bifurcation case (section 3.3) are presented. These results are summarized and discussed in section 4. The appendix contains details on the proofs of the theorems and on the numerical computations.

The model

We consider two layers of homogeneous fluids, each with a different and constant density ρ 1 and ρ 2 and with equilibrium layer thicknesses H 1 and H 2 , on a mid-latitude β-plane with Coriolis parameter is f = f 0 + β 0 y. The lighter fluid in layer 1 is assumed to lie on top of the heavier one in layer 2 so that the stratification is statically stable, i.e., ρ 1 < ρ 2 ; bottom topography is neglected.

This flow can be modeled by the two-layer QG model [START_REF] Pedlosky | Finite-Amplitude Baroclinic Waves[END_REF] using the geostrophic stream function ψ i and the vertical component of the relative vorticity ζ i in each layer (i = 1, 2). The quantities ψ i and ζ i are nondimensionalised by U L y and U/L y , respectively, wind stress with τ 0 , length with L y , and time with L y /U , where U is a characteristic horizontal velocity.

By choosing U = τ 0 /(ρ 0 β 0 L y H 1 ), where ρ 0 is a reference density, the dimensionless equations on the domain (0, 2/a) × (-1, 1) become (2.1)

∂ ∂t + {ψ 1 , •} (∆ψ 1 + F 1 (ψ 2 -ψ 1 ) + βy) = F 1 -τ β sin kπy ∂ ∂t + {ψ 2 , •} (∆ψ 2 + F 2 (ψ 1 -ψ 2 ) + βy) = -r 2 ∆ψ 2
where {f, g} = f x g y -f y g x is the usual Jacobian operator and F 1 represents the damping of upper layer vorticity due to frictional processes (to be specified below). In the bottom layer, we include a linear (Ekman) friction term -r 2 ∆ψ 2 ; in both layers, Laplacian friction terms are neglected due to the absence of continental boundary layers making such terms much smaller than the other ones. The expressions for the dimensional and dimensionless parameters, with their standard values at a latitude 45 • N, are given in Table 1.

For the boundary conditions, we assume periodicity in the x-direction and free-slip boundaries in the ydirection. Hence, the conditions are (2.2)

ψ i | x=0 = ψ i | x=2/a , i = 1, 2. ψ i | y=±1 = ∂ 2 ψ i ∂y 2 | y=±1 = 0, i = 1, 2.
In actual ocean basins, a steady zonal jet is generated by the applied wind stress through Ekman pumping, a Sverdrup balance and a western boundary layer flow [START_REF] Pedlosky | Geophysical Fluid Dynamics[END_REF]. Due to the periodic boundary conditions used here, such a flow cannot be captured in this model. However, the equations will allow a steady state of the form (2.3) ψ s 1 = Ψ sin kπy, ψ s 2 = 0 which relates to the wind stress field, when F 1 -τ β sin kπy = 0. In this paper, we will assume that the wind-stress vorticity input is balanced by vorticity decay due to the linear friction term F 1 = -r 1 ∆ψ 1 , being aware that a larger friction coefficient r 1 is needed than can be justified from existing dissipative processes in the ocean. In this case, it follows that

(2.4) Ψ = τ β (kπ) 2 r 1 .
The parameter Ψ appearing in (2.3) can then be chosen as the control parameter as is the case in this study, instead of τ . By considering the perturbation ψ i = ψ i -ψ s i , i = 1, 2, we can write the system in the following operator form.

(2.5)

M∂ t ψ = N ψ + G(ψ), ψ = (ψ 1 , ψ 2 )
where M and N are the linear operators defined as

(2.6) Mψ =   ∆ψ 1 + F 1 (ψ 2 -ψ 1 ) ∆ψ 2 + F 2 (ψ 1 -ψ 2 )   (2.7) N ψ =   Ψkπ cos kπy (kπ) 2 ∂ψ1 ∂x + F 1 ∂ψ2 ∂x + ∂∆ψ1 ∂x -β ∂ψ1 ∂x -r 1 ∆ψ 1 -ΨkπF 2 ∂ψ2 ∂x -β ∂ψ2 ∂x -r 2 ∆ψ 2  
Lastly, the bilinear nonlinearity is given explicitly by

(2.8) G(ψ) =   -{ψ 1 , ∆ψ 1 + F 1 (ψ 2 -ψ 1 )} -{ψ 2 , ∆ψ 2 + F 2 (ψ 1 -ψ 2 )}  
In terms of function spaces, the operators G and N are the mappings, G :

H 1 → H -1 and N : H 0 → H -1
, where

H 1 = {ψ = (ψ 1 , ψ 2 ) ∈ H 4 (Ω) × H 4 (Ω)| ψ satisfies (2.2)}, H 0 = {ψ = (ψ 1 , ψ 2 ) ∈ H 2 (Ω) × H 2 (Ω)| ψ satisfies (2.2)}, H -1 = L 2 (Ω) 2 .
Here Ω = (0, 2/a) × (-1, 1) and H 4 (Ω), H 2 (Ω), L 2 (Ω) are the usual Sobolev and Lebesgue function spaces endowed with their natural inner products.

These space account for spatial regularity of the solution ψ for which H p (Ω) denotes the space of square-integrable functions that possess pth-order derivatives (in the distribution sense) that are themselves square-integrable; see e.g. [START_REF] Brézis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF].

Results

In this section, we first present the linear stability analysis of the basic shear flow and then we move on to describe the first transitions due to the instabilities, covering both the Hopf and double Hopf bifurcations.

3.1. Linear Stability Analysis. We first investigate the linear stability of the basic solution. For this purpose, we denote the eigenmodes of the linear problem by

ψ m,j (x, y) = e iαmx Y j (y), j ∈ N, m ∈ Z, α m := amπ.
with eigenvalues σ m,j , i.e.

(3.1) σ m,j Mψ m,j = N ψ m,j .

Since the linear operators M and N are real, we have This eigenvalue problem is solved numerically by means of a standard Legendre-Galerkin method; see Appendix A. A typical picture of the spectrum near the criticality is given in Figure 3. This figure shows that many eigenvalues are clustered near the imaginary axis at the critical value of Ψ c as defined below.

σ m,j = σ -m,j , ψ m,j = ψ -m,j , ∀m ∈ N.
We assume (as will be confirmed by the numerical results) that the eigenvalues are ordered so that for each m ∈ Z, σ m,1 has the largest real part among σ m,j , j ∈ N.

For each m ∈ N, we define Ψ m , if it exists, to be the value of Ψ for which the eigenvalue σ m,1 crosses the imaginary axis, that is

(3.2) Re(σ m,1 ) = Re(σ -m,1 ) =      < 0, if Ψ < Ψ m = 0, if Ψ = Ψ m > 0, if Ψ > Ψ m
Hence Ψ = Ψ m defines a neutral stability curve in the a -Ψ plane. In Figure 1, these neutral curves are plotted for zonal wave numbers m = 1, 2, 3, 4.

Our numerical analysis suggests that for m in N, Ψ m exists only for aspect ratios of the basin characteristic lengths smaller than a threshold a m , that is for a < a m . The threshold a m is defined by the vertical asymptote condition, lim

a→am- Ψ m = ∞. Moreover, ∞ > a 1 > a 2 > • • •
We define the critical maximal amplitude Ψ c of the steady state given in (2.

m c = argmin m∈N Ψ m .
The typical structure of the spectrum at the critical parameter Ψ = Ψ c is shown in Figure 3 where a pair of complex conjugate eigenvalues cross the imaginary axis. The eigenvalues on the real axis belong to wavenumber m = 0 and are always stable although in Figure 3 it looks like as if there is an additional critical real eigenvalue. To describe the solutions near the onset of transition Ψ = Ψ c , we define the spatio-temporal function

(3.3) f m (x, y, t) = 2 Re e i Im(σm,1)t ψ m,1 (x, y) , m ∈ Z,
where σ m,1 is the first eigenvalue and ψ m,1 is its associated eigenfunction. The spatial structure of the eigenmodes ψ mc,1 is shown in Figure 2, revealing the well-known 'banana-shaped' patterns characteristics of baroclinic instability.

The values of Ψ c with respect to the aspect ratio a for k = 1, 2, 3 is shown in Figure 5. By the previous remarks, (3.4) lim a→a1-

Ψ c = lim a→a1- Ψ 1 = ∞.
By (3.4), for a > a 1 , the system is linearly stable. As is expected, the neutral stability curves (Figure 1) approach the asymptote Ψ c → ∞ as a increases to the critical aspect ratio a 1 over which the system is linearly stable for all Ψ. The value of Ψ c ≈ 0.09/k for small a for each k = 1, 2, 3, see Figure 5. This value of critical maximal shear velocity corresponds to an upper layer friction, (2.4) which is approximately,

r 1 ≈ τ β k 2 π 2 Ψ c ≈ 1000.0 k ,
which is indeed much larger than can be justified from dissipative processes in the ocean but, as explained in section 2, is needed here to connect the zonal jet top the wind-stress field. The friction term in the lower layer however is physical (Ekman friction) and for this study, it is fixed at r 2 = 5.0, see Table 1. Also, from Figure 1, we see that for small aspect ratios, many modes become unstable as Ψ c is exceeded.

For a < a 1 , the system has a first transition at Ψ = Ψ c and exactly one of the following two principal of exchange of stability (PES) condition holds:

(3.5) Re(σ m,1 ) = Re(σ -m,1 )      < 0, if Ψ < Ψ c = 0, if Ψ = Ψ c > 0, if Ψ > Ψ c if m = m c Re(σ m,1 ) = Re(σ -m,1 ) < 0 if m = m c (3.6) Re(σ m,1 ) = Re(σ -m,1 )      < 0, if Ψ < Ψ c = 0, if Ψ = Ψ c > 0, if Ψ > Ψ c if m ∈ {m c , m c + 1} Re(σ m,1 ) = Re(σ -m,1 ) < 0 if m / ∈ {m c , m c + 1}.
According to (3.5) and (3.6) either, two or four eigenvalues become unstable as Ψ crosses Ψ c . The case of two critical eigenvalues is generic and results in a Hopf bifurcation. The case of four critical eigenvalues results in a double Hopf bifurcation and requires the fine-tuning of the aspect ratio a = a DH so that Ψ c = Ψ mc = Ψ mc+1 . The values of a DH where double Hopf transition occurs is given in Table 2. Although the double Hopf transition is not generic, its analysis gives an insight into the moderate Ψ regime where multiple eigenvalues are unstable. In Figure 4, we present the dominant part of the spectrum of the linearized operator at a critical aspect ratio a DH .

-1 -0.5 0 Our numerical results in Figure 1 and Figure 5 verify the PES condition for the chosen parameter values. The PES condition (3.5) has been rigorously verified for Kolmogorov flows in [START_REF] Meshalkin | Investigation of the stability of a stationary solution of a system of equations for the plane movement of an incompressible viscous liquid[END_REF] via a continued fraction method. This method has later been extended for the single layer QG model for the k = 1 case in [START_REF] Simonnet | Hopf bifurcation in quasi-geostrophic channel flow[END_REF] and for k ≥ 2 in [START_REF] Lu | On the spectral instability and bifurcation of the 2D-quasigeostrophic potential vorticity equation with a generalized Kolmogorov forcing[END_REF] where k is the forcing frequency in (2.1). It is still an open problem to prove this claim for the current problem. Lemma 3.1. Assume that the first critical eigenvalue is complex simple so that the PES condition (3.5) holds. Then the transition and stability of the steady state solution (2.3) of the equation (2.5) in the vicinity of the critical maximal shear velocity Ψ = Ψ c and for any sufficiently small initial condition are equivalent to the stability of the zero solution of the equation

•10 -2 -20 0 20 Re(σ m,j ) Im(σ m,j ) m = -2 m = -1 m = 0 m = 1 m = 2
(3.7) dz dt = σ mc,1 z + P z|z| 2 + o(|z| 3 ).
where z(t) which lies in C denotes the amplitude of the projection of the solution onto the first eigenfunction ψ mc,1 and the complex number P denotes the transition number defined in Eq. (B.8) of Appendix B.

The analysis of Lemma 3.1 yields the following theorem.

Theorem 3.1. Assume that the first critical eigenvalue has simple complex multiplicity, that is the PES condition (3.5) is satisfied. Then the following assertions hold true.

(1) If Re(P ) < 0, the system (2.5) undergoes a continuous transition accompanied by a supercritical Hopf bifurcation on Ψ > Ψ c . In particular, the steady-state solution bifurcates to a stable periodic solution ψ on Ψ > Ψ c , satisfying ψ → 0 as Ψ → Ψ c and has the following approximation

(3.8) ψ(x, y, t) = -Re(σ mc,1 ) Re(P ) 1 2 f mc (x, y, t) + o |Ψ -Ψ c | 1 2
.

The spatial structure of the time periodic solution ψ is shown at t = 0 for different aspect ratios a in Figure 2. (2) If Re(P ) > 0, the system (2.5) undergoes a jump transition on Ψ < Ψ c accompanied by a subcritical Hopf bifurcation. In particular, an unstable periodic orbit ψ given by (3.8) bifurcates on Ψ < Ψ c and there is no periodic solution bifurcating from 0 on Ψ > Ψ c . Moreover, there is a singularity separation at some Ψ s < Ψ c generating an attractor and an unstable periodic orbit ψ.

When the PES condition (3.5) holds, the system exhibits a Hopf bifurcation as described by Theorem 3.1. The type of transition boils down to the determination of the transition number P given in (B.8) of Appendix B. For the practical calculation of this number, we refer to Appendix C. Our analysis given in Figure 6 shows that for low forcing frequencies k = 1, 2, 3, generally Re(P ) < 0 and as a result only continuous transition (supercritical Hopf bifurcation) is possible for the parameter regime we have selected. In Figure 6, we also display the critical wavenumber m c . In the range 4 ≤ a ≤ 20, the critical wavenumber is found to 1 ≤ m c ≤ 4.

T dim (in days) k = 1 k = 2 k = 3

Ly

U compared to the (dimensional) length of the channel L dim = 2L y /a of the both stable and unstable bifurcated time periodic solution (3.8) in the Hopf bifurcation case. Here L y and U are the characteristic scales defined in Table 1. The jumps in the derivative of the time period of the bifurcated solution is due to the change of the imaginary part Im σ mc,1 of the critical eigenvalue at the double Hopf aspect ratios a DH .

There are discontinuities in P vs a plot in Figure 6 of the transition number which are due to the changes in the critical zonal wavenumber m c . These discontinuities take place at double Hopf bifurcation aspect ratios where two consecutive zonal wavenumbers become critical simultaneously which is investigated in the next section. However, there are also discontinuities in Figure 6, k = 1 case (for example near a = 16) whose origin is mysterious.

As detailed in Appendix B, the transition number P accounts for two types of nonlinear interactions between the eigenmodes, and is written

P = P 0 + P 2 ,
where P 0 accounts for nonlinear interactions between the critical modes and the zonally homogeneous modes m = 0 (see (B.9) below), while P 2 accounts for the interactions between the critical modes and the modes having wavenumbers twice that of the critical modes (see (B.10) below). A comparison of typical numerical values of P 0 and P 2 shows that P 2 is several orders of magnitudes smaller than P 0 ; see Figure 8. We refer to [3, Theorem III.1] for a transition number diagnosing also the type of Hopf bifurcations arising, generically, in delay differential equations, and whose nature is also characterized by the interactions of linearized modes through the model's non-linear terms. where P 0,j measures the contribution of the j-th mode with zero wavenumber. Figure 9 shows that the decay of P 0,j as j → ∞ is essentially linear. We believe the results in Figure 8 and Figure 9 may help when choosing the modes to include in a simulation when the maximal shear velocity is well above the criticality. We also compare the dimensional time period of the bifurcated solution (3.8) to the (dimensional) length of the channel in Figure 3.2. With the default parameters as chosen in Table 1, our simulations yield a solution with time period of 180-380 days depending on the channel length of 100-700kms.

Double Hopf Bifurcation.

In this section we are interested in the transitions that take place at the critical aspect ratios a DH where four modes with consecutive wavenumbers m c , m c + 1 become unstable as given by the PES condition (3.6).

We first present the reduced equations in this case (for proofs, see Appendix D).

Lemma 3.2. Assume that the first critical eigenvalues have complex multiplicity 2 so that the assumption (3.6) holds. Then, the transition and stability of the steady state solution (2.3) of the equations (2.5) in the vicinity of the critical Reynolds number Ψ = Ψ c and for any sufficiently small initial condition are equivalent to the stability of the zero solution of the following equations (3.9)

dz 1 dt = σ mc,1 z 1 + z 1 (A|z 1 | 2 + B|z 2 | 2 ) + o(|z|) 3 ), dz 2 dt = σ mc+1,1 z 2 + z 2 (C|z 1 | 2 + D|z 2 | 2 ) + o(|z|) 3 ),
where z 1 (t), z 2 (t) ∈ C denote the amplitudes of the projection of the solution onto the first eigenfunctions ψ mc,1 , ψ mc+1,1 and the transition numbers A, B, C, D are determined by the nonlinear interactions of the first eigenfunction with higher modes given by (D.2). More precisely, the terms A and D account for the self-interactions among the critical modes, while the terms B and C account for the cross-interactions between the critical modes with the higher modes.

T 1 T 2 λ 2 λ 1 I II III IV V VI (a) θ > 0, δ > 0, θδ > 1 T 2 T 1 λ 2 λ 1 I VII III IV V VI (b) θ > 0, δ > 0, θδ < 1 T 2 T 1 λ 2 λ 1 VII VIII VI IX IV V (c) θ > 0, δ < 0
It is known that the equation (3.9) exhibit a zoo of dynamical behaviors. We refer to [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF] for a detailed analysis of all possible cases. Here, we restrict our attention to the case (3.10) Re(A) < 0, Re(B) < 0, Re(B) + Re(C) < 0, Re(D) < 0, which is the only case we observe in our numerical experiments, see Table 2. Under these conditions it is known that the transition is continuous (see Theorem 2.3 in [START_REF] Li | Dynamic Transitions of the Swift-Hohenberg Equation with Third-Order Dispersion[END_REF]). For the next theorem, let us define the numbers (3.11) λ 1 = Re(σ mc,1 ), λ 2 = Re(σ mc+1,1 ).

(3.12)

δ = Re(C) Re(A) , θ = Re(B) Re(D) η 1 = λ 1 -θλ 2 ) (θδ -1) Re(A) 1 2 , η 2 = λ 2 -δλ 1 (θδ -1) Re(D) 1 2
.

Recalling f mc defined in (3.3) (with m = m c ), we define the following spatio-temporal profiles (3.13)

ψ mc p (x, y, t) = - λ 1 Re(A) 1 2 f mc (x, y, t) + o |Ψ -Ψ c | 1 2 , ψ mc+1 p (x, y, t) = λ 2 Re(D) 1 2 f mc+1 (x, y, t) + o |Ψ -Ψ c | 1 2 , ψ qp (x, y, t) = η 1 f mc (x, y, t) + η 2 f mc+1 (x, y, t) + o |Ψ -Ψ c | 1 2
Theorem 3.2. Assume that the conditions of Lemma 3.2 as well as the condition (3.10) hold. Then the equations (2.5) undergo a continuous transition at Ψ = Ψ c , and an S 3 attractor Σ bifurcates on Ψ > Ψ c , which converges to 0 as Ψ ↓ Ψ c . Depending on the values of θ and δ, there are three transition scenarios as shown in Figure 10. In each scenario, near the onset of transition (λ 1 , λ 2 ) = (0, 0), λ 1 -λ 2 plane is dissected into several regions with distinct topological structures for the attractor Σ as given in Figure 11. (2) We note that the features of the spatial structures of upper vs lower layer of the bifurcated periodic solutions in Figure 2 and the quasi-periodic solution in Figure 12 do not alter much. We expect that the situation would be different if bottom topography is included.

The transition scenario of double Hopf transition is given by Theorem 3.2 by Figure 10 and Figure 11. We find that near the onset of transition, depending on the fluctuations the basic state transitions either to a time periodic solution or a quasi periodic solution. Our results in Table 2 show that all three of the scenarios presented in Figure 10 are realizable.

In particular, near a double Hopf transition point, one of the following three possibilities must occur post transition, Ψ > Ψ c :

• there is only a single stable limit cycle,

• there are two distinct stable limit cycles, and an unstable quasi-periodic solution • there is a stable quasi-periodic solution and either one or two unstable limit cycles. For the double Hopf transition, Theorem 3.2 basically tells that all of the above local structures, the time periodic solution and the 2D torus, if they exist, reside in a local attractor homeomorphic to the three dimensional sphere. The existence of this attracting 3D sphere is guaranteed by the attractor bifurcation theorem, Theorem 6.1 in [START_REF] Ma | Bifurcation Theory and Applications[END_REF].

Summary and Discussion

In this paper, we investigated the stability of a parallel zonal jet forced by a Kolmogorov-type wind stress in a periodic zonal channel, using a two-layer quasi-geostrophic (QG) model. This problem is, in terms of complexity, situated between the horizontally unbounded problem [START_REF] Phillips | A simple three-dimensional model for the study of large-scale extratropical flow patterns[END_REF]27] and the fully bounded gyre problem [START_REF] Dijkstra | Temporal variability of the Wind-Driven Quasi-geostrophic Double Gyre Ocean Circulation: Basic Bifurcation Diagrams[END_REF]. More precisely, the effect of boundaries is captured by the interactions of only a few modes (solutions of the linear stability problem), while still keeping a parallel flow which is connected to the wind-stress field.

Our numerical results show that, as expected based on earlier one ( [START_REF] Kieu | On the Hopf (double Hopf) bifurcations and transitions of twolayer western boundary currents[END_REF][START_REF] Cai | Baroclinic Instability and Transitions in a Two-Layer Quasigeostrophic Channel Model[END_REF][START_REF] Dijkstra | Dynamic Transitions of Quasi-geostrophic Channel Flow[END_REF][START_REF] Lu | On the spectral instability and bifurcation of the 2D-quasigeostrophic potential vorticity equation with a generalized Kolmogorov forcing[END_REF]), the zonal shear flow is linearly stable if its maximal amplitude Ψ (or equivalently the maximal amplitude τ of the wind-stress curl) is below a critical threshold Ψ c . Moreover, as this critical threshold is exceeded, generically a Hopf bifurcation occurs. We approach the problem using dynamic transition theory determining all the attractors near a transition. Under characteristic values of parameters which describe the midlatitude ocean we find that a continuous transition in the form of supercritical Hopf bifurcation occurs. The numerical results also show that instead of infinitely many modes which effect the type of transition (supercritical vs subcritical Hopf), the transition is indeed determined by the interaction of the first two critical modes with only the first few zonally homogeneous (m = 0) modes. We also investigated the double Hopf bifurcation scenario which takes place at critical length scales where four modes with consecutive wavenumbers become critical. By a rigorous center manifold analysis, we obtain the coefficients of the 4D-ODE system. Our results show that for the parameters we have considered, there exists a quasi-periodic solution which is a linear combination of two periodic solutions and may be stable depending on the parameters. From a transition point of view, in the double Hopf transition, an attractor homeomorphic to 3D sphere bifurcates. This attractor contains stable/unstable limit cycles and an stable/unstable invariant torus (a quasi-periodic solution).

The results add more detail to the nonlinear development of baroclinic instabilities on a non-constant parallel zonal jet, in that the periodic orbits can become unstable to torus bifurcations and give rise to quasi-periodic behavior. Such a scenario was also found for the barotropic double gyre flow [28], but only in a weakly nonlinear framework using a set of (reduced) amplitude equations. The transition scenario found for the zonally periodic zonal jet is likely to be more relevant for the ocean circulation than the sideband instabilities in the zonally unbounded zonal jet case which require a nearby band of wavenumbers to be unstable.

In our set-up, the linear friction coefficient in the upper layer is relatively large and as explained needed to balance the vorticity input by the wind stress for generating the zonal jet. When this friction is decreased, more modes will become unstable near the critical point and their interaction will allow to give a detailed view on the eddy formation process due to baroclinic instabilities. In this way, the dynamic transition theory approach, possibly with extensions as in [START_REF] Mickaël D Chekroun | Variational Approach to Closure of Nonlinear Dynamical Systems: Autonomous Case[END_REF], provides a way forward to develop a mathematical theory of such ocean-eddy formation processes. Here the linear operators M and N (y) are defined as

(A.3) M =   ∆ m -F 1 F 1 F 2 ∆ m -F 2   , N (y)
N 11 = c 1 cos kπy (kπ) 2 + ∆ m -βiα m -r 1 ∆ m N 12 = c 1 F 1 cos kπy N 21 = 0 N 22 = -c 1 F 2 cos kπy -βiα m -r 2 ∆ m ∆ m = D 2 -α 2 m , D = ∂ ∂y c 1 = Ψkπiα m .
The eigenvalue problem (A.2) is supplemented with the following boundary conditions

(A.5) Y i (±1) = D 2 Y i (±1) = 0, i = 1, 2.
We use Legendre-Galerkin method to discretize and solve the (A.2) with boundary conditions (A.5). We refer to [START_REF] Shen | Spectral Methods[END_REF] for the details of the Legendre-Galerkin method and to [START_REF] Dijkstra | Dynamic Transitions of Quasi-geostrophic Channel Flow[END_REF] for its use in dynamical transition problems.

Let {L j } be the Legendre polynomials and consider compact combinations of the Legendre polynomials

f j (y) = L j (y) + 4 k=1 c jk L j+k (y)
with c jk chosen so that f j satisfy the boundary conditions (A.5), i.e.

f j (±1) = D 2 f j (±1) = 0.
To discretize the eigenvalue problem, we plug

(A.6) Y Ny i (y) = Ny-1 j=0 y (i) j f j (y), Y i = [y (i) 0 , . . . , y (i)
Ny-1 ] T , i = 1, 2. into (A.2) to obtain (A.7) σ   ∆ m -F 1 A 3 F 1 A 3 F 2 A 3 ∆ m -F 2 A 3     Y 1 Y 2   =   N 11 N 12 N 21 N 22     Y 1 Y 2   (A.8) N 11 = c 1 (kπ) 2 A 5 + c 1 A T 4 -α 2 m A 5 -βiα m A 3 -r 1 ∆ m N 12 = c 1 F 1 A 5 N 21 = 0 N 22 = -c 1 F 2 A 5 -βiα m A 3 -r 2 ∆ m
Here (A.9)

A 1 = (D 4 f j , f k ), A 2 = (D 2 f j , f k ), A 3 = (f j , f k ), A 4 = (cos kπyD 2 f j , f k ), A 5 = (cos kπyf j , f k ), ∆ m = A 2 -α 2 m A 3 , with (f, g) = 1
-1 f (y)g(y) dy. The explicit expression of the matrices A i , i = 1, . . . , 5 can be found in [START_REF] Dijkstra | Dynamic Transitions of Quasi-geostrophic Channel Flow[END_REF].

Appendix B. Proof of Lemma 3.1 and Theorem 3.1

We first proceed with the proof of Lemma 3.1 For this, we denote the adjoint modes by ψ * m,j = e iαmx Y * m (y). We denote the critical eigenmode and the critical eigenvalue by

ψ c = ψ mc,1 , σ c = σ mc,1 .
We denote the bilinear operator G as

G(u) = G 2 (u, u) where G 2 (u, v) is linear in each component. Let us define now (B.1) G s (u, v) = G 2 (u, v) + G 2 (v, u).
The center part of the solution is

(B.2) u c = z(t)ψ c + c.c.
where c.c. stands for complex conjugate of the terms before. The evolution of z(t) near the onset of transition is obtained by the projection onto the critical mode ψ c .

(B.3) ż = σ c z + 1 Mψ c , ψ * c G(u c + Φ), ψ * c .
where Φ is the center manifold function. We will obtain its quadratic approximation Φ 2 given by

Φ = Φ 2 (z, z) + o(2)
Here o(n) = o (|(z, z)| n ) denotes higher than n-th order terms in z, z.

Using the notation (B.1), the reduced equation (B.3) can be written

(B.4) ż = σ c z + 1 Mψ c , ψ * c G s (u c , Φ 2 ), ψ * c + o(3).
To obtain a closed system, we need to approximate the center manifold function. The approximation of the center manifold in this case reads, see [START_REF] Taylan | Dynamic Transitions and Baroclinic Instability for 3D Continuously Stratified Boussinesq Flows[END_REF],

(B.5) Φ 2 = (2σ c -L) -1 Π s G 2 (zψ c , zψ c ) + (σ c + σ c -L) -1 Π s G 2 (zψ c , zψ c ) + c.c.
where L = Π s M -1 N and Π s is the projection on the stable space. Using the formula (B.5), we obtain the expansion of the center manifold as

(B.6) Φ 2 = z 2 j≥1 g 2mc,j ψ 2mc,j + |z| 2 j≥1 g 0,j ψ 0,j + c.c.
Here (B.7)

g 0,j = 1 (σ c + σ c -σ 0,j ) Mψ 0,j , ψ * 0,j G 2 (ψ c , ψ c ), ψ * 0,j g 2mc,j = 1 (2σ c -σ 2mc,j ) Mψ 2mc,j , ψ * 2mc,j G 2 (ψ c , ψ c ), ψ * 2mc,j ,
are the coefficients of the center manifold function.

We write (B.4) as (3.7), that is ż = σ c z + P z|z| 2 + o(3). which finishes the proof of Lemma 3.1.

Recalling the definition of G s given in (B.1), the transition number P can then be written as (B.8)

P = P 0 + P 2 ,
where (B.9)

P 0 = j≥1 P 0,j , P 0,j = 1 Mψ c , ψ * c (g 0,j + c.c.) G s (ψ c , ψ 0,j ), ψ * c ,
denotes the contribution of m = 0 modes ψ 0,j while (B.10)

P 2 = j≥1 P 2,j , P 2,j = 1 Mψ c , ψ * c g 2mc,j G s (ψ c , ψ 2mc,j ), ψ * c ,
denotes the contribution of 2m c modes ψ 2mc,j on the transition number respectively. The transition depends on the real part of the transition number P . The proof of Theorem 3.1 follows from the standard Hopf bifurcation analysis of the reduced equation.

Appendix C. Practical aspects for the calculation of the transition number

The practical calculation of the P 0 -term in (B.9) and the P 2 -term in (B.10), boils down to the efficient calculation of the inner and trilinear products involved therein. In that respect, we provide here explicit expressions of the latter. They are given by

Mψ m,j , ψ * m,j = iα m 1 -1 (D 2 -α 2 m -F 1 )Y 1 m,j + F 1 DY 2 m,j Y * 1 m,j dy +iα m 1 -1 (D 2 -α 2 m -F 2 )Y 2 m,j + F 2 DY 1 m,j Y * 2 m,j dy, and
G 2 (ψ m,j , ψ n,k ), ψ * p,l = -δ m+n-p 1 -1 G 1 2 Y * ,1 p,l + G 2 2 Y * ,2 p,l dy, G 1 2 = iα m Y 1 m,j D(D 2 -α 2 n -F 1 )Y 1 n,k + F 1 DY 2 n,k -iα n (D 2 -α 2 n -F 1 )Y 1 n,k + F 1 Y 2 n,k DY 1 m,j G 2 2 = iα m Y 2 m,j D(D 2 -α 2 n -F 2 )Y 2 n,k + F 2 DY 1 n,k -iα n (D 2 -α 2 n -F 2 )Y 2 n,k + F 2 Y 1 n,k DY 2 m,j .
In practice, the integrals can be evaluated by any commonly used quadrature rules in which the values of the integrand are evaluated at quadrature points. In our calculations, we use As the reduction in the case of (3.6) is similar to the case of (3.5) given in the previous section, we will only mention the differences between these two cases. Under the assumption (3.6), we write the center part of the solution as u c = z 1 (t)ψ 1 + z 2 (t)ψ 2 + c.c. where the first two critical modes are

ψ 1 = ψ mc,1 , ψ 2 = ψ mc+1,1 , ψ -1 = ψ -mc,1 , ψ -2 = ψ -mc-1,1 .
with corresponding eigenvalues

σ 1 = σ mc,1 , σ 2 = σ mc+1,1 , σ -1 = σ -mc,1 , σ -2 = σ -mc-1,1
The equation (B.4) becomes the system

(D.1) żj = σ j z j + 1 Mψ j , ψ * j G(u c + Φ 2 ), ψ * j + o(3), j = 1, 2.
and the center manifold function (B.5) is replaced by

Φ 2 = 2 |j|,|k|=1 z j z k Ψ j,k + o(2), Ψ j,k = (σ j + σ k -L) -1 Π s G 2 (ψ j , ψ k ).
Now, the equations (D.1) become (3.9) with the coefficients defined as below.

(D.2)

A = g 1,1,-1,1 + g 1,-1,1,1 + g -1,1,1,1 B = g 1,
g i,j,k,l = 1 Mψ l , ψ * l G s (ψ i , Ψ j,k , ψ * l ) , Ψ j,k = (σ j + σ k -L) -1 Π s G 2 (ψ j , ψ k )
We note that the above coefficients contain only g i,j,k,l for which i + j + k = l. The expansion of the center manifold coefficients can be written more explicitly as

Ψ j,k = ∞ i=1 G 2 (ψ mj , ψ m k ), ψ * mj +m k ,i Mψ mj +m k ,i , ψ * mj +m k ,i (σ j + σ k -σ mj +m k ) -1 ψ mj +m k ,i
Now we analyze the equations (3.9) by first putting them in polar form .

Since the Jacobian matrix of the right hand side of (D.3) at the quasi-periodic solution has determinant -4(σ 2 -δσ 1 )(σ 1 -θσ 2 ) θδ -1 .

With this information, the transition scenarios summarized in Figure 10 and Figure 11 can be obtained by a standard analysis. To prove the claim on the bifurcation of an S 3 -attractor, we need to prove that (ρ 1 , ρ 2 ) = (0, 0) is locally stable equilibrium of (D.

Figure 1 .

 1 Figure 1. Neutral stability curves Ψ m , m = 1, 2, 3, 4 defined in (3.2), for wind-stress profiles defined by k = 1, 2, 3. For values of Ψ > Ψ m , the shear flow Ψ s becomes unstable to a perturbation pattern with wavenumber α m .

3 )Figure 2 .

 32 Figure 2. Real part of the upper and lower layers of the time-periodic solution f mc (or equivalently of the dominant eigenmode ψ mc,1 ) at t = 0 for k = 1 and a = 10 and a = 5, respectively.

Figure 3 .

 3 Figure 3. The first 240 eigenvalues at the critical parameter when m c = 2, a = 10, k = 1 and N y = 240.

Figure 4 .

 4 Figure 4. The spectrum near the double Hopf aspect ratio a DH = 11.263 for k = 1. The first four critical eigenvalues σ m,1 , m ∈ {-2, -1, 1, 2}, can be seen on the imaginary axis.

Figure 5 .

 5 Figure 5. The critical maximal shear steady state velocity Ψ c with respect to channel aspect ratio a for k = 1, 2, 3 where the PES condition (3.5) holds.

3. 2 .

 2 Hopf Bifurcation. We first investigate the generic Hopf transition scenario based on the attractor bifurcation theorem [15, Theorem 5.2] and the dynamical transition theorem[START_REF] Ma | Phase Transition Dynamics[END_REF] Theorem 2.1.3]. For proofs of the following lemma and theorem, see Appendix B).

3 m c = 1 m c = 2 m c = 3 m c = 4 Figure 6 .L

 346 Figure 6. The real part of the transition number P compared to the channel aspect ratio a normalized by the largest absolute value of Re P . The parameters are as set in Table1.

Figure 7 .

 7 Figure 7. The (dimensional) time period T dim = 2π Im(σm c ,1)

Figure 8 . 2 Figure 9 .

 829 Figure 8. The effect of P 0 (of m = 0 modes) compared to the effect P 2 (of 2m c modes) on the transition number.

Figure 10 .

 10 Figure 10. The regions in the λ 1 -λ 2 plane with different dynamical behaviors. In region V, the basic steady state is locally asymptotically stable. In regions IV and VI, the system undergoes a supercritical Hopf bifurcation. The dynamics in regions I, II and III is the double Hopf bifurcation scenario and the details are given in Figure11. The lines T 1 and T 2 in the Figure have slopes 1/θ and δ (as defined in (3.12)) respectively.

Remark 3. 1 . ( 1 )Figure 11 .

 1111 Figure 11. The dynamics in the regions given in the first quadrant of Figure 10. ψ mc p , ψ mc+1 p are time-periodic with zonal wavenumbers m c and m c + 1 respectively. ψ qp is the quasiperiodic solution given in (3.13).

Figure 12 .

 12 Figure 12. Real part of the upper and lower layers of the time quasi-periodic solution ψ qp given in Theorem 3.2 at t = 0 for k = 1 and a = 5.722 where the first two critical wavenumbers are m c = 1 and m c = 2.

Appendix A .

 . Numerical treatment of the linear stability problem To solve the eigenvalue problem numerically, we first plugin the ansatz (A.1) ψ(x, y) = e iαmx Y j (y), j ∈ N, m ∈ Z, α m := amπ. into the eigenvalue problem σMψ(x, y) = N (y)ψ(x, y), to obtain (A.2) σ MY (y) = N (y)Y (y), Y (y) = (Y 1 (y), Y 2 (y))

1 - 1 f

 11 (y) dy = Ny n=0 f (y n )ω n , where y n and ω n are Legendre-Gauss-Lobatto quadrature points and weights respectively. Appendix D. Proof of Lemma 3.2 and Theorem 3.2

z j = ρ i e iγj , j = 1 1 + 1 +κ 1 = -2σ 1 , κ 2 = σ 2 -δσ 1 κ 1 = -2σ 2 , κ 2 = σ 1 -θσ 2

 11111211222 Re(σ 1 )ρ 1 + ρ 1 (Re(A)ρ 2 Re(B)ρ 2 2 ) + h.o.t. ρ2 = Re(σ 2 )ρ 2 + ρ 2 (Re(C)ρ 2 Re(D)ρ 2 2 ) + h.o.t. and γ1 = Im(σ 1 ) + h.o.t. γ2 = Im(σ 2 ) + h.o.t.For the specific case of (3.10), the equations (D.3) always admit the solutions which represent the periodic solutions (ρ 1 , ρ 2 ) = -Also, the equations (D.3) admit the following solution which represents a quasi-periodic solution(ρ 1 , ρ 2 ) = σ 1 -θσ 2 Re(A)(θδ -1), σ 2 -δσ 1 Re(D)(θδ -1)

H 2 1 f- 2 ρ 0 reference density 10 3 kgm - 3 F 1 0 r 2

 2123102 3) at Ψ = Ψ c , that is when Re(σ 1 ) = Re(σ 2 ) = 0. In this case by assumption (3.10), from (D.3), we can obtain d dt upper layer depth 750m∆ρ density difference ρ 2 -ρ 1 = 1 kgm -3 U characteristic velocity τ 0 /(ρ 0 H 1 β 0 L y ) = 0.02 m s -1 τ 0 characteristiczonal wind stress 0.1 Pa β 0 planetary vorticity gradient 2 × 10 -11 (ms) -upper layer Froude number f 2 0 L 2 y /(g H 1 ) = 4000 F 2 lower layer Froude number f 2 0 L 2 y /(g H 2 ) = 4000/3 β planetary vorticity factor β 0 L 2 y /U = 1000 τ wind-stress parameter 1.lower layer linear friction coefficient r 2 = 0 L y /U = 5.0

  2,-2,1 + g 1,-2,2,1 + g 2,1,-2,1 + g 2,-2,1,1 + g -2,1,2,1 + g -2,2,1,1 C = g 2,1,-1,2 + g 2,-1,1,2 + g 1,2,-1,2 + g 1,-1,2,2 + g -1,2,1,2 + g -1,1,2,2 D = g 2,2,-2,2 + g 2,-2,2,2 + g -2,2,2,2

Table 1 .

 1 Model parameters. k m c a DH Ψ c Re(A) Re(B) Re(C) Re(D) δ θ δθ 1 1 11.26 0.10 -0.25 -1.00 0.01 -0.33 -0.05 3.08 -0.15 1 2 7.30 0.10 -0.37 -1.00 -0.13 -0.40 0.36 2.51 0.89

	1 3 5.38 0.10 -0.95 -1.00 -0.22 -0.61 0.23 1.65 0.37
	1 4 4.23 0.09 -0.81 -1.00 -0.29 -0.61 0.36 1.63 0.58
	2 1 11.18 0.05 -0.15 -1.00	0.12	-0.23 -0.77 4.39 -3.38
	2 2 7.20 0.05 -0.34 -1.00 -0.04 -0.46 0.13 2.20 0.28
	2 3 5.27 0.05 -0.37 -1.00 -0.14 -0.45 0.38 2.24 0.86
	2 4 4.14 0.05 -0.38 -1.00 -0.21 -0.44 0.55 2.25 1.25

Table 2 .

 2 The double Hopf transition numbers. Here k = 1, 2 is the wavenumber of the forcing, m c and m c + 1 are the wavenumbers of the first two critical modes which become unstable simultaneously at the critical aspect ratio a DH . A, B, C, D are the coefficients of the double Hopf transition normalized with respect to the maximum absolute value of those coefficients. The parameters δ and θ are defined by (3.12).

Acknowledgments: This work has been partially supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 810370). This study was also supported by a Ben May Center grant.