Gregory Font

Delphine Sinoquet

Hoël Langouët

Michel Castagne

Sébastien

Grégory Font

Michel Castagné

Sébastien Magand

Derivative free optimization method and physical simulations coupled with statistical models for transient engine calibration

Introduction

Future regulations enforcement concerning emissions, customers requirements and economic competition lead to higher constraints on internal combustion engines development. It paves the way for increasing complexity in terms of design and control. Diesel engines, with many sophisticated devices for air loop as well as for injection system, have now potentially thousands of control parameters to be calibrated. Fast and reliable tuning of these parameters are key issues for engine manufacturers. The first stage of the calibration process consists in basic engine maps tuning for steady-state conditions. The conventional method for engine map optimization is based on statistical engine models depending on control parameters obtained from optimally designed experiments at test bench in steady-state conditions. The cumulative emissions on the driving cycle computed with these models are optimized ( [START_REF] Röpke | DoE -Design of Experiments, Methods and applications in engine development[END_REF], [START_REF] Castagné | Engine calibration : towards an integrated approach[END_REF]). In this approach, one assumes a quasi static engine behaviour on engine speed and load trajectory, indeed the transient effects during the driving cycle are not taken into account.

Engine physical model coupled with statistical and control models is a promising solution to deal with the transient effects optimisation. Eventually the real air and fuel trajectories during the driving cycle could be used as inputs to optimise the control parameters to fulfil noise and emissions targets. Powertrain system simulation, with phenomenological 0D combustion and air loop models, appears as a good compromise between accuracy and low CPU time for engine physical modelling. One of the main challenges concerns the coupling process in order to limit CPU time and loss of accuracy. Usual optimisation algorithms require derivatives of optimised functions which are not available for such models. Moreover, the large size of the optimization problem (number of parameters for engine map description) prevents from approximating these derivatives by perturbation techniques (finite differences). Hence a derivative free optimization method whose efficiency has been highlighted in several application areas ( [START_REF] Langouët | A derivative free optimization method for reservoir characterization inverse problem[END_REF]), is proposed to tackle this topic. Thus this paper focuses on the calibration process (design of experiments, statistical modelling of engine responses), engine physical models, to couple those models and to design of adequate optimization tools .

Experimental data and statistical modelling 2.1. Description of the test case

The engine used during this study is a 1.6 litre four-cylinder Diesel engine (DV6 ATED 4 from PSA Peugeot-Citroën). This engine has a maximum power output of 66 kW. The bore and the stroke of this engine are respectively of 75 mm and 88.3 mm. The compression ratio is 18:1. It is a turbocharged engine with a high pressure cooled exhaust gas recycling circuit. It has a second generation injection system using Bosch CRI 2.2 injectors with injection pressure as high as 1600 bar. The exhaust line has a conventional Euro4 Diesel oxidation catalyst (DOC) and is not equipped with a Diesel particulate filter (DPF).

Moreover, this study is part of a project concerning the use of biofuels to deal with the issues of energy diversification and greenhouse gas emissions. Therefore an Ethanol based Diesel fuel was chosen; the fuel formulation is deeply described in the reference [START_REF] Magand | Optimization of a Euro 5 Vehicle Powered by an Ethanol Based Diesel Fuel[END_REF]. The vehicle and the powertrain were not modified to use this blend but the engine calibration was updated as explained in previous publications [START_REF] Magand | Optimization of a Euro 5 Vehicle Powered by an Ethanol Based Diesel Fuel[END_REF], [START_REF] Chaudoye | Advanced calibration methods applied with an innovative ethanol/diesel fuel formulation[END_REF].

The most significant point concerning the facilities is the use of a dynamic test bench driven by Morphée for the reproduction of normalised European driving cycle (NEDC). Such equipment and procedure have been detailed in a previous paper [START_REF] Magand | Optimization of a Euro 5 Vehicle Powered by an Ethanol Based Diesel Fuel[END_REF]. The test cell provides all the equipment required for transient tests (fuel consumption, opacimeter for smoke measurement) and the design of the fluid circuits has been optimised to ensure good vehicle behaviour reproduction.

Experimental data

Different approaches to tune engine maps are proposed to deal with the high number of parameters to optimise and to reduce this time-consuming task during the development process. These methods (local, mixed and global) based on mathematical models and Design of Experiments (DoE) were compared in previous publication [START_REF] Castagné | Comparison of engine calibration methods based on DoE[END_REF]. Global modelling approach for calibration, considering load and speed as parameters, reduces time spent for whole calibration workflow because an unique DoE is necessary to model engine responses for entire domain (or several DoE if this domain is divided in several zones with common parameters) [START_REF] Chaudoye | Modelling engine operating space for DoE calibration methods[END_REF]. In fact the global method considers the trace of the driving cycle (load and speed) as inputs of the models of the engine responses. Figure 1 shows the operating area chosen for the calibration in the speed/load area. The next stage consists in defining the operating space of each engine control parameter. Concerning injection parameters, previous experiences [START_REF] Chaudoye | Modelling engine operating space for DoE calibration methods[END_REF] show that fixed boundaries or simple linear limits compared to the original maps could be a suitable answer to ensure a sufficient potential of optimization.

The main issue involves the air loops limits. It is required to have a wide enough operating space close to the physical domain (to ensure a high potential of optimization) and to make it easy to model. Four air limit parameters should be defined: two for the boost pressure (maximum and minimum limits) as well as two for the air mass flow. Three limits are physical ones and could be straightforwardly obtained. To determine the last one (minimum air mass flow), different criteria could be chosen (air fuel ratio, stability, etc.). Then each limit should be modelled, in this case with low degree polynomial. The accuracy of the model is illustrated in Figure 2 for the minimum air flow limit by comparing experiences used for the modelling (blue crosses) and the validation points (red circles). The coefficient of determination R² as well as the root mean squared error (RMSE) show the good estimation obtained with the model. Some of the main challenges when using the global modelling approach concern complex constraints for description of parameter variation domain (e.g. non-linear) and increased number of parameters to manage. Thus dedicated statistical tool and method should be used to take into account these issues. This process was described in previous reference [START_REF] Chaudoye | Advanced calibration methods applied with an innovative ethanol/diesel fuel formulation[END_REF]. After defining the speed and load boundaries as explained previously, the values of these parameters are fixed; this permits to determine acceptable values for injection parameters and then air loop limits depending on other parameters could be given. The most valuable experience is determined by using a Maximin criterion (maximization of minimal distance between experimental points). Finally, the global DoE includes also the corners of the operating space, some mono parametric variations and a space filling.

Global modelling of the engine responses

Modelling the complex engine response surfaces requires to use advanced statistical models instead of standard polynomial models. In this case the choice of kriging for modelling was the most relevant since it is flexible and takes into account the experimental repeatability and the spatial dependency providing therefore consistent values for close operating points. A total of 1200 experiences, added to 100 points for repeatability and 200 points for validation of the models, were performed. The engine responses that are modelled are : HC, NOx, CO, CO 2 and PM emissions (g/h), combustion noise (dB) and fuel consumption (kg/h). The statistical criteria are satisfying for all models except HC and PM emissions. Concerning PM emissions, the oxygenated fuel used for these tests implies near-zero soot emissions and there is a lack of accuracy of the measurement device for these levels. Concerning HC emissions the bad repeatability is the main explanation of the poor quality of the model. Some examples of validation of the estimated value from the model are given in Figure 3 and in

Coupling physical engine simulator with statistical models

IFPEN is strongly involved in powertrain modelling and simulation to support technological development. It develops specific libraries (IFP-Engine, IFP-Drive and IFP-Exhaust) in the LMS.IMAGINE.Lab AMESim environment which includes a large set of commercial libraries used for various applications, in particular for control strategy development, validation and pre-calibration as it is presented in [START_REF] Albrecht | Towards a stronger simulation support for engine control design : a methodological point of view[END_REF].

The purpose of this section is the coupling of a physical engine simulator with statistical models and control algorithm model in order to have both well representative transient effects and accurate pollutants model. The complete engine simulator is based on the LMS.IMAGINE.Lab AMESim platform with in particular the use of the IFP-Engine library. Its description and validation are given in the following sections.

Complete engine simulator creation and validation in steady state

Combustion model

IFP-Barba phenomenological combustion model

In this study, a combustion model based on the Barba approach [START_REF] Lafossas | Development and Application of a 0D D.I. Diesel combustion model for Emissions Prediction[END_REF] [10] [START_REF] Lebas | A Phenomenological Approach to Model Diesel Engine Combustion and In-Cylinder Pollutant Emissions Adapted to Control Strategy[END_REF] is chosen to compute the combustion heat release rate. In this approach, the combustion process is divided into two parts. In a first step, the fuel is burnt using a pre-mixed model with the hypothesis of flame propagation in the pre-mixed zone. When the pre-mixed zone is burnt, the remaining fuel is oxidized using a mixing controlled combustion model. Using a simple example with a pilot and a main injection, different steps can be isolated during the combustion. At the beginning of the pilot injection, the combustion model is initialised and the different variables are computed. During injection, the auto-ignition delay is computed and the injected fuel is introduced in the premixed zone.

The auto-ignition delay is defined at the beginning of the injection and is computed using a simple Arrhenius law. When the auto-ignition delay is reached at the end of the phase, the pre-mixed combustion starts. In a first step, the ascendant part of the combustion heat release rate (burning mode 1) is generated by the flame propagation in a turbulent field. Then, because of the multiplication of the auto-ignition sites, the hypothesis of the propagation of a single flame in the pre-mixed zone is no longer valid and the flame interactions must be modelled (burning mode 2), leading to a reduction of the total flame surface. When the main injection starts, the model creates a new pre-mixed zone for the main injection and computes the auto-ignition delay. The definition of the pre-mixed zone is similar to that of the pilot injection.

As for the first injection, the auto-ignition delay of the main injection is computed from the beginning of the injection, leading to the beginning of the combustion. The autoignition model is the same as for the pilot injection, but thanks to a higher pressure and temperature, the auto-ignition delay is smaller and reached before the end of the injection.

As the auto-ignition delay is reached before the end of the injection, the fraction of fuel available in the premixed zone is lower than for the pilot injection. The combustion process is the same as for the first injection.

In parallel to the pre-mixed combustion mode, the remaining injected fuel starts to burn progressively with a mixing controlled combustion model. In this burning mode, the combustion is piloted by the mixture speed between the injected fuel and the surrounding air. Generally, the mass of fuel burnt in this mode is defined with a simple equation depending on the turbulent kinetic energy. In order to enlarge the range of applications of this combustion model, there is the possibility to compute more injections and combustions.

Combustion model set up

A single cylinder simulator with closed valves is used to calibrate the combustion chamber. The calibration is made between Intake Valves Closure (IVC) to Exhaust Valve Opening (EVO). The initial in-cylinder thermodynamic conditions are defined using the experimental results for the steady state points. The advantage of this approach is to lower the error generated by the air path on the in-cylinder conditions at the IVC.

Concerning the thermal transfers, a Woschni model is used [START_REF] Woschni | Universally Applicable Equation for the Instantaneous Heat Transfer Coefficient in the Internal Combustion Engine[END_REF]. In order to fit correctly the combustion, the wall temperatures are functions of the BMEP with the following law :

19 100 450 ) ( * 19 100 ) ( - + = bar BMEP K T wall (1)
To set up correctly this combustion model, 9 points are chosen among the 218 available in the engine map (cf. Figure 5 -red surrounded crosses). The calibration process consists in fitting the simulated cylinder pressure on the experimental one for each operating condition. Then a validation on the whole operating points is done. Only one set of parameters is used for all the engine map.

Figure 5: Engine map. 9 points (red circles) are chosen among the 218 available to calibrate the combustion model

The cylinder pressure curves shown in this paper (cf. Figure 6) are three out of the nine points used for the calibration of the combustion model (cf. Figure 5 -red surrounded crosses) . The four experimental cylinder pressures are in blue and the mean of the four is in dashed black. The simulated cylinder pressure is in red. In this paper, all the comparisons are made between the mean experimental cylinder pressure and the simulated one. The cylinder pressures of the model are visually closed to the cylinder pressures of the bench. Moreover as described in the Figure 7 the percentage of points matching the criteria are high. It must be reminded that these results are obtained with only one set of parameters for the whole engine map.

Steady state engine model

Figure 8 presents the complete engine simulator. The whole air path with the turbocharger, the charge air cooler, the exhaust line and the EGR loop is modelled. After its tuning the combustion model is embedded in this complete engine simulator.

Figure 8 : Four-cylinder simulator in AMESim environment

In order to validate the behaviour of the complete engine simulator, it is coupled with a simple control implemented in Matlab Simulink. Two PID controllers regulate the engine : one for the intake manifold pressure with the control of the VNT position and one for the air mass flow rate with the control of the EGR valve position.

In order to appreciate the quality of the results, some criteria are defined :

-Intake manifold pressure (P2S) : +/-0.05 bar -Intake manifold temperature (T2S) : +/-10 K -Exhaust manifold pressure (P3) : +/-0. Figure 10 gives the percentage of the points which respect the criteria on IMEP, air mass flow rate and intake and exhaust manifold temperatures and pressures. According to these criteria, results are consistent with test bench data and the complete engine simulator can be considered as set up.

Coupling with statistical models : presentation and validation on NEDC cycle

Figure 11 presents a way to estimate pollutant emissions during a transient thanks to a coupling between the AMESim simulator described in previously section, the engine control and the statistical models described in section 2.3. Inputs of these statistical models are air mass flow rate, intake manifold pressure, start of main injection, rail pressure, injected mass in the pilot injection, gap between pilot and main injection, BMEP and engine speed. The coupling is made with Matlab Simulink.

In this study, the real control software was not available. So, a basic open-loop control based on maps is used. Thus maps function of engine speed and total injected mass are used for the pilot and main start of injection, injected mass for the pilot injection, rail pressure, VNT and EGR valve positions.

A NEDC cycle was performed on this platform. As the optimization described in the next section is run on the EUDC part, only the results concerning EUDC are shown. Figure 12 presents comparison between virtual engine (in red line) and experimental data (in blue dashed line) for intake and exhaust manifold pressures (P2S and P3), air mass flow rate (Air_mfr), fuel air equivalence ratio (FAR_sonde), turbo speed (TurboSpeed), BMEP and CO 2 , CO, NO x and HC mass flow rates (CO2GH, COGH, NOxGH and HCGH). A comparison on cumulative emissions (NOx, CO 2 , CO and HC) is also presented in Table 1.

Table 1 Comparison of the cumulative emissions between model and experimental data on EUDC

Keeping in mind that engine control is different for experimental tests and simulation, one can conclude that the behaviour of the engine during the EUDC is well reproduced.

This virtual engine is sufficiently predictive to be used for an optimization process.

Engine map optimization with coupled models 4.1. Principle of engine map optimization

The engine map optimization consists in directly optimizing the cumulative pollutant emissions over the cycle via distortions of the engine maps instead of optimizing individually some selected operating points and building afterwards the engine maps by a smoothing step ( [START_REF] Castagné | Comparison of engine calibration methods based on DoE[END_REF], [START_REF] Sinoquet | Multi-objective constrained optimization of engine maps, 5th IAV conference: Design of Experiments (DoE) in Engine Development[END_REF]). Consequently, the objectives to be minimized (or constrained) are the engine responses cumulated on the considered driving cycle. The engine responses can be modelled by a global statistical model (as described in section 2) or by a physical simulator that models transient behaviour of the engine (as described in section 3).

The map optimization problem is formulated as

≤ ≤ ≤ ∈ j P P T P j p P P T P i R m S dt t c t r m t c t r m t c t r m t c t r F c r u c r Am c r l to subject dt t c t r m t c t r m t c t r m t c t r F Np Np p N p 0 0 (2)
where : o (r(t); c(t)) indicate the trace of the considered driving cycle within the engine speed-load domain, o F i is the global model of the engine response i depending on the engine control parameters but also on the speed and load, o m

Pi are the 2D engine maps of the control parameters in the engine (speed, load) operating domain.

Additional smoothing constraints such as global smoothing constraints (to preserve the regularity of the original maps), robustness constraints (to take into account parameter dispersions) and local constraints (for example limits on the gradients of the maps, or limits on engine noise response) are also introduced. This formulation requires an adapted parameterisation of the engine maps, m Pi (r,c), in order to limit the total number of unknowns in the optimization process. For this application, an adapted discretization of the engine speed/engine load domain with a bilinear approximation in cells is implemented (see Figure 13):

rc c r c r m i P 3 2 1 0 ) , ( α α α α + + + =
within a given cell with

Derivative Free Optimization (DFO)

Problem ( 2) is a nonlinear constrained optimization problem which requires dedicated optimization methods to be solved. For engine responses modelled by statistical models, a classical SQP (Sequential Quadratic Programming) method is applied: first derivatives of objective function and constraints are available. When the physical model described in section 3 is used, a Derivative Free Optimization method is implemented. In this case, no derivative of modelled engine responses is available and moreover simulations require a high computational effort (a couple of hours of calculation for one simulation). Therefore, a dedicated optimization method is proposed, the SQA method (Sequential Quadratic Approximation), already applied successfully on several industrial applications [START_REF] Langouët | A derivative free optimization method for reservoir characterization inverse problem[END_REF]. This method is an extension of the DFO method proposed by Powell in [START_REF] Powell | The NEWUOA software for unconstrained optimization without derivatives[END_REF] to nonlinear constrained problems: in order to save simulations, quadratic interpolating models are used as surrogate of the simulator responses. These quadratic models are minimized and updated thanks to additional simulations performed along the optimization process iterations. The main steps of the SQA algorithm are described in Table 2.

SQA Algorithm

Step

Application on EUDC driving cycle

First, a quasi-static optimization (QSO) is applied. It is based on cumulative engine responses on the EUDC cycle (1 Hz discretization) modelled by global statistical models presented in section 2.3 (conventional optimization method is used, derivatives of objective and constraint functions being available for surface response models). Four engine maps controlling injection are optimized: main injection timing, fuel injection pressure, pilot fuel injection quantity, pilot injection timing whereas mass air flow and boost pressure maps remain fixed to initial values. Bilinear modelling of maps allows to limit the number of map parameters to be optimized: 80 parameters. The minimized objective is the cumulative NOx emissions along the cycle, CO2 emissions are constrained to remain smaller or equal to the initial value, whereas an increase of 50% for HC and CO emissions is allowed. The optimization process is then applied but this time with physical AMESim model coupled with pollutant statistical models described in section 3.2 with DFO optimization method presented in Table 2. In order to save simulation computational time, the smallest number of interpolation points for a quadratic model: m=n+2 is chosen for this application. This incomplete quadratic model is completed during optimization process with a constraint on model variations along iterations. Additional terms of the model are then determined. Therefore, as displayed in Figure 14, the initialization phase costs 82 simulations used to build the initial quadratic model. After this phase, the optimization work starts with a rapid decreasing of NOx emissions.

Results are presented in Table 3 and compared to results obtained with quasi-static optimization. The optimization coupled with the physical AMESim model, compared to the quasi-static optimization, enables to drop off the NOx emissions with a slight increase of HC and CO emissions. But they are still largely under the defined thresholds and even under initial maps results. CO 2 emissions remain close to the initial value. 2) to build the initial quadratic model.

Physical coupled models

NOx emissions CO2 emissions CO emissions HC emissions

Conclusions and outlook

Finally an innovative method has been presented to deal with transient effects on engine map calibration process on a Diesel engine, while performing only tests in steady-state conditions. It completes or replaces the standard quasi-static approach.

The main features of the method consist in the different stages listed below:

-Obtain statistical models taking into account engine speed and load as parameters by applying advanced methodologies for domain definition, design of experiments and engine response modelling; -Tune phenomenological 0D combustion and air loop models on a limited data set, then validate on the engine operating domain; -Couple of statistical and physical engine models and engine control model to manage the effects of transient operations (especially regarding air loop parameters) on emissions; -Adapt the formulation of the map optimization and design a derivative-free optimization algorithm to optimize engine responses on a driving cycle. This methodology has been successfully applied to minimize NOx cumulative emissions on a EUDC cycle of a 1.6 litre engine with biofuel while maintaining initial CO 2 emissions. Nevertheless the results show that the potential of emission reduction remains low on EUDC cycle, because of the smoothness of such a cycle. Transient operations should become a key issue in engine calibration process since driving cycle might change from the NEDC to the World Light Duty Procedure (WLTP). This method could also be helpful to optimize the engine control transient corrections (for example correction of injection parameters to compensate turbo lag). One of key issues for system simulation emergence concerning optimization on driving cycles is to have both reduced CPU time and accurate models. Therefore multimodel multi-core simulation platforms ( [START_REF] Ben Gaid | Heterogeneous Model Integration and Virtual Experimentation Using xMOD: Application to Hybrid Powertrain Design and Validation[END_REF]) which allow to reach real-time computation for a virtual engine in keeping the same model accuracy will be investigated in further work.

Figure 1 :

 1 Figure 1: Engine speed and load area for global method

Figure 2 :

 2 Figure 2: Comparison of estimated and measured minimum air mass flow limit. Blue crosses correspond to experimental data used for modelling whereas red circles are validation points.

Figure 4 .

 4 First NOx emissions model is validated with experimental specific tests. CO 2 emissions estimation during monoparametric variations are repre-RMSE = 11 mg/cp R² = 0.989 sentative of test bench results in Figure 4. In fact, for the different emissions, the models show good trends and consistent estimated levels.

Figure 3 :

 3 Figure 3: NOx emissions model validation. Blue crosses correspond to experimental data used for modelling whereas red circles are validation points.

Figure 4 :

 4 Figure 4: Comparison of estimated and measured CO 2 emissions on an engine speed variation

Figure 6 :Figure 7 :

 67 Figure 6 : Cylinder pressure -2500 RPM -BMEP = 1, 10 and 16.9 bar. The four experimental cylinder pressures are in blue and the mean of the four is in dashed black. The simulated cylinder pressure is in red.

Figure 9 :Figure 10 :

 910 Figure 9: Steady state engine simulator results. Left figures : blue crosses represent experimental data and red squares modelled data results for each operating point. Right figures : red lines represent criteria to respect.

Figure 11 :Figure 12 :

 1112 Figure 11 : Description of the coupling between engine Control, engine AMESim model and emission statistical models

  corners. Other map representations based on LoLiMoT (Local Linear Model Tree) are proposed in[START_REF] Sinoquet | Multi-objective constrained optimization of engine maps, 5th IAV conference: Design of Experiments (DoE) in Engine Development[END_REF].

Figure 13 :

 13 Figure 13 : Map discretization (for bilinear parameterization) in the engine (speed, load) operating domain and EUDC cycle trace.

Figure 15

 15 compares engine emission traces on EUDC cycle for initial engine maps, QSO maps and maps obtained by optimization with coupled physical models. On 70-100km/h acceleration of EUDC cycle, one notices the gain of optimization considering transient effects compared to quasi-static optimization for NOx reduction. Nevertheless, globally, the NOx reduction is low due to the smoothness of EUDC cycle. Larger gains are expected on more realistic cycles.The physical AMESim model allows to take into account the transient effects of the air loop. It is then possible to take care of boost pressure lag for example to optimize the maps and not only to work on steady-state conditions.

Figure 14 :

 14 Figure 14 : Objective function decreasing and constraint satisfactions along SQA iterations (Derivative Free Optimization method). Dotted lines indicate upper constraint thresholds. The 82 first simulations belong to the SQA initialization phase (see algorithm in Table2) to build the initial quadratic model.

Figure 15 :

 15 Figure 15 : Pollutant emissions (physical AMESim models coupled with statistical models) along EUDC cycle (1Hz discretization) for initial, Quasi-Static optimized maps, optimized maps (with physical models). (Bottom) Zoom on the 70-100km/h acceleration of EUDC cycle.

Table 2

 2 0. INITIALIZATION A quadratic model (QM) of each optimized engine response is constructed in a neighbourhood of the initial guess, namely the trust region. They interpo-Main steps of SQA algorithm: a DFO (Derivative Free Optimization) method for nonlinear constrained optimization.

	late m points (n+2 ≤ m ≤ (n+1)(n+2)/2 with n the number of
	parameters), obtained by evaluation of the engine simula-
	tor.
	Step 1. MINIMIZATION At each iteration, we substitute in prob-
	lem (2) engine responses with their associated QM and solve
	it within the trust region (TR).
	Step 2. EVALUATION and UPDATE of TR We evaluate the engine re-
	sponses at the solution of previous problem and update QM
	with this new point. The size of TR is also updated accord-
	ing to the predictivity of QM at this point :
	i.	if predictions by QM are very closed to simulation
		results, TR is enlarged. Go to step 1.
	ii.	if only the trend is correct, TR is not modified. Go
		to step 1.
	iii.	else reduce TR and go to 3.
	Step 3. IMPROVEMENT OF QM Choose and simulate a new point
	(evaluation of engine responses) in order to improve the
	accuracy of QM, update QM and go to 1.

Table 3

 3 Map optimization results with physical AMESim models coupled with statistical models

1 2 1

Acknowledgements

The authors kindly acknowledge Fabien Chaudoye for his significant contribution to this paper.

))) (


Yann Poirette

Martin Guiton

Guillaume Huwart

Delphine Sinoquet

Jean Marc Leroy

Jean Leroy An Optimization

An Optimization method for the Configuration of Inter Array Cables for Floating Offshore

INTRODUCTION

The first commercial offshore wind farms appeared in the early nineties giving then new opportunity for the development of renewable energy sources. Today, fixed offshore wind turbines largely contribute to renewable energy produced from ocean (in 2015, 12 GW of offshore wind capacity [START_REF]Annual Report[END_REF] vs 0.53 GW for other renewable marine energy sources [3]). Nevertheless, the deployment of fixed bottom structures for offshore wind is limited in water depth mainly due to economic reasons. The floating offshore wind turbines are currently developed to extend this water depth limitation. The first floating offshore wind turbines have been connected to grid in 2009 (Hywind, in Norway) and 2011 (WindFloat, in Portugal), see [START_REF]Deep Water. The next step for offshore wind energy[END_REF]. The first tender for floating offshore wind pilots has been launched in 2015 in France [START_REF]Appel à projets -Fermes pilotes éoliennes flottantes[END_REF]. Nevertheless, the costs of such projects remain still high and need to be optimized, to accelerate the deployment of floating wind turbines. In the cost breakdown of an offshore wind farm, the inter array cable cost accounts for about 10 to 15% and represents a significant amount of the capital cost [START_REF] Blanco | The economics of wind energy[END_REF]. This motivates the optimization of the configuration of inter array cables.

Automated optimization remains still new in the design process of offshore structures. Nevertheless, the global energy economic context and the emergence of new technologies for renewable marine energy sources lead the engineers to develop and use optimization tools, see for example the pioneering work of [START_REF] Clauss | Hydrodynamic Shape Optimization of Large Offshore Structures[END_REF] on the hydrodynamic shape optimization of large offshore floating structures for oil and gas production, and more recently [START_REF] Martins | Compliant vertical access riser assessment: DOE analysis and dynamic response optimization[END_REF] and [START_REF] Yang | Multi-objective optimization for deepwater dynamic umbilical installation analysis[END_REF] on the optimization of oil and gas risers and umbilicals.

In their review about design optimization of wind turbine structures, Muskulus et al [START_REF] Muskulus | Analytical gradientbased optimization of offshore wind turbine substructures under fatigue and extreme loads[END_REF] classify first the papers according to the kind of analysis: static analysis, frequency domain analysis or time domain analysis, analysis which integrate increasing complexity in the model. Most of the studies they referenced aim to save structure weight. They also mentioned the specificity of floating wind turbine optimization with a more non-linear dynamic problem. At last, they deliver some recommendations for optimization studies. In particular, they underline the crucial importance to correctly estimate the gradient when using a gradient based optimization method, recommending to use as much as possible analytical approaches (e.g. [START_REF] Kennedy | An Adjoint-based Derivative Method for Time-Dependent Aeroelastic Optimization of Flexible Aircraft[END_REF] for an aeroelastic optimization of flexible aircraft with an adjoint approach). However, Muskulus et al. also explain in the same review paper that users of the dedicated aero-hydrodynamic software generally do not have access to the source code or to a proper documentation required for analytical gradient estimation.

In their recent works, Chew et al. [START_REF] Chew | Analytical gradient-based optimization of offshore wind turbine substructures under fatigue and extreme loads[END_REF] propose an optimization methodology of fixed foundation for offshore wind turbine with analytical gradients. The simulations are performed in the time domain using a linear hydro-elastic solver with aero-servo loads calculated by a specific model. The constraints include sizing, eigenfrequency, extreme and fatigue loads. They highlight strong limitations in the Finite Difference gradient approximation when compared to analytical gradients. Especially, the choice of the perturbation step size of Finite Differences is critical. They also indicate that the decoupled aerodynamic and hydrodynamic approach could be insufficiently accurate for more flexible non-linear structures.

Brommundt et al. [START_REF] Brommundt | Mooring system optimization for floating wind turbines using frequency domain analysis[END_REF] use the simulation consuming Nelder-Mead simplex algorithm for a mooring system optimization of floating wind turbines, with constraints approximately represented by penalty functions. Simulations are performed in frequency domain. The authors underline the limitations of such calculations for the case study of semisubmersible floater that they investigate.

In [START_REF] Fylling | WINDOPT-An optimization tool for floating support structures for deep water wind turbines[END_REF], the authors propose to use an integrated optimization tool to minimize the cost of the floating support, mooring system and power cable for offshore wind turbine. The computations are performed in the frequency domain. The specific optimization tool is extended from a project for optimization of mooring and riser in deep water and uses a Sequential Quadratic Programming implementation (NLPQL, [START_REF] Schittkowski | NLPQL: A fortran subroutine solving constrained nonlinear programming problems[END_REF]), based on a gradient search method. A limitation of this optimization application is, as already pointed out, the use of Finite Difference approximation of derivatives.

Other studies are dedicated to optimization of the inter array layout design of the wind farm in order to maximize the energy yield and minimize the wake losses. The earliest studies on wind farm optimization used genetic algorithms [START_REF] Grady | Placement of wind turbines using genetic algorithms[END_REF] or evolutionary strategy algorithms [START_REF] Kusiak | Optimization of wind turbine energy and power factor with an evolutionary computation algorithm[END_REF] which are expensive in terms of simulations.

In the present paper, we propose to implement a methodology to optimize the configuration of the inter array cable of a floating offshore wind farm based on the derivativefree Sequential Quadratic Approximation (SQA) algorithm developed by IFPEN, which allows to integrate non-linear constraints ( [START_REF] Langouët | Optimisation sans dérivées sous contraintes : deux applications industrielles en ingénierie de réservoir et en calibration des moteurs[END_REF], [START_REF] Font | Derivative free optimization method and physical simulations coupled with statistical models for transient engine calibration[END_REF], [START_REF] Langouët | A derivative free optimization method for reservoir characterization inverse problem[END_REF]). This algorithm relies on a local quadratic interpolation within a trust region which is updated during the optimization ( [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF], [START_REF] Powell | The NEWUOA software for unconstrained optimization without derivatives[END_REF]). It provides an efficient way to avoid the time and error limitations of Finite Difference gradient computations, or for the general case, when the aerohydrodynamic simulation does not provide any gradient information. The constraints are in that case labeled as "blackbox". Furthermore, we propose a methodology for choosing the relevant set of parameters to be optimized via a sensitivity analysis and for the choice of the initial points of a multi-start optimization approach. Finally, Finite Element Method (FEM) computations in time domain performed in this paper allow to achieve a high accuracy even for highly non-linear extreme loading cases. The potential of this methodology is illustrated with a case study of a lazy-wave power cable configuration. We consider 28 Ultimate Limit State loading cases from a Mediterranean environment, and an IFPEN floater design for a 3.6 MW wind turbine.

CASE STUDY

The power cable considered in this study is a classical three-core inter array cable of offshore wind farm in a lazy wave configuration, shown in Figure 1. To define the extreme load cases on the power cable, a floater derived from the concept described in [START_REF] Poirette | A new ballasted floating support for offshore wind turbine[END_REF] has been used. The environmental conditions considered in this study correspond to typical mean values of various metocean data from the Mediterranean sea. Two extreme design load cases were then considered, depending on whether the wind turbine is in production or not (Table 1). Note that the parked cases associating 50 years marginal return period for wave, current and wind, are particularly severe. From these data, 28 load cases have been defined based on wave, current or wind heading combinations.

The dynamic behavior of the floater is calculated using Deeplines TM [START_REF] Heurtier | Fully coupled dynamic analysis of rigid lines and floaters behaviours in deep water[END_REF]. This software is dedicated to offshore structure analysis with FEM and it is co-developed by IFPEN and Principia. It computes hydrodynamic and aerodynamic loads applied on fixed or moored offshore structures, and internal stress in the structure or mooring lines. The FEM simulation first computes a static state, and then a dynamic evolution around it.

The wind turbine is modeled as a mass point where mean aerodynamic loads depending on wind speed are applied. The tower is modeled using beam elements. The floater is modeled as a rigid body with hydrodynamic database calculated using Diodore TM [START_REF] Berhault | Diodore: a Numerical Tool for Frequency and Time Domain Analysis of the Behaviour of Moored or Towed Floating Structures[END_REF].

The electrical cable is modeled using beam elements with finite rotations and small strain kinematics.

The environment is modeled as follow:

 the waves are defined as JONSWAP spectrum using Hs, Tp shown in Table 1;  the current, if any, is defined by a constant piecewise linear speed along the depth,  the modeling of the wind effect depends on the load case considered: in production case, the wind is not modeled but the thrust and torque are calculated depending on constant wind speed considered in Table 1 and applied at the hub level. In parked conditions the wind is defined by a constant vertical profile. Note that the current and wind loads are taken into account for the static state calculation, while the wave loads are computed during the dynamic simulation.

Considering the low mass of the power cable compared to the floating support, no interaction is considered between the floater behavior and the cable behavior. Then, the movements of the floater have been calculated beforehand for each of the 28 loading cases in time domain calculation. The time series of displacement for the six degrees of freedom are applied on the top of every studied power cable.

The set of parameters related to the power cable for the configuration optimization process are the following:

 L 1c : length of the head part of the cable before the buoyancy modules,  L b : length of the buoyancy part of the cable,  L c : total cable length,  D c : cable diameter,  D b : buoyancy modules diameter,  L s : length of the stiffener  D sb : diameter of the base of the stiffener.

LARGE SCALE MECHANICAL PROPERTIES OF THE POWER CABLE

A classical cross section of the inter-array power cable is shown in Figure 2. The geometry of the internal core which contains in particular three insulated conductors will be assumed to be fixed in this study, just as the external sheath diameter. The cable structure design focuses here on the wire armor dimensions. The function of the steel armors is to withstand the mechanical tension and to limit the curvature due to the submarine movements imposed by waves and currents. . Each wire is made of steel with elasticity modulus E = 210 GPa and shear modulus G =80 GPa. The average number of wires for the two armoring layers is denoted 𝑛; the radius of an armor is r, the lay angle of each armor (i.e. angle between the axis and the armoring wire) is  fixed here to 10° and the cable diameter is denoted D c . For each value of D c during the optimization process, r and 𝑎 ̅ are then updated according to linear relationships in D c .

After classical mechanical considerations, the axial stiffness EA and the flexural stiffness EI of the cable is written as:

(𝐸𝐴) = 𝑛 ̅ ⋅ cos 3 (𝛼) ⋅ 𝜋𝑟 2 ⋅ 𝐸 + (𝐸𝐴) 𝑐𝑜𝑟𝑒 + (𝐸𝐴) 𝑠ℎ𝑒𝑎𝑡ℎ Eq. 1 (𝐸𝐼) = 𝑛 ̅ ⋅ 2 ⋅ 𝑐𝑜𝑠(𝛼)

SMALL SCALE STRESSES ON THE ARMOR WIRES

In the following, a simplified analysis, based on the helix geometry of armors and assuming a constant radius is used to assess the stresses in the armors. A more realistic approach could be used, such as the numerical model described in [START_REF] Leroy | Assessing mechanical stresses in dynamic power cables for floating offshore wind farms[END_REF].

Considering that armor layers withstand the axial force, the mean axial stress within an armor wire is defined as:

𝜎 ̅ = 𝐹 𝑡𝑜𝑡𝑎𝑙 𝑛 ̅ cos(𝛼) 𝐴 Eq. 4
where 𝐹 𝑡𝑜𝑡𝑎𝑙 is the resulting axial force for all the wires and A the wire cross section area.

The Figure 4 represents an elementary part of an armor wire with the tangent and normal unit vector denoted 𝑡 ⃗ and 𝑛 ⃗⃗, respectively. The tangential forces applied at the left and right boundaries are equilibrated by a linear friction force f F n , with f=0.15 the friction coefficient and F n the linear contact force.

Reminding the normal curvature of an helix to be: The equilibrium projected on 𝑡 ⃗ provides the Coulomb criterion:

𝐴|𝑑𝜎| ≤ 𝑓𝐹 𝑛 𝑑𝑠.

Eq. 7

Let us define a circumferential angle  along the armor helix during the bending set to 0° at the outer-arc and 180° at the inner arc. When the sliding is not triggered, axial stress reads:

𝜎 𝑠𝑡𝑖𝑐𝑘 = 𝐸 𝑎 𝑅 𝑐𝑜𝑠 2 𝛼 cos 𝜃.
Eq. 8

Differentiating Eq.8, Eq.7 becomes:

𝐸 1 𝑅 𝑐𝑜𝑠 2 𝛼 sin 𝜃 sin 𝛼 ≤ 𝑓𝜎 𝑠𝑖𝑛 2 𝛼 𝑎 .
Eq. 9

The sliding is triggered first at the neutral fibers (𝜃 = 𝜋 2

). Integrating the Coulomb criterion, one obtains the increment of axial stress due to bending:

∆𝜎(𝜃) = 𝑓 𝐴 𝐹 𝑡𝑜𝑡𝑎𝑙 𝑛 ̅ cos (𝛼) sin(𝛼) . ( 𝜋 2 -𝜃),
Eq. 10 so that the total axial stress is finally computed as:

𝜎(𝜃) = 𝜎 ̅ [1 + 𝑓 sin(𝛼) . ( 𝜋 2 -𝜃)] 𝑓𝑜𝑟 𝜃 ∈ [0; 𝜋]. Eq. 11

POWER CABLE CONFIGURATION, OPTIMIZATION PROBLEM

The admissible domains of the seven configuration design parameters introduced in the previous section are defined as: ) ⋅ 𝐿 𝑠 ⋅ 𝐶 𝑣𝑠 Eq. 12 with 𝐶 𝑣𝑐 , 𝐶 𝑣𝑏 and 𝐶 𝑣𝑠 the volumetric costs of the cable, of the buoyancy modules and of the stiffener respectively and 𝐷 𝑠𝑒 = 𝐷 𝑐 + 10 𝑚𝑚, the diameter at the end of the stiffener. The first, second and third terms represent the cost of the armors, the buoyancy modules and the stiffener respectively.

 L 1c in [60m, 90m]  L b in [30m, 60m]  L c in [520m, 550m]  D c in [120mm, 135mm]  D b in [2 D c ,4 D c ]  L s in [3m, 6m]  D sb
The particularity of our optimization problem lies in the introduction of constraints to make the design solution acceptable. Indeed, the optimization solution has to respect the following six black-box constraints which are nonlinear with respect to the optimization parameters: The max and min operators refer to the 28 loading cases (Table 1), the curvilinear abscissa along the cable s, and the time t. C1 prevents from any compression in order to avoid a disorganization of the cable components. C2 corresponds to the armor limit state in tension, as 66% of the yield strength 𝑅𝑒=1000 MPa as for the normal safety class in [START_REF]DNV-RP-F401 -Electrical Power Cables in Subsea Applications[END_REF]. To be conservative, we choose the worst case which is at the outerarc. The curvature limit of C3 comes from manufacturer's data. C4 elongation threshold is the same as [START_REF] Nasution | Fatigue analysis of copper conductor for offshore wind turbines by experimental and FE method[END_REF]. Finally, C5 and C6 aim to prevent the cable from touching the sea floor, and the buoyancy modules from emerging from the water, respectively. In the following, these nonlinear constraints are normalized (to be close to [-100,100]) and reformulated so that the feasible domain is ℝ - * .

 C1: max 𝐷𝐿𝐶,𝑠,𝑡 𝜎 > 0  C2: min 𝐷𝐿𝐶,𝑠,𝑡 σ < 2 3 𝑅𝑒 = 660𝑀𝑃a
Copyright © 20xx by ASME Finally, this problem is solved with a workflow chaining different software: the optimization loop is controlled with an IFPEN in house general optimization library and communicates through Python scripts to Deeplines TM FEM solver which is run through multi-threading parallelization.

In this study, the SQA algorithm developed by IFPEN ( [START_REF] Langouët | Optimisation sans dérivées sous contraintes : deux applications industrielles en ingénierie de réservoir et en calibration des moteurs[END_REF], [START_REF] Font | Derivative free optimization method and physical simulations coupled with statistical models for transient engine calibration[END_REF], [START_REF] Langouët | A derivative free optimization method for reservoir characterization inverse problem[END_REF]) is used. This is a local derivative free optimization method based on quadratic interpolation models controlled in a trust region, generalizing the NEWUOA algorithm [START_REF] Powell | The NEWUOA software for unconstrained optimization without derivatives[END_REF] to the case of non-linear black-box constrained problems.

RESULTS

Firstly, a global sensitivity analysis has been conducted to check the influence of the selected design parameters on the simulation outputs thanks to a Latin Hypercube Sampling (LHS) of 600 simulated points. To reduce computational time, only 6 loading cases over the 28 are selected to be particularly severe, all from parked DLC.

The presented results are the Sobol's indices based on the Sobol decomposition of the variance of the outputs with respect to the inputs [START_REF] Saltelli | Sensitivity Analysis in Practice. A Guide to Assessing Scientific Models[END_REF]. The total and the first order indices for each design parameters are displayed in percentiles of the total variance on Figure 5, Figure 6 and Figure 7 for the most active constraints : C1, C3 and C5. The first order index for a given parameter measures the part of the response variance explained by the effect of this parameter and effect of functions of only this parameter, whereas the total index measures the part of the response variance explained by all the effects in which this parameter plays a role (including interactions with other parameters). Note that a parameter is not displayed if its influence is smaller than 2% of the total variance. Constraint C1 (no compression) is mainly influenced by both the buoys with 𝐷𝑏, and the cable with 𝐿𝑐 and 𝐷𝑐, but not so much by the stiffener parameters.

For the maximum curvature constraint C3, the buoys are crucial as 𝐷𝑏 and 𝐿𝑏 are the first and second larger primary effects and the first and fourth larger total effects. The stiffener parameters 𝐷𝑠𝑏 and 𝐿𝑠 are, as expected, also essential. Finally, the constraint on the cable with respect to the sea floor (C5) is largely influenced by the diameter of the buoys 𝐷𝑏. The length of the buoys 𝐿𝑏 has a less significant role because it contributes less to the volume. Also, change of cable geometry with 𝐿𝑐 or 𝐿1𝑐 is prone to influence the risk of touching the sea floor.

The optimization phase can now be implemented. Because of the six nonlinear constraints, the feasible set is not convex in the parameter space, and the algorithm must faces the difficulty of local minima. To overcome this difficulty one can use a global optimization algorithm like in [START_REF] Martins | Compliant vertical access riser assessment: DOE analysis and dynamic response optimization[END_REF] and [START_REF] Yang | Multi-objective optimization for deepwater dynamic umbilical installation analysis[END_REF]. It is however difficult to control the stopping criteria and thus the computational time of such an algorithm. Alternately, when using an algorithm like SQA that converges to a local optimum, the choice of the initial points strongly determines the location of this solution.

Copyright © 20xx by ASME To find good diversified candidates for the initial points, a non parametric kriging response surface of each constraint is first built from the 600 points of the LHS used for the sensitivity analysis. The optimization problem defined in Eq.12 and Eq.13 is then solved fastly using these response surfaces. It is then possible to have multiple starting points (multi-start method see e.g. [START_REF] Saltelli | Sensitivity Analysis in Practice. A Guide to Assessing Scientific Models[END_REF]) to have a better chance to determine the global minimum. Considering the possible error of the response surface based on the limited LHS, we select not only the best solution, but two solutions as initial points of the accurate optimization based on FEM simulations and the 28 loading cases. The selection of these two initial points also considered the constraint activation with real simulations considering the 28 loading cases. The results of the simulation-based optimization for initial points defined in Table 2 are first presented. The total cost at optimum point is composed of 45% for the armors, 20% for the buoys and 35% for the stiffener.

The evolution of the constraints along the simulation based optimization process is illustrated on Figure 8. The maximum curvature (C3) particularly drives the optimization evolution, limiting the exploration of the parameter space with C1 (no compression) and C5 (risk of touching the sea floor). At the optimal solution, C1 and C3 are slightly positive, meaning that these constraints are not strictly satisfied.

The evolution of the parameters and the objective function are displayed on Figure 9. The first 15 iterations show peaks around the initial values because simulations around the initial points are performed in order to build first quadratic interpolation models of the objective function and of the constraints . Then, these quadratic models are optimized and updated with the simulations results at the current solutions. Finally, in the last iterations, we observed not converged simulations that make the SQA algorithm stop because of inaccurate quadratic models due to the lack of converged FEM simulations. They correspond to load cases with lateral current which require more static increments to guarantee the convergence. The results show that we are actually at the limit of capturing this non linearity. The amplitude of the parameter variations can be very disparate, up to two thirds of the range for the length of the cable 𝐿𝑐 to less than 1% for the ratio The total cost at optimum point is composed of 38% for the armors, 20% for the buoys and 42% for the stiffener.

The evolution of the constraints is illustrated on Figure 10 and again the optimization is driven by C3, C1 and C5. However this optimal solution strictly satisfies all the constraints.

CONCLUSION

This paper presents a new methodology for the design optimization of offshore wind farm inter array cable exploiting FEM hydrodynamic simulations based on the derivative free SQA optimization algorithm developed by IFPEN. This methodology is an efficient way to avoid limitations of Finite Difference gradient computations when the user has no access to the sources of the aero-hydrodynamic simulation solver.

To illustrate the potential of this approach, a case study is considered with a lazy wave cable configuration for 28 ULS extreme load cases. We compute solutions for seven parameters describing the geometry (length and diameter) of the cable, buoyancy modules and stiffener. The optimization problem which aims to minimize the material cost, is completed by six non-linear black-box constraints to ensure the cable integrity, computed from the simulation outputs. This case study could easily be extended to fatigue loads and constraints.

Analytical expressions are provided to compute the global scale mechanical properties of the cable when changing the design parameters, to set the hydrodynamic simulation inputs.

Copyright © 20xx by ASME Also, analytical expressions are given for the local stress on armor wires to compute the risk of tensile rupture.

At least because of non-convexity of the feasible set induced by the constraints, the optimization has to face the difficulty of local minima. To overcome this difficulty a multistart strategy based on a pre-optimization computed with a surrogate model is used.

The methodology involves first a computation of a kriging response surface with many fast optimizations for different initial points sampled in the design parameter space, for 6 loading cases. From these results, it is possible to select among these first solutions some candidates to initialize a more precise optimization, with all the loading cases and the FEM solver. This response surface model is used also to perform a global sensitivity analysis to rank the parameter effects on the variance of the optimized functions. It allows for instance to detect parameters that produce weak effects on the responses and thus to check the relevance of the chosen parameterization.

Two different optimized solutions obtained from different initial points are documented. The solutions provide a similar total cost but a very different cost distribution among the cable, the buoyancy module and the stiffener individual costs.

As a perspective to this work, we consider enriching the methodology by defining more accurate for the initial point selection. Another different approach could be to replace the multi-start by a first step with a global optimization technique based on adaptive surrogate response models [START_REF] György | Efficient multi-start strategies for local search algorithms[END_REF], [START_REF] Parr | Improvement Criteria for Constraint Handling and Multiobjective Optimization[END_REF] and a second step with a fine local optimization with our local derivative free optimization method SQA.

Figure 1 :

 1 Figure 1: initial cable geometry. The pink part of the curve represents the cable with buoyancy modules.

Figure 2 :

 2 Figure 2: typical cross section of MVAC cable.

Figure 3

 3 Figure 3 illustrates the cable geometry.

Figure 3 :

 3 Figure 3 : cable geometry for the optimization. The cable is composed of two armor layers with a mean radius 𝑎 ̅ = 𝑎 1 +𝑎 2 2

  found in the literature, axial stiffness and flexural stiffness are fixed as: (𝐸𝐴) 𝑐𝑜𝑟𝑒 + (𝐸𝐴) 𝑠ℎ𝑒𝑎𝑡ℎ = 90 ⋅ 10 6 𝑁 (𝐸𝐼) 𝑐𝑜𝑟𝑒 + (𝐸𝐼) 𝑠ℎ𝑒𝑎𝑡ℎ = 6.8 𝑘𝑁. 𝑚 2 . These properties are used for the large scale hydrodynamic momentum balance solved by Deeplines TM which also requires to define the linear mass of the cable : 𝑚𝑙 𝑐𝑎𝑏𝑙𝑒 = 𝑛 ̅⋅𝜋𝑟 2 cos(𝛼) ⋅ 𝜌 + 𝑚𝑙 𝑐𝑜𝑟𝑒 + 𝑚𝑙 𝑠ℎ𝑒𝑎𝑡ℎ Eq. 3 with 𝜌 = 7850 𝑘𝑔. 𝑚 -3 ; 𝑚𝑙 𝑐𝑜𝑟𝑒 + 𝑚𝑙 𝑠ℎ𝑒𝑎𝑡ℎ = 15 𝑘𝑔. 𝑚 -1 .

Figure 4 :

 4 Figure 4: mechanical equilibrium on an elementary wire part.

  in [2.5 D c , 3.5 D c ]. The optimization problem aims to find a set of parameter values in their feasible domain which minimizes an objective material cost function defined as: 𝐶(𝐿 𝑐 , 𝐿 𝑏 , 𝐷 𝑐 , 𝐷 𝑏 , 𝐿 𝑠 , 𝐷 𝑠𝑏 ) = 𝑛 ⋅ cos(𝛼) ⋅ π ⋅ r 2 𝐿 𝑐 • 𝐶 𝑣𝑐 +𝐿 𝑏 • 𝜋 • (

Figure 5 :

 5 Figure 5 : First order (light blue) and total (blue) effects of the design parameters in percent, for C1 constraint variance.

Figure 6 :

 6 Figure 6 : First order (light blue) and total (blue) effects of the design parameters in percent on C3 constraint variance.

Figure 7 :

 7 Figure 7 : First order (light blue) and total (blue) effects of the design parameters in percent, on C5 constraint variance

Figure 8 :

 8 Figure 8 : evolution of the normalized constraint for solution 1.

Figure 9 : parameter evolution during the optimization loop for solution 1 .

 91 Figure 9: parameter evolution during the optimization loop for solution 1.

Figure 10 : evolution of the normalized constraint for solution 2 .

 102 Figure 10 : evolution of the normalized constraint for solution 2.The evolution of the parameters and the objective function are displayed on Figure11. The initial parameters were already close to their final values with the notable exception of 𝐿1𝑐 which was reduced because of active constraints.

Figure 11 :

 11 Figure 11: parameter evolution during the optimization loop for solution 2.

Table 1 : Design Load Cases (ULS).

 1

			Wave	Current	Wind
		Hs	Tp	Speed	Speed
	Production case	H s-1year H s-1year H s-1year	Associated T p -2s Associated T p Associated T p +2s	no no no	V cut-out V cut-out V cut-out
	Parked case	H s-50year	Associated T p -2s	Cur. 50years	V 50year
		H s-50year	Associated T p	Cur. 50years	V 50year

Table 2 : Parameters value for the solution 1.

 2

		L 1c	L b (m) L c (m) D c	𝐷 𝑏	𝐿𝑠	𝐷 𝑠𝑏
		(m)			(mm)	𝐷 𝑐	(m)	𝐷 𝑐
	Initial	66.2	34.1	534.8	133	2.1	3.3	2.8
	Optim. 65.9	33.4	541.1	131	2.0	3.3	2.8

ACKNOWLEDGMENTS

The authors would like to thank IFP Energies Nouvelles which has funded this work, and have benefited from the help of peoples from the Applied Mechanics Department and the Mechatronics Departments. Finally, the continuous support of Principia on Deeplines TM has been greatly appreciated.