Stanislav Sokolenko
email: stanislav.sokolenko@dal.ca

Tangi Jézéquel

Ghina Hajjar

Jonathan Farjon

Serge Akoka

Patrick Giraudeau

Robust 1D NMR lineshape fitting using real and imaginary data in the frequency domain

Keywords: NMR, lineshape, fitting, quantitative, frequency, imaginary

Quantitative NMR is intrinsically dependent on precise, accurate, and robust peak area calculation. In this work, we demonstrate how the use of complexvalued peak descriptions can improve peak fitting in the frequency domainincorporating phase and baseline correction as well as apodization while working with commonly used Fourier-transformed data. The method has been implemented in an open source R package called rnmrfit that is available for download on GitHub (https://github.com/ssokolen/rnmrfit). Application to real data suggests that this approach can also result in dramatically higher precision than can be achieved with existing software. Simulation data indicates that coefficients of variation below 0.1% can be readily achieved at signal to noise (SNR) ratios of approximately 100. The use of complex-valued data in the frequency domain is demonstrated as a relatively simple and effective means of improving peak fitting for quantitative NMR analysis.

Introduction

The calculation of NMR resonance peak areas is a key aspect of NMR analysis and forms the basis of so-called "quantitative" NMR (typically referred to as qNMR). This quantification relies on the basic principle that the area of a resonance is proportional to the number of nuclei that corresponds to that resonance [START_REF] Bharti | Quantitative 1H NMR spectroscopy[END_REF], allowing the calculation of relative compound concentrations from ratios of peak areas (and the calculation of absolute concentrations with the addition of an internal standard). qNMR has been used across a wide array of applications that include metabolomics [START_REF] Wishart | Quantitative metabolomics using NMR[END_REF], natural product analysis [START_REF] Pauli | Quantitative 1H NMR: development and potential of a method for natural products analysis[END_REF], drug discovery and development [START_REF] Singh | The application of absolute quantitative (1)H NMR spectroscopy in drug discovery and development[END_REF], as well as isotopic analysis (the most demanding application as it requires coefficients of variation less than 0.1%) [START_REF] Jézéquel | The new face of isotopic NMR at natural abundance[END_REF]. And whereas there are considerable differences across these different fields, one thing they all have in common is the need for precise and accurate peak area calculation.

The basic method of area quantification is the integration of the NMR lineshape between two chemical shifts. Malz and Jancke [START_REF] Malz | Validation of quantitative NMR[END_REF] have shown that the simple integration of isolated peaks can achieve 95% confidence intervals as low as 1.5% of the compound concentrations. However, particularly in the analysis of natural products and biofluids typical of metabolomics, isolated peaks may be few and far between. Furthermore, the same researchers [START_REF] Malz | Validation of quantitative NMR[END_REF] have also commented that the integration of even well isolated peaks can result in relative uncertainties of 11% with even "slightly wrong" phase and baseline corrections [START_REF] Malz | Validation of quantitative NMR[END_REF]. The general alternative to direct integration of the spectrum is some form of lineshape fitting with the overall goal of separating the contribution of different chemical species from each other and confounding components such as baseline and phase error. Roughly speaking, this process can be divided into two steps -the identification of the chemical species (or residues) that make up the observed resonance and their subsequent quantification. In many cases, the compounds or residues to be quantified are known in advance -making the final step of quantification the more general of the two. And in many applications of qNMR, the basic challenge is how to achieve the most precise quantification of resonance peaks [START_REF] Jézéquel | The new face of isotopic NMR at natural abundance[END_REF][START_REF] Giraudeau | Challenges and perspectives in quantitative NMR[END_REF].

From a mathematical perspective, one of the basic divisions of existing lineshape fitting approaches is whether the NMR lineshape is modelled in the time or frequency domain. In the time domain, an NMR resonance is typically modelled as a decaying sinusoidal function (more specifically, a complex-valued exponential) with constraints added to account for multiplet relations. The use of Bayesian inference to estimate the number of resonances as well as resonance frequency, decay rate, and intensity has been described as early as 1990 [START_REF] Bretthorst | Bayesian analysis. I. Parameter estimation using quadrature NMR models[END_REF][START_REF] Bretthorst | Bayesian analysis. II. Signal detection and model selection[END_REF][START_REF] Bretthorst | Bayesian analysis. III. Applications to NMR signal detection, model selection, and parameter estimation[END_REF]. Some of the more recent advances include the improvement of the Markov Chain Monte Carlo technique used to draw samples [START_REF] Rubtsov | Time-domain Bayesian detection and estimation of noisy damped sinusoidal signals applied to NMR spectroscopy[END_REF][START_REF] Rubtsov | Application of a Bayesian deconvolution approach for high-resolution (1)H NMR spectra to assessing the metabolic effects of acute phenobarbital exposure in liver tissue[END_REF], the implementation of downsampling for isolating specific regions of interest [START_REF] Krishnamurthy | CRAFT (complete reduction to amplitude frequency table) -robust and time-efficient Bayesian approach for quantitative mixture analysis by NMR[END_REF], and the introduction of more generalized lineshapes to account for distortions [START_REF] Matviychuk | An experimental validation of a Bayesian model for quantification in NMR spectroscopy[END_REF]. In the frequency domain (following the Fourier transform of the FID), lineshapes are typically fit as Lorentz peaks using real-valued data following phase correction. Most approaches rely on a database of compounds, whether generated in-house or available online [START_REF] Mercier | Towards automatic metabolomic profiling of high-resolution one-dimensional proton NMR spectra[END_REF][START_REF] Zheng | Identification and quantification of metabolites in (1)H NMR spectra by Bayesian model selection[END_REF][START_REF] Hao | BATMAN -an R package for the automated quantification of metabolites from nuclear magnetic resonance spectra using a Bayesian model[END_REF][START_REF] Ravanbakhsh | Accurate, fully-automated NMR spectral profiling for metabolomics[END_REF], and many of these are also based on Bayesian inference [START_REF] Zheng | Identification and quantification of metabolites in (1)H NMR spectra by Bayesian model selection[END_REF][START_REF] Hao | BATMAN -an R package for the automated quantification of metabolites from nuclear magnetic resonance spectra using a Bayesian model[END_REF][START_REF] Ravanbakhsh | Accurate, fully-automated NMR spectral profiling for metabolomics[END_REF]. The use of libraries can be very effective for standardized samples and pulse sequences (such as those used in metabolomics), but their application is limited to specific combinations of compounds, biological matrices, and pulse sequences on which the database is built. However, individual peaks (with constraints to define multiplet relations) have also been fit using classical non-linear optimization techniques based on the least-squares fit [START_REF] Mihaleva | Automated quantum mechanical total line shape fitting model for quantitative NMR-based profiling of human serum metabolites[END_REF][START_REF] Bernstein | Optimization and automation of quantitative NMR data extraction[END_REF]. Some general trade-offs can be identified when considering the various approaches. Lineshape fitting in the time domain has the advantage of being more general, with the use of complex-valued data allowing for built-in phase correction. However, even a relatively crowded NMR spectrum contains large areas of relatively signal-free regions in the frequency domain that only contribute noise to an overall fit in the time domain. Although it is possible to downsample [START_REF] Krishnamurthy | CRAFT (complete reduction to amplitude frequency table) -robust and time-efficient Bayesian approach for quantitative mixture analysis by NMR[END_REF],

this adds extra processing steps and potential for artefacts, whereas the Fourier transform is already a well accepted means of isolating frequencies of interest and disregarding more problematic ones. Fitting in the frequency domain can be used to target specific peaks of interest while greatly speeding up the fit process (by simply considering fewer points). There is a similar trade-off between Bayesian and classical non-linear optimization techniques. Bayesian approaches are more general as they seek out a global solution whereas optimization based on gradient descent or related techniques require good initial guesses to avoid getting trapped in local solutions that are not globally optimal. On the other hand, Bayesian fitting methods are also considerably slower for high precision calculations.

The method proposed in this work bridges some of these trade-offs by considering a more generalized form of fitting in the frequency domain. One way that is accomplished is by describing NMR peaks as complex-valued functions.

There are a number of practical and theoretical advantages to this approach.

First, phase correction can be implemented as part of the fitting process in the frequency domain, thus avoiding user-dependent phasing that affects the accuracy and precision of peak area determination. Second, baseline distortions can be better separated from the peak shapes as one set of peak parameters is used to describe two different shapes (as they appear in the real and imaginary components). Third, the broader shape of typical Lorentz/Voigt peaks in the imaginary domain allows for more robust gradient-descent style optimization -initial guesses at peak position can be further away from their real positions without resulting in a spurious fit. It should be noted that the use of complex-valued data in the frequency domain has been suggested in the past (e.g. [START_REF] Dunkel | Identification of organic molecules from a structure database using proton and carbon NMR analysis results[END_REF]), but to the best of our knowledge, the approach has not been described or explored indepth. On top of using both real and imaginary data, the proposed method also adds a convolution term that can be used to incorporate changes in peaks shape that result from truncation or different forms of apodization. The proposed approach is implemented in an open source R package called rnmrfit that is available for download on GitHub (https://github.com/ssokolen/rnmrfit).

Here, we present the theory underlying this approach, as well as an evaluation of its performance on both simulated and experimental data.

Theory

Complex-valued Lorentz function

Using the general notation of Keeler [START_REF] Keeler | Understanding NMR spectroscopy[END_REF], a simple description of a magnetic resonance singlet in the time domain is that of a decaying complex exponential expressed as:

S(t) = S 0 exp(iΩt) exp(-Rt) (1)
where the signal (S) is a function of time (t), relative intensity (S o ), resonance frequency (Ω), and relaxation rate constant (R, which can also be expressed as the reciprocal of the relaxation time 1/T 2 ). The Fourier transform of Equation 1is the complex-valued Lorentz function:

F T [S(t)] = S 0 (A(ω) + iD(ω)) (2)
= S 0 R 2 + iR(ω -Ω) R 2 + (ω -Ω) 2 (3)
where ω is frequency and A and D are the absorption and dispersion lineshapes.

Written explicitly:

A(ω) = R 2 R 2 + (ω -Ω) 2 (4)
D(ω) = R(ω -Ω) R 2 + (ω -Ω) 2 (5)
S 0 is generally taken to have arbitrary units while ω, Ω, and R are all in the same units of frequency, so a variable transformation can be used to make the peak function dimensionless:

f (z) = S 0 1 + iz 1 + z 2 (6)
z = ω -Ω R ( 7
)
As z is dimensionless, the units of Ω and R will be ignored for the rest of this section (with Hz and ppm used interchangeably).

Parameter estimation

Observed Fourier transformed data can be considered as a collection of (x, y) pairs, where y has both a real (y) and imaginary (y) component. Minimizing the sum of squared deviation of both real and imaginary data offers a convenient optimization condition. Considering the fit of a single Lorentz peak to j points, the deviation of the fit from the observed data is expressed as:

z = x -Ω R (8)
f (z) = S 0 1 + iz 1 + z 2 (9)
ε = y -f (z) (10) SS = (ε) 2 + (ε) 2 (11)
where x and y are vectors of data points and SS is the total sum of squares across all points. The calculation of analytical derivatives allows for the convenient use of gradient-descent optimization algorithms:

∂ε ∂S 0 = - 1 + iz 1 + z 2 (12)
∂ε ∂z = -S 0 -2z + i(1 -z 2 ) (1 + z 2 ) 2 (13)
∂ε ∂Ω = ∂ε ∂z ∂z ∂Ω (14)
∂ε ∂R = ∂ε ∂z ∂z ∂R (15)
Thus, the derivative of SS with respect to any parameter θ can be calculated as:

∂SS ∂θ = 2 (ε) ∂ε ∂θ + 2 (ε) ∂ε ∂θ (16)

Gauss and Voigt functions

The same principles can be expanded to consider both Gauss and Voigt functions. The Voigt function can be seen as an extension of the Lorentz, where the signal decay rate is not constant, but takes on normally distributed values -essentially leading to convolution of the Lorentz and Gauss peaks. Although it is possible to approximate the Voigt function using a sum of Lorentz and Gaussian terms, performing the optimization with complex values allows for a more concise representation. The optimization relies on the same general equations as above [START_REF] Bretthorst | Bayesian analysis. III. Applications to NMR signal detection, model selection, and parameter estimation[END_REF][START_REF] Rubtsov | Time-domain Bayesian detection and estimation of noisy damped sinusoidal signals applied to NMR spectroscopy[END_REF][START_REF] Zheng | Identification and quantification of metabolites in (1)H NMR spectra by Bayesian model selection[END_REF], with the specific parameters outlined in Table 1.

It should be noted that the direct evaluation of exp(-z 2 )erfc(-iz) is numerically unstable, so it is necessary to use a numerically optimized implementation of the Voigt function for the calculation (often referred to as the Faddeeva function) [START_REF] Auguie | RcppFaddeeva: 'Rcpp' Bindings for the 'Faddeeva' Package, R package version 0[END_REF].

Lorentz

Gauss Voigt

f (z) S 0 1 + iz 1 + z 2 S 0 exp(-z 2 )erfc(-iz) S 0 exp(-z 2 )erfc(-iz) z x -Ω R x -Ω √ 2R x -Ω + iR i √ 2R G ∂f (z) ∂z -2z + i(1 -z 2 ) (1 + z 2 ) 2 -2zf (z) + i2 √ π -2zf (z) + i2 √ π 2.

Baseline and phase correction

Apart from considering Lorentz, Gauss, and Voigt peak functions, f (z), it is also possible to add baseline terms f b (z). In rnmrfit, f b (z) is modelled as a complex-valued spline polynomial with a variable order and a variable number of equidistant knots. The result is that the baseline is expressed as a piecewise polynomial that is capable of modelling most observed baseline distortions.

Baseline smoothness can be indirectly controlled by both polynomial order and knot number -with 2nd or 3rd order polynomials and 1 or 2 interior knots sufficient for correcting minor distortions. In principle, the baseline for the real and imaginary domains should be the same. In practice, however, differences can arise from a number of possible sources. For example, baseline distortions caused by macromolecules will naturally have different real and imaginary components.

Furthermore, peaks are considerably wider in the imaginary domain than in the real and therefore feature a greater degree of overlap, so fitting a subsection of the spectrum may feature contributions from outside of the selected region that only appear in the imaginary data. To deal with these discrepancies, the default option in rnmrfit is to fit separate real and imaginary baselines, but extra precision can be gained by forcing a single common baseline where appropriate.

Despite the use of two different baselines, peaks are estimated using a single set of parameters, so more information is nonetheless captured than if the imaginary domain is entirely ignored.

The use of complex-valued data also allows the inclusion of phase correction directly in the fit. A phase transform in the frequency dimension can be expressed as y , with the real and imaginary components calculated as:

(y ) = (y) cos(φ) + (y) sin(φ) (17)
(y ) = (y) cos(φ) -(y) sin(φ) (18)
with the corresponding change in ε to ε :

(ε ) = (y ) -(f (z)) ( 19
) (ε ) = (y ) -(f (z)) (20)
Apart from the difference in nomenclature, this change does not impact any of the previously developed equations. But it does add another set of derivatives to fit the local phase angle φ:

∂ (ε ) ∂φ = -(y) sin(φ) + (y) cos(φ) (21)
∂ (ε ) ∂φ = -(y) sin(φ) -(y) cos(φ) (22)

Apodization

Apodization and truncation can be expressed as a simple product in the time domain or a convolution in the frequency domain. Although the time domain calculation is naturally more efficient, the difference in computation time is less pronounced when considering that the frequency domain fit typically considers only a fraction of the data points. The addition of a convolution term can be directly included in the model fit, transforming Equations 10 and 13 into:

ε = y -f (z) * g ( 23
)
∂ε ∂z = - ∂f (z) ∂z * g ( 24
)
where g is a convolution vector with a Fourier transform. Examples of g include the Fourier transforms of a step function to represent truncation or common apodization functions such as exponentials and sinusoidals. The result is that the change in shape of the peak resulting from truncation (the formation of "sinc wiggles") or apodization (line sharpening/broadening) can be incorporated directly in the fit. The apodization can be combined with any of the peak functions as well as baseline and phase correction. Furthermore, g can take on

any function that has a Fourier transform -allowing the application of reference deconvolution or other forms of signal correction in addition to apodization.

Multiple peaks

Although the above derivation considers the fit of only one peak at a time, multiple peaks can be fit by including a sum of all peak contributions, with minimal changes to the derived equations. The definition of multiplets requires the use of constrained optimization algorithms, where the distance between multiple peaks and intensity ratios can be specified. The rnmrfit package allows for both hard (linear) and soft (inequality) constraints to allow for approximate multiplet definitions when the coupling constant may not be known exactly.

Materials and methods

Software implementation

The proposed algorithm has been implemented in pure R in the rnmrfit package by leveraging the nloptr [START_REF] Ypma | nloptr: R interface to NLopt, R package version 1[END_REF] and RcppFaddeeva [START_REF] Auguie | RcppFaddeeva: 'Rcpp' Bindings for the 'Faddeeva' Package, R package version 0[END_REF] packages (which provide R wrappers around the C/C++ NLopt and Faddeeva code developed by the Ab Initio group). Of the various optimization algorithms available through NLopt, the sequential quadratic programming (SQP) algorithm for nonlinearly constrained gradient-based optimization (SLSQP) was chosen for its ability to support both inequality and equality constraints. The rnmrfit package provides the user with a relatively simple API for defining, fitting, and plotting NMR multiplets within 3-5 lines of code. An in-depth tutorial for the package is available on GitHub (https://github.com/ssokolen/rnmrfit).

Simulated NMR data

Sample NMR data was generated in the frequency domain using the Lorentz and Voigt lineshapes (with the Gauss component of the Voigt lineshape set to 50% of the Lorentz). Three peak combinations were chosen to compare algorithm performance with various levels of overlap: a singlet, a singlet and a doublet, and two doublets. Each dataset was generated with 256 real and 256 imaginary data points, with the width of each peak at half height corresponding to approximately 10-20 points. The ideal peak shapes were distorted using different levels of noise as well as phase and baseline errors. Gaussian noise was applied to generate signal to noise (SNR) ratios of 5, 10, 50, and 100, where SNR is calculated as the variance of the pure peak data divided by the variance of the noise. Baseline errors were generated as quadratic polynomials with the position of the maximum/minimum randomly selected to fall within the domain of the data and the magnitude set to 20% of the maximum peak intensity. Phase errors were added at three different levels: 0 • , 15 • , and 30 • , with the sign of the phase error chosen at random.

Cholesterol NMR data

For the basic illustration of the fitting approach, a 13 C-NMR spectrum of cholesterol was recorded using the adiabatic INEPT pulse sequence reported in [START_REF] Bussy | Isotopic finger-printing of active pharmaceutical ingredients by 13C NMR and polarization transfer techniques as a tool to fight against counterfeiting[END_REF]. 90 mg of cholesterol was dissolved in 600 µL of CDCl 3 . The spectrum was recorded on a 500 MHz Bruker Avance-III spectrometer equipped with a 5 mm dual 13 C/ 1 H cryoprobe (tuned to the 13 C recording frequency of 125.76

MHz). The temperature of the probe was set to 288 K. The 90 • 1 H and 13 C pulse widths were calibrated to 10 µs and 11 µs, respectively, with an acquisition time of 0.8 s, recovery time of 10.8 s, 4 dummy scans, and 16 scans. Both inversion and 1 H decoupling were performed with an adiabatic full-passage pulse [START_REF] Tenailleau | Adiabatic 1H decoupling scheme for very accurate intensity measurements in 13C NMR[END_REF].

Refocusing was achieved with a composite adiabatic pulse.

Vanillin NMR data

1 H-NMR spectra of vanillin were used to evaluate the accuracy and precision of the algorithm on real data (previously reported in [START_REF] Jézéquel | Optimized sliceselective H NMR experiments combined with highly accurate quantitative C NMR using an internal reference method[END_REF]). Briefly, 77 mg of 3.5. Synthetic mixture and plasma NMR data

Results and discussion

Basic fit

A basic demonstration of the fit is shown in Figure 1 using 13 C-NMR analysis of a cholesterol sample. Although a similar high quality fit of the basic lineshape can be expected using existing approaches, Figure 1 demonstrates the versatility of the proposed algorithm in correcting phase and baseline errors as well as the incorporation of apodization directly in the fit. The incorporation of resolutionenhancing apodization is likely to be useful for fitting highly overlapping spectra or dealing with highly truncated spectra without the need for significant line broadening. The overall fit including phase and baseline correction is achieved in a fraction of a second on a mid-range laptop with an i7-2640M processor (with the apodization adding a couple of extra seconds). Although the overall calculation time scales with increasing numbers of peaks, there is potential for speed improvement by translating the main optimization code from R into C++ using the Rcpp package [START_REF] Eddelbuettel | Rcpp: Seamless R and C++ integration[END_REF] -ensuring rapid calculation even at larger scales.

Varying noise and error

The impact of signal to noise ratio (SNR) as well as baseline and phase errors was explored using simulated data. Three different resonance combinations were considered as simple test cases -a singlet, a singlet and doublet, as well as Although Figure 2 presents the data for only the Lorentz peak shape, the results were practically equivalent for the more general Voigt shape (see Figure S1 in the Supplementary Information). It was found that the extra flexibility offered by the Voigt shape (in modulating the relative width of a peak at different peak heights) led to an increased possibility of optimization divergence. Effectively, a dramatic change in the baseline and phase terms led to a sub-optimal (global) convergence before the peak width could be correctly fit. However, this effect was eliminated by fitting Lorentz peaks first and using the results as an initial guess to fit the Voigt peaks. Overall, the best results were obtained by fitting Lorentz peaks without phase or baseline correction, refitting with the addition of a phase term, and then changing the Lorentz peak type to Voigt for a final fit. It should be noted that the rates of divergence were rarely higher than 1%, meaning that a direct Voigt fit should be sufficient for most general cases.

Varying initial parameter values

Despite the potential for a high degree of precision and accuracy, iterative optimization algorithms can be very sensitive to the initial parameter values used to begin the optimization. However, the use of both real and imaginary data has made the proposed algorithm more robust to poor initial values. As it is currently implemented, the algorithm requires the user to supply an estimate of only the multiplet chemical shift (as well as the coupling pattern) -the peak heights are estimated as the intensity value corresponding to the chemical shift of the singlets that make up the multiplet and the initial peak width is taken as 1 Hz (which was found to be sufficient for all real data tested, although there is an option to override the default value). The ability of the algorithm to converge onto a global optimum was tested on a simulated singlet with an SNR of 5, with and without phase and baseline error correction. The results are summarized in Figure 3. Despite the significant noise and phase or baseline errors, 100% convergence was achieved for any initial chemical shift that fell within approximately 25% of the maximum peak height. With no baseline or phase errors, the initial chemical shift could be as far away as 5% or 10% of maximum peak height while ensuring convergence. Such estimates could be made by eye or with the aid of separate peak-picking software. The initial peak width was found to have minimal impact (data not shown) and was set at 3× narrower than the true width for the simulation. Some complex combinations of multiplets may suffer from worse performance, but it is unlikely that these can be identified a priori. Performing an initial fit without phase or baseline correction to improve on the initial chemical shift values was sufficient to ensure convergence in all practical cases tested (data not shown). When present, the phase error was set to ±30 • and the baseline error was set to a quadratic function with a magnitude of up to 20% of the maximum peak height.

Application to high precision data

The performance of the proposed approach was assessed in the context of high precision quantitative 1 H-NMR using recently published data [START_REF] Jézéquel | Optimized sliceselective H NMR experiments combined with highly accurate quantitative C NMR using an internal reference method[END_REF]. 1 H-NMR spectra were obtained with the DWET, MWET and PWET pulse sequences capable of yielding a high degree of precision on peak areas from successively recorded spectra, as required in the context of isotopic analysis for which they have been developed. Briefly, these pulse sequences combine several spatial encoding elements to remove the effect of radiation damping in concentrated samples with the goal of determining the relative amount of internal standard at a higher level of accuracy that would be possible with gravimetric analysis. Seven 1 H spectra of vanillin for each pulse sequence were fit using the proposed algorithm and the resulting precision was compared to that achieved by PERCH [START_REF] Laatikainen | A computational strategy for the deconvolution of NMR spectra with multiplet structures and constraints: analysis of overlapping 13C-2H multiplets of 13C enriched metabolites from cell suspensions incubated in deuterated media[END_REF] and Mnova software. PERCH is commonly used for NMR peak fitting and is recognized by the qNMR community for its good performanceparticularly in the field of isotopic NMR where it acts as a reference, whereas Mnova offers the added benefit of automated peak-picking on top of spectral deconvolution. The results of the comparison are summarized in Figure 4 with the complete area quantification data available as a Microsoft Excel file in the Supplementary Information (the Mnova CV values are excluded from the figure as they were considerably higher than both of the other methods). Overall, the proposed approach was found to be more precise across practically all spectral regions and all three of the pulse sequences tested, in both absolute terms (Fig- ure 4C) and when normalized by the area of the DMSO 2 peak (Figure 4D).

The magnitude of the difference, however, varied considerably across the different regions. The absolute quantification of DMSO 2 and methoxy peak areas using the proposed approach was observed to have a coefficient of variation in the range of 0.02%-0.05%, approximately five times more precise than with PERCH. As shown in Figure 4B, the signal to noise ratio for these two peaks was visibly larger than for the other spectral regions and suggests that the proposed approach can take advantage of strong signal data -at smaller signal to noise ratios, the fitting algorithm is unlikely to have as much impact. However, it should be noted that the primary strength of the proposed algorithm is its robustness and generality, with the increase in precision serving as a useful side benefit.

Application to more complex data

Whereas the analysis above has focused on relatively simple spectra, the proposed approach can also be used for more complex data, featuring significant spectral crowding or baseline distortions. The following examples are not in-tended to be comprehensive but offer a glimpse of what the approach is capable of. Moving beyond basic singlets, a mixture of 17 peaks with some light overlap (stemming from glutamate, methionine, creatine, histidine, and phenylalanine in the 3.70-4.05 ppm range) could be fit within a couple of seconds and achieve a coefficient of variation of approximately 0.5% for the peak area corresponding to each compound (see Figure 5 for a visualization of the fit). The peaks were fit as four double doublets and one singlet, with the J-coupling frequency obtained from the Human Metabolome Database (HMDB, [START_REF] Wishart | HMDB 4.0: the human metabolome database for 2018[END_REF]) -rnmrfit can also make use of inequality constraints to allow for cases where the J-coupling frequency estimates may not be exact. An example from a human plasma sample featuring a complex baseline is presented in Figure 6. Modelling the complex baseline required the number of interior baseline knots to be increased from 2 or 3 (used in the previous examples) up to 25, but the resulting fit appears to be quite good -the baseline estimate follows the general contour of the data while leaving a sufficient gap at the base of the peaks. Despite the increased sample complexity, the median coefficient of variation for the peak area was found to be approximately 2%, with the fit taking several seconds per spectrum.

Although rnmrfit is primarily intended to be used for fitting spectra where the identity of many compounds are known ahead of time (and their corresponding NMR spectra can be defined as a series of multiplets), the fit algorithm itself can also be used in the context of global spectrum decomposition, where a complex spectrum is broken down into a series of singlets (or other simple spectral features). An example from a human plasma sample featuring complex spectral crowding is presented in Figure 7, where the proposed approach demonstrates a very good fit of the data using 86 singlets and a baseline with 20 interior knots. Despite the good fit, it is also important to note a number of key limitations. First, the initial position of the singlets was identified by eye. Automated peak selection is a challenging task that was deemed to be outside the scope of the current implementation, which focuses squarely on the fit process itself 1 .

1 Future versions of the software may include an implementation of an existing peak-picking Second, achieving a reasonable fit of complex data required liberal use of optimization constraints. Once picked, peak positions were not allowed to vary more than 1% of the subspectrum ppm range (or 0.004 ppm for the 0.4 ppm range of data being fit) and peak widths were constrained to be in the range of 0.5-2 Hz to better separate peak and baseline contributions. Furthermore, singlets were constrained from "passing" or completely overlapping each other during optimization -thereby limiting the possibility of the optimization converging on a non-useful local minimum. Phase correction was also constrained to within approximately 10 • as it was expected that dedicated phase correction software method (e.g. [START_REF] Abbas | Automatic peak selection by a Benjamini-Hochberg-based algorithm[END_REF][START_REF] Cobas | Automatic assignment of 1h-NMR spectra of small molecules[END_REF][START_REF] Tikole | Peak picking NMR spectral data using non-negative matrix factorization[END_REF]) or present a novel approach should be capable of getting relatively close to the true value. Although rnmrfit assumes rational defaults for all of these constraints, the default values may not be suitable in all cases. Third, increasing the number of peaks in a single fit considerably slowed down the overall optimization time. The time it takes to achieve a reasonable overall fit increases approximately linearly with each peak and baseline knot -so 50-100 peaks can be fit within a couple of minutes.

However, a small fraction of peaks require considerably more time than the rest of the peaks to achieve a good fit, slowing down computation to more than an hour to ensure high levels of precision for all peaks. One solution to this problem is computational -implementing the minimization function in Fortran sped up computation by a factor of 10-100 but the resulting code had to be dropped as it was not easily portable across all platforms (a C++ version is planned for the near future). Another solution is to segment the overall spectrum into smaller regions prior to fitting, thereby preventing one or more problematic peaks from slowing down the fit.

Conclusion

The imaginary component of NMR data is not typically considered when fitting peaks in the frequency domain. However, adding this component is relatively straightforward and enables the incorporation of both phase and baseline

DMSO 2

 2 and 250 mg of vanillin were dissolved in 510 µL of acetone-d 6 (with 11.76 mM Cr(acac) 3 relaxing agent). 1 H-NMR spectra were recorded using the recently published pulse sequences for high-precision 1 H quantitative NMR (DWET, MWET and PWET)[START_REF] Jézéquel | Optimized sliceselective H NMR experiments combined with highly accurate quantitative C NMR using an internal reference method[END_REF] on a 400 MHz Bruker Avance-I spectrometer equipped with a 5 mm dual13 C/ 1 H probe. The temperature of the probe was set to 303 K. The 90 • 1 H pulse width was calibrated to 10 µs, with an acquisition time of 4s, delay time of 5s, 2 dummy scans, and 4 scans.Seven spectra were collected for each pulse sequence using a single sample to estimate precision and accuracy. Peak areas were calculated using PERCH Software (Perch solutions Ltd, Kuopio, Finland), global spectral deconvolution (GSD) implemented in Mnova 12.0 (Mestrelab Research, S.L., Santiago de Compostela, Spain), as well as the proposed algorithm. The overall spectra were divided into 5 regions, corresponding to the methoxy, aromatic, hydroxy, and aldehyde groups of vanillin along with the DMSO 2 peak. All but the aromatic regions were fit with one singlet each, while the aromatic region required four singlets.

Figure 1 :

 1 Figure 1: Demonstration of proposed algorithm using the C2, C7, and C8 resonances of cholesterol (see Figure S2 in the Supplementary Information for reference). The data fit in the first row is the original data with 2 Hz line broadening, the data fit in the second row had simulated baseline and phase errors added to it, and the data fit in the third row has had sine bell apodization applied.

30 DFigure 2 :D

 302 Figure 2: Summary of Lorentz peak fit accuracy and precision as a function of SNR as well as phase and baseline error. 100 peaks were randomly generated for each condition -combining an SNR of 5, 10, 50, or 100 with a phase error of 0 • , 15 • , or 30 • . A quadratic baseline with a magnitude of up to 20% of the maximum peak height was added at each condition. A Example spectrum at an SNR of 100 with 0 • phase error. B Example spectrum at an SNR of 5 with 30 • phase error. C Median error of calculated area as a percent of the true area. D Coefficient of variation of the calculated area (standard deviation divided by the average area).

Figure 3 :

 3 Figure 3: Impact of initial chemical shift on global convergence with and without baseline and phase error correction. A Example singlet with an SNR of 5 and no phase or baseline error. B Percent convergence based on 1000 iterations as a function of initial chemical shift.

Figure 4 :

 4 Figure 4: Comparison of peak area quantification precision between rnmrfit and PERCH Software based on the analysis of 7 1 H-NMR spectra of vanillin (see Materials and methods for more detailed sample description). A Hydrogens corresponding to each spectral region highlighted with parentheses. B Example spectra corresponding to each spectral region fit with the software. C The coefficient of variation (in percent of mean value) from the quantification of area for each spectral region. D The coefficient of variation (in percent of mean value) from the quantification of area for each spectral region, normalized by the area of the DMSO 2 peak.

Figure 6 :

 6 Figure 6: Demonstration of proposed algorithm on a subspectrum of a human plasma sample (see Materials and methods for more details). A Overlap of the overall fit (black lines) and estimated baseline (dashed line) over the observed data (gray fill). B Breakdown of the overall fit into constituent peaks.

Figure 7 :

 7 Figure 7: Demonstration of proposed algorithm on a subspectrum of a human plasma sample (see Materials and methods for more details). A Overlap of the overall fit (black lines) and estimated baseline (dashed line) over the observed data (gray fill). B Breakdown of the overall fit into constituent peaks.

Table 1 :

 1

Summary of Lorentz, Gauss, and Voigt functions, where S 0 is relative intensity of the peak (its relative height), Ω is the resonance frequency (relative peak position), and R is the relaxation rate constant (relative peak width). The R G term in the Voigt function relates to the distribution of possible R values (it can be seen as the Gaussian component of the overall peak).

correction as well as apodization directly in the peak fitting process -resulting in the robust fit of small or overlapping peaks with the added computational efficiency of working in the frequency domain. Application to real data suggests that this approach can also result in dramatically higher precision than can be achieved with commonly used PERCH software (although the increase in precision is not uniform across different spectral regions or different spectra).

The proposed algorithm has been implemented in the rnmrfit package for the R programming language (which has been made available on GitHub), but it is not specifically tied to this particular software. The implementation of this approach in other software is likely to be of benefit for a broad community of qNMR users.

Acknowledgment

The authors would like to thank the Natural Sciences and Engineering Research Council of Canada (NSERC) for the NSERC Postdoctoral Fellowship awarded to SS as well as the Lebanese CNRS and the research council of the Saint Joseph University of Beirut for their financial support of GH.

The authors would also like to thank two anonymous reviewers whose input helped improve the manuscript.

Conflict of Interest

The authors declare that they have no conflicts of interest.


Diego Matos
email: diego.matos@uerj.br

Silva Lopes

Americo Cunha
email: americo.cunha@uerj.br

On the physical consistency of evolution laws obtained with sparse regression

Keywords: nonlinear dynamics, data-driven dynamical system, sparse regression, SINDy, physics-informed machine learning

This work employs the technique known as sparse identification of nonlinear dynamics (SINDy) to infer, from a set of information provided by a given time-series, the evolution law of a dynamical system of interest, accessing the physical consistency of the obtained dynamic model. The Duffing oscillator is used as a benchmark due to the variety and richness of its dynamical behavior. The numerical experiments attempt to identify whether the method is capable of recognizing the correct evolution law and respecting basic principles of physics such as the balance of momentum and energy. The numerical results illustrate the method's ability to obtain approximations to the system evolution law that provides physically consistent behavior for short time intervals.

Introduction

Most classical dynamic systems, such as a harmonic oscillator, a simple pendulum have their evolution law inferred through the first principles [START_REF] Corben | Classical mechanics: 2nd edition[END_REF][START_REF] Meirovitch | Methods of analytical dynamics[END_REF]. However, some modern dynamical systems such as those that appear in areas like epidemiological modeling [START_REF] Dantas | Calibration of a SEIR-SEI epidemic model to describe the Zika virus outbreak in Brazil[END_REF][START_REF] Dantas | An uncertainty quantification framework for a Zika virus epidemic model[END_REF], neuroscience [START_REF] Marblestone | Toward an integration of deep learning and neuroscience[END_REF][START_REF] Glaser | The roles of supervised machine learning in systems neuroscience[END_REF][START_REF] Richards | A deep learning framework for neuroscience[END_REF], have the basic evolution laws not well understood or even unknown so that the use of first principles to obtain the governing equations is unfeasible. In other applications, structural health monitoring [START_REF] Villani | Damage detection in nonlinear systems based on stochastic Volterra series[END_REF][START_REF] Villani | Damage detection in an uncertain nonlinear beam based on stochastic Volterra series: An experimental application[END_REF], for instance, it may be possible that the basic principles of dynamics are well understood, but some peculiarities not, which make it difficult to construct high-fidelity predictive models.

In this context, discover the evolution law for this kind of dynamical system in an analytical way becomes an almost impossible task. However, with the emergence of the information age, where data is exponentially generated and stored, it is natural to try to use information embedded in datasets associated with the dynamical system of interest to infer its underlying evolution law.

Neural networks (NN) [START_REF] Hasssanibesheli | Reconstructing complex system dynamics from time series: a method comparison[END_REF][START_REF] Lecun | Deep learning[END_REF][START_REF] Zhang | Learning in modal space: solving timedependent stochastic PDEs using physics-informed neural networks[END_REF] have a high capacity to learn complex patterns, which makes them at first glance natural candidates to be used in the search for mathematical expressions that represent the desired evolution law. However, despite having incredible power of prediction (interpolation) in the domain spanned by the training data, several NN architectures lack in providing interpretable expressions for the identified dynamic model, which makes it difficult to generalize (extrapolate) system behavior for unknown conditions.

Alternatively, regression techniques are more proper for providing the interpretability that NN does not have, which makes them an attractive tool to extract dynamic evolution laws from datasets [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF]. In this scenario, stands out the technique known as sparse identification of nonlinear dynamical systems (SINDy) [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Kaiser | Sparse identification of nonlinear dynamics for model predictive control in the low-data limit[END_REF][START_REF] Bhadriraju | Operable adaptive sparse identification of systems: Application to chemical processes[END_REF][START_REF] Zhang | On the convergence of the SINDy algorithm[END_REF], a sparse regression method to identify evolution equations from data that proven to be efficient in different areas and problems, and stands out for three aspects: (i) interpretability of the obtained equation; (ii) excellent generalization (extrapolation) capability; (iii) computational efficiency. As several dynamic systems have evolution laws with only a few terms, SINDy looks for a sparse and parsimonious differential equation that best fits the known data.

If the method is applied in a system where the first principles are unknown, it is interesting to have some guarantee that the identified evolution law respects fundamental physical laws since their violation can translate into a dynamic behavior that is divergent from the original one. For instance, if the balance of momentum is not respected, mass or velocity are affected, potentially changing the inertial effects. When the energy balance is violated, the contribution from kinetic/potential energy to the mechanical energy can be in a proportion far from the original system, which can imply changes in the displacement/velocity. Even when a qualitatively good result is observed, small differences in the system parameters can significantly change quantitatively the long-term behavior.

This work aims to verify the physical consistency of the evolution law obtained by SINDy, checking if the identified dynamical system respects physical principles such as the balance of momentum and energy, and if the underlying balance of these quantities is an accurate approximation for those associated with the original system. For this purpose, the Duffing oscillator is used as a benchmark, where the underlying evolution equation is obtained from synthetic data that emulate experimental measurements for the displacement time-series of this dynamical system.

The rest of this manuscript is organized as follows. In section 2, it is presented a brief explanation of the Duffing oscillator and the methodology used to identify the dynamic evolution laws. Some examples and test results are presented in section 3. Finally, in section 4, the final remarks are highlighted.

Methodology

The dynamical system chosen to test the physical consistency of the identified dynamic laws obtained with SINDy is the Duffing oscillator, once this kind of oscillator has many possible applications, such as structural dynamics and energy harvesting [START_REF] Cottone | Nonlinear energy harvesting[END_REF][START_REF] Lopes | Numerical study of parameters influence over the dynamics of a piezo-magneto-elastic energy haversting device[END_REF][START_REF] Lopes | Nonlinear characterization of a bistable energy harvester dynamical system[END_REF][START_REF] Peterson | Maximization of the electricalpower generated by a piezo-magneto-elastic energy harvesting device[END_REF][START_REF] Guyomar | Energy Harvesting using Non-linear Techniques[END_REF][22], and a very rich and well known dynamic behavior. Fig. 1 shows a schematic of a vibratory system that behaves like a Duffing oscillator, in which the dynamic behavior evolves according to

ẍ + δ ẋ + α x + β x 3 = γ cos(ωt) , (1)
where the function x = x(t) is the displacement of the beam tip, and the ẋ and ẍ are respectively the first and second derivative of x, i.e., the velocity and acceleration [START_REF] Brennan | Examples of physical systems described by the Duffing equation[END_REF]. The parameters in the equation of motion are the damping coefficient δ; the linear stiffness α; the nonlinear stiffness β; the external excitation amplitude γ; and the external excitation frequency ω.

Sparse identification of nonlinear dynamics

The SINDy method can handle first order dynamical systems of form

ẋ = f (x) , (2)
where x is the state vector, ẋ is the time-derivative of the state vector, and f (x) is the unknown evolution law. Higher-order systems may be written as a first-order system by increasing the state dimension. For instance, the Duffing oscillator in Eq.( 1) can be rewritten as

ẋ1 = φ x 2 , ẋ2 = -δ x 2 -α x 1 -β x 1 3 + γ cos x 3 , (3) ẋ3
= ω.
To determine the dynamic evolution law f it is necessary to have observations (time-series) of the system state x and velocity ẋ, which are organized as follows

X =      x T (t 1 ) x T (t 2 ) . . . x T (t m )      =      x 1 (t 1 ) x 2 (t 1 ) . . . x n (t 1 ) x 1 (t 2 ) x 2 (t 2 ) . . . x n (t 2 ) . . . . . . . . . . . . x 1 (t m ) x 2 (t m ) . . . x n (t m )      , (4)
Ẋ =      ẋT (t 1 ) ẋT (t 2 ) . . . ẋT (t m )      =      ẋ1 (t 1 ) ẋ2 (t 1 ) . . . ẋn (t 1 ) ẋ1 (t 2 ) ẋ2 (t 2 ) . . . ẋn (t 2 ) . . . . . . . . . . . . ẋ1 (t m ) ẋ2 (t m ) . . . ẋn (t m )      . ( 5
)
After having the data arranged in a matrix form, it is possible to build a library (dictionary) of candidate functions Θ(X) to construct the evolution law terms. In general, this library includes polynomials, trigonometric, exponential, logarithmic functions (among other possibilities), i.e.,

Θ(X) =   1 X X P2 X P3 X P4 . . . X P k sin(X) cos(X) e X    , (6)
where each column in this matrix represents a candidate function, the X Pn notation indicates all possible n-order polynomials formed by combining the state variables. What SINDy does is to solve the following regression problem

Ẋ ≈ Θ(X)Ξ, (7)
where the coefficient matrix Ξ is determined by minimizing a misfit function that measures the discrepancy between Θ(X)Ξ and Ẋ, i.e.,

Ξ * = arg min Ξ || Ẋ -Θ(X) Ξ|| 2 + λ ||Ξ|| 1 , (8)
where || • || 2 and || • || 1 denote the norm-2 and norm-1, respectively. Note that, to produce a sparse solution, a user-controlled λ parameter is introduced, which acts as a threshold. A first least-squares regression results in a vector with estimated values for Ξ. Coefficients that represent the exact evolution law have more meaningful value than those that are not. After that, the absolute values of the coefficient vectors Ξ, which is smaller than a λ, have their values changed to zero. This process, known as sequential thresholded leastsquares, is repeated until a parsimonious evolution law is discovered [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF]. For systems where the best parameter λ is unknown, machine learning techniques for model validation are necessary to ensure a good result, e.g., cross-validation, among others [START_REF] Bishop | Pattern recognition and machine learning[END_REF].

Training dataset

In a typical scenario, the dataset used to infer the evolution law of the dynamic system are observations obtained with some type of sensor, in the field or in the laboratory. To emulate a set of displacement and velocity measurements for the Duffing oscillator, synthetic data are considered in this work. The corresponding time-series are generated with the aid of the numerical integrator ODE45 from MATLAB, fixing the time step as ∆t = 0.01. Concretely, a dataset with 151 equally spaced samples between 0 and 15 seconds is obtained with the numerical integrator and polluted with a zero-mean Gaussian white noise (with variance η 2 ) to emulate the fluctuations in experimental measurements.

Balance of energy and momentum

The total energy, for both the original and the identified system, is the sum of the kinetic and potential energy

E = 1 2 x 2 2 + 1 2 α x 2 1 + 1 4 β x 4 1 . (9)
Although this scalar quantity is not conserved, as the considered dynamics is dissipative, the balance between the mechanical energy of the system and the work done by external and dissipation forces must be respected, i.e.,

1 2 x 2 2 + 1 2 α x 2 1 + 1 4 β x 4 1 t 0 = t 0 -δ x 2 2 (τ ) + γ cos (ωτ ) x 2 (τ ) dτ. (10)
Similarly, the linear momentum variation has to balance the impulse of the forces acting on the system

x 2 (t) -x 2 (0) = t 0 -δ x 2 (τ ) -α x 1 (τ ) -β x 3 1 (τ ) + γ cos (ωτ ) dτ. (11)
Therefore, if SINDy identifies the correct candidate functions, the balance of energy and momentum is respected by construction. However, even in this case, a question persists: Compared to the original dynamics, how representative (in quantitative terms) are these balance principles?

Results and discussion

Two parameter settings are used to evaluate the balance of energy and momentum in the SINDy method, nominally, the unforced conservative case, and the forced dissipative dynamics. For all the numerical experiments done and the database used, SINDy obtained enough parsimonious results to activate just the right candidate functions.

In the first case, the unforced conservative dynamics, 151 equally spaced time-samples are used in the training step, with a library of polynomial functions up to third-order without trigonometrical terms, with λ = 0.2, and initial conditions x 0 = 2 and v 0 = -2. Table 1 shows the values that SINDy inferred for two different noise intensities. The more noise, the less similar is the identified dynamics to their original counterpart.

In Fig. 2 the reader can see a verification test that estimates if the identified dynamics present an energy/momentum balance that is compatible with the invariances intrinsic to the original dynamics. The right side of Eqs. [START_REF] Hasssanibesheli | Reconstructing complex system dynamics from time series: a method comparison[END_REF] and [START_REF] Lecun | Deep learning[END_REF] are computed for both the identified dynamics and the original dynamics, and the differences between the respective quantities are computed, to estimate the balance error due to deviations in the identified coefficients. It is possible to notice that the error in the momentum balance is directly proportional to the noise intensity, emphasizing the importance of obtaining good quality data.

Table 1. Parameters of the original and identified dynamics in the conservative unforced case, with noise intensities equal to η = 0.001 and η = 0.010. Fig. 2. The difference in the balance of energy and momentum between the original and identified dynamics for the unforced conservative case, with η = 0.001 (left) and η = 0.010 (right). The two insets below each graph show a zoom from the upper rectangles, where it is possible to see the symmetrical format of the curve. The velocity is plotted together with the momentum error, to show that the error in the momentum is a local maximum whenever the oscillator velocity is zero.

Therefore, despite identifying candidate functions correctly and with coefficients close to the original, the difference in both balances is a sinusoidal growing function over time. For the two noise intensities, the coefficients identified for the linear stiffness are lower than the original, while the nonlinear stiffness is slightly higher. Both systems after a long time become out of phase with the original one. For the systems with η = 0.010, this is achieved early than the case where η = 0.001, with time for that bigger than 300.

In Fig. 3 a case without noise is considered, to assess the fundamental limit of noiseless measures. Note that, for different initial conditions, the error in the momentum balance is comparable to floating point noise, so that for any practical effect, SINDy can be considered to respect the linear momentum balance. Now the dissipative forced Duffing oscillator is considered. Two different noise intensities, η = 0.001 and η = 0.010, are employed once again, as well as the library of polynomial functions up to third-order with trigonometrical terms, 151 equally spaced time-samples, and λ = 0.02, initial conditions x 0 = 2 and v 0 = -2. The SINDy is able to identify the correct terms with great accuracy, despite the small amount of data, as can be seen in Table 2, which shows the system parameters for the original and the identified dynamics.

The balance of energy and momentum are plotted for this second case in Fig. 4, where it is possible to notice a richer oscillatory pattern, for both the momentum and energy errors, when compared to the first case. Looking at the energy error, initially, the local pikes and valleys are coinciding with those from the momentum error curve, but after some time, they become out of phase.

Table 2. Parameters of the original and identified dynamics in the dissipative forced case, with noise intensities equal to η = 0.001 and η = 0.010. The difference in the balance of momentum and energy have both the same asymptotic behavior as in the conservative case, but now it is possible to observe an oscillatory pattern that is completely different from the conservative unforced case. A transient is noticed that is quickly dampened giving rise to a periodic stationary response, with an amplitude that increases (very slowly) with time.

Conclusion

This work presented a study on the physical consistency of the evolution law inferred from a given dataset employing the SINDy method. Two configurations of the Duffing oscillator are used to investigate the difference in energy and momentum balance with respect to the original dynamics. In all cases, SINDy identified the right candidate functions, therefore the balance of energy and momentum is respected. However, analyzing the difference in the energy and momentum balance identified in relation to the original, that function is growing rapidly over time. This implies that for long time periods, the identified dynamic will be out of phase to the original.

Fig. 1 .

 1 Fig. 1. Schematic of Duffing oscillator driven by a sinusoidal external excitation.

9 Fig. 3 .

 93 Fig.3. The difference in the balance of energy and momentum between the original and identified dynamics with noise-free data, for different initial conditions. On the top, the initial conditions are x0 = 2 and v0 = -2, on the middle x0 = 3 and v0 = 2, and for last, on the bottom x0 = 1 and v0 = 2.5.

1 Fig. 4 .

 14 Fig.[START_REF] Dantas | An uncertainty quantification framework for a Zika virus epidemic model[END_REF]. The difference in the balance of energy and momentum between the original and identified dynamics for the dissipative forced case, with η = 0.001 (left) and η = 0.010 (right).

Acknowledgments

The authors would like to thank the financial support given to this research by the Brazilian agencies Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES) -Finance Code 001, and the Carlos Chagas Filho Research Foundation of Rio de Janeiro State (FAPERJ) under the following grants: 211.304/2015, 210.021/2018, 210.167/2019, and 211.037/2019.


Pham Tran

Anh Quang

Yassine Hadjadj-Aoul
email: yassine.hadjadj-aoul@irisa.fr

Abdelkader Outtagarts
email: abdelkader.outtagarts@nokia-bell-labs.com

On Using Deep Reinforcement Learning for VNF Forwarding Graphs Placement

The placement of Virtual Network Function -Forwarding Graphs (VNF-FGs) is one of the fundamental operations in the next generation networks. However, embedding a VNF-FG is a complicated process which has been proved as an NP-hard problem. In this demo, we leverage one of the most advanced approaches in Deep Reinforcement Learning (DRL) in order to efficiently solve the VNF-FG embedding problem. In particular, we have proposed an extension of this approach in order to better explore the space of solutions while making it safer. The experimental platform we have developed, which is based on Mininet and containers, allows us to clearly show the advantage of our approach over existing approaches.

I. INTRODUCTION

With the advent of Network Function Virtualization (NFV), today's network infrastructures are becoming more flexible, allowing them to deploy services almost in real time and ondemand [START_REF] Stampa | A deep-reinforcement learning approach for software-defined networking routing optimization[END_REF], making the concept of network slicing a reality [START_REF] Pham | Deep Reinforcement Learning Based QoS-Aware Routing in Knowledge-Defined Networking[END_REF].

Resolving the placement of complex services, in the form of a Virtual Network Function-Forwarding Graph (VNF-FG), remains a very difficult problem that can hardly be scaled up [START_REF] Dolati | Deepvine: Virtual network embedding with deep reinforcement learning[END_REF]. Indeed, the problem is known to be NP-hard, which limits the exact resolution of this kind of problem to very small service instances. Therefore, several heuristics have been developed in the literature [START_REF] Li | Virtual network function placement and resource optimization in nfv and edge computing enabled networks[END_REF]. However, although very fast, these approaches often get stuck in local minima and adapt very poorly to the addition of constraints or change of objective. The use of metaheuristics does not solve the problem either. Indeed, the latter do not benefit from past experiences, hence the benefit of using learning strategies [START_REF] Quang | A deep reinforcement learning approach for vnf forwarding graph embedding[END_REF].

Reinforcement learning is certainly the most appropriate approach among the existing learning techniques, in view of the conditions prevailing in the networks (i.e. variable bitrate traffic, presence of congestion, unpredictable behavior of certain services, etc.). However, only few approaches exist in this category [START_REF] Quang | A deep reinforcement learning approach for vnf forwarding graph embedding[END_REF].

In this paper we present a demonstrator of a reinforcement learning technique that combines deep learning and a heuristic, which makes it safer [START_REF] Quang | A deep reinforcement learning approach for vnf forwarding graph embedding[END_REF]. Indeed, the proposed solution, named EDDPG, not only improves the exploration mechanism suggested for the Deep Deterministic Policy Gradients (DDPG) technique [START_REF] Lillicrap | Continuous control with deep reinforcement learning[END_REF], it allows, by combining it with a heuristic, to have at worst the performance of the latter.

The rest of this paper is structured as follows. In Section II, we introduce the architecture of the proposed solution. In Section III, we introduce the experimentation carried out and discuss the obtained results. We conclude our paper in Section IV.

II. THE PROPOSED SYSTEM ARCHITECTURE

A. Overview

The emergence of deep reinforcement learning has unveiled effective solutions to complex problems, in particular in networking. In [START_REF] Stampa | A deep-reinforcement learning approach for software-defined networking routing optimization[END_REF], [START_REF] Pham | Deep Reinforcement Learning Based QoS-Aware Routing in Knowledge-Defined Networking[END_REF], DRL has been adopted to solve routing problems. Furthermore, it can also be exploited to solve more complicated problem such as virtual network embedding [START_REF] Dolati | Deepvine: Virtual network embedding with deep reinforcement learning[END_REF]. However, the high dimensional, in both state space and action space of VNF-FG embedding problems, poses a challenge to DRL. Thanks to recent advances in DRL, e.g., DDPG, we can tackle the high dimensional issues. In the previous study [START_REF] Quang | A deep reinforcement learning approach for vnf forwarding graph embedding[END_REF], we confirmed the capabilities of DRL in coping with VNF-FG embedding problems via simulations. This paper is the continuation in which we implement our solution in an emulation environment based on Mininet1 , with a componentbased SDN controller framework (RYU2 ). The DRL agent that we proposed is integrated within an orchestrator designed by the authors.

B. Architecture description

Fig. 1 presents the architecture of the demo. The MininetController loads the substrate network stored in the shared memory. Then, it can initialize the substrate network via MininetRunner . Within the Orchestrator , the LoadSN (Load Substrate Network) function loads the substrate network from the shared memory and feed substrate network information to the algorithm. For each iteration, the GenerateVNFFGs function generates a set of VNF-FGs which will be loaded into the algorithm via the LoadVNFFGs function. The algorithm computes the solution based on the substrate networks and the virtual network request. Then, it writes the solution to the shared memory and it sends a command to MininetController to execute the simulation via the RestClient . The • MininetRunner : Defined in the file mininet runner.py to enable controlling the mininet emulator.

• RCG : Defined in routing config generator.py to convert a solution calculated by the algorithm into routing configurations to be sent to the SDN controller.

III. EXPERIMENTS

All the modules of our demonstrator have been installed on a laptop, although we can consider putting each module in a separate element. We considered three placement strategies: the First-Fit heuristic, the DDPG algorithm and the proposed EDDPG solution. The First-Fit heuristic consists in placing the virtual nodes, each time in the first available physical node, and then placing the paths using the Dijkstra algorithm.

To assess the performance of the proposed approach, we consider, for the first experiments, the network topology of BtEurope [START_REF]The internet topology zoo[END_REF] with 24 nodes and 37 full-duplex links. We considered the placement of 10 VNF-FGs, comprising between 2 and 10 VNFs each. For the physical nodes, we considered normalized CPU, RAM and Storage resources in the interval [0.5, 2.0]. For physical links we have considered latencies in the interval [1.0, 10.0]ms3 , and a loss rates in the interval [1%, 4%]. For the virtual nodes, we considered resources in the interval [0.1, 1.0]. For the virtual links we considered latencies in the interval [10, 100]ms and loss rates in the interval [1%, 10%].

Figures 2, 3 and4 show the results obtained with the different strategies: First-Fit, DDPG and EDDPG, respectively. The reward, in the different figures, represents the percentage of VNF-FGs that have been placed.

Note that a successful placement is when you manage to place all the nodes of a request but also all its links. Note as well that each point in the plots represents a request for placement of a batch of VNF-FGs, for which the difficulty of placement may vary, hence the variability of the results obtained.

We can easily notice that the heuristic has an average success rate of 20.6% while DDPG has a success rate of only 17.1%, which reveals the limits of applying a reinforcement learning approach directly, without adapting it to the placement problem. On the other hand, the strategy we proposed has an average success rate of 28.75%, which represents an improvement of 68.1% over DDPG and 39.6% over the heuristic. 4 also show the placement of virtual nodes and links on the physical network. It is quite easy to see that the number of nodes placed varies from one approach to another. As mentioned above, the number of placed VNFs does not reflect the quality of the placement since a VNF-FG is placed only when all its nodes are placed as well as all its links.

When performing experiments over a larger number of episodes (i.e., 1000 episodes), we observed that the results remained in the same order. In other words, learning converges fairly quickly thanks to the technique we proposed. It is, indeed, easier to learn how to assign weights (i.e. how to make a particular shuffle when assigning virtual nodes), which is then exploited by the heuristic we have proposed, than to learn the general model for the placement. Thus, the methodology we have proposed not only ensures that we have at least the performance of the heuristic right from the start, but also allows to converge towards optimal performance in a very short period of time.

IV. CONCLUSIONS

Service placement combined with deep reinforcement learning techniques are certainly part of the technological building blocks to converge towards a completely autonomous network, known as Zero-touch network.

Based on a learning strategy that we have previously proposed, we provide in this paper a demonstrator the potential of this type of solution in dealing with service placement. This demonstrator also shows that such a solution can be considered for the control of a real network infrastructure. Indeed, the approach that we have proposed provides guarantees in the quality of the placement, since in the worst case, and from the beginning of the learning process, it allows to obtain the performance of the heuristic.

In our future work, we will consider the use of an even more realistic environment by integrating an orchestrator such as Kubernetes while considering the placement of real network services.

Fig. 1 :

 1 Fig. 1: Demo architecture

Fig. 2 :

 2 Fig. 2: First-Fit strategy

  Figures 2, 3 and 4 also show the placement of virtual nodes and links on the physical network. It is quite easy to see that the number of nodes placed varies from one approach

  MininetRunner module. This module loads the substrate network defined in ../runs/nodes.txt and ../runs/edges.txt initialize: Set up the MininetRunner module and initialize RCG. In the setup function of MininetRunner , the network is created with IP addresses and links. Then, the switches and the connections of the substrate networks are entered as arguments to initialize the RoutingConfigGenerator module. load data: Load VNFs and VLs configuration to RoutingConfigGenerator. setup routes: Send routes of the VLs to the SDN controller. send traffic: Running the Distributed Internet Traffic Generator (D-ITG). del route: Delete all the routes at the SDN controller. del data: Delete all the QoS measures. Defined in mininet controller.py. It comprises the definitions of the REST interface to control the MininetController.2) SDN Controller: The main file is reset router.py, which enables to configure the controller via a REST API. It could be deployed by running the following command line: It loads the substrate network to the algorithm.• RestClient : It is in charge of communicating with the MininetController. • RetrieveQoS : This function role is to load the QoS results from the shared memory. This allows computing the rewards.

					• Algorithm (defined in orchestrator.py, main function): It
					defines the available algorithms to compute the VNR
					embedding. ALG variable can be 0, 1, 2 corresponding to
					FirstFit-Dijkstra algorithm (hueristic), DDPG, and EED-
					DPG, respectively. We can define a new algorithm in this
					section.
					• Constants (defined in const.py): It defines the following
					metrics
					-The number of episodes for the training
					(EPISODE COUNT).
					-The number of steps in each episode
					(MAX STEPS).
					-The minimum and maximum request bandwidth
					(MIN REQ BW, MAX REQ BW).
					-Latency,	and	lossrate	of	VLs
					(MIN REQ LATENCY, MAX REQ LATENCY,
					MIN REQ LOSS, and MAX REQ LOSS).
					-The number of VNF-FGs (NUM VNFFGS), and the
					minimum and maximum VNFs in each VNF-FG
	• MininetController	:	Defined	in	(MIN VNFS and MAX VNFS)
	mininet controller.py.
	-init : Intialize ryu-manager rest\_router.py
	3) Orchestrator: All functions of the Orchestrator are
	defined in orchestrator.py

• Action :

• GenerateVNFFGs : It generates the VNF-FGs.

• LoadVNFFGs : It loads VNF-FGs to the algorithm.

• LoadSN :

http://mininet.org/

https://ryu-sdn.org/

The effective latency results from the placement, it is therefore calculated dynamically.


Julio Cesar De Castro Basilio

José Geraldo

Telles Ribeiro

Americo Cunha

Tiago Roux Oliveira

Julio Cesar Basilio
email: basilio.julio@posgraduacao.uerj.brjose.ribeiro@uerj.bramerico.cunha@uerj.brtiagoroux@uerj.br

An optimal fractional LQR-based control approach applied to a cart-pendulum system

Keywords: fractional controller, cross-entropy method, nonlinear dynamics, genetic algorihm

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The idea of the fractional calculus was introduced in 1695, when L'Hôpital in a letter to Leibniz questions what would be the mathematical interpretation for the notation of the derivative d n f /dx n (created by Leibniz) if it has a non-integer order (n = 1/2, for example). The search to find this interpretation involved many mathematicians like Euler, Fourier, Laplace, among others. Its first practical application was made in 1823 by Niels Henrik Abel, to solve problem of tautochrone curve, for which the time of descent of a body abandoned on it and subject to the action of gravity is the same regardless of the starting point where the body is abandoned [START_REF] Monje | Fractional-order Systems and Controls -Fundamentals and Applications[END_REF][START_REF] Ortigueira | What is a fractional derivative[END_REF].

During the following centuries, many pure and applied mathematicians contributed to the development of the theory of fractional calculus and many different fractional operators were proposed, including the fractional derivatives of Grünwald-Letnikov, Riemann-Liouville and Caputo [START_REF] Li | Remarks on fractional derivatives[END_REF].

The Riemann-Liouville fractional integral [START_REF] Li | Remarks on fractional derivatives[END_REF] of order α, where α ∈ IR and α > 0, is defined in terms of a convolution type operation between the real-valued function y(t) and the kernel t α-1

I α a,t y(t) = 1 Γ (α) t a (t -τ ) α-1 y(τ ) dτ , t > a (1)
being Γ the Gamma function, and the Riemann-Liouville fractional derivative [START_REF] Li | Remarks on fractional derivatives[END_REF] is defined in terms of the classical derivative of order n ∈ Z Z + of this fractionalorder integral

RL D α a,t y(t) = 1 Γ (n -α) d n dx n t a (t -τ ) n-α-1 y(τ ) dτ , (2)
which is a global operator (not local as in classical calculus) that presents "memory", being t > a and n -1 < α < n.

The development of new definitions is an active field of research in mathematics [START_REF] Caputo | A new definition of fractional derivative without singular kernel[END_REF][START_REF] Zheng | A new definition of fractional derivative[END_REF], including unified definitions [START_REF] D'elia | A Unified Theory of Fractional, Nonlocal, and Weighted Nonlocal Vector Calculus[END_REF]. On the other hand, the use in applications was possible mainly due to the development of numerical methods to simulate fractional systems [START_REF] Tepljakov | FOMCOM: a MATLAB toolbox for fractional-order system identification and control[END_REF], which also allowed the use of fractional operators in engineering analysis [START_REF] Lin | Eigenvalue and eigenvector derivatives of fractional vibration systems[END_REF] and in control theory [START_REF] Tepljakov | Fractional-order Modeling and Control of Dynamic Systems[END_REF], especially when dealing with delays [START_REF] Martelli | Stability of PID-controlled second-order time-delay feedback systems[END_REF] or chaotic nonlinear systems [START_REF] Niu | Analysis of Duffing oscillator with time-delayed fractional-order PID controller[END_REF].

In addition to the favorable feature of fractional-order operators in providing a natural framework for describing phenomena with memory, the fractional exponent in control theory also offers a type of additional degree of freedom to tune a controller, opening up opportunities for further performance improvements in the controller project. In this context, this work proposes a LQR-based control system to stabilize an inverted cart-pendulum system and analyzes the possibility of improving the performance of the controller using fractional integrators, comparing this with the integer-order. Performance is evaluated through indexes related to the control error, the settling time of the system output and the control effort. To find the parameters that result in the best performance for the controller, including the possibility of fractional-order, and find the most efficient method to optimize systems that use fractional control, two methods of global optimization are used, the genetic algorithm [START_REF] Holland | Genetic Algorithms[END_REF] and the cross-entropy method [START_REF] Rubinstein | The Cross-Entropy Method: A Unified Approach to Combinatorial Optimization, Monte-Carlo Simulation and Machine Learning[END_REF][START_REF] Cunha | Enhancing the performance of a bistable energy harvesting device via the cross-entropy method[END_REF]. The inverted pendulum is chosen as a reference because it is a classic control problem, a nonlinear system widely studied using integer-order controllers [START_REF] Kharola | A comparison study for control and stabilisation of inverted pendulum on inclined surface (IPIS) using PID and fuzzy controllers[END_REF], which it has also started to be tested in the fractional control literature [START_REF] Mousa | Optimal fractional-order proportional-integral-differential controller for inverted pendulum with reduced order linear quadratic regulator[END_REF].

Nonlinear dynamic system

An inverted pendulum has its center of mass above its pivot point, so when in the vertical position is in its unstable equilibrium position, and a small disturbance can bring down the pendulum. Maintaining this equilibrium position or at a desired reference angle is done by changing the cart position. The inverted cartpendulum used in this work is shown in Figure 1 (left), where M is the cart mass; m is the inverted pendulum mass; L is the distance from the center of pendulum's mass to the fixation point; x is the cart's horizontal displacement; θ is the angle between the pendulum and a perpendicular axis through the cart's centroid; and u is the force applied to the cart to control the system. The mass of the rod is neglected, the surface is considered to be frictionless and the mechanical joint is considered smooth.

The dynamics of the system evolve according to the nonlinear model

(J + m L 2 ) θ -m g L sin θ = -m L ẍ cos θ (3)
(m + M ) ẍ + m L θ cos θ -m L θ2 sin θ = u ( 4
)
where J is pendulum moment of inertia. As can be seen, the dynamics of the plant is nonlinear and a suitable controller is needed to keep the pendulum in a certain position. In this study, are considered the values of m equal to 0.1 kg, M equal to 2 kg, L is the length of 0.5 m and the moment of inertia J is equal to 0.006 kg.m 2 . The upper dot is an abbreviation for time-derivative.

3 Control strategy In this proposed controller, r is the reference input signal; e is the tracking error; u is the control signal; K 1 , K 2 , K 3 , K 4 , K l are the control gains; I α1 , I α2 , I α3 , I α l are the integrators, if α ∈ Z Z (integer-order) and if α / ∈ Z Z (fractionalorder). In this plant the feedback is done through two state variables, angular speed θ and cart speed ẋ, chosen to be the observed states. Therefore, to obtain the angular position of the pendulum and the position of the car, these variables are integrated by I α1 and I α3 , respectively, as defined by the fractional integral of Riemann-Liouville in Eq.( 2). In addition, a feedback in the input through the integral I α2 and I α l of the car position is performed.

e I α l K l + - u θ θ x ẋ I α1 K 1 K 2 K 3 I α3 K 4 I α2
To know the performance of an integer-order controller, used as a reference in the comparisons of this work, the Pole Placement method was used to obtain the controller gains K i . This is a classic and widely used method, developed for linear systems using integer-order controllers, where the location of the poles is chosen based on the desired performance of the closed-loop system. Since this method is based on the linear system hypothesis, the inverted cart-pendulum model is linearized. Using the Matlab to calculate the gains, the following values are obtained:

K 1 = -200.6; K 2 = -50.3; K 3 = -70.1; K 4 = -46.8; K l = -63.8.

Controller performance

Performance indexes of control systems are indicators of the quality of the elaborated project, being possible to evaluate quantitatively the performance of different controllers. Mathematically, a performance index is a function of system states, inputs and time. In addition, performance indices can also be used as an objective function for obtaining controllers. In this work, the performance of the controllers is evaluated using the following criteria:

-Integrated square error:

ISE = t 0 e(t) 2 dt -Settling time (ST) -Integrated square control signal: ISU = t 0 u(t) 2 dt
The first two indices refer to the state variables (θ, θ, x and ẋ). The ISE index, is related to the area displaced by the car or the pendulum during the evaluated period until reaching the reference. The settling time, as the name says, is the time it takes the signal to reach the reference. The latter, the ISU index, quantifies the control effort during the evaluated period.

Controller optimization problem

LQR-based fitness function

The objective of this work is to find the set of controller parameters that maximize the performance of the controller. Thus, the design variables chosen for the optimization problem are the gains and the order of the integrators, especially the latter that define whether the controller it will be fractional or integer.

The lower and upper bounds of the variables must be defined according to the computational capacity of the optimization method and in a way that allows the search for the global without one or more variables being close to the stipulated limit. This problem with the bounds of the variables is one of the reasons for using more than one optimization method.

As explained in section 3.2, the performance evaluation is made through the indexes chosen here, however for the optimization process, evaluating the indexes separately would make a multi-objective optimization of high computational cost, since the number of design variables is relatively high. Therefore, it is necessary to develop an objective function that takes into account all these performance indices, including the one that assesses the control effort. Control theory has an optimized controller design technique well known in the literature called linear quadratic regulator (LQR). This technique is used in linear systems and has a quadratic cost function that includes the state variables and the control signal, in order to find the optimized gains for the control project [START_REF] Anderson | Optimal Control: Linear Quadratic Methods[END_REF].

A fitness function based on the LQR was developed mainly because it also takes into account the state variables and the control signal in the formulation

F = τ 0 w 1 x 2 (t) ISE x I + w 2 θ 2 (t) ISE θ I + w 3 u 2 (t) ISU u I + w 4 ẋ2 (t) ISE ẋI + w 5 θ2 (t) ISE θI dt ( 5
)
where t is the time variable, τ is the final evaluation time and the w's are the weights for each component of the objective function. It is worth mentioning that as the reference for the control is zero, the state variables used directly in the function quantify the error, one of the performance indices evaluated in this article. State variables are divided by the respective performance value (ISE and ISU ) in integer-order control, minimizing the difference in order of magnitude of each component, trying to allow the weight to be defined as much as possible by w 1 = w 2 = w 3 = 0.3 and w 4 = w 5 = 0.05. Thus, the value of the fitness function for the classical controller is equal to 1. Therefore, to obtain a better performance F < 1 must be obtained, that is, the objective function must be minimized. To turn this minimization problem in a maximization problem, as is treated in the optimization methods applied in the article, it is considered S = -F. The constraints are incorporated into the objective function through a penalty function

P = 100 (max(0, ST θ -ST * θ )) 2 + 100 (max(0, ST x -ST * x )) 2 (6)
where ST θ and ST x are the settling time of the angular position and the position of the car, ST * θ and ST * x are the respective values in the classic controller.

Genetic algorithms

Genetic algorithms (GA) are heuristic search approaches applicable to a wide range of optimization problems. Invented by Holland (1975[START_REF] Holland | Genetic Algorithms[END_REF], this stochastic optimization and global search technique are successful methods for use in problems with difficult solution spaces, for example, if no derivatives are available and if the fitness scenario suffers from poorly conditioned parts [START_REF] Holland | Genetic Algorithms[END_REF]. Using the principle of genetics and natural selection, evolution is the basis of genetic algorithms, together with their main genetic operators, which are crossover, mutation and selection. Crossover operators combine the genomes of two or more solutions. The mutation operator adds randomness to the solutions. And finally, the selection operator chooses the best solutions in a population for survival [START_REF] Kramer | Genetic Algorithm Essentials[END_REF].

Cross-Entropy method

Global search methods like genetic algorithms are quite effective in overcoming local minimums or maximums, unlike gradient-based methods, however its control parameters are not intuitive and most of the time its tunning is done through trial and error, which can lead to loss of performance and accuracy. In an attempt to overcome such difficulties and succeed in the search for the controller with the best performance, this work uses, in addition to the GA method, another global search algorithm known as the cross-entropy (CE) method. More robust and simple, it was proposed by R. Rubinstein in 1997 [START_REF] Rubinstein | The Cross-Entropy Method: A Unified Approach to Combinatorial Optimization, Monte-Carlo Simulation and Machine Learning[END_REF] initially for simulation of rare events, and then its effectiveness in application in optimization problems was observed [START_REF] Kroese | Cross-Entropy Method[END_REF][START_REF] Rubinstein | Simulation and the Monte Carlo Method[END_REF]. The idea of the method is to transform a non-convex optimization problem into a rare event estimation problem, that can be solved by Monte Carlo sampling technique. The process consists of: (i) sample the feasible region according to a given probability distribution; (ii) evaluate the objective function in each of these samples; (iii) identify the samples that produced the highest values for the objective function, this subset being defined as the elite sample set; (iv) update the parameters of the probability distribution based on the mean and standard deviation of the elite sample set, modifying the distribution in order to try to make it as close as possible to the global optimum; (v) repeat steps ( 2) to ( 4) as long as a stop criterion is not met. This iterative process can be classified into two stages: sampling (i and ii) and learning (iii to v). Furthermore, this process only requires the user to define the number of samples N s , the number of elite samples N e <N s , a tolerance of convergence tol and the maximum number of levels (iterations) l max [START_REF] Cunha | Enhancing the performance of a bistable energy harvesting device via the cross-entropy method[END_REF].

Numerical experiments

The numerical experiments presented here, consist of the execution of optimizations using GA and CE methods through numerical simulations in Matlab, where several evaluations of the objective function are performed and the ideal parameters of the controller that result in its best performance are found. For the operations of fractional integrators the FOMCON toolbox is employed [START_REF] Tepljakov | FOMCOM: a MATLAB toolbox for fractional-order system identification and control[END_REF].

The initial conditions for the numerical simulations are defined: x 0 = 0 m; ẋ0 = 0 m/s; θ 0 = 10×π/180 rad; θ0 = 0 rad/s. Note that the initial condition of the angular position other than 0 (zero) will force the control system to act to find stabilization. The dynamics are integrated over the time interval [0; τ ] = [0; 10] s.

The optimizations will be divided into two stages according to the choice of design variables to be analyzed. In the first round of optimization, it seeks to find the optimal values of α 1 and α 3 , which results in the best performance for the controller, keeping the other parameters fixed. The strategy is to use the controller's gains, obtained by the linear Pole Placement method and consider the order of the integers α 2 and α l equal to 1. Then, optimizations are performed with all system parameters (gains and integrators) to find the one that results in the best control performance.

The α 1 and α 3 optimization

Initially, the GA and CE optimizations are made with the design variables being α 1 and α 3 , keeping the other controller parameters fixed and equal to 1. The results obtained in each of these optimizations are presented in Table 1. The first highlight can be made on the efficiency of the CE method compared to GA to obtain the same result more quickly and with less computational cost. Another observation is that the use of fractional-order integrators can improve the performance of the controller, with an optimization close to 10% of the overall performance (based on the fitness function) of the controller. The ISE of the angular position had a reduction of 20% compared to the integer-order. The ISE of the angular velocity, the settling time of x and the control effort were slightly higher, however, compared to the improvement obtained in the others, they are irrelevant. Some of these results are highlighted in Figure 2, with the time series of the angular position, car position and control effort.

Optimization of all parameters

Based on the results of section 5.1, the following question is asked: "if the gains are also optimized, will the integrators also have a fractional order?". To answer this question, new optimizations are developed with the GA and CE methods, now with all the parameters of the controller as design variables. The results are compared with the integer-order controller, as done previously. Table 2 shows the configuration of the methods, the bounds of all design variables and the performance of each method for optimization evaluated with all parameters. Unlike the previous optimization, the optimizations performed with all parameters required a higher computational cost. However, as shown in Table 3, the results of the performance indices and the fitness function of the optimized parameters show that the use of fractional-order integrators together with the optimal gains can further improve the performance of the controller designed through the method of Pole Placement. The optimization methods GA and CE allowed results in general (fitness function) with 27% and 37% better than the controller used as reference, respectively. This difference in methods is linked to the better efficiency of the CE in relation to the operational cost (time/function evaluates). Because in GA, due to the CPU time and the number of variables in the problem, it was necessary to establish lower bound limits of each parameter/design variable for the method to run. This problem can be seen in the values of K 1 and K 3 , which were very close to the lower bounds established for them. Another important highlight is how errors can be reduced by up to 50% with fractional integrators, yet requiring 9% less effort from the controller. The time series of the angular position, the car position and the control effort are shown in Figure 3, with the comparison of the result of each method with the integer controller. This work used the stability problem of an inverted cart-pendulum to develop a study on the performance of fractional controllers compared to those with integer-order integrators. So there is a nonlinear non-convex optimization problem that aims to maximize the performance of the control system, taking into account the error, the settling time and the control effort. For this, two different optimization methods were used: genetic algorithms and cross entropy. The latter being faster and more robust for optimization that involves fractional-order integrators, i.e. the CE method presents better performance for problems with many variables, which is the case of systems with fractional control.

As in [START_REF] Tepljakov | Fractional-order Modeling and Control of Dynamic Systems[END_REF], there is an improvement in the performance of the control circuit that has fractional-order operators, but in this work for a LQR-based controller. The results show that the use of controllers with fractional integrators instead of integer-order controllers makes the control more efficient and faster, approximately 20% faster to reach the benchmark (θ = 0 rad), with less error and requiring less effort of control. With the reduction of all control parameters, integrators and gains, the fitness function can be reduced by more than a third, allowing to find a better answer from the results of the Pole Placement method.

3. 1 Figure 1 (

 11 Figure 1 (right) shows the LQR-based controller proposed to control the inverted cart-pendulum system.

Fig. 1 .

 1 Fig. 1. Schematic illustration of the cart-pendulum system (left), and the proposed LQR-based controller (right)

Fig. 2 .

 2 Fig. 2. Time series of the angular position (top-left), car position (top-right) and the control signal (bottom) for results firsts optimization

Fig. 3 .

 3 Fig. 3. Time series of the angular position (top-left), car position (top-right) and the control signal (bottom) for results optimization of the all parameters

Table 1 .

 1 Optimal values of the design variables (α1 and α3) with their respective performance indexes and fitness function value

			integer	fractional	fractional
			controller	controller GA	controller CE
		α1	1.00		1.03		1.03
		α2	1.00		1.00		1.00
		α3	1.00		0.93		0.93
		α l	1.00		1.00		1.00
	parameters	K1	-200.60		-200.60		-200.60
		K2	-50.30		-50.30		-50.30
		K3	-70.10		-70.10		-70.10
		K4	-46.80		-46.80		-46.80
		K l	-63.80		-63.80		-63.80
		ISE θ	0.021	0.017	-20%	0.017	-20%
		ISE θ	0.190	0.204	7%	0.204	7%
		ISE x	0.257	0.230	-11%	0.230	-11%
	performance ISE ẋ	0.660	0.509	-23%	0.509	-23%
	index	ISU	23.659	24.412	3%	24.411	3%
		ST θ	4.873	4.088	-16%	4.088	-16%
		ST x	5.073	5.100	1%	5.100	1%
	fitness function	1.000	0.909	-9%	0.909	-9%

Table 2 .

 2 Bounds and performance of GA and CE algorithm for optimization evaluated

	method	α 1,2,3,l	bounds K1	K 2,3,4,l	levels	function evalution
	GA	[0.8 1.1] [-250 -150] [-90 -30]	327 generations 29520
	CE	[0.5 1.5] [-350 -1]	[-150 -1]	88 iterations	4400

Table 3 .

 3 Controller configuration with the optimal values of the design variables (all parameters) and their respective performance indexes and fitness function value

			integer	fractional	fractional
			controller	controller GA	controller CE
		α1	1.00		0.96		0.72
		α2	1.00		0.94		0.79
		α3	1.00		0.91		1.22
		α l	1.00		0.97		0.80
	parameters	K1	-200.60		-249.98		-205.68
		K2	-50.30		-60.08		-86.97
		K3	-70.10		-89.96		-67.54
		K4	-46.80		-70.35		-124.60
		K l	-63.80		-83.01		-65.56
		ISE θ	0.021	0.012	-45%	0.009	-55%
		ISE θ	0.190	0.169	-11%	0.157	-17%
		ISE x	0.257	0.182	-29%	0.136	-47%
	performance ISE ẋ	0.660	0.357	-46%	0.283	-57%
	index	ISU u	23.659	22.158	-6%	21.497	-9%
		ST θ	4.873	3.952	-19%	3.949	-19%
		ST x	5.073	4.899	-3%	5.100	1%
	fitness function	1.000	0.730	-27%	0.630	-37%

Acknowledgements

The authors would like to thank the financial support given to this research by the Brazilian agencies Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES) -Finance Code 001, the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and the Carlos Chagas Filho Research Foundation of Rio de Janeiro State (FAPERJ).