Anh Quoc

Bastien Chevalier Pierre Tran

Breul

Spectral analysis of the response of coarse granular material to dynamic penetration test modelled with DEM

Keywords: dynamic cone penetration, discrete element method, granular material, particle size distribution, Discrete Fourier Transform

Dynamic penetration tests are often used to determine the strength properties of surface soils.

The paper presents a study on the use of spectral analysis on dynamic cone penetration tests results, modelled with Discrete Element Method. This method is applied to assess the effect of the variation of the grain size distribution of the soil on test results. A two-dimensional discrete model is used to reproduce cone penetration tests in dynamic conditions: the tip of the penetrometer is driven in the material by successive impacts of a hammer on the penetrometer.

For each impact of the hammer, a curve of the load applied by the tip on the soil is obtained versus the penetration distance of the tip. The curves of the load vs. penetration traditionally used to calculate the tip resistance of the soil are analyzed with Discrete Fourier transform in order to investigate curve's shape. The effect of the variation of the grain size distribution of the soil on these curves is investigated, i.e. average particle diameter and span of particle size distribution. It was found out that the grain size distribution influences tip resistance but also the shape and oscillation modes of the curve of the stress-penetration curve. Based on these indicators, the exploitation of the load-displacement curve obtained with dynamic penetration tests could be enlarged to determine other properties of the soils.

Introduction

In the practice of geotechnical engineering, the use of in situ testing is widespread. Among existing testing techniques, lightweight dynamic penetration testing device such as Panda penetrometer is used to characterize mechanical properties of surface soils [START_REF] Chaigneau | Caractérisation des milieux granulaires de surface à l'aide d'un pénétromètre[END_REF][START_REF] Breul | Penetration Test Modelling in a Coarse Granular Medium[END_REF]. Based on recent technological improvements, this device is able to record the curve of the tip stress or tip force versus the tip penetration distance for each impact of the hammer on the penetrometer [START_REF] Benz Navarrete | Mesures dynamiques lors du battage du pénétromètre Panda 2[END_REF][START_REF] Escobar | Dynamic Characterization of the Supporting Layers in Railway Tracks Using the Dynamic Penetrometer Panda 3®[END_REF]. This curve, referred to as load-penetration curve, provides information on dynamic tip resistance but also on additional mechanical parameters involved during the driving of the tip [START_REF] Benz Navarrete | Mesures dynamiques lors du battage du pénétromètre Panda 2[END_REF][START_REF] Benz | Dynamic measurements of the penetration test -Determination of the tip's dynamic load penetration curve[END_REF][START_REF] Valencia | Mise au point et exploitation d'une nouvelle technique pour la reconnaisance des sols: le PANDA 3[END_REF]. Cone penetration test is a blind testing technique, because there is no sampling of the different layers of soils that are cut across. It would be very interesting to be able to get information about the nature of the soils, from the actual measurements recorded by the device.

The present paper proposes to study the effect of the particle size distribution on the load vs. penetration distance measured at the tip. Numerical modelling using Discrete Element Method was chosen so that we have a strict control on the parameters of the particle size distribution (PSD).

The numerical model of penetration tests using Discrete Element Method (DEM) used to reproduce the penetration tests in dynamic conditions is the same model as the one presented in [START_REF] Tran | Discrete modeling of penetration tests in constant velocity and impact conditions[END_REF]. Many authors have modeled cone penetration tests before with DEM in 2D [START_REF] Huang | An analytical study of cone penetration test in granular material[END_REF][START_REF] Huang | Advanced calibration chambers for cone penetration testing in cohesionless soils[END_REF][START_REF] Calvetti | Micro-macro relationships from DEM simulated element and in-situ tests[END_REF][START_REF] Jiang | Discrete element modeling of deep penetration in granular soils[END_REF][START_REF] Jiang | Future continuum models for granular materials in penetration analyses[END_REF][START_REF] Jiang | Investigating mechanism of inclined CPT in granular ground using DEM[END_REF][START_REF] Tran | A numerical study of the penetration test at constant rod velocity[END_REF][START_REF] Janda | DEM modeling of cone penetration and unconfined compression in cohesive solids[END_REF] and in 3D [START_REF] Benz Navarrete | Mesures dynamiques lors du battage du pénétromètre Panda 2[END_REF][START_REF] Arroyo | Cone penetration tests in a virtual calibration chamber[END_REF][START_REF] Butlanska | Mapping deformation during CPT in a virtual calibration chamber[END_REF][START_REF] Mcdowell | A particle refinement method for simulating DEM of cone penetration testing in granular materials[END_REF][START_REF] Quezada | Penetration test in coarse granular material using Contact Dynamics Method[END_REF][START_REF] Ciantia | DEM modelling of cone penetration tests in a double-porosity crushable granular material[END_REF]. However, most of them modelled the penetration in static conditions, i.e. with a constant velocity of the tip, and very few focused on penetration tests in dynamic conditions, i.e. with impacts [START_REF] Valencia | Mise au point et exploitation d'une nouvelle technique pour la reconnaisance des sols: le PANDA 3[END_REF][START_REF] Tran | Discrete modeling of penetration tests in constant velocity and impact conditions[END_REF][START_REF] Quezada | Penetration test in coarse granular material using Contact Dynamics Method[END_REF].

At the macroscopic scale, the tip resistance quantifies the mechanical response of the granular material to the driving of the tip (Fig. 1). At the scale of the contacts or of the particles, the shapes of the load-penetration curves are analyzed in terms of both frequency and amplitude of signal oscillation by using Discrete Fourier Transform (DFT). Figure 2 shows an example of the load penetration curve. The objective of the study is to evaluate the effect of the PSD on the load-penetration curves, based on DFT analysis. Two parameters were tested: the average particle diameter, D50 and the span of the PSD given as the ratio of maximal particle diameter Dmax to minimal diameter Dmin.

At first, we will present the numerical model and parameters used to reproduce the dynamic penetration test as well as the method used to analyze the load-penetration curve with DFT.

Finally, the influence of the particle size distribution is discussed.

Numerical model

Discrete Element Method in two dimensions is used with Itasca software PFC 2D [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. Dense assemblies of disks are generated without gravity and without friction in a container of 0.60 m in width and of 0.45 m in height (Fig. 3). The lateral walls of the container are fixed. A study of boundary conditions with different sample sizes on this numerical model was conducted by Tran et al. [START_REF] Tran | A numerical study of the penetration test at constant rod velocity[END_REF] and revealed that for a container width larger than 0.60 m, there is no effect of the lateral walls on the tip resistance anymore. A linear contact model is used along with a Coulomb friction criterion after the generation process. The normal contact stiffness of 1.25 x 10 8 N/m is chosen in order to assess the assumption of rigid particles during penetration tests [START_REF] Combe | Mécanique des matériaux granulaires et origines microscopiques de la déformation[END_REF][START_REF] Roux | Discrete numerical simulation and the mechanical behavior of granular materials[END_REF][START_REF] Tran | Modélisation numérique du comportement des milieux granulaires à partir de signaux pénétrométriques: approche micromécanique par la méthode des éléments discrets[END_REF]. The tangential contact stiffness has been set to 0.75 of normal contact stiffness.

The samples are generated without friction and without gravity in order to reach a random close packing volume fraction, corresponding to the minimal void ratio of the considered PSD. After sample generation, gravity is applied to the system as well as a confining vertical stress of 40 kPa. The confining stress is applied on the top face of the container to simulate an overlaying layer of material and to prevent the effects of free surface to be observed [START_REF] Chaigneau | Caractérisation des milieux granulaires de surface à l'aide d'un pénétromètre[END_REF].

The penetration tests were performed with a frictionless rod of width Drod =14 mm linked to a tip of width Dtip =16 mm at its bottom edge which has a friction coefficient µtip of 0.3 [START_REF] Quezada | Penetration test in coarse granular material using Contact Dynamics Method[END_REF]. First, the rod is driven with constant velocity until a depth of 0.15 m is reached inside the granular material. Then, the rod is released and stabilized under its own weight. Finally, series of five successive impacts are produced in the sample with a hammer represented by an additional disk hitting the top of the rod (Fig. 3). The mass of the impacting cylinder is equal to the rod mass.

The impacting velocity at the impact is equal to 1,25 m.s -1 in order to obtain an average penetration distance that is representative of experimental tests. The description of the model and the effect of the impact velocity are addressed by Tran et al. in [START_REF] Tran | Discrete modeling of penetration tests in constant velocity and impact conditions[END_REF]. Table 1 summarizes the main parameters of the model. Table 2 summarizes the parameters of the five granular materials studied, with two types of PSD. The first distribution type (I) keeps constant the ratio between maximal and minimal particle diameters Dmax/Dmin = 2, and particle number varies from 10 000 to 160 000 making the average diameter D50 varying (Fig. 4 a). The second distribution type (II) keeps constant maximal diameter Dmax = 7.02 mm and the ratio between maximal and minimal particle diameters Dmax/Dmin varies from 2 to 10 (Fig. 4 b). All the samples were tested with their maximal volume fraction (φmax) and with a particle friction coefficient equal to 1.0.

For each type of granular material, the testing process was repeated on three different samples corresponding to the same sample conditions but with different initial particle arrangements.

The following section explains how the load-penetration curves obtained from dynamic cone penetration tests were analyzed in the frequency domain by using DFT. At the end of the driving process, the tip reaches a final position corresponding to a residual penetration distance sres. We can notice that there is a difference between the maximal penetration distance smax and the final residual penetration distance sres (Fig. 5). However, the mechanical work of the tip force between these two positions is negligible due to the value of the tip force between these two positions. Consequently, the dynamic tip resistance Rd of the granular material for one impact was calculated as the average tip force Fd for penetration distance between 0 and smax:

Investigation of load -penetration curve

R = 1 t ds t (1)
with: t the time; tsmax the time when penetration distance is maximal and equal to smax. Then, <Rd> is the average value of dynamic tip resistances obtained for the five impacts and for the three samples.

Frequency analysis using Discrete Fourier Transform

Each load-penetration curve presents variations -i.e. oscillations, peaks -observed between the time of the impact and the time when the system stabilizes again (Fig. 5). As the properties of the granular material change, the shape and size of these variations change too. Consequently, the load-penetration curve can provide not only the tip resistance of granular media but also information on the granular material properties.

Fourier transformation provides a powerful way for study discrete data acquisition in the frequency domain. It allows to decompose a signal into the frequencies that make it up. Discrete Fourier Transform (DFT) changes the N "temporal" points yn(x) in to N "frequency" points Yk and inverse DFT the other way around by using the following equations:

Y = !" #/% % & y ( = ) # !" #/% % & # (2)
with: yn the original function, Yk the transformed function and N-points DFT.

The DFT bins (fk) represent frequencies in the discrete Fourier transform that are spaced at intervals of ∆fk = Fs/N, where Fs is the sample rate equal to 1/∆x, with ∆x the recording interval.

All of load-penetration curve have been analyzed in frequency domain with Fs = 10 5 [m -1 ] and frequency resolution ∆fk = Fs/N = 0.1. The equation of inverse DFT can be computed by the following formula:

y ( = * # +, 2./ # 0 + 2 % & # (3)
with: fk the k th frequency, Ak the associated amplitude and ϕ the phase of Yk with

= tan & Img )8 /Real )8 (4)
In order to obtain a stabilized and accurate amplitude spectrum, it is necessary to ensure that the number of DFT points is sufficient. Thus, the original signal of tip force versus penetration distance is padded with trailing zeros to increases its length, before computing the DFT. Figure 6 presents one example of the amplitude spectrum when the N-points of DFT is equal to 1, 4 and 16 times respectively the number of data of the signal detected between s = 0 and s = smax for one impact in the material A. We found that the amplitude spectrum becomes stable when N-points DFT increases.

Figure 7 presents two examples of amplitude spectrum related to the materials A and C corresponding respectively to maximal and minimal average particle diameters. The first part of the spectrum i.e. fk < 100 [m -1 ] corresponds to a transition zone where amplitude spectrum decreases rapidly. Thus, we introduce one parameter called (ftrans) which is the frequency corresponding to the first local minimum value of the amplitude spectrum (Fig. 6).

In order to analyze the oscillations of the load-penetration curves, we studied the range of frequencies that provides the most significant information to rebuild the signal with inverse DFT and to filter the signal noise. For that purpose, the relative error RE between the reconstructed signal of the k first components (yk) based on Eq. 3 and the original signal (y) descripted by the following equation has been computed.

;< = ∑ ∥ " -# ∶" " ∥ % & × 0 " -0 " & ∑ " × 0 " -0 " & % " (5)
Figure 8 presents the relative error (RE) for one impact of the material A. The range of frequencies [0, flimit] that gives an RE smaller than 10% was chosen for rebuilding the signal and for the load-penetration curves analysis. Note that, the more RE decreases, the more the reconstructed signal is accurate in comparison to the raw signal. Figure 9 presents one example of signal reconstruction for the material A. We found that the signal reconstructed with the frequency range [0, ftrans] provides the general trend (baseline) of the load-penetration curve. The signal reconstructed with the frequency range [0, flimit] provides a reliable reconstruction of the raw signal without noise. Thus, we define the band-pass of frequency (BF range) from ftrans to flimit, for obtaining the full oscillation information of loadpenetration curve. The reconstructed signal with [ftrans, flimit] frequency range presents the major oscillations of load-penetration curve.

Effect of PSD variation

The effect of PSD on the load-penetration curve in terms of tip resistance <Rd> and in frequency domain by using the frequencies detected in band-pass range [ftrans, flimit] is presented in this section. The PSD types I and II are studied in this section (Tab.2).

Figure 10 presents the load-penetration curves of 5 impacts for materials A, B and C. These materials have the same ratio Dmax/Dmin = 2 and D50 of respectively 7.02 mm, 3.51 mm and 1.76 mm. Figure 11 presents the load-penetration curves of 5 successive impacts for the materials A, D and E. These materials have a ratio Dmax/Dmin of respectively 2, 5 and 10, and the same maximal diameter Dmax of 7.02 mm. We found that the final penetration distance and the signals oscillations amplitude increases as D50 decreases and Dmax/Dmin ratio increases.

Concerning the macroscopic response, Figure 12 shows <Rd> as function of D50 for the 5 different psd A to E. In general, <Rd> increases when D50 increases and the standard deviation decreases for smaller particles sizes. Thus, the response is more variable for coarser material.

In addition, as Dmax/Dmin increases, we found that <Rd> sightly decreases and that the standard deviation decreases. Thus, it seems that <Rd> decreases faster as function of D50 than for the effect of spreading. As the quantity of smaller particles increases, the response given by the load-penetration curve is more repeatable.

The oscillations of the load-penetration curve vary as function of the combination of their frequencies fk and the associated amplitudes Ak. We can observe on the load-penetration curves of Fig. 10 and 11 than the amplitudes of oscillations decrease and the frequency increases as the material becomes finer. It is then more difficult to detect the frequency features for finer materials i.e C or E.

In order to quantify the oscillations of the load-penetration curve for each material and in order to compare them to each other, we defined the coefficient Af as the average value of the product of the frequency fk and the associated amplitude Ak detected in the band-pass range [ftrans, flimit].

*/ = ∑ * # × / # # & B (5)
with n the number of frequencies detected in band-pass range [ftrans, flimit].

Figure 13 presents the coefficient Af averaged on 15 impacts for each material A to E as function of D50 and of PSD spreading. We found that <Af> and its standard deviation decrease as the average diameter D50 decreases. The coefficient <Af> was also found to decrease as the spreading of PSD increases. It means that the oscillations of load-penetration curve, quantified by a coefficient <Af>, are dependent on the PSD of the material tested.

Conclusions

In this paper, we presented a method to analyze the dynamic penetration test results, based on the Discrete Fourier Transform. This method was used on results of dynamic penetration test of Panda modelled in 2D with DEM. The method was applied to the evaluation of the variation of particle size distribution of a coarse granular material, on and more particularly the influence of average particle diameter and particle size distribution span. Five different materials were simulated under dynamic penetration test.

The discrete Fourier transform was applied on the load vs. penetration curves so that:

• the oscillations of the penetration curve can be reconstructed • a band-pass range that characterizes the material in terms of particle size distribution can be identified • influence of mean diameter D50 on the both global response and signal oscillations can be quantified with a scalar called here amplitude coefficient <Af> calculated from the frequencies and associated amplitudes in the characteristic band-pass range. This coefficient captures only the nature of the oscillations and is influenced by particle size distribution of the material.

Concerning the macroscopic response, we observed that • average tip resistance decreases when the content of smaller particles increases: when D50 decreases and when the spreading of the particle size distribution increases • the DFT applied on signal obtained with finer granular materials shows a distribution of frequencies and associated amplitude that is more flat, in the characteristic band-pass range, making it more difficult to detect the characteristic frequency and amplitude of the oscillations for such materials.

The method presented in this paper based on spectral analysis will be useful to extract more information from the penetrometer test results, in terms of physical properties of the material, and beyond the classical tip resistance. Further investigations using three-dimensional modeling and experimental validation will have to be performed in order to confirm and extend the 228 benefits of this method.

3. 1 .

 1 Macroscopic exploitation for the tip resistanceTip force Fd is defined as the vertical component of the force applied by the granular material on the tip as it drives in the granular with a penetration distance s. Figure5shows examples of load-penetration curves Fd = f(s) obtained for 3 impacts with the numerical model on the material B. The response obtained with the model is similar to the one classically obtained experimentally[START_REF] Escobar | Dynamic Characterization of the Supporting Layers in Railway Tracks Using the Dynamic Penetrometer Panda 3®[END_REF].

Figure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 8 .Figure 9 .Figure 10 .Figure 11 .Figure 12 .Figure 13 .Figure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 8 .Figure 9 .Figure 10 .Figure 11 .Figure 12 .Figure 13 .

 1234567891011121312345678910111213 Figure 1. Experimental result of dynamic penetration test [12] Figure 2. Example of experimental load-penetration curve obtained in a dynamic penetration test for one impact [3] Figure 3. View of the numerical model of the tip and of the container for the penetration test in dynamic conditions (Tran et al., 2015) Figure 4. Particle size distribution (a) for different average particle diameters (materials A, B and C) and (b) for different particle size distribution shapes (materials A, D and E) Figure 5. Example of load-penetration curves obtained in a dynamic cone penetration test for 3 impacts in material A Figure 6. Example of the amplitude spectrum of one impact of the material A with different numbers of N-points DFT Figure 7. Example of the amplitude spectrum of five impacts for two materials A and C Figure 8. Relative error for one impact of the material A Figure 9. (Left) Comparison of the raw load-penetration curve for the material A with signal reconstructed with frequency ranges [0, ftrans] and [0, flimited]. (Right) Signal reconstructed with frequency range [ftrans, flimited] showing the major oscillations of load-penetration curve Figure 10. Tip force as a function of penetration distance for 5 successive dynamic penetration tests for materials A, B and C Figure 11. Tip force as a function of penetration distance for 5 successive dynamic penetration tests for materials A, D and E Figure 12. (Left) <Rd> as function of D50 with five PSD A to E. (Right) <Rd> as function of PSD shape for three materials A, D and E (height of vertical bars is twice the standard deviation of Rd) Figure 13. (Left) <Af> as fonction of D50 for five particle size distributions A to E. (Right) <Af> as function of psd shape for three particle size distributions A, D and E (height of vertical bars represent twice the standard deviation of <Af>)

Table 1 .

 1 Parameters of the model. 318

	Parameter	Symbol	Value	Unit
	Width box	L	0.60	m
	Height box	H	0.45	m
	Particle density	ρ	2 700	kg.m -3
	Normal contact stiffness	kn	1.25 x 10 8	N.m -1
	Tangential contact stiffness	ks	9.375 x 10 7	N.m -1
	Rod friction coefficient	µrod	0.0	-
	Tip friction coefficient	µtip	0.3	-

Table 2 .

 2 Material characteristics used in penetration tests with different particle size distribution.

	Ref	Particle number	φmax [-] Dmax [mm] D50 [mm] Dmax/Dmin [-] Type
		Np
	A	10 000	0.844	7.02	5.80	2	I;II
	B	40 000	0.843	3.51	2.91	2	I
	C	160 000	0.843	1.76	1.45	2	I
	D	30 000	0.856	7.02	3.98	5	II
	E	97 831	0.869	7.02	2.61	10	II


D Collet
email: david.collet@ifpen.fr

M Alamir
email: mazen.alamir@grenoble-inp.fr

D Di Domenico

G Sabiron

Non quadratic smooth model of fatigue for optimal fatigue-oriented individual pitch control

Operation and maintenance cost of wind turbine is intimately related to fatigue damage. However, considering fatigue directly in an optimization needs to be carefully done because its faithful model does not fit standard forms. Recent results have shown that the fatigue damage of a signal estimated using fatigue theory is non-linearly related to its corresponding variance. Therefore, this relationship is used to design a convex cost function approximating fatigue. In this paper, the sensitivity of open-loop optimizations using this fatigue-oriented nonquadratic cost function to the time horizon considered is studied, in terms of both fatigue and computational cost. Moreover, in order to alleviate the undesirable effects of the computational cost, an efficient fixed-point iteration-based way to optimize this fatigue-oriented cost function is presented. Results show that the fatigue-oriented cost function allows to obtain sensitive fatigue reduction compared to optimizations using a family of classical quadratic cost functions, especially for long optimization horizons. It is eventually, shown that the prediction horizon length required for the non-quadratic criterion to be efficient in an MPC makes its use prohibitive in such framework. However, prospects of reaping the benefits of the fatigue-oriented optimization using imitation learning are given.

Introduction

The economic profit that an Horizontal Axis Wind Turbine (HAWT) owner makes is mainly driven by the difference between gains, e.g. energy sold, and losses, e.g. Operation and Maintenance (O&M) cost [START_REF] Loew | Proc. European Control Conference[END_REF]. O&M cost can be approximated as a weighted sum of the prices of replacement and mechanical fatigue damage of given HAWT components [START_REF] Collet | American Control Conference (ACC)[END_REF], as fatigue damage is equal to 1 when a component is broken. This weighted sum is also called fatigue trade-off cost. Therefore, in order to maximize HAWT economic profit, it is important to optimize this fatigue trade-off. Dedicated control of HAWT blade pitch angle can contribute to this cost reduction.

The main objectives of HAWT blade pitch control are to regulate output power and rotor speed while reducing mechanical fatigue. In the early work on blade pitch control, wind was assumed to be uniformly distributed over the rotor area. Therefore, all the blades were pitched to the same angle, this technique is called Collective Pitch Control (CPC). With recent increase in rotor diameter, this assumption becomes questionable. Aerodynamic forces on the blades fluctuate with the azimuth angle while the blade pitch angle remains constant [START_REF] Hansen | Aerodynamics of wind turbines[END_REF]. Therefore, by varying each blade pitch angle individually depending on its azimuth, blades fatigue loads

The Science of Making Torque from Wind (TORQUE 2020) Journal of Physics: Conference Series 1618 (2020) 022004 IOP Publishing doi:10.1088/1742-6596/1618/2/022004 2 can be alleviated. This technique is called Individual Pitch Control (IPC) [START_REF] Bossanyi | [END_REF]. IPC is usually divided in two stages:

(i) The CPC stage whose objective is to regulate the rotor rotational speed and power, while alleviating the fatigue loads on the tower. (ii) The IPC stage that gives a differential pitch angle on each blade, added to the collective pitch angle, in order to reduce the unbalanced loads on the rotor, contributing to the rotating components damage (i.e. blades, rotor bearing, blade bearings)

Therefore, the goal of the IPC stage is essentially to optimize the fatigue damage trade-off.

Optimal control is a good solution to achieve this objective. However, expressing this objective as an economic cost function for optimal control is a challenging task due to the fatigue damage estimation [5,6]. Fatigue damage is quantified using the Palmgrem-Miner fatigue theory, where it is expressed as a sum of damages caused by hysteresis load cycles [7]. These cycles are counted using a RainFlow Counting (RFC) algorithm [8] which cannot be expressed as a simple algebraic function. Consideration of fatigue as an explicit control objective remains an open topic [9]. To our best knowledge, it is possible to find in the literature five methods considering directly fatigue as an objective cost function for HAWT control. In [5], the spectral properties of the stress history are used to design an optimal feedback controller. In [START_REF] Collet | American Control Conference (ACC)[END_REF], a data driven surrogate model is designed to select the best controller among several candidates. In [10], a quadratic cost function parameters are varied to match an on-line estimated economic fatigue cost function. In [START_REF] Loew | Proc. European Control Conference[END_REF], a sequential optimization is made where the cost function parameters are iteratively modified between gradient steps, based on RFC. Only the two last methods are used to design an MPC controller, as their approach is temporal. In [START_REF] Collet | A fatigue-oriented cost function for optimal individual pitch control of wind turbines[END_REF], several aspects of the previous methods are bridged together. A data-driven surrogate model relating quadratic norms (variance) and fatigue of given signals is used in order to approximate the fatigue trade-off cost with a non-quadratic cost function. This approximated objective function is used to optimize the system input trajectory during an open-loop optimization.

This paper is organized as follows:

• In Section 2, the equations governing the rotor dynamics of an HAWT are presented.

• In Section 3 the parameterization of a quadratic cost function used as benchmark, which was firstly presented in [START_REF] Collet | A fatigue-oriented cost function for optimal individual pitch control of wind turbines[END_REF], is depicted. • In Section 4, the novel fatigue-oriented cost function presented in [START_REF] Collet | A fatigue-oriented cost function for optimal individual pitch control of wind turbines[END_REF] is introduced.

• In Section 5, an effective optimization method for optimal control using the presented fatigue-oriented cost function is detailed, involving a fixed-point method and a decomposition of the optimization horizon. • In Section 6, performance of an optimization using the fatigue-oriented cost function is compared to one using a parameterized quadratic cost function, along with suggestions suggestions for closed-loop implementation. • In Section 7, a conclusion and an outlook on the undergoing work are given.

Presentation of the system: An HAWT rotor in MBC coordinates

This section presents the system under interest, which is a multi-body HAWT disturbed by the wind, whose different bodies vibrations are coupled. It is a multi-inputs/multi-outputs system, that can be simulated using high fidelity nonlinear aero-elastic HAWT simulators, such as FAST [START_REF] Jonkman | Definition of a 5-MW reference wind turbine for offshore system development[END_REF]. It can be noted that FAST allows to linearize the nonlinear equations governing the HAWT behavior for control purposes. In this study, the focus is set on the IPC stage of an HAWT blade pitch controller. As mentioned previously, the IPC stage objective is to reduce the fatigue of multiple HAWT components by giving differential pitch angles on each blade, to be added to the collective pitch angle. These differential pitch angles can be obtained from a non-rotating reference frame by using the Coleman [START_REF] Coleman | National Advisory Committee for Aerounautics[END_REF] or Multi-Blade Coordinate (MBC) transform [START_REF] Bir | Proc. 46th AIAA aerospace sciences meeting and exhibit[END_REF].

FAST and its MBC module allows to obtain a linearization of the dynamics around an operating point, relating the hub-height wind velocity, denoted by v and, the yawing and tilting pitch angles, denoted respectively by θ yaw and θ tilt , to the yawing and tilting out-of-plane blade root bending moments, denoted respectively by M yaw and M tilt . The operating point is defined by its steady wind velocity v 0 , yawing and tilting steady out-of-plane blade root bending moments, denoted respectively by M 0 yaw and

M 0 tilt . Let δM = [M yaw -M 0 yaw , M tilt -M 0 tilt ]
T be the differential yawing and tilting out-of-plane blade root bending moments, moreover Θ = [θ yaw , θ tilt ] T , x ∈ R 6 is the state of the system corresponding to the blades first flapwise mode, and δv = v -v 0 is the differential wind speed. By augmenting the system obtained from FAST linearization and MBC modules with the actuators dynamics modeled by a first order low-pass filter, similar to the one presented in [START_REF] Collet | A fatigue-oriented cost function for optimal individual pitch control of wind turbines[END_REF], the following Linear Time Invariant (LTI) system can be obtained:

ẋ = Ax + Bu + B d δv (1a) y = Cx + Du + D d δv (1b) where y = [δM T , Θ T ] T , u = [θ sp yaw θ sp tilt ]
T , θ sp yaw and θ sp tilt denote respectively the yawing and tilting blade pitch angles to be given to the actuators via the inverse MBC transform. A, B, B d , C, D and D d are matrices of appropriate dimensions allowing to model the dynamics of the system. It should be noticed that in this paper, it is assumed that the system is not restricted by any constraints whatsoever, or that it cannot reach its constraints in the operational range it is working on. This can be forced by appropriate choice of the weighting coefficient and checked a posteriori. In the sequel, in order to alleviate notations δv := v.

Parameterized quadratic cost function benchmark

Quadratic cost functions are widely used in the formulation of fatigue-oriented open-loop optimal control problems. This is the reason why such cost functions are used as a benchmark for comparison with other open-loop formulations. The tuning of a quadratic cost functions can be computationally expensive and might be sub-optimal. However, as far as the tuning of a quadratic structure is concerned, for the specific case of HAWT rotors, it is possible to use this knowledge in order to derive a low-dimensional parameterization of the quadratic cost function weights.

Quadratic cost function parameterization

In [START_REF] Collet | A fatigue-oriented cost function for optimal individual pitch control of wind turbines[END_REF], a parameterization of a quadratic cost function for the specific problem of HAWT fatigue alleviation with IPC is described. The quadratic cost function is denoted by J quad and parameterized by ρ:

J quad (u|v, x 0 , ρ) = N k=1 y k (u, v, x 0 ) T Q(ρ)y k (u, v, x 0 ) + u T k R(ρ)u k (2)
where y k (u, v, x 0 ), u and v are respectively the output vector at instant k, inputs and disturbance trajectories over the optimization horizon. The index k represents the k th time instant, x 0 is the initial state of the open-loop optimization and N is the number of prediction steps considered in the optimization horizon. Q ∈ R 4 and R ∈ R 2 are respectively semi-definite and definite positive matrices, parameterized by ρ > 0, such that that y is therefore a function of u, v and x 0 , driven by the dynamic system (1). Moreover, this parameterization is designed such that the variations of blade pitch angles set-points do not impact the fatigue trade-off.

Q(ρ) = diag ([ρ, ρ, 1, 1]) and R(ρ) = min(ρ,

Parameterized quadratic cost function tuning procedure

The parameterized quadratic cost function is thus used to solve the following unconstrained open-loop optimal control problem:

min u J quad (u|v, x 0 ) (3)
for given v and x 0 . This open-loop optimization yields an output trajectory, which can be fed to the fatigue cost function denoted by J fat (.) in order to estimate the output trajectory fatigue cost. Therefore, each triplet (ρ, x 0 , v) corresponds to a fatigue cost, denoted by J fat (ρ, x 0 , v).

The tuning procedure consists in finding the parameter ρ that minimizes the fatigue cost for given x 0 and v:

min ρ J fat (ρ, x 0 , v) (4)
where x 0 and v are randomly drawn from a given relevant distribution to be carefully chosen.

The novel fatigue-oriented cost function

In [START_REF] Collet | American Control Conference (ACC)[END_REF], the fatigue cost function J fat is defined to be a weighed sum of the damages and prices of replacements, denoted respectively by D i and π i for the i th component, out of n in total:

J fat (y) = n i=i π i D i (y i ) (5)
where y i is the i th output trajectory.

In [START_REF] Collet | A fatigue-oriented cost function for optimal individual pitch control of wind turbines[END_REF], it was shown that it is possible to derive a data-driven relationship between a variance such as variance of a signal and fatigue damage, thanks to a linear regression between the logarithms of variance and fatigue damage of a given signal. This allows to approximate the fatigue damage of the i th component considered in the fatigue cost estimation with a function denoted by Di :

Di (y i ) = e b i (L ult i ) m i Var(y i ) a i (6)
where L ult i , m i , a i and b i are respectively the ultimate load, Wöhler exponent, linear coefficient and bias of the regression corresponding to the i th component. It should be noticed that the ultimate load is the maximum load that the component can bear and the Wöhler exponent depends on the material (4 for steel and 10 for glass fiber). Therefore, J fat can be approximated by the fatigue-oriented cost function, denoted by Jfat , for a given output trajectory y:

Jfat (y) = n i=1 π i e b i (L ult i ) m i (Var(y i )) a i (7)
Hence, Jfat allows to approximate J fat with a differentiable cost function based on variance. This cost function can be efficiently implemented in an optimization, in order to better reduce J fat in optimal control problems.

The fatigue-oriented optimization

The cost function Jfat is particularly interesting for optimal control problems because it is convex provided that a i ∈] -∞, 0]∪]1, +∞] ∀i ∈ {1, . . . , n}, differentiable and expressed algebraically. Therefore, there is no loss of convexity nor differentiability in using it instead of standard quadratic function of time series. The resulting problem can be solved directly using standard NonLinear Programming (NLP) solvers:

min u Jfat (y(x 0 , u, v)) (8)
where x 0 and v are given. Moreover, as Jfat is based on variance, it is possible to efficiently solve the problem with sequential Quadratic Programming (QP) in a fixed-point problem, by approximating Jfat with a first order Taylor expansion. Besides, it is possible to decompose the horizon of the successive QP free of inequality constraints, in order to significantly reduce the computational time. In this section are presented the derivation of the fixed-point problem from the NLP and the decomposition of the horizon, which together significantly reduce the computational time needed for the resolution of the NLP.

Fixed-point problem formulation

It is a common approach to approximate NLP with QP which are easier to solve than NLP and convex, in order to efficiently approach a suboptimal solution of the NLP. The approximated fatigue cost function Jfat has the advantage to be based on a quadratic cost function which is variance. Therefore, a first order Taylor expansion of Jfat , defined by (7), around an output trajectory y (σ) is sufficient to make appear a quadratic form:

J(σ) fat (y) = Jfat y (σ) + n i=1 π i e bi a i (L ult i ) mi Var(y (σ) i ) ai-1 αi(y (σ) )
Var(y i ) -Var(y

(σ) i ) (9)

Relationship between fatigue damage and variance

The experience described in Section 4 is conducted and the results are plotted in Figure 2. The relation can be effectively approximated by a linear regression in Logarithmic scale and yield couples of parameters (a, b) of value (2.6, 4.6) and (5.6, 15.6), for Wöhler coefficients of 4 and 10 respectively.

Fatigue reduction and computing time-sensitivity analysis

Concerning the sensitivity study of optimal control problem solutions using the parameterized quadratic and non-quadratic fatigue-oriented cost function to the horizon's length, a sampling time of 0.1 second is considered.

For both optimizations, the optimization horizon's temporal decomposition scheme described above is implemented. The parameter M takes the values {5, 5, 5, 5, 5, 10, 10, 10, 10, 10, 20, 20, 20, 20, 20, 20} together with the horizon lengths {2, 3, 4, 5, 7.5, 10, 20, 30, 40, 50, 100, 200, 300, 600} seconds respectively. The parameterized quadratic cost function, defined in Section 3, is optimized for a disturbance whose filter parameters are randomly drawn, in the same distribution that the one used for disturbance generation. In this presented case, the static gain and characteristic time drawn are 6950 and 2.86 seconds respectively. The optimization of ρ is conducted for each horizon considered. The fatigue cost of the 1000 optimizations solutions using the non-quadratic cost function, defined by (7), is estimated and denoted by J NQ fat . The fatigue cost of the optimization solutions using the finely tuned parameterized quadratic cost function, defined by ( 3) is also estimated and denoted by J Q fat . For the fatigue estimations, parameters are carefully chosen in order to yield a realistic fatigue trade-off. The values of L ult , π and m for the yawing and tilting blade root bending moments are 3 × 10 3 , 10 3 , 10 respectively, while for the yawing and tilting pitch angles they are respectively 6.98 × 10 -1 , 1 and 4. In Figure 3, the mean, the median and the 90 th percentile of the ratio between the fatigue costs open-loop optimal trajectories using the non-quadratic and quadratic cost functions, defined by

J NQ fat J Q fat are
shown. The thick red line represents the border between cases where non-quadratic optimization induces higher fatigue cost compared to the quadratic one (above), and cases where it induces lower fatigue costs (below). It can be seen that the longer the time horizon is, the smaller the relative fatigue. More precisely, the median relative fatigue cost is lower than 1 for time horizon length longer than 4 seconds, while for the mean and 90 th percentile of the relative fatigue cost it is respectively 7.5 and 50 seconds. Moreover, beyond 20 seconds, the fatigue reduction ratio induced by the non quadratic cost is as high as 50%.

In Figure 4, the mean and the range of the CPU times needed for the 1000 quadratic and non-quadratic open-loop optimizations are plotted as a function of the horizon time length, for optimization run on a desktop PC equipped with a 2.60GHz processor and 16Gb RAM. It can be seen that the optimization using the non-quadratic cost function requires a much longer CPU time than the quadratic one. This results was expected as the fixed-point method requires to run a quadratic optimization per iteration.

In Table 1 are summarized the CPU times needed to solve the optimization problem defined by (8) for N = 100 and 1000, using the NLP solver CasADi [17] and the fixed-point method with and without a decomposed horizon. It can be seen that, compared to the NLP solver, the fixed-point method with a decomposed horizon is much less time consuming, especially for long optimization horizons. Moreover, the fixed-point method avoids the compilation time, which even though it is performed off-line, can be highly time consuming for long optimization horizons and mandatory for the NLP solver. It also emphasizes the efficiency of the horizon decomposition, which allows to decrease the CPU time from one to three order of magnitudes compared to the fixed-point method without decomposition. Even though, they might have possibilities to decompose the optimization horizon of the NLP as well, the optimization problem managing the initial states of the intervals might be very challenging to express for the NLP case, while the fixed-point method allows to consider only QPs and to relatively easily decompose the optimization horizon.

Prospects of closed-loop implementation

The non-quadratic criterion is used in an MPC and simulated in closed-loop for several prediction horizons ranging from 0.2 to 2 seconds. This MPC is compared to another MPC using the parameterized quadratic cost function on 100 closed-loop simulations of 200 seconds long. In Figure 5, the fatigue averaged on the 100 closed-loop simulations and normalized by the average fatigue of the quadratic MPC for a 0.2 seconds prediction horizon is plotted in function of the prediction horizon. It can be seen that for relatively short prediction horizon, compared to the ones used for the open-loop optimization, the MPC using the non-quadratic cost function is less efficient in reducing fatigue than the quadratic MPC. This coincides with the previous observation, stating that the open-loop optimization using the non-quadratic criterion becomes efficient compared to the parameterized quadratic one for horizons longer than 4 to 50 seconds. This is due to the fact that considering fatigue on such short prediction horizons is very limited as very few hysteresis cycles can be counted, for the dynamics considered in this application and therefore the statistical relation observed between variance and damage might be weak. However, the CPU time needed for the open-loop optimizations using the non-quadratic cost function makes its use in an MPC prohibitive. Nevertheless, it might be possible to learn from the open-loop optimization using the non-quadratic criterion the behavior of a non-linear feedback controller, using information on the future disturbance on an acceptable receding horizon.

Conclusion and perspectives

In this paper, the sensitivity analysis of fatigue-oriented open-loop optimal cost and its computational burden with respect to the length of the prediction horizon is conducted for different choices of cost functions. The analysis shows that the use of a non quadratic cost induces significant benefit over standard quadratic (even finely tuned) cost functions provided that prediction horizons longer than 50 seconds are involved. However, the computational burden of the fatigue-oriented optimization becomes rapidly intractable to be incorporated in MPC design. Undergoing works focus on the extraction of fitted static nonlinear feedback from the off-line computation of open-loop formulations with long prediction horizons and non quadratic cost leading to solutions that inherit the benefit from the newly introduced cost while being real-time implementable. Preliminary results on this behavioral cloning solution are encouraging and will hopefully be presented in a forthcoming paper.

10 - 2

 2 10 -1 10 0 10 1 10 2 10 3

Figure 2 .

 2 Figure 2. Plot of the relation between Var and (L ult ) m D in Logarithmic scale, for two Wöhler coefficients.

Figure 3 .

 3 Figure 3. Plot of the influence of prediction horizon time length on the relative system fatigue cost between the non-quadratic and quadratic cost functions.

Figure 4 .

 4 Figure 4. Influence of time horizon length on the CPU time. Thick lines correspond to the mean and the shaded areas to the range of values.

Table 1 .

 1 Summary of the compilation and solving time for the NLP solver, fixed-point and fixed-point method with a decomposed horizon.

	Relative Fatigue	Q MPC NQ MPC NQ OL
		Predicition horizon (s)
	Figure 5. Plot showing that the prediction
	horizon typically used in common MPC
	design is far too short to be relevant.

where σ is the iteration number of the fixed-point problem.

Let Q(y (σ) ) = diag [α 1 (y (σ) ), . . . , α n (y (σ) )] and R(y (σ) ) = I 2 min [1, α 1 (y (σ) )] × 10 -3 be the weighting matrices of an QP that approximates the solution of (8) at the σ th iteration:

where R avoids singularities in the optimization problem. It should be noticed that the optimization (10) gives an optimal solution u (σ+1) , which integrated with the system dynamic equation [START_REF] Loew | Proc. European Control Conference[END_REF] for given x 0 and v yield an output trajectory y (σ+1) . This y (σ+1) gives a new optimal control problem P(x 0 , v, y (σ+1) ), defined by (10), resulting in a u (σ+2) different from u (σ+1) , so on and so forth. There is thus a fixed-point problem consisting in finding the right u (σ) and y (σ) . The following algorithm is proposed to solve this fixed-point problem:

(i) σ ← 0, Give an initial guess for y (0) and u (0) (ii) Compute the weighting matrices Q(y (σ) ) and R(y (σ) ) (iii) u (σ+1) ←-Solve the optimal control problem P(x 0 , v, y (σ) ) defined by (10) (iv) y (σ+1) ←-Integrate the dynamic equations, defined by (1), for u (σ+1) , x 0 and v (v) If y (σ) and u (σ) have not converged, increment σ and go to step (ii), otherwise stop This fixed-point method is used in the sequel instead of NLP solvers. However, for a high number of decision variables, the QP is computationally and time expensive because of the inversion of high dimensional matrices. That is the reason why a time decomposition is introduced in the next section in order to define a change in the decision variable's definition and significantly reduce the RAM and time needed for the QP resolution.

Optimization horizon decomposition scheme

The solution to the QP defined by ( 10) is a linear combination of x 0 and v, which requires a matrix inversion of dimension [N × 2, N × 2]. The longer the prediction horizon is, the more RAM and CPU time are required, moreover the QP must be solved at every iterations of the fixed-point method. Therefore, reducing significantly the time needed for the QP resolution can also significantly reduce the time to solve the fixed-point method and the NLP defined by (8).

The optimization horizon of an QP can be temporally decomposed as in [START_REF] Chang | A hierarchical, temporal decomposition approach to long horizon optimal control problems[END_REF][START_REF] Xu | [END_REF], allowing to significantly alleviate the computational cost. The temporal decomposition of an QP optimization horizon consists in decomposing the problem into a series of smaller QPs with shorter optimization horizons, constrained by their initial and final states, provided that the cost function on the interval is independent from the others. Then, the solution to (10) is recovered by solving another QP that manages the initial states of the smaller QPs. However, for the QP defined by (10), the cost function on the smaller intervals is not independent from the others, as µ, which is the average of the outputs on the whole horizon has an influence on the cost function. Therefore, µ is first considered as an exogenous parameter in the smaller QPs and retrieved later with the QP managing the initial states. A schematization of the optimization horizon decomposition is proposed in Figure 1. Let M be the number of intervals the horizon is decomposed into, which is also the number of smaller QPs. The QP on the j th interval, for j ∈ {1, . . . , M }, at the σ th iteration, denoted by QPj, is parameterized by the initial condition and disturbance on the interval, denoted respectively by x (j) 0 and v (j) , and the initial condition on the next interval x (j+1) 0

. All the QPj are using the same cost function parameterized by y (σ) and µ which are considered as external parameters:

The where u (j) is the input trajectory on the j th interval, u (j) k is the input value at the k th time instant in u (j) , k 0 (j) = (j -1) N M + 1 and k f (j) = j N M are respectively the index of the first and last time instant of the j th interval. y (j) k and x f are affine functions giving respectively the output value at the k th time instant and the j th interval final state for given x 0 , v (j) and u (j) . It should be noticed that that x (1) 0 = x 0 , which is the initial state of the QP defined by (10). The solution to [START_REF] Collet | A fatigue-oriented cost function for optimal individual pitch control of wind turbines[END_REF], denoted by u (j,σ) , is an affine function of x (j) 0 , x (j+1) 0 , v (j) and µ, besides each resolution of the QPj needs a matrix inversion of dimension [ N M ×2, N M ×2]. Therefore, the output trajectory, denoted by y (j) , which is a concatenation of the y

is an affine function of x (j) 0 , u (j,σ) and v (j) , moreover by composition it is an affine function of x (j) 0 , x (j+1) 0 , v (j) and µ. Hence, J j can be expressed as a quadratic function of x

] be the vector of the interval initial conditions. The cost function J j can thus be expressed as quadratic cost function of X 0 and µ, given the disturbance trajectory v. Therefore the cost function J(σ) fat (y) of the QP defined by ( 10) can be expressed as the sum of the J j ∀j ∈ {1, . . . , M }. Moreover, µ is the average of y which is a concatenation of the y (j) ∀j ∈ {1, . . . , M }, which makes thus µ an affine function, denoted by F, of X 0 and µ, given x 0 and v. The master QP than manages the vector of initial conditions X 0 , such that the concatenation of y (j,σ) ∀j ∈ {1, . . . , M } matches the solution of ( 10) is the following:

where the implicit equation defined by (12b) is actually a linear equality constraint allowing to retrieve µ. The solution of ( 12) is thus a linear combination of x 0 and v, requiring a matrix inversion of dimension [6 × (M + 1), 6 × (M + 1)] as the state dimension is 6. Eventually, the optimal input trajectory u (σ) can be recovered by concatenating the solutions of the QPj defined by [START_REF] Collet | A fatigue-oriented cost function for optimal individual pitch control of wind turbines[END_REF], denoted by u (j,σ) with the x (j) 0 and x (j+1) 0 obtained from the solution of [START_REF] Jonkman | Definition of a 5-MW reference wind turbine for offshore system development[END_REF].

From the above discussion, it comes out that M determines the computation burden involved in the solution of [START_REF] Collet | A fatigue-oriented cost function for optimal individual pitch control of wind turbines[END_REF] and [START_REF] Jonkman | Definition of a 5-MW reference wind turbine for offshore system development[END_REF]. The computational complexity of the QP problem, for the considered system (1) without temporal decomposition is O (2N ) 2 (2N + 6 + N ) , while using the temporal decomposition it becomes O M 6(M -1) 2 (6 + 1) + (2N ) 2 M ( N M (2 + 1) + 12 + m) . This means that if N takes values of 100, 1000 or 10000, the M giving the lowest computational complexity are respectively 5, 16 and 49, allowing to divide the computational complexity of the global QP problem by 12, 120 and 1200 respectively.

Results

In this section, an assessment regarding the relationship between variance and fatigue damage is first presented. Then the sensitivity of the optimal solution to horizon's length is analyzed in terms of both mechanical fatigue and CPU time. Furthermore, comparison with optimization using a finely tuned quadratic cost function is proposed. Finally, in a prospect of using the openloop optimization with the fatigue-oriented cost function in closed-loop, a sensitivity analysis to the prediction horizon of a linear quadratic MPC in terms of fatigue is conducted. Statistics on each results are made using 1000 disturbances which are Gaussian noises filtered through a family of first order low pass filters. The filters are parameterized by their static gain and characteristic time, taking values in uniform distributions ranging respectively on ]0, 10 4 ] and ]0, 10]. In the sequel, all simulations are considering an initial state located at the origin.

The Science of Making Torque from Wind (TORQUE 2020) Journal of Physics: Conference Series 1618 (2020) 022004


Léa Blaise
email: lblaise@localsolver.com

Christian Artigues

Thierry Benoist

Solution Repair by Inequality Network Propagation in LocalSolver

Keywords: constraint propagation, inequality networks, local search, solution repair, solver, disjunctive scheduling

Introduction and context

This paper introduces a solution repair algorithm based on constraint propagation, designed to overcome the difficulties met by small neighborhood search on a family of tightly constrained problems.

LocalSolver is a mathematical programming solver, whose goal is to offer a model-andrun approach to optimization problems, including combinatorial, continuous, and mixed problems, and to offer high quality solutions in short running times, even on large instances. It allows OR practitioners to focus on the modeling of the problem using a simple formalism, and then to defer its actual resolution to a solver based on efficient and reliable optimization techniques, including local search (but also linear, non-linear, and constraint programming). The local search algorithms implemented in LocalSolver are described in [START_REF] Gardi | Mathematical Programming Solver Based on Local Search[END_REF].

This paper focuses on problems whose constraints include a network of two-variable inequalities. More precisely, it focuses on problems whose constraints comprise either linear inequalities between two variables or disjunctions of such inequalities. Many such problems arise in the field of scheduling: for example the Job Shop Problem [START_REF] Fisher | Probabilistic learning combinations of local job-shop scheduling rules[END_REF] and the Bridge Building Problem [START_REF] Bartusch | Optimierung von Netsplanen mit Anordnungsbeiziehungen bei knappen Betriebsmitteln[END_REF], which are both characterized by generalized precedences and disjunctive resource constraints. However, this kind of structure is also typical of packing, layout, and mining problems.

These problems are highly constrained: in a good solution of a Job Shop or Bridge instance, the precedences and disjunctive resource constraints are often very tight. Because of that, moving from a solution of makespan x to a solution of makespan x -1 requires a lot of small changes on many integer variables -the start times of the tasks. Being able to move from a good feasible solution to another using random small neighborhoods is then very unlikely: one would have to randomly target the right set of integer variables and to randomly shift them all by the right amount. For these reasons, the algorithms described in [START_REF] Gardi | Mathematical Programming Solver Based on Local Search[END_REF] encounter serious difficulties on these problems. In the vast literature on job-shop scheduling by local search, these difficulties are overcome by exploiting higher level dedicated solution representations, such as the disjunctive graph (Vaessens et. al. 1994). In this work, we aim at keeping the modeling elements simple and we wish to target other problems as well. Hence, we focus on the direct integer variable representation.

Solution repair

The solution that we envisioned and implemented in LocalSolver to tackle this problem is a kind of constraint propagation: a promising but infeasible solution is gradually repaired one constraint at a time. This gradual procedure might remind one of the min-conflicts heuristic, used to repair an inconsistent assignment for the variables of a CSP, or of ejection chains algorithms, both recently studied in the context of scheduling problems, respectively in [START_REF] Ahmeti | Min-conflicts Heuristic for Multi-mode Resource-constrained Projects Scheduling[END_REF] [START_REF] Ding | Parallel machine scheduling with completion-time-based criteria and sequence-dependent deterioration[END_REF].

Before the search, specific constraints (such as generalized precedences and disjunctive resource constraints) are detected in the model and made into a constraint network. Those constraints will be the ones on which the repairing algorithm will apply when necessary. Let's consider any iteration of the local search, and let's assume that the current solution S 0 is feasible. A local transformation is applied to the solution (for example, modification of a few start times), rendering it infeasible (S). In the simplest case, the repair phase is similar to the classic constraint propagation of Constraint Programming [START_REF] Rossi | Handbook of Constraint Programming (Foundations of Artificial Intelligence)[END_REF]), but has a few singularities. First, while constraint propagation generally aims at reducing the variables' domains, here the constraints are only propagated when they are violated and need repairing. Indeed, our only concern here is to build a feasible solution S as close to S as possible. Then, we impose that, during the local move and along the successive constraint repairs, the variables must always be shifted in the same direction: if a variable's value was increased, it can still be increased further when repairing a constraint, but it can never be decreased, and reciprocally. This ensures that all the decisions taken during the current iteration (during the successive constraint repairs, and more importantly during the local move) are respected. This way, the successive repairs "follow the move's direction": they lead to finding a feasible solution as close as possible to S, by amplifying the move rather than by cancelling it. At the end of the propagation, the algorithm was either able to repair all the constraints, and thus found a feasible solution S , or found a constraint that could not be repaired. In the latter case, all the changes are cancelled and the whole procedure starts over: the current solution is set back to S 0 , and a new move is applied.

Two-variable linear inequalities

We consider inequalities of the form

aX + bY ≤ c
where X and Y are integer or floating point variables, and a, b, and c are any constants. When a = 1 and b = -1, these inequalities correspond to the generalized precedence constraints encountered in scheduling problems, but the following algorithm is not limited to this specific case.

Let's assume that the inequality aX + bY ≤ c is violated. Since the initial solution was feasible, at least X or Y has already been shifted in the "wrong" direction, and therefore cannot be shifted now in the repair direction. At most one variable (X by symmetry) can then be shifted in the repair direction: there exists only one way to repair the constraint, which consists in shifting X just enough (X ← c-bY a ). This is a necessary decision, as in every feasible solution following the move's choices, X is at least (resp. at most) c-bY a . When the only repairable constraints in the model are inequalities, the repair phase is equivalent to a particular kind of constraint propagation, which will be referred to as "half bound consistency" in the remainder of the paper. The reduction of a variable X's domain is propagated only if it excludes its current support x (only if the constraint is violated and needs repairing). After repairing the constraint, the new support of X is written x . When propagating the reduction of X's domain, only one of its bounds is modified: either x > x (then X's value can no longer be decreased) and x is its new lower bound, or x < x and x is its new upper bound. Since each variable must always be shifted in the same direction, one of its bounds at most can be modified throughout the propagation.

An interesting property when the only repairable constraints in the model are inequalities is that, if there exists a feasible solution that respects the decisions of the local move, the algorithm is guaranteed to find one at the end of the propagation.

Disjunctions of two-variable linear inequalities

We consider disjunctions of inequalities of the form

i (a i X i + b i Y i ≤ c i )
where the X i and Y i are integer or floating point variables, and the a i , b i , and c i are any constants. When a i = 1 and b i = -1 ∀i, these disjunctions correspond to disjunctive resource constraints (disjunctions of size 2), or packing constraints in higher dimensions. However, the algorithm is not limited to this specific case here either.

Let's assume that a disjunction is violated. We will try to repair it in a non deterministic way. Since a priori none of the generalized precedences of the disjunction should prevail over the others, the algorithm chooses one at random and tries to repair it. If it cannot be repaired, it tries to repair the following one, and so forth. If none of them could be repaired, the propagation fails.

Let aX +bY ≤ c be the inequality that was randomly chosen for repair in the disjunction. If only one of its variables can be shifted in the repair direction, then the constraint is repaired as described in 2.1. It is also possible that both variables can be shifted in the repair direction, since the chosen inequality may not have been the one that was respected in the initial feasible solution S 0 . If so, the algorithm randomly chooses how to shift them. Let ∆ be the distance to feasibility: ∆ = aX + bY -c, and let δ X and δ Y be the shares of the repair respectively attributed to X and Y , verifying δ X + δ Y = ∆. There are four possible ways for the algorithm to repair the constraint: either X repairs it alone (δ X = ∆), or Y repairs it alone (δ Y = ∆), or X and Y equitably share the repair (δ X = δ Y = ∆ 2 ), or X and Y randomly share the repair (δ X = random(1, ∆ -1) and δ Y = ∆ -δ X ).

Since the repair procedure of a disjunction is non deterministic, the previous properties, of half bound consistency and guarantee to find a feasible solution, do not hold anymore when such constraints are detected among the repairable constraints in the model. However, if there exists a solution that respects the decisions of the move, there is always a non-zero probability that the propagation leads to finding one, depending on whether the algorithm always takes the "right" random decisions when repairing a disjunction.

Application to scheduling problems

This method of solution repair by constraint propagation dramatically improves the performances of LocalSolver on problems with a network of two-variable linear inequalities.

Results on the Bridge Building Problem

The optimum value of the Bridge Building Problem is 104. Without our repair mechanism, within ten seconds of search, LocalSolver 9.0 is only able to find solutions of value 115 on average, and virtually never finds an optimum solution. But with the integrated repair mechanism, it always finds a solution of value 104 within four seconds of search, and very often finds one in less than one second.

Results on the Job Shop Problem

We compared the performances of LocalSolver 9.0 with and without this repair mechanism on three classic Job Shop instance classes: the FT class by [START_REF] Fisher | Probabilistic learning combinations of local job-shop scheduling rules[END_REF], the LA class by [START_REF] Lawrence | Resource constrained project scheduling: an experimental investigation of heuristic scheduling techniques[END_REF], and the ORB class by [START_REF] Applegate | A computational study of the job-shop scheduling problem[END_REF]. As shown above in table 1, our repair mechanism enables LocalSolver not only to find very good solutions of many Job Shop instances, but also to find good solutions very quickly: less than 10% from the optimum within 10 seconds of search.

Conclusion

In this paper, we considered a family of problems, whose constraints comprise a network of two-variable linear inequalities. Although small neighborhood search algorithms may encounter serious difficulties on these problems, we introduced a solution repair algorithm based on constraint propagation, overcoming the difficulties met by small neighbourhood search. The two main specificities of our propagation algorithm are that a domain reduction is only propagated if it excludes the current support of the variable, and that each variable must always be shifted in the same direction. Its integration into LocalSolver dramatically improves its performances on the targeted problems, not only on classic scheduling problems such as the Job Shop Problem, but also on some 3D packing and mining industrial instances of our test base.

Table 1 .

 1 Evolution of the gap between the average solution found by LocalSolver and the optimum solution, in 10 seconds and in 60 seconds, on different Job Shop instance classes

		Gap 10s		Gap 60s
	Instance class No repairs With repairs No repairs With repairs
	FT	73%	7%	15%	3%
	LA	246%	10%	91%	4%
	ORB	120%	6%	22%	4%


Benny Reinmart
email: benny.reinmart@gmail.com

Séverine Le Roy
email: severine.leroy@laplace.univ-tlse.fr

Gilbert Teyssedre
email: gilbert.teyssedre@laplace.univ-tlse.fr

Ngapuli I Sinisuka
email: n_sinisuka@yahoo.com

Effect of Cumyl Alcohol in Multilayer Dielectric on Space Charge Build Up

Keywords: space charge, crosslinking byproducts, cumyl alcohol, low density polyethylene

One of the issues of using cross-linked Polyethylene (XLPE) in DC insulation systems is the propensity of crosslinking byproducts to promote charge buildup inside the insulation, leading potentially to insulation failure. Field-induced ionization of the molecules, like in insulating liquids, is thought to be one of the processes of charge generation. In this paper, space charge measurement is realized to probe the effect of cumyl alcohol as one of the crosslinking byproducts. Low Density Polyethylene (LDPE) is used as polymer matrix, and soaked in cumyl alcohol. We show that polyethylene naphthalate (PEN) layers constitute efficient barriers to the evaporation of cumyl alcohol from LDPE films. The space charge behavior of soaked LDPE is compared to that of reference LDPE with PEN. Significant charge build-up occurs at the LDPE/PEN interfaces due to the conductivity gradient. PEN appears not suited for such analysis as the field is concentrated in the layer and is weak in the bulk LDPE.

I. INTRODUCTION

Crosslinked Polyethylene (XLPE) is used for long as insulation in High Voltage AC systems, particularly cables. XLPE has better thermomechanical stability, lasts longer, and shows less moisture sensitivity compared to LDPE [START_REF] Peake | The history of high voltage direct current transmission[END_REF]. Despite those advantages, there is some weakness in XLPE due to the effect of crosslinking byproducts and these effects are detrimental to the use of XLPE for HVDC systems.

During the cross-linking process, peroxides decompose to create free radicals that can initiate chemical bonding between polymeric chains. This process forms byproducts such as dimethyl benzyl alcohol (cumyl alcohol), acetophenone, and α-methylstyrene. These byproducts remain in the insulation and migrate out slowly during the time and affect the electrical properties of the material. Conductivity is generally increased in presence of byproducts, but the underlying mechanisms are not agreed upon. Three processes can be envisaged, being injectionassisted process for electronic carriers [START_REF] Li | Influence on spatial charge distribution of cross-linking agent residues in XLPE[END_REF], transport-assisted process for these charges [START_REF] Hussin | The effect of acetophenone and alphamethylstyrene on the space charge properties of low density polyethylene[END_REF], and promotion of ionic transport either by direct ionization under field or indirectly through ion solvation effects [START_REF] Kobayashi | Effect of acetophenone on the space charge evolution in LDPE and LLDPE[END_REF].

Byproducts can also cause deep chemical traps in XLPE, i.e. centers that stabilize electrical charges as it is probed with photoluminescence techniques [START_REF] Tardieu | Role of additives as recombination centres in polyethylene materials as probed by luminescence techniques[END_REF]. Hirai et al. show that byproducts increase the conductivity, but only cumyl alcohol creates carrier traps [START_REF] Hirai | Effect of byproducts of dicumyl peroxide on space charge formation in low-density polyethylene[END_REF].

Our purpose in this work is to investigate possible ionization effects associated with crosslinking byproducts. Though the matrix is a solid, this has relation to dielectric liquids, as ionization processes appear to be much better understood in liquids than in solids, with as common driving parameter the ionization potential of molecules [START_REF] Smalø | Field dependence on the molecular ionization potential and excitation energies compared to conductivity models for insulation materials at high electric fields[END_REF] to explain streamer propagations in liquids [START_REF] Davari | Field-dependent molecular ionization and excitation energies: Implications for electrically insulating liquids[END_REF] and resistance to treeing processes in solids [START_REF] Ashcraft | Laboratory studies of treeing in solid dielectrics and voltage stabilization of polyethylene[END_REF][START_REF] Englund | High efficiency voltage stabilizers for XLPE cable insulation[END_REF] for example. The other link to dielectric liquids is the need to maintain the liquid additive into the matrix and control the product uptake.

In this work, cumyl alcohol is chosen as investigated byproduct. It is incorporated into Low Density Polyethylene (LDPE) as a base material of XLPE. The tests consist of space charge measurement using Pulse Electro-Acoustic (PEA) method and conduction current measurements. To keep the cumyl alcohol amount as much controlled as possible inside the specimen for several hours during the electrical characterization, films of Polyethylene Naphthalate (PEN) are intercalated at both faces of LDPE samples. PEN was chosen since available in relatively thin form and for its good barrier properties. This constitutes an originality of the work carried out. Though the purpose of the measurement is not to investigate in particular PEN properties, its impact on the results must be considered.

II. SAMPLES PREPARATION

A. Materials

The LDPE samples used in these measurements have an average thickness of 340 μm and diameter of 8 cm. They were obtained by press-molding technique at 150°C from LDPE pellets. Polyimide layers were used to ease the sample removal after the molding process is completed. Polyimide layer is chosen because it does not affect the LDPE sample with decomposition products when it is processed at high temperature. To approach the desired thickness, an epoxy spacer with 350 μm thickness is selected. PEN films, 12 or 25µm thick, were cut into circular shape following the shape of LDPE sample. Gold metallization is done on one face of each PEN layer. The purpose of the metallization is to have same electrode conditions on each side of the sample. Gold electrodes were deposited by cold sputtering with diameter of 2 cm and a thickness of 50 nm. Thus, the PEN film will be placed in upper and bottom side of the LDPE sample. This structure is used in each measurement as multi-layered dielectric. The illustration of the sample can be seen in Fig. 1. time. The evaporation rate was measured after dipping for 5h.

B. Cumyl Alcohol uptake

The first measurement is to observe the PEN barrier effect by comparing samples soaked with and without PEN for up to 5 hours at 70°C. The absorption rate is taken every 30 minutes, and before the sample is put in balance to weight, it is rinsed with acetone, ethanol, and deionised water to remove cumyl alcohol from the surface of the sample. Fig. 2a shows that after 5 hours both samples have a same cumyl alcohol percentage (5.34% and 5.35%).

The evaporation rate is an important characteristic as it controls the amount of chemical during electrical measurements. For the evaporation rate, the weighing process is repeated in intervals of 20 min for the first 4 hours after the dipping process, and the data is continually recorded until 18 hours after the dipping process. Both samples are weighed at room temperature. Sample 1 is weighed without PEN layers, while sample 2 is weighed with 12 µm thick PEN layers on the entire surface.

To make sure that the cumyl alcohol is not absorbed by PEN layers, the later were weighed before and after the measurement process and no weight variation was detected. Based on Fig. 2.b, which shows that after 18 hours, the percentage of cumyl alcohol in sample 2 is substantially higher than in sample 1, quantitatively it can be stated that the effect of PEN layers as a barrier to cumyl alcohol evaporation is effective. Thus, this weight measurement will be done after each space charge and conduction current measurement. The same protocol, depicted in Fig. 3, was applied for every sample configuration with 30 min polarization and 30 min depolarization. The applied field is from 5 to 50 kV/mm for positive polarity and then -5 to -50 kV/mm for negative polarity. Considering the presence of cumyl alcohol, all the space charge measurements were taken at 25°C to prevent any evaporation. The temperature is controlled automatically in the oven chamber. Charge profiles were recorded all along the protocol in 30 min intervals.

III. ELECTRICAL MEASUREMENTS

B. Space charge measurement results

Fig. 4 shows space charge results obtained for the different materials. Data for 25 µm thick PEN layers are not shown for sake of space saving.

a) Untreated LDPE: The space charge pattern obtained in LDPE alone is plotted in Fig. 4a. It constitutes the reference for the results. LDPE was not provided with gold electrodes. From the cartography, positive charge dominates in the insulation, consistently with previous results obtained on the same material. The charge density inside the bulk increases with the increase of applied voltage on the sample. After polarity reversal a significant amount of negative charges is detected. Presumably, the positive charge stored in the first run increases the field at the electrodes and triggers negative charge injection. b) Untreated LDPE with PEN layers. Fig. 4b shows the charge pattern when 12 µm thick gold-metallized PEN layers are interposed at both sides of the LDPE. LDPE is not soaked with cumyl alcohol in this case, to see the effect of PEN layers on space charge behavior. At time 4-5 h a positive heterocharge builds up from 30 kV/mm near the ground electrode (bottom) and stays there until the applied field is -30 kV/mm. The charge distribution inside the bulk is having a periodic repetition of plans of charges with alternated polarity. The (positive) charge build up occurs at the region 0-25μm which corresponds to the region of PEN layers in the sample. It can be assumed that the charges build up is due to the existence of PEN layers on the sample. In positive polarity especially at low field, there is no heterocharges detected. To distinguish the effect of the PEN layer, space charge measurements were achieved with 25 µm thick layers. The pattern (not shown) is qualitatively similar to that of Fig. 4b but the charge amount is substantially less. The positive heterocharges build up appear at an applied field of 30 kV/mm. The charge inside the bulk which comes periodically with different orientation polarization is still present but this time in the region 0-50 μm. The period of oscillations is also multiplied by 2.

The space charges found in this configuration is rather complex and may have several origins: i) orientation polarisation of PEN, which is a polar material [START_REF] Hoang | Electrical characterization of PEN films using TSDC and PEA measurements[END_REF]; ii) charges migrating into LDPE and being blocked at the interface with PEN; iii) charges injected in the PEN layer; iv) interface charge build-up due to the difference in permittivity of the two materials. Finally, because of mismatch of acoustic properties of the materials, multiple reflection in the PEN layers may produce apparent acoustic signal in the bulk. The periodic nature of the charge found in the buk is certainly associated with this later process. The period of oscillations is about twice for the thicker PEN because the waves need more time to be reflected. Charge injection from PEN is not obvious since it would provide homocharge rather than heterocharge. An interfacial charge (Maxwell-Wagner-Sillars effect) [START_REF] Vu | Correlating conductivity and space charge measurements in multi-dielectrics under various electrical and thermal stresses[END_REF] is certainly present, but it does not explain all the feaures: the charge seems clearly stabilized at the LDPE/PEN interface and needs strong reversed field to be modified. A similar process seems to be at play at the high voltage (top) electrode.

c) LDPE soaked in cumyl alcohol with PEN layers.

The space charge pattern for soaked LDPE with 12 µm thick PEN layers is shown in Fig. 4.c. The measurement takes 15 hours time to be completed. The sample was weighed before soaked in cumyl alcohol, after soaked and after doing space charge measurement. After 15 hours measurement, there is still about 3.7% cumyl alcohol inside the sample, vs. 5.7% initially, which is consistent with Fig. 2b.

Comparatively to the case of Fig 4 .b there is clearly much more charges. Already after the first step at 5 kV/mm, heterocharges settle near both electrodes. For the highest fields applied, there are two prominent peaks at each electrodes and it is really difficult to localize the dielectric/dielectric interface. In the 3 rd to 5 th voltage steps, there seems to be negative charge injection combined to heterocharge build up. The bulk of the insulation shows periodic charge patterns, dominated by positive charges in the first half of the cycle and by negative charges in the end of the second half: charges seem to be a repetition of the heterocharge.

C. Conduction current measurement

For conduction current measurement, new samples with the same thickness (340 μm) and the same diameter are used. The sample configuration for conduction current measurement is the same as for space charge measurement (PEN layers with 2 cm diameter gold electrodes). The temperature for this measurement is 25°C with the applied electrical field from 5 kV/mm to 40 kV/mm. Each measurement took 30 minutes for polarisation and for depolarisation at each of the 5 field steps. The current was acquired every 2 seconds all over the cycle. The results from current measurements are shown in Fig. 5. From Fig. 5a, it is obvious that the sample soaked in cumyl alcohol has much more current than the reference LDPE. The current difference is almost two decades for each applied field. For 40 kV/mm, the current in the sample seems to increase along with time in the two cases. This phenomenon has already been observed in LDPE and XLPE especially when stressing under high field [START_REF] Teyssedre | Charge distribution and electroluminescence in cross-linked polyethylene under dc field[END_REF]. Fig. 5.b shows the apparent conductivity vs. electrical field in the sample. Apparent means that we do not take into account the field distribution in the different layers; the field is supposed homogeneous. The difference in apparent conductivity between soaked and unsoaked sample is by two orders of magnitude. This obviously affects the electrical field distribution under DC voltage. From Figure 5.b, gross value of the conductivity of LDPE without cumyl alcohol is 10 -15 S/m. For the sample with cumyl alcohol it is 10 -13 S/m as obtained for a field of 30 kV/mm. Considering these values and value of the conductivity for PEN of 2.6 x 10 -17 S/m [START_REF] Hoang | Electrical characterization of PEN films using TSDC and PEA measurements[END_REF], it can be anticipated that the field is contained essentially in the PEN layers, and is very high (300 kV/mm or more) even for the unsoaked LDPE. This is consistent with the large interface charge found at the PEN/LDPE interface. The use of a less resistive layer material than PEN is necessary to investigate charging effects in the bulk.

IV. CONCLUSION

The PEN layers tend to prevent charge injection in LDPE, but produces significant interfacial charges between dielectric layers due to differences in conductivity. Periodic charge peaks appearing inside the bulk is assumed as originating from multiple reflections of acoustic waves in the PEN layer due to acoustic impedance mismatch with LDPE. PEN is preventing the cumyl alcohol to evaporate and permits to make measurements for several hours with keeping a reasonable amount of product. Cumyl alcohol has an impact in terms of space charge build up. The presence of cumyl alcohol in the sample creates a strong heterocharge build-up near both electrodes. The LDPE conductivity is substantially enhanced when introducing cumyl alcohol, which explains the increase of interfacial charge between PEN and LDPE. Space charge results show that this charge tends to be persistent and high field is necessary to reverse it: it means that ions migrate and/or deep traps are present at the interface. Ions generated by the cumyl alcohol into the sample could produce charge accumulation and make the sample more conductive. Because PEN is a high resistivity material, the field tends to be concentrated in this layer and the configuration is not really suited to investigate field induced dissociation in the bulk material.

Fig. 2 .

 2 Cumyl alcohol absorption rate (a) and evaporation rate (b) vs.

Fig. 3 .Fig. 1 .

 31 Fig. 3. Space charge measurement protocol

  A. Protocol Space charge measurements have been performed in different configurations of sample such as:  Untreated LDPE only  Untreated LDPE with 12µm PEN layers  Untreated LDPE with 25µm PEN layers  LDPE soaked in cumyl alcohol for 4h with 12µm PEN layers

Fig. 4 .

 4 Fig. 4. Space charge distributions obtained in the different samples. The color bar represents charge density scale in C/m³.

Fig. 5 .

 5 Fig. 5. Current measurements for cumyl alcohol soaked (S) and unsoaked (US) LDPE. Both are with 12 µm thick PEN layers.