Daniel Wagner
email: wagneda1@fel.cvut.cz

Didier Henrion
email: henrion@laas.fr

Martin Hromčík
email: hromcik@fel.cvut.cz

† D Henrion

Measures and LMIs for Validation of an Aircraft with MRAC and Uncertain Actuator Dynamics

It is well known that actuator dynamics with degraded performance can lead to poor stability and tracking for adaptive control. In this paper, we use occupation measures and LMI relaxations (called the moment sums of squares or Lasserre hierarchy) for verification and validation of a longitudinal F-16 aircraft model with MRAC and higher order actuator dynamics with uncertain parameters and deflection saturation. These uncertain parameters are represented explicitly in the space of occupation measures. To improve numerical scalability, we exploit sparsity for ordinary differential equations (ODEs) using parsimony to partition the dynamics. We finally compare these main results to those obtained using Monte-Carlo methods. I. Nomenclature α(t) = angle of attack [rad] q(t) = pitchrate [rad/s] v(t) = elevator actuator output [deg] x p (t) = short-period longitudinal dynamics x c (t) = third order actuator dynamics

II. Introduction

A challenge in the design of adaptive control laws is addressing system stability and acceptable command following in the presence of uncertain, unmodeled actuator dynamics. In particular, degraded actuator performance can lead to poor tracking and even instability [START_REF] Drydek | Adaptive Control and the NASA X-15-3 Flight Revisited[END_REF]. In real world applications, actuators seldom follow a first order model and their dynamics cannot be safely neglected as they can limit achievable stability of model reference adaptive control (MRAC). Therefore, it remains imperative that there exists a way to ensure safe interaction between the higher order actuator and the MRAC in the design phase.

Current literature attempts to address this by using MRAC control modifications. See [START_REF] Johnson | Pseudo-Control Hedging : a New Method for Adaptive Control[END_REF][START_REF] Johnson | Limited authority adaptive flight control for reusable launch vehicles[END_REF] where the authors consider a pseudo control hedging modification to include actuator dynamics in the design phase of MRAC. The authors of [START_REF] Gruenwald | An LMI-Based Hedging Approach to Model Reference Adaptive Control with Actuator Dynamics[END_REF][START_REF] Gruenwald | Computing Stability Limits for Adaptive Control Laws with High-Order Actuator Dynamics[END_REF] revisit pseudo-control hedging with a novel Lyapunov linear matrix inequality (LMI) method to ensure ideal bounded reference trajectories for a range of admissible bandwidths, natural frequency, and damping inherent in the actuator dynamics.

As discussed mainly in [START_REF] Gruenwald | Computing Stability Limits for Adaptive Control Laws with High-Order Actuator Dynamics[END_REF], the relationship between the damping and natural frequency can be bilinear. In other words, any variation in the damping can adversely affect the range of admissible natural frequencies required to maintain closed-loop stability for higher order actuators and vice-versa. There is also a relationship between the MRAC learning rate and the closed-loop stability. Consequently, traditional Monte-Carlo methods can become costly for systems with these parameter uncertainties inherent in their dynamics. Every combination of these parameters must be simulated, and there remains risk of state-space explosion or unexplored, unsafe trajectories for medium and large scale systems.

Moment sum of squares (SOS) is a state-of-the-art method for verification and validation (VV) in aerospace. See [START_REF] Chakraborty | Nonlinear Robustness Analysis Tools for Flight Control Law Validation and Verification[END_REF][START_REF] Chakraborty | Susceptibility of F/A-18 Flight Controllers to the Falling-Leaf Mode: Nonlinear Analysis[END_REF] where the authors focus on polynomial dynamical models and polynomial SOS Lyapunov candidate function. Occupation measures and LMI relaxations (called the moment sums of squares or Lasserre hierarchy) have also been used in VV of aircraft models with MRAC [START_REF] Wagner | Measures and LMIs for Adaptive Control Validation[END_REF]. See also more recent development in using this framework where the authors [START_REF] Wagner | Measures and LMIs for Lateral F-16 MRAC Validation[END_REF][START_REF] Tacchi | Approximating regions of attraction of a sparse polynomial differential system[END_REF] exploit sparsity for ODEs using parsimony. Specifically, scaling issues brought about by optimizing directly over medium to large scale systems can be addressed by partitioning the dynamics.

The procedure remains similar to [START_REF] Tacchi | Approximating regions of attraction of a sparse polynomial differential system[END_REF][START_REF] Henrion | Measures and LMI for space launcher robust control validation[END_REF][START_REF]Optimization on linear matrix inequalities for polynomial systems control[END_REF]. First, our control law validation problem can be written as a polynomial dynamical optimization problem similar to that in optimal control. This problem is then written as its equivalent infinite dimensional linear programming (LP) problem by introducing occupation measures supported over admissible trajectories. Lastly, we relax the infinite dimensional LP problem of measures to a finite dimensional moment LMIs of truncated sequences.

The main contributions are the following:

• We validate an F-16 longitudinal model with MRAC and uncertain third order elevator dynamics with deflection saturation. The elevator contains bounded uncertainties inherent in its damping and natural frequency. Since uncertain parameters have an explicit representation in the space of occupation measures, we can identify numerically the worst case behavior caused by inherent uncertainties in damping and natural frequencies coming from the actuator. To avoid scaling problems, we exploit sparsity for ODEs to partition out the dynamics into a class of smaller systems that serve as control inputs for each other. In total, we are able to validate medium and large scale systems (up to 11 states) using this approach.

• Two examples are considered where the MRAC tolerates uncertainty in the natural frequency and damping, and one example where it obviously fails. Then, we compare our main results to those achieved using Monte-Carlo simulations and a search algorithm. The search algorithm locates the worst case damping and natural frequency behavior from the actuator. Unlike the results of [START_REF] Johnson | Pseudo-Control Hedging : a New Method for Adaptive Control[END_REF][START_REF] Gruenwald | An LMI-Based Hedging Approach to Model Reference Adaptive Control with Actuator Dynamics[END_REF][START_REF] Gruenwald | Computing Stability Limits for Adaptive Control Laws with High-Order Actuator Dynamics[END_REF], we do not use the the pseudo control hedging modification to address actuator bandwidth uncertainty. Furthermore, none of these results take into account deflection saturation of the actuator. The results in [START_REF] Chakraborty | Nonlinear Robustness Analysis Tools for Flight Control Law Validation and Verification[END_REF][START_REF] Chakraborty | Susceptibility of F/A-18 Flight Controllers to the Falling-Leaf Mode: Nonlinear Analysis[END_REF] provide only ultimate stability for aircraft models with linear feedback. Our previous work [START_REF] Wagner | Measures and LMIs for Adaptive Control Validation[END_REF][START_REF] Wagner | Measures and LMIs for Lateral F-16 MRAC Validation[END_REF] do not take into account uncertainties with a polynomial structure. In other words, our framework guarantees finite-time convergence for non-autonomous systems with piecewise dynamics. This is solution better than the solution provided by Barbalat's lemma [START_REF] Haddad | Nonlinear dynamical systems and control: A Lyapunov-based approach[END_REF] provided by the readily available Lyapunov certificate.

The organization of this work is as follows: Section III contains some necessary mathematical preliminaries, section IV discusses the short period of the F-16 with third order actuator dynamics, section V discusses our main results, and section VI gives our main conclusions and suggestions for future research.

III. Mathematical Preliminaries

We begin by very briefly stating some standard notation used for occupations measures. If X is a compact subset of R n , C (X) denotes the space of continuous functions on X and M (X) (resp., M + (X)) denotes the cone of (resp., non-negative) measures. We also assume that there is no relaxation gap (J = J ∞ ) between the polynomial dynamical optimization problem and its representation in the infinite dimensional space of measures. For the proceeding sections, we rely on the main theoretical results discussed extensively in [START_REF] Wagner | Measures and LMIs for Adaptive Control Validation[END_REF][START_REF] Wagner | Measures and LMIs for Lateral F-16 MRAC Validation[END_REF]. In particular, consider the polynomial dynamical optimization with piecewise polynomial different constraints

J = inf h T , h h T (T, x(t)) + ∫ T 0 h(t, x(t))dt s.t. x(t) = f j (t, x(t)), x(t) ∈ X j , j = 1, . . . , N x(0) ∈ X 0 , x(t) ∈ X T , t ∈ [0, T], (1)
with given costs h, h T ∈ R[t, x], polynomial dynamics f j ∈ R[t, x], j = 1, . . . , N, state trajectory x(t) constrained in compact basic semialgebraic sets X j , and initial (resp. terminal) given compact basic semialgebraic sets X 0 (resp., X T ). We also assume that the state space partitioning sets, or cells X j , are such that all of their respective intersections have zero Lebesgue measure and all belong to a given compact semialgebraic set X.

The nonconvex optimization problem eq. ( 1) can be approximated as a generic convex infinite dimensional LP problem of measures using parsimony

J ∞ = inf ∫ h T (T, x(t))dµ T + ∫ h(t, x(t))dµ s.t. ∂ µ ∂t + divf j µ j + µ T = µ 0 ∂ν ∂t + divf k ν k + ν T = ν 0 N j=1 π t,y# µ j = M k=1 π t,y# ν k ∫ µ 0 = 1, ∫ ν 0 = 1, ( 2
)
where div is the divergence operator and the infimum is with respect to the occupation measure µ, ν

∈ M + ([0, T] × X), initial measure µ 0 , ν 0 ∈ M + ({0} × X 0 ), terminal measure µ T , ν T ∈ M + ({T } × X T ), terminal time T > 0.
Each measure µ j (resp. ν k ) supported on their respective cell X j (resp., V k ) so that the global occupation measure becomes

µ = N j=1 µ j , ν = M k=1 ν k ,
with marginal π t,y# µ (resp., π t,y# ν) of measure µ (resp., ν) with respect to variables t, y. As noted in [START_REF] Tacchi | Approximating regions of attraction of a sparse polynomial differential system[END_REF], this approach allows some flexibility with the computational limit caused by the largest moment SDP block. With appropriate partitioning strategies, we can solve problems with medium to large number of states.

IV. F-16 Short Period with Actuator Dynamics

The short-period mode of an F-16 aircraft can be expressed by the dynamics

x p (t) = Ax p (t) + Bv(t), x p (0) = x p0 (3)
where x p (t) = α(t) q(t) are the short period dynamics, α(t) ∈ R is the angle of attack, q(t) ∈ R is the pitchrate,

A = -1.0189 +0.9051 +0.8223 -1.0774 , B = -0.0022 -0.1756
are the open loop short period dynamics. Elevator actuator output v(t) ∈ R is described the third order unmodeled actuator dynamics

x c (t) = (F + L(k 1 , k 2 )) x c (t) + Gu(t) (4)
v(t) = Hx c (t), v(0) = v 0 ( 5
)
where u(t) ∈ R is a measurable control input,

F =        0 1 0 0 0 1 -aω 2 n ω 2 n + 2ζω n a; 2ζω n + a        , G =        0 0 1        , H =        aω 2 n 0 0        , ζ
, ω n , a ∈ R + are the damping, natural frequency, and gain respectively,. Matrix

L =        0 1 0 0 0 1 -ak 2 1 k 2 1 + 2k 2 k 1 a; 2k 2 k 1 + a        contains inherent
uncertainties within damping and natural frequency of the actuator dynamics described by bounded parameters k 1 , k 2 ∈ -k max k max , k max ∈ R + . We assume that resulting matrix (F + L) is Hurwitz, the DC gain between the control input u(t) and actuator output v(t) is unity, and pair (A, B) is controllable.

A. Model Reference Adaptive Control

Given a measurable command signal c(t), we want to design a control law

u(t) = u n (t) + u a (t) (6)
such that the combined nominal u n (t) and adaptive control u a (t) laws allow eq. ( 3) to asymptotically track reference model x r (t) = A r x r (t) + B r c(t), x r (0) = x r0 [START_REF] Chakraborty | Susceptibility of F/A-18 Flight Controllers to the Falling-Leaf Mode: Nonlinear Analysis[END_REF] where x r (t) = α r q r and the error dynamics

e(t) = x p (t) -x r (t) (8)
satisfies lim t→∞ e(t) = 0. Matrices A r = (A -BK 1 ) and B r = BK 2 are from the nominal controller

u(t) = K 1 x p (t) + K 2 c(t) (9)
with feedback/feedforward gains

K 1 = +4.7432 -2.3163 , K 2 = -4.3396
obtained using the LQR method [START_REF] Fravolini | A Model Reference Adaptive Control Approach for Uncertain Dynamical Systems with Strict Component-wise Performance Guarantees[END_REF][START_REF] Stevens | Aircraft control and simulation: Dynamics, controls design, and autonomous systems: Third edition[END_REF]. We assume that there exists a feedback matrix K 1 such that matrix A r is Hurwtiz.

The adaptive controller

u a (t) = -ŴT (t)Φ x p (t) , (10)
with the given basis function Φ i (x p (t)) = (1 + e x pi ) -1 , i = 1, 2 and weight estimate Ŵ(t), satisfies the update law

Ŵ(t) = ΓΦ x p (t) e T (t)PB p , Ŵ(0) = Ŵ0 (11)
where Γ = diag 10 , << 1 is the learning rate and positive definite matrix P is generated by the Lyapunov equality 0 = A T r P + P A r + I 2×2 .

Theorems that highlight the boundedness and long range stability of MRAC are discussed extensively in [START_REF] Haddad | Nonlinear dynamical systems and control: A Lyapunov-based approach[END_REF][START_REF] Lavretsky | Robust and Adaptive Control With Aerospace Applications[END_REF]. For some discussion on control modifications, refer to [START_REF] Narendra | A New Adaptive Law for Robust Adaptation Without Persistent Excitation[END_REF]. For papers that discuss using these control modifications, refer to [START_REF] Wagner | Measures and LMIs for Adaptive Control Validation[END_REF][START_REF] Wagner | Measures and LMIs for Lateral F-16 MRAC Validation[END_REF].

V. Main Results

We wish to validate our existing closed loop polynomial aircraft model eq. ( 3) with actuator dynamics eq. ( 5) by finding the initial state that maximizes of the norm of the concave cost function J =ê(T) 2 2 where

ê(t) = α(t) -α r (t) q(t) -q ss (t) ≈ e(t)
approximates the error dynamics and q ss = lim t→∞ q r (t) . We also have given terminal time T = 25 s, c(t) = 0.1 and actuator constants

a = 5, ω n = 20, ζ = 0.1,
which are all guaranteed feasible solutions from [START_REF] Gruenwald | Computing Stability Limits for Adaptive Control Laws with High-Order Actuator Dynamics[END_REF]. If we can show that for every chosen initial state

x p (0) x c (0) x r (0) Ŵ(0) x(0) ∈ X 0 -0.2 0.2 2 × -1 × 10 -6 1 × 10 -6
6 that all trajectories remain bounded in the box

x p (t) x c (t) x r (t) Ŵ(t) x(t) ∈ X -2 2 7 × -3 3
until they reach the terminal cost target ê(T) ∈ {J ≤ 4 × 10 -3 }, then our control law is validated. For our main results, we consider a 1% and 10% increase in the uncertainty the elevator dynamics by increasing k max . We also include an example where the MRAC fails due to severe degraded actuator performance invoked by a 5000% increase in the dynamic uncertainty.

To include saturation in the elevator deflection, we partition the dynamics using locally affine functions in three cells X j , j = 1, . . . , 3 corresponding to linear, upper saturation, lower regimes

X 1 {x(t) ∈ R 8 : |v(t)| ≤ θ max }, v(t) = Hx c (t) (13)
X 2 {x(t) ∈ R 8 : v(t) ≥ θ max }, v(t) = θ max (14)
X 3 {x(t) ∈ R 8 : v(t) ≤ -θ max }, v(t) = -θ max (15)
where θ max = 25π 180 . With the combined short period dynamics eq. ( 3), actuator dynamics, in-state actuator uncertainties k 1 and k 2 , and controller dynamics all partitioned into the saturation regimes, we can write eq. ( 1) as

J = inf ê(T ) -ê(T) 2 2 s.t. x p (t) = f p (t, x p (t), v(t)), x r (t) = f r (t, x r (t)), x c (t) Ŵ(t) = f j (t, x c (t), x p (t), Ŵ(t), α r , k 1 , k 2 ), x(t) ∈ X j , k 1 , k 2 ∈ -k max k max , j = 1, . . . , 3 x(0) ∈ X 0 , x(t) ∈ X T , t ∈ [0, T], (16)
and the measure LP of parsimony eq. ( 2) as

J ∞ = inf µ T - ∫ ê(T) 2 2 dµ T s.t. ∂ µ ∂t + divf p µ(t, x p , v) + µ T = µ 0 ∂ν ∂t + divf r ν(t, x r ) + ν T = ν 0 ∂ξ ∂t + divf j ξ(t, x c , Ŵ, x p , α r , k 1 , k 2 ) + ξ T = ξ 0 π t,x p ,v# µ = π t,x p ,v# ν π t,x r # ν = 3 j=1 π t,x r # ξ ∫ µ 0 = 1, ∫ ν 0 = 1, ∫ ξ 0 = 1 ( 17
)
with respect to each occupation measure µ, ν, ξ. Marginal π t,x p ,v# µ (respectively, π t,x p ,v# ν) of measure µ (respectively, ν) with respect to variables t, x p , and v. Marginal π t,x r # µ (respectively, π t,x r # ξ) of measure µ (respectively, ξ) with respect to variables t, x r . As discussed in section III we use this strategy to address scaling issues inherited from the problem. With this approach, we can reduce the overall size of the problem by 3 variables. As discussed in [START_REF] Wagner | Measures and LMIs for Adaptive Control Validation[END_REF][START_REF] Wagner | Measures and LMIs for Lateral F-16 MRAC Validation[END_REF], we also employ some basic scaling strategies so that all dynamics are normalized on the interval -1 1 .

Table 1 Gloptipoly 3 + MOSEK Upper Bounds for section V

The Monte-Carlo use a fixed time step and evenly distributed initial conditions to simulate the trajectories. For every simulation, we also use a search algorithm to find the matrix L(k 1 , k 2 ) that maximizes the terminal cost J. The worst case Monte-Carlo simulations and upper bounds are given in Figs. 1 to 3 and Table 2, respectively,. Worst case parameters for natural frequency and damping are also shown in the captions.

As shown in Table 1, the upper bounds obtained using our framework (Gloptipoly 3 + MOSEK) show the MRAC can tolerate a 0.1% and 10% uncertainty in damping and natural frequency. At k max = 50, the actuator performance degrades significantly and tracking is lost. We can also extract the time spent in the saturation regime, by extracting the mass of each occupation measure [START_REF] Henrion | Measures and LMI for space launcher robust control validation[END_REF].

VI. Conclusion

We validated a short period F-16 model with MRAC and a higher order elevator dynamics. The elevator contained parameter uncertainties and deflection saturation. We compared the results using our framework with those obtained using Monte-Carlo and a search algorithm. We also showed an example where the MRAC obviously fails due to poor performance from the elevator. For future results, we wish to find new ways to partition the dynamical problem to solve even larger scale systems and use region of attraction framework to find the numerical upper bounds for the uncertain parameters.

Fig. 1 F

 1 Fig. 1 F-16 Monte-Carlo (k max = 0.01, k 1 = -0.01, k 2 = 0.01)

1 3. 13 ×

 113 10 -5 7.96 × 10 1 50 2.27 × 10 -2 7.43 × 10 1

Fig. 2 F

 2 Fig. 2 F-16 Monte-Carlo (kmax = 0.1, k 1 = -0.1, k 2 = 0.1)

Fig. 3 F

 3 Fig. 3 F-16 Monte-Carlo (k max = 50, k 1 = -50, k 2 = 0)

Table 2 Monte-Carlo Upper Bounds for section V

 2


Chloë Bonnenfant
email: chloe.bonnenfant@supagro.fr

Nathalie Gontard

Chahinez Aouf

Active packaging : Incorporation of polyphenols in polyhydroxyalkanoate (PHA) : Thermal stabilization and antioxidant properties

Conclusion:

We expect that this PhD research will lead to a new formulation of PHA derivatives with enhanced thermal stability. Polyphenols

Introduction : Pollution of environment by plastic is a worldwide issue. Today, packaging is mainly produced from conventional petroleum-based plastics especially in food industry. Plastic is everywhere, that is why ongoing research are attempting to develop biobased and biodegradable polymers. Polyhydroxyalkanoates (PHA) seem to be good candidates. The aim of this study is to enhance their thermal properties and confer antioxidant properties to PHA by adding polyphenols. Evolution of Mw with the number of extrusion process Thermogravimetric analysis of neat PHB (solid line), PHB/TA5 (dotted line), PHB/TA15 (broken line) [2] TA5 = Tannic acid 5% TA15=Tannic acid 15%

  Zainuddin, S., et al. (2019). "Mechanical, Fire Retardant, Water Absorption and Soil Biodegradation Properties of Poly(3-hydroxy-butyrate-co-3valerate) Nanofilms." Journal of Polymers and the Environment 27(10): 2292-2304. [2] Auriemma, M., et al. (2015). "Blending poly(3-hydroxybutyrate) with tannic acid: Influence of a polyphenolic natural additive on the rheological and thermal behavior." European Polymer Journal 63: 123-131.


Valentin Antuori
email: vantuori@laas.fr

Emmanuel Hebrard
email: hebrard@laas.fr

Marie-José Huguet
email: huguet@laas.fr

Siham Essodaigui
email: siham.essodaigui@renault.com

Alain Nguyen
email: alain.nguyen@renault.com

Leveraging Reinforcement Learning, Constraint Programming and Local Search: A Case Study in Car Manufacturing

Keywords: Constraint Programming, Reinforcement Learning, Local Search, Scheduling, Traveling Salesman Problem

The problem of transporting vehicle components in a car manufacturer workshop can be seen as a large scale single vehicle pickup and delivery problem with periodic time windows. Our experimental evaluation indicates that a relatively simple constraint model shows some promise and in particular outperforms the local search method currently employed at Renault on industrial data over long time horizon. Interestingly, with an adequate heuristic, constraint propagation is often sufficient to guide the solver toward a solution in a few backtracks on these instances. We therefore propose to learn efficient heuristic policies via reinforcement learning and to leverage this technique in several approaches: rapid-restarts, limited discrepancy search and multi-start local search. Our methods outperform both the current local search approach and the classical CP models on industrial instances as well as on synthetic data.

Introduction

Improving the production line is a constant concern in the industry. The car manufacturer Renault has long been interested in models and techniques from Operations Research and Constraint Programming to tackle the various routing and scheduling problems arising from the production process.

Recent advances in Artificial Intelligence and in particular in Machine Learning (ML) open up many new perspectives for solving large scale combinatorial optimization problems with promising results popularized by the success of Al-phaGo and AlphaZero [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]. In particular, several approaches combining reinforcement learning and deep neural networks to guide the reward strategy have been proposed for solving the traveling salesman problem (TSP) [START_REF] Bello | Neural Combinatorial Optimization with Reinforcement Learning[END_REF][START_REF] Dai | Learning Combinatorial Optimization Algorithms over Graphs[END_REF]. Moreover, the combination of ML and classical combinatorial techniques seems very promising. For instance, the integration of standard TSP heuristics with a neural network heuristic policy outperforms the pure ML approaches [START_REF] Michel Deudon | Learning Heuristics for the TSP by Policy Gradient[END_REF].

In Section 2 we introduce the problem of planning the flow of vehicle components across the assembly lines. More precisely, given fixed production cycles, logistics operators are in charge of collecting components from, and delivering them to, working stations in such a way that there is no break in the manufacturing process. This problem is in many ways similar to the Torpedo Scheduling Problem [START_REF] Martin | Solving the Torpedo Scheduling Problem[END_REF], from the 2016 ACP challenges [START_REF] Schaus | The Torpedo Scheduling Problem[END_REF]: pickup and delivery operations are to be scheduled, and these operations are repeated in time because the commodity is being produced constantly at a fixed rate. However, in our problem, components are carried over using trolleys that can be stacked into a "train" that should not exceed a given maximal length. It should be noted that in this paper we consider the problem associated to a single operator whose route is to be optimized. However, the more general problem for Renault is to assign components (or equivalently working stations) to operators as well as to plan the individual routes, and the longer term objective is to proactively design the layout of the assembly line to reduce the cost of logistics operations.

We first evaluate in Section 3 two basic constraint programming (CP) models in Choco [START_REF] Prud'homme | Choco Documentation[END_REF] and compare them to the current method used by Renault: a local search (LS) method implemented in LocalSolver4 . From this preliminary study, we observe that although the problem can be hard for both approaches, CP shows promising results compare to LS. Moreover, if solving the problem via backtracking search seems very unlikely given its size, and if stronger filtering techniques seem to be ineffective, greedy "dives" are often surprisingly successful.

We therefore propose in Section 4 to learn effective stochastic heuristic policies via reinforcement learning. Then, we show how these policies can be used within different tailored approaches: constraint programming with rapid restarts, limited discrepancy search, and multi-start local search approach.

Finally, since industrial benchmarks are easily solved by all new methods introduced in this paper, we generated a synthetic dataset designed to be more challenging. In Section 6 we report experiments on this dataset that further demonstrates the efficiency of the proposed methods.

Problem Definition

The Renault 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 a specific trolley. When a trolley is filled at the production point for that component, an operator must bring it to its consumption point. Symmetrically, when a trolley is emptied it must be brought to the corresponding production point. A production cycle is the time c i taken to produce (resp. consume) component i, that is, to fill (resp. empty) a trolley. 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. There are two pickups and two deliveries at every cycle k of each component i: the pickup pe k i and delivery de k i of the empty trolley from consumption to production 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 a is denoted p a and the travel time between operations a and b is denoted D a,b .

Let P be the set of pickup operations and D the set of delivery operations:

P = m i=1 ni k=1 {pe k i , pf k i } , D = m i=1 ni k=1 {de k i , df k i } .
The problem is to compute the bijection σ (let ρ = σ -1 be its inverse) between the |A| = n first positive integers to the operations A = P ∪ D satisfying the following constraints.

s σ,j = r σ(j) if j = 1 max(r σ(j) , s σ,j-1 + p σ(j-1) + D σ(j-1),σ(j) ) otherwise
Then, the completion time e σ,j = s σ,j + p σ(j) of operation σ(j) must be lower than its due date:

∀j ∈ [1, n], e σ,j ≤ d σ(j) (1)
Precedences. Pickups must precede deliveries.

ρ(pf k i ) < ρ(df k i ) ∧ ρ(pe k i ) < ρ(de k i ) ∀i ∈ [1, m] ∀k ∈ [1, n i ] (2)
Notice that there are only two possible orderings for the four operations of a production cycle. Indeed, since the first delivery (of either the full or the empty trolley) and the second pickup take place at the same location, doing the second pickup before the first delivery is dominated (w.r.t. the train length and the time windows). Hence Equation 3 is valid, though not necessary:

ρ(df k i ) < ρ(pe k i ) ∨ ρ(de k i ) < ρ(pf k i ) ∀i ∈ [1, m] ∀k ∈ [1, n i ] (3)
Train length. The operator may assemble trolleys into a train5 , so a pickup need not be directly followed by its delivery. However, the total length of the train of trolleys must not exceed a length T max . Let t i be the length of the trolley for component i, and let t a k i = t i if a k i ∈ P and t a k i = -t i otherwise, then:

∀j ∈ [1, n], j l=1 t σ(l) ≤ T max (4)
This is a particular case of the single vehicle pickup and delivery problem with capacity and time windows constraints. However, there is no objective function, and instead, feasibility is hard. Moreover, the production-consumption cycles entail a very particular structure: the four operations of each component must take place in the same time windows and this is repeated for every cycle. As a result, one of the best method, Large Neighborhood Search [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF], is severely hampered since it relies on the objective to evaluate the moves and the insertion of relaxed requests is often very constrained by the specific precedence structure.

Baseline Models

We designed two CP models: a variant of a single resource scheduling problem and a variant of a TSP. Then, we describe the current LocalSolver model.

Scheduling-based model. The importance of time constraints in the problem studied, suggests that a CP model based on scheduling would be relevant [START_REF] Beck | Vehicle Routing and Job Shop Scheduling: What's the Difference[END_REF]. The problem is a single resource (the operator) scheduling problem with four types of non overlapping operations (pickup and delivery of empty and full trolleys). For each operation a ∈ A, we define the variable s a ∈ [r a , d a ] as the starting date of operation a. Moreover, for each pair of operations a, b ∈ A, we introduce a Boolean variable x ab standing for their relative ordering. In practice, we need much fewer than n 2 Boolean variables as the time windows and Constraint (2) entails many precedences which can be either ignored or directly posted.

x ab = 1 ⇔ s b ≥ s a + p a + D a,b 0 ⇔ s a ≥ s b + p b + D b,a ∀(a, b) ∈ A (5)
x df k i pe k i ∨ x de k i pf k i ∀i ∈ [1, m] ∀k ∈ [1, n i ] (6)
Constraint (5) chanel the two sets of variables, and constraint (6) encodes Equation (3). Finally, Constraint (4) can be seen as a reservoir resource with limited capacity, which is filled by pickups, and emptied by deliveries. We use the algorithm from [START_REF] Laborie | Algorithms for propagating resource constraints in AI planning and scheduling: Existing approaches and new results[END_REF] to propagate it on starting date variables using the precedence graph implied by Boolean variables and precedences from Constraint [START_REF] Bello | Neural Combinatorial Optimization with Reinforcement Learning[END_REF].

TSP-based model. The second model is an extension of the first one, to wich we add extra variables and constraints from the model for TSP with time windows proposed in [START_REF] Sylvain Ducomman | Alternative Filtering for the Weighted Circuit Constraint: Comparing Lower Bounds for the TSP and Solving TSPTW[END_REF]. We need two fake operations, 0 and n + 1, for the start and the end of the route. For each operation a ∈ A ∪ {0, n + 1}, there is a variable next a that indicates which operation directly follows a. Also, a variable pos a indicates the position of the operation a in the sequence of operations. We need another variable train a ∈ [0, T max ] which represents the length of the train before the operation a. The following equations express the constraints of the problem:

train 0 = 0 ∧ train nexta = train a + t a ∀a ∈ A ∪ {0} (7)
pos pf k i < pos df k i ∧ pos pe k i < pos de k i ∀i ∈ [1, m] ∀k ∈ [1, n i ] ( 8
)
s 0 = 0 ∧ s nexta ≥ s a + p a + D a,nexta ∀a ∈ A ∪ {0} (9)
Constraint [START_REF] Martin | Solving the Torpedo Scheduling Problem[END_REF] encodes the train length at every point of the sequence using the Element constraint and constraint [START_REF] Harvey | Limited Discrepancy Search[END_REF] ensures that pickups precede their deliveries. Constraint ( 9) ensure the accumulation of the time along the tour. Moreover, we use the Circuit constraint to enforce the variables next to form an Hamiltonian circuit. Additional redundant constraints are used to make the channeling between variables and improve filtering in addition to the constraint from the first model:

x ab =⇒ next b = a, pos b > pos a + 1 =⇒ next a = b, pos b > pos a ⇔ x ab , pos a = b∈A∪{0,n+1}
x ab and AllDifferent(pos).

Search strategy. Preliminary experiments revealed that branching on variables in an ordering "consistent" with the sequence of operations was key to solving this problem via CP. In the TSP model, we simply branch on the variables next in ascending order, and choose first the operation with least earliest start time.

In the scheduling model, we compute the set of pairs of operations a, b such that {s a , s b } is Pareto-minimal, draw one pair uniformly at random within this set and assign x ab so that the operation with least release date comes first. In conjunction with constraint propagation, this strategy acts as the "nearest city" heuristic in TSP. Indeed, since the sequence of past operations is known, the earliest start time of an operation depends primarily on the distance from the latest operation in the partial sequence (it also depends on the release date).

LocalSolver. The LocalSolver (LS) model is similar to the TSP model. It is based on a variable of type list seq, a special decision variable type that represent the complete tour: seq j = a means operation a is performed at position j. This variable is channeled with another list variable pos with pos a = j ⇔ seq j = a.

We need two other list variables: train and s to represent respectively the length of the train and the start time of the operation at a given position. 14)

s 1 = 0 ∧ s j = max(r seqj , s j-1 + p seqj-1 + D seqj-1,seqj ) ∀j ∈ [2, n] (10)
train 1 = 0 ∧ train j = train j-1 + t seqj ∀j ∈ [2, n] (11)
s j + p seqj ≤ d seqj ∀j ∈ [1, n] (12)
pos pf k i < pos df k i ∧ pos pe k i < pos de k i ∀i ∈ [1, m] ∀k ∈ [1, n i ] ( 13
)
pos df k i < pos pe k i ∨ pos de k i < pos pf k i ∀i ∈ [1, m] ∀k ∈ [1, n i ] (
Count(seq) = 0 (15)
The last constraint (15) acts like the global constraint All Different and therefore ensures that seq is a permutation. Surprisingly, relaxing the due dates and using the maximal tardiness as objective tends to degrade the performance of LocalSolver, hence we kept the satisfaction version.

Preliminary Experiments

The industrial data we have collected fall into three categories. In the industrial assembly line, each category is associated to one logistic operator who has the charge of a given set of components. In practice, these datasets are relatively underconstrained, with potentially quite large production cycles (time windows) for each component. For each category, denoted by S, L and R, we consider three time horizons: 43 500, 130 500 and 783 000 hundredths of a minutes which corresponds to a shift of an operator (7 hours and 15 minutes), a day of work (made up of three shifts) and a week (6 days) respectively. We then have 9 industrial instances from 400 to more than 10 000 operations. The main differences between those three categories are the number of components, and the synchronicity of the different production cycles. For S instances, there are only 5 components, and their production cycles are almost the same and very short. The other two instances have more than 30 components, with various production cycles (some cycles are short and others are very long).

The CP models were implemented using the constraint solver Choco 4.10 [START_REF] Prud'homme | Choco Documentation[END_REF], and the LocalSolver version was 9.0. All experiments of this section were run on Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz with a timeout of 1 hour.

The results are given in Table 1. The first column (Cl) denotes the category of the instance and second column (H) denotes the temporal horizon. For the two CP models (Scheduling and TSP), two indicators are given, the CPU time (in seconds), and the numbers of fails for solved instances. For LocalSolver, we give the CPU time when the instances are solved before the time limit.

The CP scheduling-based model solves all of the industrial instances without any fail except for instances R and L on the whole week. Notice, however, that when using generic heuristics such as Domain over Weighted Degree [START_REF] Boussemart | Boosting systematic search by weighting constraints[END_REF], only 3 instances could be solved by the same model in less than an hour. Moreover, LocalSolver cannot solve the largest instances, and requires much more CPU time in general. Although industrial instances are clearly very underconstrained, they are not trivial. Moreover, even on underconstrained instances, wrong early choices may lead to infeasible subtrees from which we cannot easily escape. In particular, the number of fails for the largest instance of category R shows that the very deep end of the search tree is likely explored in a brute-force manner.

One key factor in solving these instances is for the variable ordering to follow the ordering of the sequence of operations being built. Indeed, the propagation is much more efficient in this way, and in particular, the earliest start times of future operations can be easily computed and, as mentioned in the description of the heuristic, it reflects the distance from the last operation in the route.

We observe that the scheduling-based model is the fastest. Moreover, the TSP-based model contains too many variables and constraints, and run out of memory for the bigger horizon. Instances without any fails show that the first branch can be very slow to explore as propagation in nodes close to the root can be very costly, in particular with the TSP-based model.

We draw two conclusions from these preliminary experiments: First, the basic CP (or LS) models cannot reliably solve hard instances 6 . It is unlikely that stronger propagation would be the answer, as the TSP model (with stronger filtering) tends to be less effective, and does not scale very well. Second, building the sequence "chronologically" helps significantly, which explains why CP models outperform local search. As a consequence, greedy runs of the CP solver are surprisingly successful. Therefore, we propose to learn efficient heuristic policies and explore methods that can take further advantage of these heuristics.

Reinforcement Learning

The search for a feasible sequence can be seen as a Markov Decision Process (MPD) whose states are partial sequences σ, and actions A(σ) are the operations that can extend σ without breaking precedence (2) nor capacity constraints (4).

In order to apply Reinforcement Learning (RL) to this MDP, it is convenient to relax the due date constraints and replace them by the minimization of the maximum tardiness: max{L(σ, j) | 1 ≤ j ≤ |σ|} where L(σ, j) = e σ,j -d σ(j) is the tardiness of operation σ(j). We also define the maximum tardiness on intervals: L(σ, j, l) = max{L(σ, q) | j ≤ q ≤ l}, and we write L(σ) for L(σ, 1, |σ|). Since in this case operations can finish later than their due dates, it is necessary to make explicit the precedence constraints due to production cycles:

max(ρ(df k-1 i ), ρ(de k-1 i )) < min(ρ(pe k i ), ρ(pf k i )) ∀i ∈ [1, m] ∀k ∈ [2, n i ] ( 16
)
Then we can further restrict the set A(σ) to operations that would not violate Constraints [START_REF] Williams | Simple statistical gradient-following algorithms for connectionist reinforcement learning[END_REF], nor (3). As a result, any state of the MDP reachable from the empty state is feasible for the relaxation. Finally, we can define the penalty R(σ, j) for the j-th decision as the marginal increase of the objective function when the j-th operation is added to σ: R(σ, j) = L(σ, 1, j) -L(σ, 1, j -1)). Now we can apply standard RL to seek for an effective stochastic heuristic policy for selecting the next operation. Moreover, as the MPD defined above is exponentially large, it is common to abstract a state σ with a small descriptor, namely, λ(σ, a) a vector of four criteria for operation a ∈ A in state σ. Let lst(a, σ) be the latest starting time of the task a in order to satisfy constraint (1) with respect to the tasks in σ and constraints (2) and (3). For each operation a ∈ A(σ), we compute λ(σ, a) as follows:

λ 1 (σ, a) = (lst(a, σ) -max(r a , e σ,|σ| -L(σ) + D σ(|σ|),a ))/ max b∈A {d b -r b } (17) λ 2 (σ, a) = max(r a -(e σ,|σ| -L(σ)), D σ(|σ|),a )/ max b,c∈A 2 {D b,c } (18)
λ 3 (σ, a) = 1 -|t a |/T max (19) λ 4 (σ, a) = 1 if a ∈ P 0 otherwise (20)
Criterion ( 17) can be seen as the operation's emergency: the distance to the due date of the task. Criterion (18) is the travel time from the last task in the sequence. For both of these criterion, we use e σ,|σ| -L(σ) instead of e σ,|σ| as the end time of the partial sequence σ to offset the impact of previous choices. Criterion ( 19) is the length of the trolley. Indeed, since all operations must eventually be done, doing the operations which have the highest consumption of the "train" resource earlier leaves more freedom for later. Finally, criterion (20) penalizes pickups as leaving a trolley in the train for too long is wasteful. We want to learn a stochastic policy π θ , governed by a set of parameters θ which gives the probability distribution over the set of actions A(σ) available at state σ. We first define a fitness function f (σ, a) as a simple linear combination of the criteria: f (σ, a) = θ λ(σ, a). Then, we use the softmax function to turn the fitness function into a probability distribution (ignore parameter β for now).

∀a ∈ A(σ) π θ (a | σ) = e (1-f (σ,a))/β b∈A(σ) e (1-f (σ,b))/β (21)

Policy Gradient

As we look for a stochastic policy, the goal is to find the value of θ minimizing the expected maximum value (i.e., maximum tardiness) J(θ) of solutions σ produced by the policy π θ . The basic idea of policy gradient methods is to minimize J(θ) by gradient descent, that is: iteratively update θ by subtracting the gradient. We resort to the REINFORCE learning rule [START_REF] Williams | Simple statistical gradient-following algorithms for connectionist reinforcement learning[END_REF] to get the gradient of J(θ):

∇ θ J(θ) = E σ∼π θ (σ) [L(σ)∇ θ log π θ (σ)] (22)
Which can then be approximated via Monte-Carlo sampling over S samples:

∇ θ J(θ) ≈ 1 S S k=1 L(σ k )∇ θ log π θ (σ k ) (23)
We can decompose Equation (23) in order to give a specific penalty to each decision (for every position j) of the generated solutions:

∇ θ J(θ) ≈ 1 S S k=1 n j=1 G(σ k , j)∇ θ log π θ (σ k (j) | σ k (1, j -1)) (24)
The penalty function takes into account that a decision can only affect the future marginal increase of tardiness: we replace the overall tardiness L(σ) by:

G(σ, j) = R(σ, j) if j = n R(σ, j) + γ * G(σ, j + 1) otherwise
The value of γ controls how much decisions impact future tardiness. For γ = 0, we only take into account the immediate penalty. Conversely γ = 1 means we consider the sum of the penalties from position j.

Learning algorithm

We learn a value of θ for the synthetic dataset (Section 6) with the REINFORCE rule. Given a set of instances I, we learn by batch of size S = q|I|, in other words, we generate q solutions for each instance. The value of θ is initialized at random, then we apply the following three steps until convergence or timeout:

1. Generate S solutions following π θ . 2. Compute ∇ θ J(θ) according to Equation (24) 3. Update the value of θ as follows: θ ← θ -α∇ θ J(θ)

We found out that using a classic softmax function did not discriminate enough, and hence acts as random policy. To circumvent this, we use the parameter β in Equation (21) to control the trade off between the quality and the diversity of generated solutions. For a low value of β the policy always chooses the best candidate, whereas a large value yields a more "balanced" policy. It turns out that β = 0.1 was a good value for learning in our case.

Then, we have evaluated the impact of γ to compute the gradient in Equation (24). Recall that γ controls how much importance we give to the j-th decision in the total penalty: the overall increase of the tardiness from j to n for γ = 1 or only the instant increase for γ = 0. Moreover, we also tried to give the penalty L(σ) uniformly for every decision j of the policy, instead of giving individual penalties G(σ, j). We denote this penalty strategy "Uniform".

For each variant, we plot in Figure (1a) the average performance L(σ) of the policy after each iteration (notice the log-scale both for X and Y ). Here we learn on all generated instances of a day horizon (40 instances), for 2000 iterations.

High values for γ are always better. For low value of γ, the (local) optimum is higher, and in some cases the method may not even converge, e.g. for γ = 0. Using uniformly the overall tardiness gives similar results to γ = 1, however it is less stable, and in Fig. (1a) we observe that it diverges after 1000 iterations. Datapoints are colored with the vector θ interpreted as a RGB value (|θ| = 4 but it can be characterized by three values after we normalize so that 4 i=1 θ i = 1) 7 . We can see for instance that γ = 0.9 finds a value of θ that is significantly different from all other methods ( 0, 0.63, 0.29, 0.07 ). Interestingly, "Uniform" and γ = 1 not only converge to the same average tardiness, but to similar θ's (respectively 0.30, 0.49, 0.16, 0.04 and 0.25, 0.56, 0.15, 0.04 ). However, γ = 1 finds values of θ that seem closer to the target ("greener") earlier, although it is not really apparent from the value of L(σ). These values of θ indicate that a good heuristic is mainly a compromise between the emergency and the travel time 8 . Moreover, the travel time tends to be more important on larger horizon, because it a longer term impact: all subsequent operations are affected.

Finally, we report in Figure (1b) the average tardiness for each instances following π θ . The gain of learning specifically for a given (class of) instance(s) is at best marginal. This is not so surprising as we abstract states with a very simple model using a few criteria. However, it means that the value of θ learnt on the full dataset is relevant to most instances.

Using the Heuristic Policy

We have implemented two types of approaches taking advantage of the stochastic policy learnt via RL: integrating it within CP as a randomized branching heuristic, and using it to generate sequences to be locally optimized via steepest descent in a multi-start local search. Here we describe the necessary modifications of the CP model, and we propose an efficient local search neighborhood.

Constraint Programming

In order to use the stochastic policy described in Section 4 in a CP solver, we need to slightly change the scheduling model. We introduce a new set of variables seq j , one for each position j, standing for the operation at that position, with the following channeling constraint:

x seqj seq l = 1 ∀j < l ∈ [1, n]
We branch only on the variables seq in lexicographic order. Therefore, the propagation for this constraint is straightforward: when branching on seq j = a, as all variables seq l ∀l < j are already instantiated, we set x ab = 1, ∀b ∈ A\{seq l | l < j}. Conversely, after assigning seq j we can remove a from the domain of seq j+1 if there exists b ∈ A\{seq l | l ≤ j} such that the domain of the variable x ab is reduced to {0}. Moreover, we can easily enforce Forward Consistency on {seq 1 , . . . , seq j } with respect to precedence and train size constraints ( 2) and ( 4) as well as Constraint (3), i.e., when {seq 1 , . . . , seq j-1 } are all instantiated, we can remove all values of seq j that cannot extend the current subsequence. Therefore, we do not need the Reservoir resources propagator anymore.

We propose two strategies based on this CP model using the learned policy.

1. Softmax policy and rapid restart. In this method we choose randomly the next operation according to the softmax policy (Eq. 21). In order to explore quickly different part of the search tree, we rely on a rapid restart strategy, following a Luby [START_REF] Luby | Optimal Speedup of Las Vegas Algorithms[END_REF] sequence with a factor 15. 2. Limited Discrepancy Search. As the key to solve those instances is to follow good heuristics, and to deviate as little as possible from them, limited discrepancy search (LDS) [START_REF] Harvey | Limited Discrepancy Search[END_REF] fits well with this approach. We run the LDS implementation of Choco, which is an iterative version: the discrepancy starts from 0, to a maximum discrepancy parameter incrementally. For this approach we use the deterministic version of the policy π(σ) = arg min a f (σ, a).

Local Search

The solutions found by the heuristic policy can often be improved by local search moves. Therefore, we also tried a multi-start local search whereby we generate sequences with the heuristic policy, and then improve them via steepest descent. Sequences are generated using the same model used for RL (i.e., with relaxed due dates). Therefore, generated sequences respect all constraints, except (1) and we consider a neighborhood that preserves all other constraints as well. Then we apply the local move that decrease the most the maximum tardiness L(σ) until no such move can be found. We use two types of moves and the time complexity of an iteration (i.e., computing, and commiting to, the best move) is in O(nm).

We recall that L(σ, j, l) = max{L(σ, q) | j ≤ q ≤ l} is the maximum tardiness among all operations between positions j and l in σ.

Swap moves. The first type of moves consists in swapping the values of σ(j) and σ(l). First, we need to make sure that the ordering of the operations within a given component remains valid, i.e., satifies constraints (2) and ( 16): a pickup (resp. delivery) operation must stay between its preceding and following deliveries (resp. pickups) for the same component. For every operation a, a valid range between the position of its predecessor pr(a) and of its successor su(a) can be computed in constant time. Then, for all j ∈ [1, n] we shall consider only the swaps between j and l for l ∈ [j + 1, su(σ(j))] and such that pr(σ(l)) ≤ j.

The second condition for the move to be valid is that the swap does not violate constraint (4), i.e., the maximum length of the train. Let τ j = j-1 l=1 t σ(l) be the length of the train before the j-th operation. After the swap we have τ j+1 = τ j + t σ(l) which must be less than T max . At all other ranks until l, the difference will be t σ(l) -t σ(j) , we only need to check the constraint for the maximum train length, that is: max{τ q | j ≤ q ≤ l} + t σ(l) -t σ(j) ≤ T max . This can be done in constant (amortized) time for all the swaps of operations a by computing the maximum train length incrementally for each l ∈ [j +1, su(σ(j))].

Then, we need to forecast the maximum tardiness of the sequence σ where the operations at positions j and l are swapped, i.e., compute the marginal cost of the swap. The tardiness of operations before position j do not change. However, we need to compute the new tardiness L(σ , j) and L(σ , l) at positions j and l, respectively. Moreover, we need to compute L(σ , j + 1, l -1) the new maximum tardiness for operations strictly between j and l and L(σ , l + 1, n) the new maximum tardiness for operations strictly after l.

The new end time e σ ,j of operation σ (j) = σ(l) and hence the tardiness at position j is L(σ , j) = e σ ,j -d σ (j) can be computed in O(1) as follows: e σ ,j = p σ (j) + max(r σ (j) , (e σ ,j-1 + D σ (j-1),σ (j) ))

Next, operations σ(j + 1), . . . , σ(l -1) remain in the same order and σ(j + 1) is shifted by a value ∆ = e σ ,j + D σ (j),σ (j+1) -e σ,j -D σ(j),σ,j+1 . However, subsequent operations may not all be equally time-shifted. Indeed, when ∆ < 0 there may exist an operation whose release date prevents a shift of ∆.

Let g j = r σ(j) -s σ(j) be the maximum left shift (negative shift of highest absolute value) for the j-th operation, and let g j,l = max{g q | j ≤ q ≤ l}. Proposition 1. If the sequence does not change between positions j and l, a time-shift ∆ < 0 at position j yields a time-shift max(∆, g j,l ) at position l.

Let L ∆ (σ, j, l) be the maximum tardiness on the interval [j, l] of sequence σ time-shifted by ∆ from position j. We can define L -∞ (σ, j, l) the maximum tardiness on the interval [j, l] for an infinite negative time-shift:

L -∞ (σ, j, l) = max{L(σ, q, l) + g j,q | j ≤ q ≤ l} (25) Proposition 2. If ∆ < 0 then L ∆ (σ, j, l) = max(∆ + L(σ, j, l), L -∞ (σ, j, l)).
Conversely, when ∆ > 0 some of the time-shift may be "absorbed" by the waiting time before an operation. However, there is little point in moving operations coming before a position j with a non-negative waiting time (i.e., where s σ(j) = r σ(j) ), as this operation and all subsequent operations would not profit from the reduction in travel time. Therefore we consider only swaps whose earliest position j is such that ∀q > j, g q < 0. As a results, if there is a positive time-shift ∆ at q > j, we know that L ∆ (σ, q, l) = ∆ + L(σ, q, l). Moreover, the values of L(σ, j, l), g j,l and L -∞ (σ, j, l) can be computed incrementally as: L(σ, j, l + 1) = max(L(σ, j, l), L(σ, l + 1)) g j,l+1 = max(g j,l , g l+1 ) L -∞ (σ, j, l + 1) = max(L -∞ (σ, j, l), g j,l+1 + L(σ, l + 1)) Therefore, when j < l-1, we can compute the new tardiness L(σ , j + 1, l -1) = L ∆ (σ, j + 1, l -1) of the operations in the interval [j +1, l-1] in constant (amortized) time since the query of Proposition 2 can be checked in O [START_REF] Beck | Vehicle Routing and Job Shop Scheduling: What's the Difference[END_REF].

The new tardiness L(σ , l) at position l is computed in a similar way as for L(σ , j) since we know the new start time of σ (l -1) from previous steps.

Finally, in order to compute the new maximum tardiness L(σ , l + 1, n) over subsequent operations, we precompute L(σ, j, n), g j,n and L -∞ (σ, j, n) for every position 1 ≤ j ≤ n once after each move in O(n). Then L(σ , l + 1, n) can be obtained in O(1) for every potential move from Proposition 2.

Therefore, we can check the validity and forecast the marginal cost of a swap in constant amortized time and perform the swap in linear time. The time complexity for an iteration is thus in O(nm) since, for a given component i, the sum of the sizes of the valid ranges for all pickups and deliveries of this component is in O(n). Indeed, let a 1 i , . . . , a 4ni i be the operations component i ordered as in σ. Then su(a k i ) = ρ(a k+1 i ) and

4ni k=1 su(a k i ) -ρ(a k i ) = ρ(a 4ni i ) -ρ(a 1 i ) ∈ Θ(n)).
Toggle moves. As observed in Section 3, there are only two dominant orderings for the four operations of the k-th production cycle of component i. The second type of moves consists in changing from one to the other of these two orderings, by swapping the values of σ(pf k i ) and σ(pe k i ) and the values of σ(df k i ) and σ(de k i ). This change leaves the size of the train constant, hence all these moves are valid.

Let j 1 = pf k i , j 2 = df k i , j 3 = pe k i , j 4 = de k i be the positions of the four operations of component i and cycle k in the current solution, and suppose, wlog, that j 1 < j 2 < j 3 < j 4 . Let σ denote the sequence obtained by applying a toggle move on component i and cycle k in σ. In order to forecast the marginal cost of the move, we need to compute the new tardiness at the positions of the four operations involved L(σ , j 1 ), L(σ , j 2 ), L(σ , j 3 ) and L(σ , j 4 ). Moreover, we need the new maximum tardiness on four time-shifted intervals:

L(σ , j 1 + 1, j 2 -1) = L ∆1 (σ, j 1 + 1, j 2 -1), L(σ , j 2 + 1, j 3 -1) = L ∆2 (σ, j 2 + 1, j 3 -1), L(σ , j 3 + 1, j 4 -1) = L ∆3 (σ, j 3 + 1, j 4 -1) and L(σ , j 4 + 1, n) = L ∆4 (σ, j 4 + 1, n).
Computing the marginal costs can be done via the same formulas as for swaps: we can first compute the new end time for the operation at position j 1 , then from it compute the value of ∆ 1 that we can use to compute L ∆1 (σ, j 1 + 1, j 2 -1) in O(j 2 -j 1 -1) time, and so forth. The difference, however, is that there is fewer possible moves (n/4), although the computation of the marginal cost cannot be amortized. The resulting time complexity is the same: O(nm).

The learning rate depends on the size of the instances: α = 2 -12 /n where n is the average number of task in the batch. The rationale is that the magnitude of the gradient depends on the tardiness L(σ) which tends to grow with the number of operations. Therefore, we use the learning rate α to offset this growth, which is key to have a stable convergence. For the two methods using the stochastic policy, we made the first run deterministic i.e. the policy becomes stochastic only after the first restart for the CP based one, and after the first iteration for the multi-start local search.

The results are presented in Table 2. We report the number of solved instances (among 10 instances for every time horizon) averaged over 10 randomized runs for each CP model in the column "#S". The synthetic dataset is more constrained than the industrial dataset and the two basic CP models fail to solve most of the instances ("-" in the table indicates a time out). However, the relative performance remains unchanged w.r.t. Table 1: the scheduling-based models shows better performance in terms of number of solved instances and CPU time while scaling better in memory. All three of the RL-based methods significantly outperform previous approaches. The results in Table 2 indicate that the rapid restarts approach dominates the others. However, it may not be as clear-cut as that: for other settings of the hyperparameters (α, β and γ) the relative efficiencies fluctuate and other methods can dominate. Moreover, one advantage of the multi-start local search method is that since due dates are relaxed, imperfect solutions can be produced, even for infeasible instances. We report the average maximum tardiness in column L max .

This global θ also works well with the industrial dataset. All instances are easily solved by all three methods, except "R" for the week horizon, which is only solved by the rapid restart approach. We learnt a dedicated policy for the industrial dataset with the same settings. It turns out that every instance was solved by the deterministic policy using the new value for θ, except the instance "L" for the week horizon. However, it is easily solved by all three methods.

Conclusion

In this paper we have applied reinforcement learning to design simple yet efficient stochastic decision policies for an industrial problem: planning the production process at Renault. Moreover, we have shown how to leverage these heuristic policies within constraint programming and within local search.

The resulting approaches significantly improve over the current local search method used at Renault. However, many instances on synthetic data remain unsolved. We plan on using richer machine learning models, such as neural networks, to represent states. Moreover, we would like to embed this heuristic in a Monte-Carlo Tree Search as it would fit well with our current approach since it relies on many rollouts. Finally, we would like to tackle the more general problem of assigning components to operators and then planning individual routes.

Fig. 1 :

 1 Fig. 1: Behavior of reinforcement learning

  Time windows. As production never stops, all four operations of the k-th cycle of component i must happen during the time window [(k -1)c i , kc i ]. Let r a k

i

(resp. d a k i ) be the release date (k -1)c i (resp. due date kc i ) of operation a k i of the k-th cycle of component i. The start time of operation σ(j) is:

Table 1 :

 1 Comparison on the industrial instances

	Cl	H	Scheduling	TSP		LocalSolver
			cpu	#fail	cpu	#fail	cpu
	S	Shift	0.275	0	3.999	0	26
	S	Day	0.840	0	48.855	0	566
	S	Week	17.747	0	memory out	timeout
	L	Shift	7.129	0	36.033	8	121
	L	Day	23.468	0	344.338	8	3539
	L	Week	318.008	6	memory out	timeout
	R	Shift	8.397	0	41.615	14	1065
	R	Day	44.215	0	417.116	15	timeout
	R	Week	timeout	773738	memory out	timeout

https://www.localsolver.com/home.html

Trolleys are designed so that they can be extracted out of the train in any order

There might be too few data points to make that claim on industrial instances, but it is clearly confirmed on the synthetic benchmark (see Table2).

To highlight the differences we also normalize the RGB values and omit γ = 0

Although the importance of a criterion also depends on the distribution of the values of λ after normalization, we are confident that the first two criteria are more important than the other two.

Avalaible at https://gitlab.laas.fr/vantuori/trolley-pb.

Experimental Results

We generated synthetic instances 9 in order to better assess the approaches. Due to the time windows constraints, it is difficult to generate certifiably feasible instances. Their feasibility has been checked on the shortest possible horizon, i.e., the duration of the longest production cycle of any component, which is about 10 000 time units depending on the instances. There are four categories of instances parameterized by the number of components (15 in category A, 20 in B, 25 in C and 30 in D). In the real dataset, several components have similar production cycles. We replicates this feature: synthetic instances have from 2 or 3 distinct production cycles in category A, to up to 7 in category D. The latter are therefore harder because there are more asynchronous productions cycles. We generated 10 random instances for each category and consider the same three horizons (shift, day and week) for each, as industrial instances.

Experiments for this section were run on a cluster made up of Xeon E5-2695 v3 @ 2.30GHz and Xeon E5-2695 v4 @ 2.10GHz. For the basic CP models, we add randomization and a restart strategy following a Luby sequence, and we ran each of the 120 instances 10 times. We could not carry on experiments on synthetic data with LocalSolver because we were not granted a license. However, from the few tests we could do, we expect LocalSolver to behave similarly as on industrial benchmarks. For each of the methods using the learnt heuristics, we normalize θ so that 4 i=1 θ i = 1 and set β to 1/150. We learn the policy by batches of size 240 formed by 6 runs on each of the 40 day-long instances, during 2000 iterations.

Table 2: Comparison of the methods on generated instances