T Pottier
email: thomas.pottier@mines-albi.frwww.institut-clement-ader.org

F Toussaint

P Vacher

Contribution of heterogeneous strain field measurements and boundary conditions modelling in inverse identification of material parameters

Keywords:

The present paper aims at applying the Finite Element Updating inverse method to several sample geometries by the means of Digital Image Correlation. The full-field data are experimentally obtained from three geometries exhibiting increasing strain fields heterogeneities. For each test, a Finite Element model is built and boundary conditions are duplicated from the measured displacements at the sample borders. Field comparisons are performed at several time steps until fracture occurs and a Levenberg-Marquardt method is used to solve the optimization problem. Six parameters of an anisotropic elastic-plastic constitutive model are identified and validated through the simulation of a deep-drawing forming operation. Results show that identification quality is improved when heterogeneous strain fields are used.

Introduction

Facing the increasing number of new materials and the difficulties of identifying constitutive models that predict properly the material behaviour during the forming process, industrials have induced in the early 90's some scientific research that challenges the usual way of identifying material models. This former approach consists in fitting analytical model with the experimental behavior (true stress/true strain curve). Basically, tests are led on standard geometries to isolate one material behaviour from the others (e.g. a tensile test is assumed to exhibit tensile behaviour only), then hypothesis concerning the strain field homogeneity are made and lead to identify the corresponding constitutive parameter. Once these constitutive parameters are determined, the Finite Element (FE) analysis is used to retrieve the stress/strain/displacement fields, assuming that both geometry and boundary conditions are known. This identification procedure exhibits two main disadvantages: the numbers of tests may increase significantly when complex models are in use and the hypothesis of homogeneous strain fields is neither always obvious nor easy to satisfy. On the other hand, the so-called inverse approach proposes to determine constitutive parameters from considering heterogeneous strain/displacement fields. In this case, geometry and boundary conditions are known and so is the resulting displacement field. Therefore, the only constitutive parameters remain unknown. The chosen constitutive equations are, of course, assumed to be relevant.

As Summarized by [START_REF] Grédiac | The use of full-field measurement methods in composite material characterization: interest and limitations[END_REF] and [START_REF] Avril | Overview of identification methods of mechanical parameters based on full-field measurements[END_REF], several different methods have been developed to obtain the solution of such inverse problems. Indeed five approaches lead to retrieve the constitutive parameters when the stain homogeneity hypothesis does not hold. Among these, the constitutive equation gap method [START_REF] Ladevèze | Constitutive relation error estimators for (visco)plastic finite element analysis with softening[END_REF][START_REF] Latourte | Elastoplastic behavior identification for heterogeneous loadings and materials[END_REF], the equilibrium gap method [START_REF] Claire | A finite element formulation to identify damage field: the equilibrium gap method[END_REF][START_REF] Crouzeix | An orthotropic variant of the equilibrium gap method applied to the analysis of a biaxial test on a composite material[END_REF], the virtual fields method [START_REF] Grédiac | Principe des travaux virtuels et identification/principle of virtual work and identification[END_REF][START_REF] Promma | Application of the virtual fields method to mechanical characterization of elastomeric materials[END_REF], the reciprocity gap method [START_REF] Ikehata | Inversion formulas for the linearized problem for an inverse boundary value problem in elastic prospection[END_REF] and the finite element model updating (FEU) technique [START_REF] Kavanagh | Finite element applications in the characterization of elastic solids[END_REF][START_REF] Molimard | Identification of the four orthotropic plate stiffnesses using a single open-hole tensile test[END_REF] are currently used for identification purpose. In the present paper the FEU approach has been chosen.

The FEU method provides constitutive parameters iteratively by comparing the measured data with the numerical ones while parameters vary. The experiment is accurately simulated in a FE model and the calculation is run with an initial guess of the parameters values. The obtained displacement fields are then compared to the measured ones, a cost function is built and leads to a new set of parameters and so forth, until convergence is reached. The inverse methods applied to parameters identification such as defined first by [START_REF] Kavanagh | Finite element applications in the characterization of elastic solids[END_REF] have known massive improvements in the last decade. Indeed, the FEU approach is now broadly used. Several authors have applied FEU to identify material parameters using so-called global measurements, such as reaction force of the whole sample, as the material response. For instance, [START_REF] Ghouati | Identification of material parameters directly from metal forming processes[END_REF] have identified 5 parameters of an elastic-plastic constitutive law of aluminum alloy using tensile tests in a first approach and a deep drawing test in a second one. Both results were very close and authors conclude that this approach is an alternative to the classical identification method. Moreover, the experiments led by [START_REF] Yoshida | Identification of material parameters in constitutive model for sheet metals from cyclic bending tests[END_REF] have successfully identified both elastic and elastic-plastic parameters using flexion tests. The torque was used as the material response and the experimental and numerical curvature/torque curves were compared for several parameters sets until fitting was reached. Both of these studies have provided satisfactory results but none of them was using full-field measurement as experimental data. Meanwhile, [START_REF] Meuwissen | An inverse method for mechanical characterization of metal[END_REF] has developed a similar approach with the use of full-field measurements. Shear-like tensile tests were performed and displacement fields were measured. A cost function was built considering global reaction force and displacement fields as the material responses. For the sake of simplicity, only the early beginning of the plastic process was taken into account and neither high strain nor necking was considered. Recently this work has been completed by [START_REF] Kajberg | Characterization of materials subjected to large strains by inverse modeling based on in-plane displacement fields[END_REF]. This latter study deals with through necking analysis of basic tensile specimens with side notches. This geometrical singularity ensures the necking to occur at the center of the sample and improves the heterogeneity of the strain fields. The identification of two elastic-plastic model is achieved but no conclusion is drawn upon the contribution of considering the whole deformation process instead of the early stages of it. Finally the work of [START_REF] Lecompte | Mixed numerical-experimental technique for orthotropic parameter identification using biaxial tensile test on cruciform specimens[END_REF] and [START_REF] Cooreman | Elasto-plastic material parameter identification by inverse methods: Calculation of the sensitivity[END_REF] also has to be noticed. As an improvement to the works presented above, it introduced the use of a significantly more heterogeneous tests thanks to bi-axial tensile tests. Cruciform samples are then used to identify elastic-orthotropic parameters.

The present paper aims to show that identification quality improves and the required number of experiments decreases as the heterogeneity of strain fields increases. For this purpose, the FEU inverse method is applied to several sample geometries, using Digital Image Correlation data. Three geometries with increasing strain field heterogeneity are tested, and provide the parameter set of an anisotropic elastic-plastic model. For the sake of comparison, the same identification procedure is used for the three geometries. Boundary conditions of the FE model are strictly applied from the experiment and allows the FE model to fit perfectly the experimental deformation process. Thus, the identification calculations are run again using virtually noised experimental data in order to exhibit and quantify the sensitivity of the identification method to measurement noise. Finally, validation tests are led on a deep-drawing operation and results are used to compare the quality of the obtained parameter sets.

Experimental Setup Material and samples geometries

All the samples under investigation are made of Commercially Pure titanium T40, the limits of the chemical composition of the material are given in Table 1. The average grain size has been estimated to 6.4 µ m. This material, broadly used in aeronautics and biomedical domains, consists of hot-rolled sheets where the 0 • direction is chosen to fit the rolling direction. As a consequence, 90 • is the transverse direction of the rolled sheet. Sheet thickness is 0.5 mm.

As shown in Fig. 1, three sample shapes were used. A basic tensile sample (no.1) as used by [START_REF] Mahnken | Aspects on the finite-element implementation of the Gurson model including parameter identification[END_REF], [START_REF] Ghouati | A finite element-based identification method for complex metallic material behaviours[END_REF] and many others, will first be investigated. Then, the sample geometry is added a hole (no.2) as suggested by [START_REF] Lecompte | Elastic and elasto-plastic material parameter identification by inverse modeling of static tests using digital image correlation[END_REF] and finally a shearlike tensile sample (no.3) will be studied. This latter sample, originally proposed by [START_REF] Meuwissen | An inverse method for mechanical characterization of metal[END_REF], is expected to exhibit both tensile and shear behaviours.

Full-field measurements

In the present study, The Digital Image Correlation method is used to assess the local displacements at the sample surface with 7D correlation software [START_REF] Vacher | Bidimensional strain measurement using digital images[END_REF]. The experimental setup is shown in Fig. 2a.

Tests are led on a 5 kN INSTRON tensile device. The samples undergo a prescribed displacement speed of 5 • 10 -5 m • s -1 (corresponding to an approximate strain rate ε ≈ 1 • 10 -3 s -1 ) until their ruin. Samples are speckled using black and white mate painting sprays. Digital images are recorded during the deformation process at a frequency of 0.25 Hz using a Nikon D200 camera (focal distance 105 mm). The image resolution is set to 2592 × 3872 grey level pixels and stored as 8 bit (Fig. 2b). Moreover, axial reaction force is recorded every 10 µs all along the test. For each test, the initial image is split in square elements that create a virtual grid upon the sample surface. The resolution of this grid (extensometric base) is set to 16×16 pixels (Fig. 2b), corresponding to 0.32 × 0.32 mm 2 . The correlation process consists in looking for the most probable deformed pattern in the neighbourhood of each node of this grid in terms of grey level. The displacement fields of each element are then assessed by the means of a bi-linear interpolation. Finally, the displacement at image n is assessed by correlating the n th image with the initial/undeformed image. The correlation parameter used in the present study are summarized in Table 2.

The assessment of DIC technique uncertainties remains a topical issue [START_REF] Bornert | Assessment of digital image correlation measurement errors: Methodology and results[END_REF], [START_REF] Tarigopula | A study of localisation in dualphase high-strength under dynamic loading using digital correlation and FE analysis[END_REF] and [START_REF] Kajberg | Characterization of materials subjected to large strains by inverse modeling based on in-plane displacement fields[END_REF]. Several parameters may have influence on the measurement noise such as lighting, strain gradients, off-plain displacements, numerical noise of the camera and of course the quality of the speckle [START_REF] Pan | Twodimensional digital image correlation for in-plane displacement and strain measurement: a review[END_REF]. However, the works of [START_REF] Bornert | Assessment of digital image correlation measurement errors: Methodology and results[END_REF] and [START_REF] Coudert | Reconstruction tridimensionnelle du volume intérieur d'une chaussure : évaluation du chaussant[END_REF] (performed with the correlation parameter presented in Table 2) state that the order of magnitude of the DIC measurement noise remains inbetween 0.01 and 0.1 pixel.

Identification Scheme

Inverse flowchart

In the present study, the FEU solving flowchart (Fig. 3) is chosen and applied to parameters identification . Experiments are led and post-processed using 7D correlation software then data are imported into Matlab. A finite element model of the sample is designed and meshed using Abaqus-Explicit. This latter FE model is multiply solved for several sets of material parameters using the in situ prescribed conditions obtained by the means of DIC. Then, the FE calculation outputs are post-processed using Matlab in order to build the cost function that leads to the solution of the global identification problem.

Cost-function

Since, the comparison between experimental and numerical fields has to be computed, the relevant point is the manner of writing this comparison down. Several authors such as [START_REF] Meuwissen | An inverse method for mechanical characterization of metal[END_REF], [START_REF] Kajberg | Characterization of materials subjected to large strains by inverse modeling based on in-plane displacement fields[END_REF] and [START_REF] Cooreman | Elasto-plastic material parameter identification by inverse methods: Calculation of the sensitivity[END_REF] are using a weighted least square formulated cost-function. The introduction of a weighting matrix has basically two main kinds of formulation. First, it could be made up of the experimental errors and make the cost function a statistical estimator as did by [START_REF] Lecompte | Elastic and elasto-plastic material parameter identification by inverse modeling of static tests using digital image correlation[END_REF] or [START_REF] Meuwissen | An inverse method for mechanical characterization of metal[END_REF]. Each nodal response is then compared with its experimental equivalent and divided by the estimated error at the considered material point. As a consequence, the more reliable a measure is, the (b) (a) more weight it has on the identification solution. Another approach consists in considering the weighting matrix as a way to give each material point the same weight. This solution consists in dividing each local residual by the measured response at this point. Such a formulation, used by [START_REF] Ponthod | A cascade optimization methodology for automatic parameter identification and shape/process optimization in metal forming[END_REF], [START_REF] Cooreman | Elasto-plastic material parameter identification by inverse methods: Calculation of the sensitivity[END_REF] and [START_REF] Khalfallah | Anisotropic parameter identification using inhomogeneous tensile test[END_REF], also allows the use, in the same cost-function, of several kind of quantities exhibiting miscellaneous orders of magnitude. This latter approach is chosen in the present study. Moreover, the residual quantity can be made of force and/or displacement data [START_REF] Kavanagh | Finite element applications in the characterization of elastic solids[END_REF]. As mentioned in [START_REF] Avril | Overview of identification methods of mechanical parameters based on full-field measurements[END_REF], considering both force and displacement in the cost function formulation lead to a better suited solution when strain concentration is involved (geometries no.2 and no.3 for instance). Therefore, in the present study, force and displacement data has been mixed using a weighting matrix in order to obtain a dimensionless formulation of the cost function. Indeed, longitudinal and transversal displacements are added to the global reaction force residual and the local difference between these three measured data and the calculated values is the function to be minimized. This cost-function can be written as follows:

f (p) = 1 2 Nt,Ns,Nn i,j,k=1   u exp x,jk (p) -u num x,jk max j,k (u exp x )   2 +   u exp y,jk (p) -u num y,jk max j,k (u exp y )   2 + N n F exp j (p) -F num j max j (F exp ) 2 i 1 2
.

(

) 1
Where p is the unknown parameter column, u x is the displacement along the x axis (longitudinal axis) ,u y the displacement along the y axis (transverse axis) and F is the global reaction force of the sample. The scalar N n is the number of nodes in the finite element mesh, N s is the number of time steps taken into account, and N t is the number of tests that are considered. Finally, max j,k (u exp x ), max j,k (u exp y ), and max j (F exp ) are respectively the nodes greatest values over one step-time of the longitudinal and transversal displacements and the global reaction force. Hence, the optimal set of parameters must satisfy three conditions: matching the two displacement fields and the global force response (which is actually not a field response but a punctual one). The next section discuss the algorithm used to solve this optimization problem.

Optimization algorithm

Various numerical algorithms are available to solve non-linear least square problems represented by a cost function such as Eq. ( 1). Here, a Levenberg-Marquardt (LM) [START_REF] Levenberg | A method for the solution of certain nonlinear problems in least-squares[END_REF] and [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF] method is used to minimize the cost-function. This method is known to have several advantages among which can be singled out that it does not require the assessment of the Hessian matrix and it's thus very easy to implement. Let's recall that the updating equation of a LM method is:

J T J + λ LM I p (k+1) -p (k) = J T r(p) (2)
Where r(p) is the residual column constituted by the three (u x , u y and F ) relative field differences at each node. Matrix J is the Jacobian matrix that is often called sensitivity matrix. This latter matrix is assessed by computing the partial derivative, among each constitutive parameter, of every single system response (see section Sensitivity analysis). Finally, λ LM is a damping parameter that has been introduced by [START_REF] Levenberg | A method for the solution of certain nonlinear problems in least-squares[END_REF] to stabilize the Gauss-Newton method in the vicinity of the optimum. The experimental errors on the measured responses (displacement fields and global reaction forces) lead to the inexistence of an optimal parameters set. In other words, updating the constitutive parameters is not sufficient to obtain numerically the exact observed behaviour. Hence, the optimization problem is ill-posed [START_REF] Lecompte | Elastic and elasto-plastic material parameter identification by inverse modeling of static tests using digital image correlation[END_REF] and a basic Gauss-Newton's method may become very unstable around the optimum. The addition of the λ LM parameter improves the robustness of the algorithm near the optimum. As proposed by [START_REF] Ma | Some research on Levenberg-Marquardt method for nonlinear equations[END_REF], the value of λ LM is assessed as follows:

λ LM = θ r(p) + (1 -θ) J T r(p) (3)
Where θ ∈ [0, 1] is a real parameter that must be set a priori and that rules both accuracy and convergence rate of the algorithm. In the present work θ = 8 • 10 6 .

Sensitivity analysis

The so-called sensitivity analysis consists in assessing the terms of the Jacobian matrix. As presented by [START_REF] Tortorelli | Design sensitivity analysis: Overview and review[END_REF], sensitivity assessment is a major issue in parameter identification. For this purpose a broadly used approach consists in computing sensitivities by the means of finite differences. Even though the calculation time is high, this approach is easy and fast to implement. Solving each iteration require only k + 1 calculations if a forward (or backward) scheme is chosen. For the sake of calculation time, a forward finite difference scheme is used in this work, therefore sensitivities are assessed by the following assumption:

J i,j = ∂r i (p) ∂p j ≈ r i (p) -r i (p + ∆p j ) 2∆p j ( 4
)
Where j is the number of constitutive parameters to be identified and i is the number of system response, ∆p j is the relative perturbation on the j th parameter. The choice of this perturbation ∆p j remains a controversial topic. Indeed, the above equation exhibits a truncation error of ∆p 2 j that tends toward zero when the perturbation decreases. On the other hand, it has to be noticed that the decreasing of the perturbation value increases the weight of the round-off error, therefore the value of ∆p j results from a compromise between these two errors. In the present study the perturbation value is set to ∆p j = 0.01 × p j .

Material and Numerical Model Numerical model

The numerical model, developed using Abaqus, duplicates the experiment in terms of geometry and boundary conditions. The geometry is measured on the samples, modeled and meshed using quadrangular shell elements (the number of nodes varies from 1458 to 2201 depending on the considered geometry). The obtained meshes fit the three geometries under investigation. As shown in Fig. 4, boundaries conditions, are measured from the experiments using DIC and the obtained displacements are prescribed at each node of the upper and lower bounds of the mesh [START_REF] Abanto-Bueno | Experimental determination of cohesive failure properties of a photodegradable copolymer[END_REF][START_REF] Witz | Mechanical properties of crimped mineral wools : Identification from digital image correlation[END_REF][START_REF] Fedele | Identification of adhesive properties in glare assemblies by digital image correlation[END_REF]. Both longitudinal and transverse displacements are imposed. However, the data measured from the DIC may not properly fit the chosen nodes of the FE mesh. In other words, the displacement data may be obtained at a material point that does not match a node of the FE mesh because DIC grid and FE mesh are different. Therefore, the DIC displacement fields are interpolated at each node of the FE mesh so that the fields comparisons (needed for the cost function assessment) can be computed at the same material points (Fig. 4). This projection is processed using bi-linear interpolation shape functions.

Constitutive equations

The presented identification scheme can basically be applied to any kind of constitutive model. In the present study, a classical anisotropic elastic-plastic model is considered. The elastic strain is here defined by the Hooke's formulation, where E and ν, respectively the Young's modulus and the Poisson's ratio are previously identified (E = 111.8 MPa and ν = 0.34). A large strain framework is used and σ is the Cauchy stress tensor, E is the Green-Lagrange strain tensor. The generalized standard material formalism provides the plastic strain, that can be written as follows:

dE p = dλ ∂f (σ) ∂σ , dλ ≥ 0 (5)
where dλ is a non-negative scalar. This latter formulation allows the introduction of a yield surface and therefore a yield criterion can be written:

f (σ) = J (σ) -R -σ y ( 6
)
Where R is the current yield stress, σ y = 368.0 MPa is the initial yield stress and J (σ) is the equivalent plastic stress. In the present study, an orthotropic material behaviour will be assumed and the Hill's 1948 [START_REF] Hill | A theory of the yielding and plastic flow of anisotropic metals[END_REF] criterion is used. This criterion is commonly used for sheet metal anisotropic behaviour and it provides an expression of the equivalent stress. In the case of plane stress, which is assumed here, the Hill's criterion gives:

J (σ) = F σ 2 22 + Gσ 2 11 + H(σ 2 11 -σ 2 22 ) + 2N σ 2 12 ( 7
)
As a consequence, the description of material anisotropy only requires four parameters to be identified. As presented by [START_REF] Hill | A user-friendly theory of orthotropic plasticity in sheet metals[END_REF], in practice the use of the Lankford's coefficients is easier. The present study will focus on the determination of these coefficients that are defined in section Standard identification. Moreover, the material hardening described by R has to be considered. Several models have been developed in the last decades in order to take into account those effects [START_REF] Gronostajski | The constitutive equations for FEM analysis[END_REF]. For the sake of simplicity, the Ludwik's flow rule is chosen, it is both easy to implement and to identify because only two parameters are needed. The writing of such a flow rule is recalled as:

R = K(E eq p ) n (8)
Where, E eq p is the equivalent plastic strain. K and n are two parameters to be identified. Finally, an anisotropic elastic-plastic model is in use and the present identification study deals with only six parameters.

Results and Validation

Standard identification

In order to offer a comparison criterion to the results obtained by the means of inverse method, a standard identification has been proceeded. As presented by [START_REF] Toussaint | Experimental and numerical analysis of the forming process of CP titanium scoliotic instrumentation[END_REF], this approach is based on several assumptions and provides a quick assessment of material constitutive parameters. The strain fields obtained using DIC are used to localize the necking area of the sample. Therefore, only the local strain data are considered and the volume conservation assumption is made to allow the assessment of the width at the necking location and the troughthickness Hencky's strain (ε p III ) [START_REF] Sell | Experimental characterization of deformation damage in solid polymers under tension, and its interrelation with necking[END_REF]. Thus, the true stress can be computed and plotted versus the plastic equivalent strain (Fig. 5a). Hence, the hardening curve is fitted using a least-square routine which leads to identify both K and n: the Ludwik's parameters.

Anisotropic behaviour is also identified and the Hill's parameters are assessed through the Lankford's coefficients. Let's recall that those coefficients are defined as follows:

r θ = dε p III dε p II ( 9
)
Where θ is the angle between the rolling and the applied tensile direction. This ratio of the transverse plastic strain rate to the thickness plastic strain rate is experimentally represented by the slope of the measured strains in the tensile tests led at 0 • , 45 • and 90 • from the rolling direction. Their respective least-square lines are also plotted in Fig. 5b. Finally the three Lankford's parameters: r 0 , r 45 and r 90 are identified and lead to assess the four Hill's parameters that are summarized in Table 3.

Inverse identification

For the sake of calculation time, the parameters identification is processed considering only a restricted number of time steps arbitrary chosen over the hardening curve. This number N s varies inbetween 5 and 7 and they can be seen as markers in Fig. 6.

For every geometry, three tests are considered (N t = 3): at 0 • ,45 • and 90 • except for geometry no.3 where two identification procedures are run: the first with N t = 3 and the second with N t = 2. This aims to show that when strain fields heterogeneity increases, the number of required tests may be decreased without impairing the identification quality.

In order to impose the necking to occur at the same location in experiments as in the FE analysis of geometry no.1, a small width default (0.05 mm on each side) is introduced on the initial mesh. No special treatment or assumption is made for the other analysis. The initial hardening parameters are arbitrarily chosen equal to K = 550 MPa and n = 0.3, the Hill's coefficients are chosen to initially match the isotropic case: F = G = H = 0.5 and N = 1.5. The identified parameters are summarized in Table 3. and compared to those identified using the standard method.

Table 3 shows that the proposed inverse method is able to decrease the number of tests needed for identification and thus the experimental work. Indeed parameters sets no.3-a and no.3-b are very similar and can be considered as equal (denoted set no.3).

As shown in Fig. 6, the four obtained sets (sets no.1, no.2, no.3 and the standard identification set) are used to predict the global reaction force of three tests with the different geometries. It can be seen that the set no.3 is the most able to provide a good force prediction regardless of geometry. The fitting of the set no.1 is very poor. Indeed this set, identified in a pure tensile load case, appears to be unable to predict the global reaction force when other kinds of strain fields are involved (e.g. test no.3). This emphasizes the improvement of parameter identification when heterogeneous strain fields are used.

Moreover, Fig. 7 shows the identified resulting displacement fields versus the experimental ones for set no.3, the residual values are also plotted for both longitudinal and transverse displacement fields. As can be seen, the optimization process leads to the fitting of experimental and numerical responses.

Experimental noise sensitivity

Since the parameter set is identified, the predicting capability of it may be investigated. As pointed in many former studies [START_REF] Gelin | Modelling and control of hydroforming processes for flanges forming[END_REF][START_REF] Harth | Identification of material parameters for inelastic constitutive models: statistical analysis and design of experiments[END_REF][START_REF] Molimard | Identification of the four orthotropic plate stiffnesses using a single open-hole tensile test[END_REF][START_REF] Lecompte | Mixed numerical-experimental technique for orthotropic parameter identification using biaxial tensile test on cruciform specimens[END_REF][START_REF] Amiot | Identification of elastic property and loading fields from full-field displacement measurements[END_REF][START_REF] Belhabib | Heterogeneous tensile test on elastoplastic metallic sheets: Comparaison between FEM simulations and fullfield strain measurements[END_REF][START_REF] Nakamura | Determination of nonlinear properties of thermal sprayed ceramic coatings via inverse analysis[END_REF][START_REF] Promma | Application of the virtual fields method to mechanical characterization of elastomeric materials[END_REF], investigating the noise sensitivity if the identification procedure is a good way to evaluate its stability. In practice, the noise of the measuring tools may be responsible for a high deviation of the obtained results. As said above, the accuracy of DIC method is limited and the measurement uncertainties assessment is a tough task. However some authors [START_REF] Coudert | Reconstruction tridimensionnelle du volume intérieur d'une chaussure : évaluation du chaussant[END_REF][START_REF] Bornert | Assessment of digital image correlation measurement errors: Methodology and results[END_REF] state that their magnitude remains inbetween 0.01 and 0.1 pixel. In the present study, displacement fields were artificially noised using a sample-wide random white noise of which magnitude varies. A random uniform noise is applied to both longitudinal and transverse displacement fields of the geometry no.3. The noise magnitude A n is set to several increasing values then the obtained parameter sets are compared. Results are summarized in Table 4.

In addition, for the three noised sets, the predicted longitudinal displacement field is compared to the noiseless field. Fig. 8 shows the relative error e(A i ) for the three noised sets versus the noiseless field such as:

e(A i ) = U num x (p(A i )) -U num x (p(A 0 )) U num x (p(A 0 )) (10)
As can be noticed, the increasing of measurement noise quickly impairs the quality of the predicted displacement fields. This analysis shows that the FEU inverse method is obviously noise sensitive but also that for the assumed noise magnitude (between 0.01 and 0.1 pixel) the results dispersion remains acceptable (inferior to 3.9%) and in the range of experimental uncertainties. Moreover, the obtained dispersion for 0.1 pixel matches other published results for inverse methods such as [START_REF] Promma | Application of the virtual fields method to mechanical characterization of elastomeric materials[END_REF] and [START_REF] Kajberg | Characterization of materials subjected to large strains by inverse modeling based on in-plane displacement fields[END_REF].

Validation test

In order to verify the quality of the identified parameters sets, validation tests are carried out. The three identified parameters sets are used to simulate the same deep-drawing operation. Results are compared to the experimental data in order to check which of the identified sets is the most able to reproduce the observed experimental behaviour.

On one hand a deep-drawing device (Fig. 9a) is used for the purpose of experimental validation. The elastic strength of both springs is previously identified and the device is operated using the same tensile machine as described above. Three circular blanks following the geometrical specification given in Fig. 9b are deep-drawn while operating force is recorded. Hence, the geometrical responses of the three obtained cups are measured. As shown on Fig. 9c the upper profile of the cups exhibit anisotropic horns and the magnitude and frequency of it are investigated and measured. A Coordinate Measuring Machine (CMM) is used to measure the height of 120 points evenly spread along the upper periphery. On the other hand, a FE model is developed to simulate the experimental deep drawing process. Solid elements are used to model the blank while die, blank-holder and punch are considered as analytical rigid surfaces. The * 2.66 GHz and RAM 0.98 Go

• -45 • -90 • 0 • -45 • -90 • 0 • -45 • -90 • 0 • -45 • -90 • 0 • -45
A 1 =0.01 pixel A 2 =0.1 pixel A 3 =1 pixel Relative error (%)
Relative error (%) Relative error (%) applied load on the blank-holder is picked from the experimental measurements. The friction coefficient is set to 0.2 in agreement with former results [START_REF] Le Port | Finite element study and sensitivity analysis of the deep-drawing formability of commercially pure titanium[END_REF]. Finally the identified parameter sets are successively used. Then, the upper profile is compared to the experimental one (Fig. 10). The simulations with inverse parameter sets exhibit very scattered anisotropic responses, indeed eventhough geometries no.2 and no.3 show a good fitting with the experimental measurements, the response of geometry no.1 is poor. This latter set predicts horns at the wrong place and with a bad magnitude. The simulation using the parameters identified using geometry no.3 predicts more accurately the magnitude of the horns than geometry no.2 and no.1. The results of the validation tests, presented in Fig. 10 show that the quality of the shape prediction increases with the strain field heterogeneity. The numerical shape of the cup has a better fitting of the experimental shape when shearlike identified parameters are used. Eventhough the standard identification set remains the best for shape prediction, among the three inverse sets, the increasing of the strain heterogeneity in the sample lead to a more accurate identification of the material models.

Conclusions

In the present study, a FEU method has been used to identify the parameters of an anisotropic elastic-plastic model using full-field measurements. Three geometries that exhibit increasing strain heterogeneities have been investigated. Boundary conditions were picked at the borders of the DIC domain and strictly duplicated in the FE model. The results show that considering two experiments can provide the same results as three. Hence, the increasing of strain fields heterogeneity allows to perform less tests than homogeneous tests. Validation tests and comparisons between the identified parameter sets have shown that the ability of predicting the real deformation process is improved when the strain heterogeneity increases. However, the standard identification method remains the most able to predict material behaviour. This highlights the necessity of more heterogeneous tests when inverse methods are used. Finally, experimental data were noised in the range of known admissible DIC uncertainties. The obtained dispersions remain in the range of other published works. The presented identification tool provides satisfying results with kinematics data but cannot identify fully coupled thermal-mechanical models so far. Thus, further works will focus on considering both kinematic and thermal fields but also on the ways to experimentally obtain heterogeneous strain fields

Table 1 :

 1 Chemical composition limits of commerciallypure grade 2

Figure 1 :

 1 Figure 1: Investigated geometries and DIC analysis zones (dimensions in mm).

Figure 2 :

 2 Figure 2: (a) Experimental setup. (b) 2592 × 3872 grey level images from Nikon D200 camera and longitudinal displacement field calculated by 7D correlation software.

Figure 3 :

 3 Figure 3: Iterative solving flowchart of the inverse problem.

Figure 4 :

 4 Figure 4: Experimental field measurements and boundary conditions were interpolated from the DIC analysis.

Figure 5 :

 5 Figure 5: Standard identification of the five investigated parameters. (a) Ludwik's parameters. (b) Lankford's coefficients.

igure 6 :Figure 7 :

 67 Figure 7: Experimental and identified displacement fields for the geometry no.3-a at 45 • (5 th time step). Residual fields show the zones of good and bad predictions.

Figure 8 :Figure 9 :

 89 Figure 8: Relative error of the longitudinal displacement fields for the three noised set.

Envelope of the 3 Figure 10 :

 310 Figure 10: Comparison of predicted anisotropic profiles.

Table 2 :

 2 Main parameters of the Digital Image Correlation

	Grid	pattern	grey level	displacement fields
	size	size	interpolation	interpolation
	16 × 16 16 × 16	bicubic	bilinear

Table 3 :

 3 Parameters identified by inverse method.

		Standard	Geom. no.1 Geom. no.2 Geom. no.3-a Geom. no.3-b
		identification
	N s	-	6	4	6	6
	N t	3	3	3	3	2
	Orientations	0

Table 4 :

 4 Parameters sensitivity to measurement noisefor shear-like geometry (no.3).

	A n (pixel)	0	0.01 0.1	1
	A n (µm)	0	0.2	2	20
	K	284 277 269 334
	n	0.50 0.48 0.49 0.48
	F	0.16 0.16 0.14 0.25
	G	0.18 0.19 0.15 0.26
	H	0.82 0.81 0.85 0.74
	N	2.02 2.01 2.01 2.01


Shanal S Kumar

Krishnil R Ram
email: krishnil.ram@usp.ac.fj

Maurizio Cirrincione
email: maurizio.cirrincione@usp.ac.fj

Ali Mohammadi
email: ali.mohamamdi@usp.ac.fj

Vincent Léchappé
email: vincent.lechappe@insa-lyon.fr

A Simplified Control Oriented Model Of an Open Cathode PEM Fuel Cell

Keywords: fuel cell, hydrogen, dynamic model, sliding mode, open cathode. I

the paper presents a simplified nonlinear model for an open cathode Proton Exchange Membrane Fuel Cell (PEMFC) and its control using two different control strategies. The model presented uses only three state variables and is therefore easy to model and control. The mass flow of oxygen, hydrogen and water were taken to be the key dynamics in the system. The unknown parameters were estimated using the experimental data of a 1.2 kW PEMFC. The model showed good agreement with experimental results. With the objective of maintaining a fixed oxygen excess ratio -a Proportional Integral (PI) Control and a Sliding Mode Control (SMC) scheme were applied to the system model. This paper compares the results of a simple PEMFC model without controller, with a PI controller and a Sliding Mode controller. The SMC performs better in terms of maintaining a fixed oxygen excess ratio of 2. The results also show that the simplified model has good use at an early design stage for systems using PEMFC and also for learning purposes.

INTRODUCTION

Fuel Cells offer an emission-free means of energy conversion in stationary power applications as well as in the transport sector. Several types of Fuel Cells such as Solid Oxide Fuel Cells (SOFC), Direct Methanol Fuel Cells (DMFC) and Polymer Electrolyte Membrane (or Proton Exchange Membrane) (PEM) Fuel Cells dominate the fuel cell market [1]. This paper focuses mainly on PEM Fuel Cells. In general, fuel cells have an anode and a cathode side. Hydrogen enters through the anode while oxygen enters through the cathode. The reactions of the two gases allow for the generation of electrical energy across the Fuel Cell terminals. With only water and heat as by-products, the PEM Fuel Cell (PEMFC) has no carbon dioxide emission. Open Cathode Fuel Cells do not use supply manifolds fed by compressors and humidifiers hence have much less patristic losses and have fewer components to maintain. This has made open cathode PEMFC a common choice, especially in transport sectors. In PEMFC, fans are used to force air through open cathodes to channel the oxygen onto the Membrane Electrode Assembly (MEA). An open cathode PEMFC has lower costs due to fewer components. The loss of auxiliary equipment requires a better control of the operating parameters for optimal operation. A key control objective to ensure longer life of PEMFC aims to avoid oxygen starvation during operation. Oxygen is normally obtained from atmospheric air, which is pulled into the PEMFC anode using axial fans. In case of oxygen starvation, the Fuel Cell suffers from degradation and would not be able to supply the required power [START_REF] Laghrouche | Load governor based on constrained extremum seeking for PEM fuel cell oxygen starvation and compressor surge protection[END_REF], [START_REF] Yang | Control challenges and methodologies in fuel cell vehicle development[END_REF]. For this reason, an excess of oxygen in the anode is necessary to avoid starvation.

Several works have proposed dynamic models to control the oxygen excess ratio in PEMFC [START_REF] Chen | Optimal Oxygen Excess Ratio Control for PEM Fuel Cells[END_REF], [START_REF] Danzer | Model-based control of cathode pressure and oxygen excess ratio of a PEM fuel cell system[END_REF]. The models vary in their focus on different dynamics in the PEMFC system. Reference [START_REF] Pukrushpan | Modeling and control for PEM fuel cell stack system[END_REF] presents PEMFC modelling for controlling a large PEMFC with a compressor used to force humidified air into the cathode side. Another paper [START_REF] Lu-Ying | An improved smallsignal model of the dynamic behavior of PEM fuel cells[END_REF] has focused on the pressure balance to effectively model the dynamic effects in PEMFC to improve internal designs and external controller designs. PEMFC emulators with dynamic models have also been proposed as a cheap alternative for hardware testing [START_REF] Marsala | A prototype of a fuel cell PEM emulator based on a buck converter[END_REF]. Detailed models are very useful to observe all the different variables in the system, but these models are difficult to implement and computationally cumbersome. Time delays in PEMFC response arise from the delays in the electrochemical reactions as well.

The present work proposes a much simplified PEMFC model which focuses mainly on the reactant behavior and power output. With the single control objective of maintaining a desired oxygen excess ratio, the present control model is a balance between accuracy and response time of the PEMFC to the demand load profile. Several control strategies can be applied to PEMFC. Classical feedback with Proportional Integral (PI) is one of the most common type of control used in fuel cells performance [START_REF] Daud | PEM fuel cell system control: A review[END_REF] . Classical PI or PID control with feedback is generally used to prevent oxygen starvation due to its simplicity and low cost [START_REF] Woo | PEM fuel cell current regulation by fuel feed control[END_REF], [START_REF] Serra | Controllability analysis of decentralised linear controllers for polymeric fuel cells[END_REF]. Sliding mode control (SMC) is a technique derived from the variable structure control and it was initially studied by Utkin [START_REF] Utkin | Variable structure systems with sliding modes[END_REF]. SMC is a simple procedure and a robust controller for linear and nonlinear systems. Additionally, it has obtained good performance in controlling chemicals processes [START_REF] Chen | A sliding mode control scheme for non-minimum phase non-linear uncertain input-delay chemical processes[END_REF]. This paper first proposes a simplified PEMFC model using three state variables and then compares the results of PI and SMC control strategies on the simplified model.

II. MATHEMATICAL MODEL

A. Stoichiometry

The stoichiometry of the combined PEMFC reaction yields the molar ratios of reactant gases and products. The chemical reaction in the PEMFC is caused by hydrogen and oxygen gases and produces electricity, heat and water. The reaction is stated as: 2H2 + O2  2H2O From the above reaction, 2 moles of hydrogen gas needs to react with 1 mole of oxygen gas to produce 2 moles of water. If molar masses are considered, then 4 g of hydrogen reacts with 32 g of oxygen gas to form 36 g of water. This gives a hydrogen to oxygen mass ratio of 1:8 for a complete reaction. The current generated during the PEMFC operation is closely linked to the moles of each gas consumed in the reaction: The stack current is given by:

= = ( 1
)
where and are the masses of the reactant gases in grams, and are the molar masses of oxygen and hydrogen gases, respectively and F is the Faraday constant given in C/mol.

B. Oxygen mass flow

At the cathode, side air is forced into the cathode channels using axial fans. The fans also provide cooling to the PEMFC. The following equation can be obtained by applying the mass conservation at the cathode side, [START_REF] Yang | Control challenges and methodologies in fuel cell vehicle development[END_REF]:

= , - , - , = (2)
By taking oxygen mass is a state variable and noting the load current as the input , then equation 2 can be rewritten as:

, = ,!" #1 - % &'( ∅ * +, -./ (3)
, = (4)
, = 0 !, (5)
= ,!" 11 - % &'( 2 * 3 4 5 0 6 - -0 !, (6)
where represents the oxygen excess ratio which is generally accepted to be around 2. The molar fraction of oxygen in the incoming air is denoted by ,!" ,while 7 " is the saturation pressure of the water at cathode side. The relative humidity is indicated by 8 ! . The inlet side pressure is assumed to be proportional to the input current linearly with the factor 0 , as the gradient. The fan conditions are represented as constant 0 6 . The fan increases its speed based on an increase in load current, and the constants are estimated using characteristic fan curves. Similarly 0 !, represents the effect of the outlet conditions. The constants were estimated empirically by comparing flow rate graphs.

C. Hydrogen mass flow

Similar to the cathode side gas dynamics, the hydrogen gas dynamics can be written by using the net flow of hydrogen where the load current is taken as the input.

= #1 - %: ' %" -0 : - -0 ", (7)
At the cathode 7; " is taken as the partial pressure of water vapor while 7< is the pressure of the anode. A valve constant 0 : is used to estimate the response of the valve to input current. The outlet flow conditions that affect the mass of hydrogen leaving the anode is taken as 0 ", .

D. Water mass flow

Most of the water forms at the cathode and this moisture passes through the layers to reach the anode side. The anode water content is essential for this simplified model. It can be assumed that is the water is less compared to the cathode side, and this can be combined directly into the cathode hydration model. The hydration of the membranes also plays an important role in power production and polarization curves.

where . A , . A , . A=, . A , . AE , . AI , . AH are constants that have been derived empirically. The flow rate of water is taken as 3

x & , ca m is the mass at cathode inlet and .

ca out m
is the mass at cathode outlet. While the hydration equation is part of the simplified model, it does not affect the control of the oxygen excess ratio.

E. Electrical model and voltage losses

The voltage (generated from the electrochemical) potential reaction is determined by the Gibbs Free energy equations [START_REF] Pukrushpan | Modeling and control for PEM fuel cell stack system[END_REF]. Since the reaction takes place across the surface area J K of the MEA, the current produced is proportional to the area of the MEA or commonly stated as the fuel cell area. The current density L is used in determination of voltage losses in the PEMFC:

L = M N OP (9)
where J K is given in [cm 2] and the stack current in [A]. The pressures of hydrogen and oxygen gases are found by using the ideal gas equation. The mass corresponds to the values found by integrating hydrogen and oxygen flow rate equations. [START_REF] Serra | Controllability analysis of decentralised linear controllers for polymeric fuel cells[END_REF] The Nernst voltage or the theoretical voltage potential the PEMFC can develop at a certain pressure of hydrogen and oxygen is given by: E = 1.229 -0.85 × 10 ]= C^-298.15D + 4.3085 × 10 ]E ^alndQ e + ln CQ Df [START_REF] Utkin | Variable structure systems with sliding modes[END_REF] While the Nernst voltage provides the theoretical maximum voltage, there are several voltage losses which occur in the PEMFC to reduce the actual output voltage for the cells:

Q = R S T *' (10) Q = R S T 'U
; 6! = g -C; "! + ; h + ; ! ! D (13)
These losses arise from three areas -activation losses, ohmic losses and concentration losses. The activation losses arise due to energy requirements for breaking the bonds and the sluggish rate of reactions at the electrode surface [20].

; "! = ; + ; " C1 -i ]! ? ) [START_REF] Ram | Model comparison and parameter estimation of polymer exchange membrane (PEM) fuel cell based on nonlinear least squares method[END_REF] where va , vo and c1 are constants that need to be estimated as shown in [START_REF] Ram | Model comparison and parameter estimation of polymer exchange membrane (PEM) fuel cell based on nonlinear least squares method[END_REF]. The resistance of electron flow for conducting electrodes and the resistance of ion flow for the membrane results in ohmic losses:

; h = L j h (15)
where j h is the internal resistance of the PEMFC. The concentration voltage losses ; ! ! arise due to concentration gradients as the reactants get used up at the cell surfaces and the slow transport of reactants to and from the reaction sites. This can be computed as:

; ! ! = L #k l'm - ! n (16)
The reactant flow dynamics in the PEMFC are used in equation 9 to estimate the total theoretical voltage that can be generated by each cell. The reactant gas flow changes in line with the demand current, which can be assumed to be the only input into the system.

F. Model Control Scheme

Two different control topologies are tested on the simplified model. The primary control objective is to ensure that the oxygen excess ratio is maintained at 2. Firstly, a PI like control scheme is applied as shown in where o = 3.65i ] and q = 100 are obtained from equation [START_REF] Pukrushpan | Modeling and control for PEM fuel cell stack system[END_REF] This transfer function has a gain equal to 2, so the oxygen excess ratio € ? € ?stu = will be equal to 2. In addition, the closed loop dynamics can be tuned by choosing proper gains . , and . . For desired closed specifications • K' and ƒ K' :

.L = " P~ … z † P~ (18) . , = … z C + { … z " P~ -qD (19)
Taking ƒ K' = 1 and • K' = 0.1 gives .L = 13.6054 and . , = -1.3605 x 10 5 . Given that the system model has been simplified, a robust control strategy needs to be adopted to ensure optimum control despite losing some details in the model. For this reason, a sliding mode control scheme is applied and then compared with a PI control.

The sliding model control equations can be derived from the state equations. Solving for all constants and assuming = ; for the auxiliary controller, thus

= ‡ ! ]!= C k + + ;D (20)
is obtained from equation ( 6). The first-order sliding mode control topology was developed with the objective of controlling the oxygen excess ratio. The following auxiliary control equation can be stated as:

; = 0 + 0 ˆ 6 + + 0 = ‰LŠ‹CiD ( 21
)
where 0 , 0 and 0 = are gains that were tuned using a Lyapunov analysis. Defining the Lyapunov function candidate as:

OE = i ( 22
)
where i = ˆ 6is the sliding surface and taking its time derivative gives: The constants used in the model are given in the Table 1 below. Parameters for the electrical model are given in [START_REF] Ram | Model comparison and parameter estimation of polymer exchange membrane (PEM) fuel cell based on nonlinear least squares method[END_REF]:

OE = i i OE = -i OE = [C1 -0 D -0 ˆ 6 -0 = ‰LŠ‹CiD]i (

III. MODEL RESULTS AND VALIDATION

Simulink model was run for 100s and the resulting polarization curve was compared to the measured experimental data. The PEMFC used to gather the experimental results had 42 cells with maximum voltage reaching 0.9 V in each cell. The PEMFC area was 62 cm 2 while its membrane thickness was 0.42mm. Parameter estimation was carried out using nonlinear least squares method. The parameter estimation details are stated in [START_REF] Ram | Model comparison and parameter estimation of polymer exchange membrane (PEM) fuel cell based on nonlinear least squares method[END_REF]. There is some delay in flow rate of hydrogen in the system when the load current changes abruptly. For oxygen the delays in following the load current is very high during load changes and even at constant current intervals -the flowrate is poorly maintained. The poor response of oxygen flow without a controller will not achieve the desired stoichiometry and results in power loss. Fig 6 shows how the oxygen excess ratio changes with time. It is noteworthy that between 30s and 40s the excess ratio drops below one. An excess ratio below one means that less oxygen is supplied than what is required to produce the load current. Not only does this result in power loss, but it also causes cell degradation. The response of oxygen flow rate to change in load current is relatively slow and takes almost 1.5s to stabilize as shown in Fig 8. Furthermore, at around 40s when the load changes to 50A, the excess ratio is seen to drop below one which is unacceptable. This would mean the PEMFC is starved of oxygen. For the SMC the gains were calculated as: . = 2000, . = -1999 and . = = 1 10 ]H . The results of the SMC with the objective of maintaining an oxygen excess ratio of 2 is shown in Fig 9 . It can be seen that the excess ratio is reasonably maintained at around 2 with perturbations of less than 0.01 from the reference. This is different from the excess ratio maintained by the PI controller where values below 1 were seen. The Sliding Mode controller performs well to ensure that the oxygen ratio never falls below 1 to avoid oxygen starvation. This control topology, therefore, achieves the objective even on a simplified PEMFC model. The proper control of the excess ratio ensures smooth power output and long life of the individual cells in the PEMFC.

Fig. 9 Oxygen excess ratio controlled using SMC.

The oxygen flow rate response is slower when using the traditional PI control with feedback, as shown in Fig 10. This is one of the reasons why the oxygen excess ratio is difficult to maintain using simple PI control. In an actual system, more delays are anticipated as the oxygen is supplied using an axial fan. Unlike large compressor driven systems, there is no storage of air in a supply manifold, and hence the response of the oxygen supply subsystem is much more critical. was increased by 60% and the SMC controller still maintained the oxygen excess ratio at 2. The SMC controller is especially suited to the nonlinear PEMFC system. Future work will involve modelling of nitrogen, humidity, and temperature and air fan dynamics to further enhance the model with the goal of eventually controlling the PEMFC system using a sliding mode controller.

CONCLUSION

The paper presents a simplified nonlinear dynamic model which uses just three state variables to estimate the changes in PEMFC system. Results show that with accurate parameter estimations, the model is in agreement with measured results for the Fuel Cell polarization curve. Hydrogen and Oxygen flow rates also correlate with the changing current in a test load profile indicating the suitability of the model for control applications. The simplification of the model allows treating some parameters as constants. While there may be loss of details from the model, it is still useful in early system design of PEMFC or for learning purposes. A Sliding Mode Control (SMC) scheme is applied with the key objective of achieving a constant oxygen excess ratio. The results of the SMC show a greater response of oxygen flow rates compared to traditional PI controllers. The study focused only on the oxygen excess ratio, given its slow dynamics compared to hydrogen flow. Future studies will focus on hydrogen gas and water flow management in the PEMFC.

  = = % *' .+ >? . !" + . A B + C. A= . B -. A D. . AE -F.+ >? .+ >G . !", . AH

  Fig 1:

Fig. 1

 1 Fig. 1 PI-like Controller scheme for maintaining oxygen excess ratio.

  From the block diagram of Fig 1, the following equation is derived:

  23)In particular choosing 0 = 1 -0 leads to:OE = C1 -0 Di -0 = ‰LŠ‹CiDi (24)Setting 0 > 1, the following inequality can be obtainedOE ≤ -0 = √2√OE (25)by which OE will reach zero in finite time or that will reach in finite time. The block diagram of the control scheme is shown in Fig.2.

Fig. 2

 2 Fig. 2 Sliding Mode Control scheme for maintaining oxygen excess ratio.

  10 ] ./ = 1.27 × 10 ]E ∅ " = 1 × 10 ] .; = 4.2 × 10 ]H ∅ ! = 0.7 G. Experimental Setup A 1.2 kW Open Cathode PEMFC was loaded using a DC electronic load as shown in Fig 3. Data was collected from the PEMFC user interface module. A single ramp current load was programmed into the electronic load, starting at 2A and ending at 60A. The PEMFC was not able to output a current lower than 2A. Voltage and current was logged in the Electronic load as well. The sampling rate in both cases are 2Hz.

Fig. 3

 3 Fig. 3 The 1.2kW Fuel Cell setup.

Fig. 4

 4 Fig. 4 Experimental and model results compared. Generally, there is good agreement between the polarization curve estimated through the model and the results obtained through experiments. It is essential for a good PEMFC model to accurately predict the flow rate of gasses for changing load currents. A load profile was applied to the model and the change in flow rates of oxygen and hydrogen are shown in Fig 5.

Fig. 5

 5 Fig. 5 Flow rate variation with respect to load current.

Fig. 6

 6 Fig.6The oxygen excess ratio without a controller.For the PI controller the gains were tuned to be . , = -1.3605e5 and .L = 13.6054. The use of a PI controller improves the systems oxygen excess ratio and oxygen flow response. The results in Fig 7clearlyshow the excess ratio stabilizes around the designed excess ratio of 2, but the response to any change in load current is still inadequate.

Fig. 7

 7 Fig. 7 Oxygen excess ratio using PI controller.

Fig. 8

 8 Fig. 8 Magnified view showing the PI response time.

Fig. 10

 10 Fig. 10 Flow rate of oxygen using PI and SMC control. Given the simplified model of the PEMFC, greater emphasis needs to be placed on the controller to account for the loss of details in model simplification. As seen in Fig 11, the SMC controller is resilient to changes in parameters as well unlike the PI controller. In Fig 11, the 0 !,was increased by 60% and the SMC controller still maintained the oxygen excess ratio at 2. The SMC controller is especially suited to the nonlinear PEMFC system. Future work will involve modelling of nitrogen, humidity, and temperature and air fan dynamics to further enhance the model with the goal of eventually controlling the PEMFC system using a sliding mode controller.

Fig. 11

 11 Fig. 11 The SMC controller maintains the excess ratio when parameters change

TABLE I CONSTANTS

 I USED IN THE MODEL


Samuele Zoboli

Vincent Andrieu

Daniele Astolfi

Giacomo Casadei

Jilles S Dibangoye

Madiha Nadri

Reinforcement Learning Policies With Local LQR Guarantees For Nonlinear Discrete-Time Systems

published or not. The documents may come

I. INTRODUCTION

In recent years Deep Reinforcement Learning (RL) received considerable attention due to its ability to autonomously learn how to solve complex high-dimensional control tasks, see, e.g., [START_REF] Silver | A general reinforcement learning algorithm that masters chess, shogi, and Go through selfplay[END_REF], [START_REF] Mnih | Playing atari with deep reinforcement learning[END_REF], [START_REF] Schulman | Proximal policy optimization algorithms[END_REF]. However, most of state-of-the-art algorithms focus on control performances and very few of them study the stability properties in a Lyapunov sense [START_REF] Khalil | Nonlynear systems[END_REF]. This led to a limited use of these methods in realworld applications. Safety requirements or strong uncertainties that cannot be learnt during the training phase eventually reduced their reliability outside of lab environment, due to the unpredictable nature of the learnt controllers.

On the flip side, Control Theory (CT) has a long history of solutions guaranteeing both performances and stability in the linear framework (see, e.g., [START_REF] Kirk | Optimal Control Theory: An Introduction[END_REF]). Unfortunately, it is hard to extend the theory to unstructured nonlinear systems. A common approach is to design a controller for the linearized model around an equilibrium point. However, this technique restricts the stability and performance guarantees to a neighbourhood of the equilibrium point, denoted as domain of attraction.

The last decades saw increasing interest in linking the fields of CT and RL (see, e.g. [START_REF] Perkins | Lyapunov design for safe reinforcement learning[END_REF], [START_REF] Randløv | Combining reinforcement learning with a local control algorithm[END_REF], [START_REF] Kumar | Improving model-based balance controllers using reinforcement learning and adaptive sampling[END_REF]), and the inclusion of Lyapunov theory in learning algorithms received increasing attention. A notable model-based solution was proposed in [START_REF] Berkenkamp | Safe model-based reinforcement learning with stability guarantees[END_REF], following the work in [START_REF] Berkenkamp | Safe learning of regions of attraction for uncertain, nonlinear systems with gaussian processes[END_REF]. However the authors assume a suitable initial safe policy and Lyapunov function are given. In [START_REF] Chow | Lyapunov-based safe policy optimization for continuous control[END_REF], the authors claim near-constraint satisfaction via policy gradient methods by projecting either the policy a Université Claude Bernard Lyon 1, CNRS, LAGEPP UMR 5007, 43 boulevard du 11 novembre 1918, F-69100, Villeurbanne, France.

b INRIA Chroma team, CITI Lab. INSA Lyon, Villeubanne, France. c Laboratoire Ampere Dpt. EEA of the École Centrale de Lyon, Université de Lyon, 69134 Ecully, France.

This research was partially supported by the French Grant ANR DELICIO (ANR-19-CE23-0006).

parameters or the action onto a feasible set described by a Lyapunov constraint. Last, [START_REF] Han | Actor-critic reinforcement learning for control with stability guarantee[END_REF] uses a Lyapunov function as critic for guaranteeing stability in the mean cost of the learnt policy.

In this paper we propose a methodology to integrate basic model knowledge in standard model-free RL algorithms. By enforcing a predefined structure to the control policy we provide theoretical guarantees concerning the local stability of the system via Lyapunov indirect theorem. We focus on the optimal control of nonlinear systems under Linear Quadratic cost functions. The approach enhances the capabilities of the linear controller by learning to operate outside its domain of attraction while maintaining its local properties.

The paper is organized as follows: Section II introduces some preliminary notions and states the problem we aim at solving. A solution is proposed in Section III and the learning behavior is discussed in Section IV. An illustration of an inverted pendulum is given in Section V. Conclusions are drawn in Section VI.

II. PRELIMINARIES A. Dynamic programming concepts

Deterministic dynamic programming deals with discretetime dynamical systems generating a sequence of states (x t ) t∈N , x ∈ R n under the influence of control inputs (u t ) t∈N , u ∈ R m . The objective is to find a sequence of control inputs, often denoted as control policy u = (u t ) t∈N , which leads to the minimum of a γ-discounted function (also known as state-value function under u), J u : R n → R, defined as J u (x) = ∞ t=0 γ t r(x t , u t ), where γ ∈ (0, 1). See, e.g., [START_REF] Bertsekas | Reinforcement Learning and Optimal Control[END_REF]. Hence, the state-value function under an optimal policy is the optimal state-value function J , i.e. the one that provides the minimum cost for each state. In the RL framework, the action-value function,

Q u : R n × R m → R, Q u (x t , u t ) = r(x t , u t )+γJ u (x t+1
), turns out to be more practical when dealing with unknown dynamics.

B. Infinite horizon LQR

Consider a deterministic discrete-time linear system

x(t + 1) = Ax(t) + Bu(t), t ∈ N,
with state x ∈ R n and control input u ∈ R m . For the sake of convenience, we denote x t := x(t) the state x at discretetime t. The same notation will adopted for all variables and inputs, obtaining the compact notation

x t+1 = Ax t + Bu t . ( 1
)
We suppose that the pair (A, B) is stabilizable. The infinitehorizon LQR problem asks for a stabilizing controller optimal with respect to a quadratic cost function. By letting γ ∈ (0, 1] be a discount factor, the discounted problem considers the following cost function

J γ (x, u) = ∞ t=0 γ t (x t Q γ x t + u t R γ u t ), (2)
where

Q γ = Q γ ≥ 0 and R γ = R γ > 0.
In [START_REF] Bertsekas | Dynamic Programming and Optimal Control 4th Edition, Volume II Chapter 4 Noncontractive Total Cost Problems[END_REF]Chapter 4.3] it is shown that, for a given discount factor γ, the optimal controller u * for (1) with cost (2) is given as

u t = K γ x t , K γ = -γ(R γ + γB P γ B) -1 B P γ A, ( 3
)
where K γ is the optimal discounted gain and P γ is the unique solution of the discounted Discrete-time Algebraic Riccati Equation (DARE)

P γ = Q γ + γA (P γ -γP γ B(R γ + γB P γ B) -1 B P γ )A. (4
) Equation ( 4) can be rewritten and solved as a standard DARE for the new matrices A = √ γA, B = √ γB. From [5, Section 3], the value functions for the discounted LQR problem under the optimal controller are

J γ (x t ) = x t P γ x t , Q γ (x t , u t ) = z t H γ z t , (5)
with z t := col(x t , u t ) ∈ R n+m and H γ given by

H γ = Q γ + γA P γ A γA P γ B γB P γ A R γ + γB P γ B .
In the discounted framework, closed-loop stability is dependent on the choice of the discount factor γ. In [17, Section IV, Corollary 3], the authors defined a conservative lower bound γ * ∈ (0, 1], depending on Q γ and R γ , such that for any γ ∈ (γ * , 1] the origin of the closed-loop (1), ( 3) is exponentially stable. Note that for γ = 1, we recover the so-called "undiscounted LQR" problem in which the cost function reads

J(x, u) := J 1 (x, u) = ∞ t=0 x t Qx t + u t Ru t , (6)
with optimal solution u given by

u t = K x t , K = -(R + B P B) -1 B P A, (7)
with P the symmetric positive definite matrix solution of the DARE

P = A P A -A P B(R + B P B) -1 B P A + Q. ( 8
)
For the undiscounted problem, stability of the closed-loop system is always guaranteed.

C. Actor-Critic algorithms

Many state-of-the-art RL algorithms use Actor-Critic architectures [START_REF] Konda | Actor-critic algorithms[END_REF] (see Figure 1). A key feature of the RL framework is that the environment (plant) returns a reward r : R n × R m → R, which embeds information about the quality of the chosen action. In Actor-Critic algorithms, this signal is used to update the parameter vector φ ∈ R q of the Fig. 1: Actor-Critic structure critic value function estimator Qφ , according to a criterion J φ : R n × R m × R q → R. The parameter vector θ ∈ R p of the actor policy π θ are subsequently updated according to a different cost function J θ : R n × R m × R p → R depending on the new critic estimation. We can then define the class of Actor-Critic algorithms A as the set of operators A := {a : R p × R q → R p × R q } with input (θ 0 , φ 0 ) and output (θ κ, φ κ) such that, for any κ ∈ [0, κ] ⊆ N, the parameters φ and θ are updated via Stochastic Gradient Descent [START_REF] Sutton | Reinforcement learning: An introduction[END_REF]Chapter 9.3]. The classic formulation imposes

θ κ+1 = θ κ -λ θ ∇ θ J θ and φ κ+1 = φ κ -λ φ ∇ φ J φ ,
where λ φ , λ θ are the positive learning rates, subscript κ identifies a parameter vector at iteration κ, ∇ θ is the gradient with respect to the policy parameters, ∇ φ is the gradient with respect to the critic value function parameters and κ ∈ N is the iteration index such that ∇ θ J θ ≤ and ∇ φ J φ ≤ for some small > 0 and for all κ ≥ κ. One of the most common objective for updating the critic parameters asks for the minimization of

J φ (x t , u t ) = ∞ t=0 r(x t , u t )+ Qφ (x t+1 , u t+1 )-Qφ (x t , u t ) 2 . (9)
Most update procedures for algorithms learning a deterministic policy exploit deterministic policy gradient theorem [START_REF] Silver | Deterministic policy gradient algorithms[END_REF]. Given a parametrized deterministic policy π θ and a set S ⊂ R n , deterministic policy gradient updates take the form

∇ θ J θ (x t , u t ) = S d π (s)∇ θ π θ (s)∇ a Qφ (s, a)ds, ( 10
)
where ∇ a the gradient with respect to the actions a = π θ (s) and d π (s) := S ∞ t=1 γ t-1 Pr(x 0 ) Pr(x 0 → s, t, π θ )ds is the discounted state distribution, being Pr(x 0 → s, t, π θ ) the probability of reaching the state s ∈ S after transitioning for t time steps with initial condition x 0 , following the policy π θ .

Remark. In the paper we will focus on algorithms learning the parameter vector θ of a deterministic parametrized policy π θ : R p → R m via [START_REF] Kirk | Optimal Control Theory: An Introduction[END_REF]. However the proposed results can be directly extended to include methods exploiting stochastic policy gradient updates since the proposed solution is independent from the choice of the RL algorithm.

D. Problem statement

Consider a deterministic discrete-time nonlinear system

x t+1 = f (x t , u t ), (11)
with state solution (x t ) t∈N taking values in R n , control input (u t ) t∈N taking values in R (m = 1), and f : R n × R → R n being a continuously differentiable function in a neighbourhood of the origin. Without loss of generality we suppose that f (0, 0) = 0 and we denote A := ∂f ∂x (0, 0) and B := ∂f ∂u (0, 0). Problem 1. Let J be an undiscounted cost function of the form [START_REF] Chow | Lyapunov-based safe policy optimization for continuous control[END_REF] and assume a linearized model of the form (1) for (11) is known, with (A, B) stabilizable. The goal is to learn an optimal parametrized control policy π θ : R n ×R p → U ⊆ R with parameters θ ∈ R p such that the origin of the closedloop system [START_REF] Konda | Actor-critic algorithms[END_REF] with u t = π θ (x t ) is locally asymptotically stable for all θ ∈ R p , namely

∂π θ ∂x (0) = K , ∀ θ ∈ R p , (12)
where K ∈ R n is the LQR optimal gain given in [START_REF] Fujimoto | Addressing function approximation error in actor-critic methods[END_REF].

Remark. Solving the optimal problem implies requirement [START_REF] Königsberger | Analysis 2[END_REF] since K is the unique optimal solution to the local problem. Moreover [START_REF] Königsberger | Analysis 2[END_REF] ensures local stability around the origin by Lyapunov indirect theorem, being K stabilizing for the linearized system.

III. ALGORITHM DESIGN

To solve the problem we rely on deep reinforcement learning algorithms. Hence, we design a suitable cost setting the learning objective. Afterwards, we enforce a specific structure on the learnt policy to guarantee that the approximation constraint (12) is satisfied.

A. Reward shaping

Rewards play a fundamental role in RL. Most RL algorithms require the value function J π to have a finite value in order to converge. However, when defined as an infinite sum over time, boundedness of J π may not be always ensured. Hence, starting from J in (6) we look for a reward function r : R n × R → R and a discount factor γ ∈ (0, 1) defining a suitable γ-discounted function

J γ,RL (x, u) = ∞ t=0 γ t r(x t , u t ) (13)
which sets the learning objective for the agent.

In the following Lemma, we show that given any undiscounted problem of the form ( 6)-( 7), we can always redefine an associated discounted problem of the form (2)-(3) so that the optimal gains (7) and (3) coincide.

Lemma 1. Consider system (1) and an associated undiscounted optimal control problem of the form [START_REF] Chow | Lyapunov-based safe policy optimization for continuous control[END_REF]. Then, for any γ ∈ (0, 1], the optimal gain K defined in [START_REF] Fujimoto | Addressing function approximation error in actor-critic methods[END_REF] is the optimal solution of the discounted problem (2) with Q γ , R γ defined as

Q γ = γQ + (1 -γ)P, R γ = γR. ( 14
)
Moreover, the state-value function (5) of the discounted problem is finite.

Proof. See Appendix A.

Based on Lemma 1, the undiscounted cost J in ( 6) can be replaced by its discounted version (2) by using the weights [START_REF] Lillicrap | Continuous control with deep reinforcement learning[END_REF]. Thus, we obtain now a suitable form of the cost J γ,RL for RL algorithms without affecting the local optimal solution to the problem. Consequently, the reward r from (13) can be set as

r(x t , u t ) = x t Q γ x t + u t R γ u t , (15)
with Q γ and R γ defined as in [START_REF] Lillicrap | Continuous control with deep reinforcement learning[END_REF]. The value of γ is a free parameter to be chosen as the most suitable one for the algorithm convergence.

B. Control policy design

Our goal is to ensure local stability of the closed-loop system [START_REF] Konda | Actor-critic algorithms[END_REF] with u t = π θ (x t ), i.e. to satisfy the local approximation constraint (12) independently from the parameter vector θ. To this end, we enforce a specific structure in our control policy. Given any k ∈ N, we define H k as the set of continuously differentiable functions h satisfying the following

H k := h ∈ C k+1 : R n → R : lim s→0 h(sx) s j |x| j = 0, j ∈ [0, k], lim s→0 h(sx) s k+1 |x| k+1 = 0, lim |s|→∞ h(sx) = 1, ∀x ∈ R n .
The proposed control law is then designed as

π θ (x t ) = π loc (x t ) + π θ learn (x t ) π loc (x t ) = K x t π θ learn (x t ) = h 1 (x t ) µ θ (x t ) -π loc (x t ) , (16)
where h 1 ∈ H 1 and µ θ : R n → R is a scalar function to be learnt by the RL agent which is parametrized by the set of parameters θ ∈ R p and satisfies the following assumption:

Assumption 1. The function µ θ is locally Lipschitz.

Remark. Since we consider deterministic policy gradient methods (e.g. DDPG [START_REF] Lillicrap | Continuous control with deep reinforcement learning[END_REF] and TD3 [START_REF] Fujimoto | Addressing function approximation error in actor-critic methods[END_REF]), the function µ θ denotes the parametrized approximator to be trained. Note that most NNs are locally Lipschitz since they are compositions of locally Lipschitz functions.

Remark. The size of the guaranteed domain of attraction for policy ( 16) can be controlled by shaping the function h 1 .

The effect of the learnt component in a neighbourhood of the origin is scaled by such a function. Hence, it is possible to strongly reduce the contribution of π θ learn (x t ) in regions where we trust the LQR controller π loc to stabilize the system.

In order to justify our choice for the control policy, in the following Proposition we prove that, by enforcing structure [START_REF] Perkins | Lyapunov design for safe reinforcement learning[END_REF] in the control policy and assuming µ θ satisfies Assumption 1, we can learn an arbitrary C 2 policy satisfying the problem constraint [START_REF] Königsberger | Analysis 2[END_REF]. Proposition 1. The following statements hold.

1) Given any K ∈ R n and any function π θ : R n → R, π θ ∈ C 2 satisfying the local approximation constraint [START_REF] Königsberger | Analysis 2[END_REF], for any function h 1 ∈ H 1 there always exists a locally Lipschitz function µ θ : R n → R satisfying the equality (16) for all x ∈ R n . 2) Let µ θ satisfy Assumption 1, h 1 ∈ H 1 and π θ : R n → R as in [START_REF] Perkins | Lyapunov design for safe reinforcement learning[END_REF]. Then (12) is satisfied.

Proof. See Appendix B.

Remark. Due to [START_REF] Perkins | Lyapunov design for safe reinforcement learning[END_REF] being always enforced, the local guaranteed properties of classical LQR are ensured even for the untrained policy.

C. Main Result

In this section, we established how to design the objective for the RL agent and how to structure the policy to be learnt. Now we present our solution to Problem 1. It can be solved via Lemma 1 and Proposition 1 by letting Assumption 1 hold.

Theorem 1. Let be given any algorithm a ∈ A. Consider the cost (13) with the reward function r(x t , u t ) given by [START_REF] Mnih | Playing atari with deep reinforcement learning[END_REF]. Then by selecting the the control policy π θ as in ( 16) and determining its parameter vector θ via a, Problem 1 is solved.

Proof. The reward function defined in [START_REF] Mnih | Playing atari with deep reinforcement learning[END_REF] allows us to use the discounted cost [START_REF] Kumar | Improving model-based balance controllers using reinforcement learning and adaptive sampling[END_REF] without affecting the local solution K thanks to Lemma 1, and without restricting our choice of γ ∈ (0, 1]. Then ( 16) ensures constraint ( 12) is always satisfied via Proposition 1 and Assumption 1. Finally a ∈ A allows the learning of the locally optimal policy parameters.

The proposed solution is independent from the choice of the deep reinforcement learning algorithm and it can be applied to a wide variety of existing model-free solutions. However we focused on Actor-Critic algorithms, since many of state-of-the art methods exploit the Actor-Critic architecture. Appendix C presents the procedure for learning the controller via an Actor-Critic algorithm.

IV. LEARNING THE POLICY

In the following section we study the behavior of the proposed solution during the training process. We also propose to structure the value function estimator in order to ensure the correctness of the estimation in the origin.

A. Training behavior

For exploring the training behavior it is sufficient to combine the update law [START_REF] Kirk | Optimal Control Theory: An Introduction[END_REF] and the policy [START_REF] Perkins | Lyapunov design for safe reinforcement learning[END_REF]. Plugging the equality [START_REF] Perkins | Lyapunov design for safe reinforcement learning[END_REF] in the deterministic policy gradient equation [START_REF] Kirk | Optimal Control Theory: An Introduction[END_REF] from Section II highlights that

∇ θ π θ (x t ) = h 1 (x t )∇ θ µ θ (x t ), (17)
being µ θ (x t ) the only term depending on the parameters. Equation [START_REF] Postoyan | Stability analysis of discrete-time infinite-horizon optimal control with discounted cost[END_REF] shows that the closer the system gets to the equilibrium point the smaller the updates become, since h 1 ∈ H 1 . Note that one can substitute the policy ( 16) in the nonlinear model ( 11) and obtain a new system under the learnt input µ θ (x t ). Hence, the exploration noise which is typically added to deterministic policies can be applied directly on µ θ (x t ). By doing so, during training π θ (x t ) is a random variable with

E π θ (x t ) = (1 -h 1 (x t ))K x t + h 1 (x t )E µ θ (x t ) , Var π θ (x t ) = h 1 (x t ) 2 Var µ θ (x t ) ,
being µ θ (x t ) the only random component. Since h 1 ∈ H 1 , the variance decreases the closer the system is to the equilibrium. Finally, the study suggests that each training episode should be concluded once the equilibrium of the state-space (i.e. an arbitrarily small neighbourhood of it) is reached. This is in accordance with the fact that the agent is actually learning only how to steer the system to the equilibrium point and not how to keep it there, being the solution to the latter problem already provided by the local controller π loc .

B. Improving learning by reshaping the value function estimation

For designing critics, model-free RL algorithms usually rely on value functions in order to drive the policy towards the optimal solution. We propose a value function estimator built on the knowledge of the policy ( 16) and of the linear framework. Suppose we are interested in the action-value function under the control policy [START_REF] Perkins | Lyapunov design for safe reinforcement learning[END_REF]. Let us denote the its estimator by Qφ π (x t , u t ). Due to the structure of π θ , equation ( 5) provides a suitable local approximation. Consequently, we impose the following local constraint

∀z t ∈ Z, Z = {z t ∈ R n+1 : z t = (0, u t ), ∀u t ∈ U} D i Qφ π (z t ) = D i Q γ (z t ), |i| ≤ 2 , (18)
with D i denoting the derivative of order i with the multiindex notation for multi-variable functions (see, Appendix D). Similarly to the policy equation ( 16), the estimated action-value function Qφ π (x t , u t ) for the nonlinear system (11) under ( 16) is modeled as

Qφ π (x t , u t ) = Q loc (x t , u t ) + Qφ learn (x t , u t ) Q loc (x t , u t ) = Q γ (x t , u t ) (19) Qφ learn (x t , u t ) = h 2 (x t )(Ω φ (x t , u t ) -Q loc (x t , u t ))
, where h 2 ∈ H 2 , Ω φ : R n × R → R is a parametrized function whose parameter vector φ ∈ R q is learnt by the RL agent and satisfying the following assumption: Assumption 2. The function Ω φ is locally Lipschitz.

As for the control policy, we justify the design choice [START_REF] Randløv | Combining reinforcement learning with a local control algorithm[END_REF]

∈ R n , ∀u ∈ U ⊆ R.
2) Let Ω φ satisfy Assumption 2, h 2 ∈ H 2 and Qφ π : R n+1 → R as in [START_REF] Randløv | Combining reinforcement learning with a local control algorithm[END_REF]. Then (18) is satisfied.

Proof. The proof follows the same steps performed in the proof of Proposition 1. Note that [START_REF] Raffin | Stable Baselines3[END_REF] holds for all points in Z instead of the single point scenario of Proposition 1. Moreover, we match the second order approximation via h 2 ∈ H 2 .

V. SIMULATION RESULTS

We run simulations in a frictionless inverted pendulum environment. The system nonlinear model is

α t+1 = α t + ω t ∆t ω t+1 = ω t - 3g 2 sin(α t + π)∆t + 3 2m 2 [sat(u t ) + d t ]∆t,
where α ∈ [-π, π) is the angle between the position of the pendulum and the top vertical one, ω ∈ R is its rate of change, sat(•) is a saturation function limiting the control input torque u t in [-2, 2], d t is wind disturbance affecting the system, ∆t = 0.05 is the discretization step, g = 10 is the gravitational acceleration, = 1 is the length of the pendulum and m = 1 is its mass. We denote x t := (α t , ω t ) . The goal is to stabilize the pendulum at the unstable equilibrium (x , u ) = (0, 0), corresponding to the top vertical position, starting from any random initial condition. The parameters of the discounted LQR problem for the linearized system are We train the proposed policy using a state-of-the-art deep reinforcement learning algorithm, namely TD3 [START_REF] Fujimoto | Addressing function approximation error in actor-critic methods[END_REF]. The nominal version of the learning algorithm comes from Stable Baselines 3 library [START_REF] Raffin | Stable Baselines3[END_REF] and it is adapted to include our policy ( 16) and value function [START_REF] Randløv | Combining reinforcement learning with a local control algorithm[END_REF], following the steps in Appendix C. Since most of RL algorithms are designed to maximize the expected reward, we simply invert the sign of r(x t , u t ) and Q loc (x t , u t ). Each training episode is stopped as soon as the state enters a small circular area of radius 10 -5 centered in the origin, or their time limit is reached. We use the functions h 1 (x) = tanh(tanh -1 (0.99) x P x c ) and h 2 (x) = tanh(tanh -1 (0.99)( x P x c )

2 ), which saturate outside the Lyapunov level set {V (x) = c}, with c ∈ R >0 . The value of c is estimated by sampling random initial conditions and testing the convergence to the equilibrium point, see Figure 2a, and we select c = 0.47.

We evaluate the performances and stability of our solution by comparing the results with the simple LQR and the nominal version of the algorithm. Standard TD3 is trained with the same reward function [START_REF] Mnih | Playing atari with deep reinforcement learning[END_REF]. In order to allow good training, episodes are not stopped as soon as a neighbourhood of the origin is reached, but at end of their time length limit. We use the same hyperparameters for both the standard and the locally asymptotically stable version of the algorithm: 3layers fully-connected NNs with ReLU activation functions, 64 units for the first hidden layer and 32 units as the second one, both for the critic and for the actor, and a learning rate of λ = 0.00371.

We first analyze the learning performances. For each training episode d t = 0, ∀t ∈ R and the initial condition is randomly sampled. From Figure 2b we can infer that classical LQR succeeds only if the initial state lies inside it's domain of attraction. Its performances strongly fluctuates since the magnitude of the return of an episode (i.e. the final cost of the episode) is very big (bad) if the initial condition lies outside of its domain of attraction or quite small (good) if the system starts close to the equilibrium. On the other hand, our solution behaves comparably to the nominal TD3 algorithm.

We also study the stability of the closed loop system. During training, we periodically evaluate the policy by running an experiment in a different environment corrupted by uncertainties. We simulate measurement noise u t = π θ (x t + w t ), w t ∼ N (0, 0.03) and mass mismatch m real = 1.2m train . Moreover, the input is perturbed by external sinusoidal wind d t = 0.36 sin( π 50 t). A "stability score" is extracted as the maximum of the norm of the error vector at steady state. We run the stability evaluation episodes twice, at first with initial condition x 0 = (0.945π, 0) , close to the position "down" of the pendulum, then with x 0 = (0, 0) , corresponding to the desired equilibrium point. Figure 3 clearly shows that the RL approach focuses solely on performances. Moreover, we note that as h 1 (x) → 1, the learnt component becomes predominant, possibly affecting stability. We presented a solution addressing the design of LQR controller for nonlinear system. By uniting basic knowledge on the system behavior around an equilibrium point and deep reinforcement learning techniques we provided a controller ensuring classical local LQR guarantees and experimental global attractiveness. The proposed method is directly implementable starting from standard RL algorithms. Even if the procedure is proposed for single input systems, it can be straightforwardly expanded to higher dimensions. Simulations proved our solution improves both the classical control via linearization and standard RL algorithms.

  γ = 0.001, γ = 0.99, K = -19.3006 -5.9918 , P = 8.088 0.4782 0.4782 0.1624 .

  (a) Estimated level set for saturated LQR (b) Episode return during training

Fig. 3 :

 3 Fig. 3: Stability of policies during training (smoothed)

  via the following Proposition, showing that (19) and Assumption 2 allow learning a generic C 3 function satisfying the local constraint (18). Given any Q loc ∈: R n+1 → R, Q loc ∈ C 3 and any function Qφ π : R n+1 → R, Qφ π ∈ C 3 satisfying the local approximation constraint (18) for any h 2 ∈ H 2 , there always exists a locally Lipschitz function Ω φ : R n+1 → R satisfying the equality (19) for all x

	Proposition 2. The following statements hold.
	1)

APPENDIX

A. Proof of Lemma 1

If γ = 1 the proof is trivial, hence we will focus on γ ∈ (0, 1). First, note that since Q, P and R are symmetric positive definite matrices, so are Q γ and R γ for any γ ∈ (0, 1). By multiplying both sides of the DARE (8) by γ and by using [START_REF] Lillicrap | Continuous control with deep reinforcement learning[END_REF], we get P = Q γ + γA P A -γ 2 A P B(R γ + γB P B) -1 B P A.

Under standard LQR problem assumptions on stabilizability and observability of [START_REF] Berkenkamp | Safe learning of regions of attraction for uncertain, nonlinear systems with gaussian processes[END_REF], the solution to the Riccati equation is unique. Hence, from the definition of the discounted DARE (4) we can conclude that P = P γ . By inserting the latter and ( 14) in the optimal gain [START_REF] Fujimoto | Addressing function approximation error in actor-critic methods[END_REF] we obtain

which shows that ( 7) is also solution to the discounted problem defined by (2) when considering weights in the form of ( 14). The state-value function for a given policy is defined as the cost-to-go starting from an initial state and following the policy. Given a quadratic cost function as in [START_REF] Chow | Lyapunov-based safe policy optimization for continuous control[END_REF], it holds that

By plugging the matrices ( 14) in [START_REF] Schulman | Proximal policy optimization algorithms[END_REF] and considering the value function ( 5) with [START_REF] Han | Actor-critic reinforcement learning for control with stability guarantee[END_REF] for the undiscounted problem under the optimal input (7), yields

Then, by solving [START_REF] Silver | A general reinforcement learning algorithm that masters chess, shogi, and Go through selfplay[END_REF] we have

where J (x 0 ) ∈ R denotes the initial condition. Since γ ∈ (0, 1) by letting t → ∞ equation ( 22) simplifies in

Being J (x 0 ) finite for the linear system (1) under the optimal controller [START_REF] Fujimoto | Addressing function approximation error in actor-critic methods[END_REF], the state-value function for the discounted problem under the undiscounted solution ( 7) is also finite and this concludes the proof.

B. Proof of Proposition 1

Let us address point 1. By keeping in mind the multi-index properties provided in Appendix D, we consider a multiindex i ∈ N n . Constraint [START_REF] Königsberger | Analysis 2[END_REF] implies

It is possible to find a suitable definition of π θ learn satisfying equation ( 16) by recalling that h 1 ∈ H 1 . This shows the first item of the proposition. We prove now the second item. If ( 16) is satisfied ∀x t ∈ R n , constraint [START_REF] Königsberger | Analysis 2[END_REF] is verified if and only if

where s ∈ R. If π θ learn is a locally Lipschitz function, then lim s→0 π θ learn (sx t ) ≤ ω π for all x t ∈ R n , being ω π ∈ R ≥0 . This implies that (25) holds if h 1 ∈ H 1 and this concludes the proof.

C. Locally Asymptotically Stable Actor-Critic

Algorithm 1 presents the procedure to implement the proposed solution with an Actor-Critic algorithm. It takes as input the undiscounted problem formulation, the linearized matrices, an Actor-Critic algorithm and a discount factor. Successively, it computes the linear optimal solution and the associated discounted problem and sets the saturation functions enforcing the local stability. Finally, it sets up the structured policy and value function(s) and runs the Actor-Critic algorithm for learning the parameters. The output is a learnt optimal locally asymptotically stable deterministic policy.

Set r as in [START_REF] Mnih | Playing atari with deep reinforcement learning[END_REF]; Set π θ as in ( 16) with parameters θ 0 ; Set Qφ π as in [START_REF] Randløv | Combining reinforcement learning with a local control algorithm[END_REF] with parameters φ 0 ; Run a(θ 0 , φ 0 ); Result: LAS control policy π θ

D. Matching approximations

In this section we will identify by x i the i th component of vector x ∈ R n , instead of identifying a time index. Firstly, let us recall some multi-index definitions. For a general multiindex α ∈ N n we denote

1 . . . x αn n , for any x ∈ R n . Then, given a function Φ : R n → R whose l-th order partial derivatives are continuous, it is possible to define its derivative of order l as

Finally, we define with o(x ν+1 ) the standard little-o notation for functions of order smaller than |x| ν . We state now the following lemma.

Lemma 2. For a given point y ∈ R n and two real-valued functions

Proof. Introduce the multi-indices i, j ∈ N n , i = (i 1 , i 2 , . . . , i n ), j = (j 1 , j 2 , . . . , j n ) . Consider an arbitrary real-valued function Φ : R n → R, Φ ∈ C l+1 , l ∈ N. Given a point y ∈ R n , by Taylor's theorem for multivariate functions [START_REF] Königsberger | Analysis 2[END_REF] it holds that

0 (1-t) |j|-1 D j Φ(y +t(x-y))dt. Then, if Ψ 1 , Ψ 2 ∈ C ν+1 they can be equivalently expressed as

).

If D i Ψ 1 (y) = D i Ψ 2 (y), |i| ≤ ν we can rearrange the last identity to obtain

Finally, by combining all previous equation we obtain Ψ 1 (x) = Ψ 2 (x) + o(x ν+2 ), thus completing the proof.