Noureddine Ait Said
email: noureddine.ait-said@univ-grenoble-alpes.fr

Mounir Benabdenbi
email: mounir.benabdenbi@univ-grenoble-alpes.fr

Katell Morin-Allory
email: katell.morin-allory@univ-grenoble-alpes.fr

Arbitrary Reduced Precision for Fine-grained Accuracy and Energy Trade-offs

Keywords: Approximate Computing, Transprecision Computing, FPU bit-width optimization, RISC-V, ASIP, ASIC

Full-precision Floating-Point Units (FPUs) can be a source of extensive hardware overhead (power consumption, area, memory footprint, etc.). As several modern applications feature an inherent tolerance to precision loss, a new computing paradigm has emerged: Transprecision Computing (TC). TC proposes several tools and techniques that trade precision for energy efficiency. However, most of these tools require developers to rewrite part or all of their existing software stacks, which is often infeasible, complex, or requires extensive development efforts. In addition to their intrusiveness, TC tools can only simulate the impact of precision loss, and they do not provide corresponding hardware designs that take advantage of the simulations.

This work proposes a non-intrusive hardware-oriented approach, requiring no modification of source code that applies approximations at the low-level in assembly. The approach can be used to approximate virtually all types of executable binaries (bare-metal applications, single-/multi-threaded user applications, OS/RTOS, etc.). We introduce AxQEMU: a software based on the well known QEMU dynamic binary translator. We demonstrate how our approach can determine the effects of FP approximations on application-level Quality of Result (QoR), and how it interfaces with other tools from the literature. A hardware-level case study on a 28-nm FD-SOI implementation is presented, demonstrating how fine-grained energy/accuracy trade-offs can be made thanks to floating-point arbitrary reduced precision (ARP). For instance, considering the well-known arclength FP application, FPU computation energy savings of up to 19.4% were achieved with an accuracy threshold of 10 significant digits, and up to 60.7% with a 4-digit accuracy when using ARP. These savings compared favorably to the limited 7.7% saving afforded by using standard variable type optimization tools.

Introduction

Approximate computing (AC) [START_REF] Xu | Approximate computing: A survey[END_REF][START_REF]A survey of techniques for approximate computing[END_REF][START_REF] Moreau | A Taxonomy of General 770 Purpose Approximate Computing Techniques[END_REF] is a paradigm based on the idea that the numerical accuracy of data and computations can be traded off against power savings, performance gains, and area/resource optimization. AC is suitable for non-critical RMS applications (Recognition, Mining, and Synthesis) that tolerate a certain reduction in their output precision (e.g. image/signal processing, computer vision, machine learning).

Current trends are pushing towards more adaptive, 10 variable, and mixed-precision computing known as Transprecision Computing (TC) [START_REF] Malossi | The transprecision computing paradigm: Concept, design, and applications[END_REF]. This paradigm targets mainly floating-point (FP) computations and storage, and the goal is to design more flexible, run-time variable precision and efficient architectures where accuracy is balanced against cost and resource savings.

In this paper, we focus on Floating-Point Units (FPUs), which are ubiquitous in modern hardware architectures, including General-Purpose Processors (GPP) and Application Specific Instruction-Set Processors (ASIP), where 20 they serve to boost the performance of computationallyintensive applications. Unfortunately, an FPU occupies a significant proportion of the CPU's core area and can cause extensive power consumption and high memory bandwidth usage. The energy consumption associated with FP arithmetic is known to be higher than that of its integer counterpart [START_REF] Matoussi | Error Analysis 780 of the Square Root Operation for the Purpose of Precision Tuning: a Case Study on K-means[END_REF], making FPU optimization a priority.

One of the techniques used to optimize FPUs is bitwidth reduction, where exponent and/or mantissa lengths are reduced to either standard bit-widths (defined in the 30 IEEE 754 standard [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF]), or custom arbitrary bit-widths. Over the years, many techniques/tools/libraries have been proposed to explore the impact of using arbitrary precision FP arithmetic [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with 790 correct rounding[END_REF][START_REF] Lefèvre | Sipe: a Mini-Library for Very Low Precision Computations with Correct Rounding[END_REF][START_REF] Tagliavini | FlexFloat: A Software Library for Transprecision Computing[END_REF][START_REF] Flegar | FloatX: A C++ Library for Customized Floating-Point Arithmetic[END_REF] in computational kernels targeting many platforms (GPPs, GPUs, FPGAs, etc.).

The study presented here contributes to this goal by introducing an extension to [START_REF] Said | FPU Bit-Width Optimization for Approximate Computing: A Non-Intrusive Approach[END_REF]. It presents a non-intrusive methodology enabling rapid design space exploration of reduced arbitrary precision FPUs in a CPU-based architecture context. The methodology proposed should help designers select the most optimized FPU configuration (exponent and mantissa bit-widths) satisfying a Quality of Result (QoR) threshold set for the application and input dataset provided by the designer, without the need to transform / rewrite / modify the source code. The approach led to the development of a software implementation which we called AxQEMU. This implementation is based on the dynamic binary translation of assembly instructions, and is built on top of the well-known functional simulator QEMU. The methodology was assessed through three case studies: 1) a direct application of AxQEMU alone, 2) using AxQEMU in conjunction with an FP variable type optimization tool from the literature, and 3) a hardware-level estimation of the energy savings provided by our methodology.

Section 2 introduces our motivations, the terminology used in the remainder of the paper, and the problem statement. Section 3 presents a formalization of the idea. Section [START_REF] Malossi | The transprecision computing paradigm: Concept, design, and applications[END_REF] shows the details of the software implementation. In Section 5, we assess the methodology developed through 60 three detailed case studies. Section 6 presents a set of related works followed by a discussion.

Background & Motivation

Terminology

A binary floating-point number can be written in the form (-1) s × (1 + m) × 2 e , where s is the sign bit, m is the mantissa (also called significand or fraction), and e is the exponent. Each of these components can be encoded either following standard IEEE 754 [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF] formats e.g., binary32 (32-bit single-precision format), binary64 (64-bit doubleprecision format), or using a custom bit-width representation. Variations in the mantissa bit-width change the precision of the number representation, whereas alterations to the exponent bit-width change its dynamic range. Custom non-standard (arbitrary) formats can be defined when some loss of precision is tolerated, or when the numbers represented have a limited range.

For this work, we focused on the software Application Binary Interfaces (ABI) that support FP arithmetic in hardware (hard-float ABI). We assumed that the ABI and the hardware FPU support at least two standard types: binary32 and binary64.

An FPU configuration is a 4-uple (E f , M f , E d , M d )
, where E f (resp. E d ) is the exponent bit-width for the single-precision (resp. double-precision) operator, and M f (resp. M d ) is the mantissa bit-width of the single-precision (resp. double-precision) computational operator. For example, an FPU that supports binary32 and binary64 is represented as [START_REF] Lefèvre | Sipe: a Mini-Library for Very Low Precision Computations with Correct Rounding[END_REF][START_REF] Waterman | I: User-Level ISA[END_REF][START_REF] Said | FPU Bit-Width Optimization for Approximate Computing: A Non-Intrusive Approach[END_REF]52). Typical hardware implementations contain an additional bit, or hidden bit, which is set to 1 by default and set to 0 when manipulating ±0.0 or denormal numbers [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF], when supported. The bit-widths M f and M d do not consider the hidden bit.

Context of the proposed approach

A typical FP algorithm implementation scenario consists mainly of three steps:

1. Algorithm design and numerical stability analysis [START_REF] Panchekha | Automatically improving accuracy for floating point expressions[END_REF][START_REF] Sanchez-Stern | Finding Root Causes of Floating Point Error with Herbgrind[END_REF][START_REF] Benz | A Dynamic Program Analysis 820 to find Floating-Point Accuracy Problems 10[END_REF]: establishing the mathematical foundations of the algorithm and their stability w.r.t the inputs. The main goal is to avoid instabilities such 100 as round-off errors. 2. Conservative (naive) implementation: a software implementation using all high precision formats. 3. Variable type optimization (VTO) [START_REF] Graillat | Auto-tuning for floating-point precision with Discrete Stochastic Arithmetic[END_REF][START_REF] Lam | Automatically adapting programs for mixed-precision floating-point computation[END_REF][START_REF] Rubio-González | Precimonious: Tuning assistant for floating-point precision[END_REF][START_REF] Ho | Efficient float-840 ing point precision tuning for approximate computing[END_REF]: the process of migrating the maximum possible number of variables from high-precision to lower precisions (e.g., changing double variables to float) in a given application source code while satisfying a Quality of Result (QoR) constraint. The process output is a (software) type configuration: a version of the 110 original application, possibly with mixed-precision, where some variables are declared as floats and others as doubles.

Although VTO is necessary to optimize the memory footprint and energy consumption of a program, it is a coarsegrained optimization process. It usually results in solutions that are over-designed for many application classes.

In our approach, we have therefore included a fourth step to minimize hardware FPU implementation through a finegrained optimization: 120 4. Fine-grained optimization using Arbitrary Reduced Precision (ARP): this approach takes advantage of non-standard reduced operators.

In this work, we demonstrate that ARP can produce numerical outputs that are as accurate as standard formats while consuming less energy.

Problem statement

Given an application that takes as an input a dataset I, and providing a numerical output result O, our objective was to identify the optimal FPU configuration (E f , M f , E d , 130 M d ) in terms of power consumption, execution time, and overall energy consumption, subject to a Quality of Result (QoR) constraint on the output O.

To reach this goal, we adopted a three-step method:

1 Starting from an application source code, steps 1 to 3 were performed using state-of-the-art tools [START_REF] Panchekha | Automatically improving accuracy for floating point expressions[END_REF][START_REF] Sanchez-Stern | Finding Root Causes of Floating Point Error with Herbgrind[END_REF][START_REF] Benz | A Dynamic Program Analysis 820 to find Floating-Point Accuracy Problems 10[END_REF][START_REF] Graillat | Auto-tuning for floating-point precision with Discrete Stochastic Arithmetic[END_REF][START_REF] Lam | Automatically adapting programs for mixed-precision floating-point computation[END_REF][START_REF] Rubio-González | Precimonious: Tuning assistant for floating-point precision[END_REF][START_REF] Ho | Efficient float-840 ing point precision tuning for approximate computing[END_REF] to generate a valid type configuration. 2 Software design space exploration: for a given type configuration, the design space must be rapidly explored to select a set of FPU bit-width configurations 140 satisfying the QoR constraint target. 3 FPU hardware assessments: for each configuration, an estimation of the overall energy savings was made to select the best final configuration.

These different steps can be automated using specific search algorithms. However, in this paper, for the sake of clarity, simulations were performed exhaustively.

Proposed Approach

In this section, we introduce some definitions along with a formalization of the approach in the form of an FSM model, the underlying details of status flag computation, and selective approximation.

Definitions & Notations

Let A be the set of FPU registers, R be the set of available rounding modes. Let X SP = {L 1 , L 2 , ...} be the set of Standard Precision formats supported (e.g., binary64, binary32), and X ARP = {l 1 , l 2 , ...} be a set of non-standard ARP floating-point formats. In the following, we denote L a format in X SP , l a format in X ARP and a format in X SP ∪ X ARP . Let F denote the set of FP numbers that can be represented in the format.

Let F be the set of FP instructions available for a given ISA. This set is partitioned in two sets F approx and F exact , where F exact represents the non-computational FP instructions (comparison, conversion, loads/stores and sign injection), whereas F approx represents the FP computational instructions: addition, subtraction, multiplication, fused multiplication-addition, square root, and division 1 .

We define the precision reduction function Reduce R L,l : F L → F l that reduces an L-bit number to an l-bit one using rounding mode R ∈ R. We also define the reverse function Extend L,l : F l → F L , that extends an l-bit number to produce an L-bit one.

Let F inst R rd, rs 1 , . . . , rs n be an FP instruction, where F inst is in F, and operates on n + 1 registers rd, rs 1 , . . . , rs n , of -bit values, using the rounding mode R. With each F inst R we associate an arithmetic operator F op R : (F ) n → F which performs the computation in the format using the rounding mode R.

To introduce selective approximation, as a means to apply approximations to a specific set of instructions, let B be the Boolean domain and do approx ∈ B an enable signal associated with F inst R which indicates whether the instruction should be approximated at run-time.

Approach formalization

The FPU can be modeled by an FSM {I, Γ, γ 0 , δ} with no output, where I = F ×R×X SP ×A n+1 ×B is the input of the FSM, representing an instruction instance characterized by its name, its rounding mode, its original FP format, its destination and source registers, and an approximation enable signal. Γ is the set of states 2 , γ 0 is the initial state where all the registers are set to zero, and δ the state-transition function, which is defined as follows:

δ : I × Γ -→ Γ (F inst R
L rd, rs1, . . . , rsn, do approx), γ -→ γ

1 Please note that in contrast to our definition, the IEEE 2008-754 Standard's [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF] computational instructions also include comparisons.

2 A state is defined by the content of the FPU registers r i ∈ A.

such that ∀r ∈ A: where,

γ (r) =          γ(r), if r = rd; Extend L,l F op R l (ṽ 1 , . . . ,
vi = γ(rsi) ∈ FL high-precision values. ṽi = Reduce R L,l (γ(rsi)) ∈ F l reduced precision values.
The assembly-level approximation is introduced by performing the computation using the ARP operator F op R l , which operates on ARP l-bit operands , with a bit-width shorter than L. Consequently, inputs should first be reduced from L to l before performing the computation, and then extended back from l to L following the computation.

Precision reduction performs a cast using a rounding 200 mode R equal to or different from R, depending on how it is implemented. In most ISAs, R is encoded either in the instruction binary code (fixed) itself, or in an FPU control status register (dynamic). Hence, an implementation can either define a fixed reduction rounding mode R , or support dynamic rounding. We chose to leave the choice to the implementation, because rounding hardware logic is generally expensive in an FPU, and in addition double rounding (i.e., rounding to reduced format plus rounding the final result) may drastically affect the QoR. When the 210 original format L and the reduced format l have similar exponent bit-widths, the Reduce R L,l function is simplified to a rounding operator. Otherwise, a complete cast operator should be implemented, which may be expensive in terms of circuit area.

Once the computation F op R l has been performed, the result should be converted back to the original format using the Extend L,l function which converts the result from l back to the L format. This operation is intended to guarantee consistency with the non-computational FP in-220 structions F exact , which are not approximated. The final result is then stored in the destination register rd.

Computation of the status register

In addition to the numerical result computed, an FPU also returns a 5-bit status register. Since the FPU is modified to take approximations into consideration, the five flags are redefined as follows:

• Inexact bit (NX): set if the final result cannot be represented precisely in the current representation.

If the instruction is executed approximately, the in-

230
exact bit is set if at least one of the reduction operations yields an inexact reduced operand ṽi , or if the result of F op R l is inexact. The extension does not affect this flag.

• Invalid bit (NV): signals an invalid FP operation (e.g., multiplying ∞ and zero). This flag is set if F op R l performed an invalid FP operation. The reduction and extension stages do not affect this flag.

• Division by zero (DZ): indicates an operation in-

volving division by zero. It is similar to the DZ flag resulting from F op R l . • Underflow (UF): indicates if at least one of the reduction results or the result of F op R l is too small to be represented in the l format.

• Overflow (OF): indicates if the result cannot be represented as a finite value in l. It has the same effect as UF.

The case of iterative operators

Some FPUs implement iterative operators F op R , particularly for division and square root calculations [START_REF] Li | Implementation of single precision floating point square root on fpgas[END_REF] [START_REF] Mach | A transprecision floating-point architecture for energy-efficient 850 embedded computing[END_REF],

to allow circuit area optimization. For such implementations, the reduction and extension stages are unnecessary since the precision of their computations can be set at runtime.

Selective Approximation (SA)

The signal do approx was added to support SA. This feature allows the programmer to apply approximations to specific parts of an application and exclude others. For instance, in a software application, the functions that compute error metrics should be executed precisely to avoid compromising the results with approximations. Hence, the developer should be able to tag such functions or instructions.

Implementation

To implement a proof of concept of the approach, a functional ISA simulator (e.g., gem5, or1k, Spike, etc.) can be used. Here, we targeted the free and open-source QEMU multi-ISA dynamic binary translator [START_REF] Bellard | Qemu, a fast and portable dynamic translator[END_REF].

QEMU overview

QEMU supports most well-known architectures (RISC-V, ARM, x86-64, MicroBlaze, etc.). We distinguish two machines:

1. the guest (or target) machine: the processor emulated, for which an application has been compiled, and 2. the host machine: the processor executing QEMU itself, which simulates the execution of the target code, even if the two processors have different architectures Figure 1 depicts QEMU's operation principle. Efficient emulation is achieved thanks to the 'Tiny Code Generator' (TCG), which is a sort of run-time compiler embedded in QEMU. This generator dynamically translates blocks of target instructions known as Translation Blocks (TB) to TCG operations (TCGops), which constitute a machineindependent intermediate representation (IR). Subsequently, the TCGops will be translated into host instructions.

When a TB is translated into its corresponding host code, the translated block is cached for later use. For example, in the case of a loop, the TB is translated only 290 once, but executed multiple times. This caching capability is one of several optimization procedures that increase the performance of this simulator when compared to others. 1) that will eventually be compiled to produce several host instructions. When translating floating-point instructions, the corresponding TCGops can either be executed directly 300 on the host hardware's FPU or emulated using C helper functions exploiting a backend FP software emulation library [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with 790 correct rounding[END_REF][START_REF] Tagliavini | FlexFloat: A Software Library for Transprecision Computing[END_REF][START_REF] Flegar | FloatX: A C++ Library for Customized Floating-Point Arithmetic[END_REF][START_REF] Hauser | SoftFloat release[END_REF]. In other words, each guest FP computational assembly instruction is associated with a custom C function which simulates the reduced precision behavior in the back-end.

Figure 1 shows an application source code that has been compiled for a RISC-V target machine. Each assembly instruction from the executable binary will be fed to the TCG to generate its corresponding intermediate repre-310 sentation (TCGop). For example the RISC-V FP addition instruction (fadd.d fa5, fa5, fa4), which should perform the operation fa5 = fa4 + fa5, is mapped to the C function helper fadd d() 3 , which operates on the values contained in emulated guest registers fa5 and fa4. The helper functions were modified to implement the behavior explained in Section 3.

Implementation of SA

SA allows developers to specify functions that should be executed precisely to facilitate integration with exist-320 ing flows. In our case, SA is supported by splitting the memory address space across two regions:

1. Non-Approximable address space: a memory region where approximations are deactivated (i.e., do approx = 0), and 2. Approximable address space: where approximations are enabled using ARP (i.e., do approx = 1).

To implement this split, the non-approximable functions requiring precise execution are annotated using a C macro. For instance, Listing 1 depicts a function 330 compute QoR, destined to compute some QoR metric using a numerical result result and a reference value ref.

It has been annotated with the PRECISE macro defined in lines [START_REF] Xu | Approximate computing: A survey[END_REF][START_REF]A survey of techniques for approximate computing[END_REF]. This macro tells the linker to place the annotated function in a particular memory section, the ".precise" section. If an assembly function belongs to the non-approximable address space (i.e., the .precise memory section) then 340 do approx is set to zero and the instruction is emulated using the standard SoftFloat [START_REF] Hauser | SoftFloat release[END_REF] FP emulation library. Otherwise, do approx is asserted, and the instruction is emulated using an ARP library [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with 790 correct rounding[END_REF][START_REF] Tagliavini | FlexFloat: A Software Library for Transprecision Computing[END_REF][START_REF] Flegar | FloatX: A C++ Library for Customized Floating-Point Arithmetic[END_REF], taking the FPU configuration specified by the user into account.

As for all similar techniques, the source code must be modified and re-compiled when using SA. However, when using our approach, the implementation is much simpler and light-weight than other current techniques, as well as being minimally-intrusive.

Key engineering decisions

For implementation, the FlexFloat library was chosen for approximation emulation, since according to [START_REF] Tagliavini | FlexFloat: A Software Library for Transprecision Computing[END_REF] it demonstrates an increase in speed of up to 2.8x compared to other mixed-precision libraries. The library provides functions for casting from full-precision types to reducedprecision types.

The reduction rounding mode R is dynamic and follows the same rounding mode R of the instruction being approximated. The RISC-V ISA [START_REF] Waterman | I: User-Level ISA[END_REF] was chosen for its 360 design simplicity. We implemented the approach for all the computational instructions of single-precision (F) and double-precision (D) extensions of the ISA.

The C helpers were modified to accept two additional arguments: the exponent and the mantissa bit-widths which define the precision of the internal computations. These arguments have been exposed so that the user can define them at launch-time. This facilitates use of the simulator when dealing with search algorithms to explore several FPU configurations. 370

Evaluation

In this section, we demonstrate the effectiveness of our methodology through three case studies. The first one is a direct application of the AxQEMU tool introduced in Section 4 on a 3D gaming application. The second demonstrates the use of AxQEMU along with a VTO tool from the literature called PROMISE [START_REF] Graillat | Auto-tuning for floating-point precision with Discrete Stochastic Arithmetic[END_REF], whereas the last one is a hardware-level trade-off study of the energy savings vs. QoR loss.

Case Study 1: Direct application of AxQEMU 380

With this first use case, we demonstrate the direct use of AxQEMU alone to examine how the FP approximation affects a computationally heavy application, Jmeint [START_REF] Yazdanbakhsh | AxBench: A Multiplatform Benchmark Suite for Approximate Computing[END_REF], in the 3D gaming domain. This application detects the intersection of two triangles in space and takes as an input 100,000 random pairs of 3D triangle coordinates. For this application, all variables are declared as doubles, and the QoR metric is the intersection detection miss rate. The source code was instrumented with SA macros as explained in Section 4.3. Figure 2 shows the QoR variation as a 390 function of the exponent E d and mantissa M d bit-widths.

These results show that the QoR is maximal and constant for M d ≥ 22, E d = 4 as well as for M d ≥ 18, E d ≥ 5. Therefore, from a software perspective, the float type is over-designed for this application. If the hardware running the application supports other formats, such as the IEEE half-precision 16-bit format [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF] which has a 5-bit exponent and a 10-bit mantissa, a miss rate of just 6.4% can be achieved alongside significant energy savings and reductions in memory footprint. Using ARP, it is possible 400 to trade accuracy against hardware overheads in a more fine-grained way.

3 3 9 1 1 1 1 1 1 1 7 1 1 1 1 1 1 1 4 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Case Study 2: AxQEMU applied to Mixed-Precision applications

In this case study, we demonstrate how AxQEMU can be used in conjunction with other VTO tools. For this experiment, we selected the tool described in [START_REF] Graillat | Auto-tuning for floating-point precision with Discrete Stochastic Arithmetic[END_REF] which is briefly introduced in Section 5.2.1. The software type configurations described in [START_REF] Graillat | Auto-tuning for floating-point precision with Discrete Stochastic Arithmetic[END_REF] were instrumented with AxQEMU to select the best hardware configurations according to the specified QoR constraint (Phase 2 ).

The PROMISE tool

PROMISE [START_REF] Graillat | Auto-tuning for floating-point precision with Discrete Stochastic Arithmetic[END_REF] 4 was selected as the primary tool for VTO (Phase 1 ). Given an application and an input dataset, it applies a delta debugging algorithm [START_REF] Zeller | Simplifying and isolating failureinducing input[END_REF] to find a configuration that minimizes the number of highprecision (double-precision) variables while nevertheless satisfying a given QoR constraint. Internally, PROMISE uses the CADNA library [START_REF] Jézéquel | CADNA: a library for estimating round-off error propagation[END_REF], which implements Discrete Stochastic Arithmetic (DSA) [START_REF] Vignes | Discrete stochastic arithmetic for validating results of numerical software[END_REF], a technique that estimates round-off error propagation and detects numerical instabilities in a program.

Benchmark QoR Metric

In the remainder of this section, we use the same QoR metric as in [START_REF] Graillat | Auto-tuning for floating-point precision with Discrete Stochastic Arithmetic[END_REF] i.e., the number of significant digits S computed as follows:

S = -log 10 Result (8,23,11,52) -Result (E f ,M f ,E d ,M d )
Result [START_REF] Lefèvre | Sipe: a Mini-Library for Very Low Precision Computations with Correct Rounding[END_REF][START_REF] Waterman | I: User-Level ISA[END_REF][START_REF] Said | FPU Bit-Width Optimization for Approximate Computing: A Non-Intrusive Approach[END_REF]52) Where

Result (E f ,M f ,E d ,M d )
is the numerical output computed when simulating the application with the (E f , M f , E d , M d ) FPU configuration on AxQEMU, and Result [START_REF] Lefèvre | Sipe: a Mini-Library for Very Low Precision Computations with Correct Rounding[END_REF][START_REF] Waterman | I: User-Level ISA[END_REF][START_REF] Said | FPU Bit-Width Optimization for Approximate Computing: A Non-Intrusive Approach[END_REF]52) is the golden reference result obtained using the standard (8, 23, 11, 52) configuration. The numerical results for these applications are non-zero.

Evaluation benchmarks

The approach was evaluated using three applications from [START_REF] Graillat | Auto-tuning for floating-point precision with Discrete Stochastic Arithmetic[END_REF]: arclength, rectangle, and squareroot. For each application, an optimized type configuration is generated by PROMISE for various QoR thresholds (Phase 1 , Section 2.3).

Table 1 lists the benchmarks studied alongside the type configurations generated. The 3 rd (resp. 4 th ) column depicts the number of double (resp. float) variables required for each type configuration. Every {type configuration, benchmark} pair achieves a maximum QoR (shaded cells in the 5 th column) when executed exclusively in full precise mode using the standard configuration [START_REF] Lefèvre | Sipe: a Mini-Library for Very Low Precision Computations with Correct Rounding[END_REF][START_REF] Waterman | I: User-Level ISA[END_REF][START_REF] Said | FPU Bit-Width Optimization for Approximate Computing: A Non-Intrusive Approach[END_REF]52). In other words, if one is limited to standard FP types, only the conservative QoR thresholds that are highlighted can be obtained. The last column depicts, for each {application, type configuration, QoR threshold} tuple, a set of candidate FPU configurations that also satisfy the given QoR constraint. Candidate FPU configurations (the 6 th column) were selected based on a design space exploration process which is explained in the next paragraph.

Design space exploration (DSE) with AxQEMU

An exhaustive DSE was performed for each type configuration and each application. Let us denote C = {c 0 , c 1 , ...} a set of FPU configurations to be studied. Our objective was to determine the most optimized FPU configuration c i ∈ C satisfying the target QoR threshold. The DSE is a two-step process: 460 Simulation. First the source code corresponding to a type configuration (e.g. the V10 configuration of arclength) is compiled. Then, the binary is executed using AxQEMU for each FPU configuration c i ∈ C. The raw QoR values (i.e., number of significant digits) for all c i are stored in an array M for processing in the next step. This simulation step can be more or less time consuming depending on the size of the input dataset, the execution time for the application itself, and the number of FPU configurations to be studied.

470

QoR normalization. Once all the simulations have been performed and the QoR array M constructed, a processing algorithm is applied to filter out anomalous local optimums caused mainly by error cancellation. Indeed, reducing the operator-level precision leads to errors in the output of an operation. However, sometimes, when multiple operations are combined, the errors may cancel each other out leading to an accidental increase in overall QoR rather than a decreased QoR. To keep results consistent, we considered the worst-case error detected (i.e., the lowest QoR) among all 480 the configurations that have an equal or greater mantissa bit-width. Thus, for each c i = (E f i , M f i , E di , M di ) we considered the minimum QoR value between the raw original value obtained for c i and the QoR values determined for any higher-precision configuration c j = (E f j , M f j , E dj , M dj ) that satisfies M dj ≥ M di or M f j ≥ M f i .

Phase (2) results

Phase 2 (from Section 2.3) consists in performing a series of simulations using AxQEMU for each benchmark, and for each type configuration (V10, V6, V4). Figure 3 shows the variation of the QoR against float mantissa bit-width M f and double mantissa bit-width M d . Thus, Phase 2 provides a set of candidate FPU configurations that produce intermediate QoR levels.

Phase (2) analysis

Figure 3 summarizes the effects of operator-level precision (M f and M d ) on application-level results. For the same application, the impact of bit-width sizing on overall accuracy differs between type configurations. Moreover, the gradient of variation of the QoR reveals some information about the sensitivity of each type configuration to low-level precision. For example, for rectangle, the output accuracy of the V10 type configuration is independent of the float mantissa bit-width (M f ). Indeed, this type configuration considers only one float variable, the impact of which on the QoR is negligible.

Figure 3 could be used by design engineers to guide their selection of the best options for application implementation. Starting from a type configuration (V10 for example), the QoR threshold (and hence the consumed energy) can be selected by varying the underlying operatorlevel precision rather than changing the variable types. This allows a more fine-grained control of the accuracy vs. energy comparison.

The last column of Table 1 shows the selected candidate FPU configurations. The most appropriate one can be selected based on a qualitative comparison of the bitwidths. However, a detailed hardware-level evaluation was performed in Case Study 3 to quantitatively determine which bit-width minimized the overall energy consumption.

Case Study 3: Hardware-level evaluation

In the previous section, we applied AxQEMU to some benchmarks to study the effects of FP approximations on application-level accuracy. In this section, we present a hardware-level evaluation on the selected candidates, which aimed to estimate the HW savings in terms of energy consumption thanks to ARP (Phase 3 ). The methodology was applied to many applications, but for the sake of concision, we only present the arclength study here.

HW evaluation methodology

To perform the HW-level evaluation, an approximateaware hardware FPU was implemented in SystemVerilog. This FPU is based on an open-source parametrized FPU [START_REF] Mach | A transprecision floating-point architecture for energy-efficient 850 embedded computing[END_REF] available online 5 . The HW design is globally parametric so that the four parameters of the FPU configuration 5 https://github.com/pulp-platform/fpnew (E f , M f , E d , M d ) are variable at design time. For the sake of concision, we have omitted the HW implementation details.

The implementation was synthesized as an ASIC, on 540 a 28-nm FD-SOI technology node, in the typical corner (Regular V t , 1.00V, 25C, No Body Biasing) for a 200-MHz frequency target. Synthesis has been performed on Synopsys Design Compiler ® with automatic clock-gating enabled and default effort levels. Post-synthesis simulations were performed using Synopsys VCS ® , and power consumption was estimated by considering both the circuit's static power as well as its dynamic power associated with the studied application using Synopsys PrimeTime ® .

Phase (3) results

550 For a fair comparison, the metric evaluation functions executed in full precision were not considered. The results presented only relate to the computational part i.e., F lo a t m a n ti s s a M f the amount of energy spent on floating-point computations. Other microprocessor and system parameters such as cache memory latency, memory bus contention, branch prediction, external peripherals, etc. are beyond the scope of this paper.

Phase (3) analysis

Figure 4 shows that the type configurations generated by PROMISE only provided 0.9%, 4.8%, and 7.7% savings on FP computation energy compared to the Ref conservative implementation when producing results with 10, 6, and 4 significant digits, respectively. Therefore, even when the accuracy is drastically reduced to just 4 digits, the energy saving does not exceed 7.7%. In contrast, by using ARP, it becomes possible to save 19.4%, 49.5%, and 60.7% energy with the pairs {V10, [START_REF] Lefèvre | Sipe: a Mini-Library for Very Low Precision Computations with Correct Rounding[END_REF][START_REF] Moreau | A Taxonomy of General 770 Purpose Approximate Computing Techniques[END_REF][START_REF] Said | FPU Bit-Width Optimization for Approximate Computing: A Non-Intrusive Approach[END_REF]44)}, {V6, [START_REF] Lefèvre | Sipe: a Mini-Library for Very Low Precision Computations with Correct Rounding[END_REF][START_REF] Hauser | SoftFloat release[END_REF][START_REF] Said | FPU Bit-Width Optimization for Approximate Computing: A Non-Intrusive Approach[END_REF]32)}, {V4, [START_REF] Lefèvre | Sipe: a Mini-Library for Very Low Precision Computations with Correct Rounding[END_REF][START_REF] Hauser | SoftFloat release[END_REF][START_REF] Said | FPU Bit-Width Optimization for Approximate Computing: A Non-Intrusive Approach[END_REF][START_REF] Chatelain | Automatic exploration of reduced floating-point representations in iterative methods[END_REF]}, respectively, without degrading QoR.

For some QoR constraints, VTO tools may not be able to identify an appropriate optimized type configuration. This was the case for all the applications presented in Case Study 2. No optimized type configuration can satisfy an 8digit QoR constraint, and the designer will therefore have to select a configuration producing a higher QoR, such as V10. This effect is due to the coarse granularity of the standard types. However, with ARP, fine-tuning can be performed until a near-threshold, or in other words "good enough", configuration is found satisfying the constraint defined.

When quality constraints are relaxed, more interesting savings are possible, by considering both software (type configurations) and hardware (FPU configurations). For example, if only 4 digits are required, the pairs {V10, [START_REF] Lefèvre | Sipe: a Mini-Library for Very Low Precision Computations with Correct Rounding[END_REF][START_REF] Moreau | A Taxonomy of General 770 Purpose Approximate Computing Techniques[END_REF][START_REF] Said | FPU Bit-Width Optimization for Approximate Computing: A Non-Intrusive Approach[END_REF][START_REF] Chatelain | Automatic exploration of reduced floating-point representations in iterative methods[END_REF]}, {V6, [START_REF] Lefèvre | Sipe: a Mini-Library for Very Low Precision Computations with Correct Rounding[END_REF][START_REF] Sanchez-Stern | Finding Root Causes of Floating Point Error with Herbgrind[END_REF][START_REF] Said | FPU Bit-Width Optimization for Approximate Computing: A Non-Intrusive Approach[END_REF][START_REF] Chatelain | Automatic exploration of reduced floating-point representations in iterative methods[END_REF]}, and {V4, [START_REF] Lefèvre | Sipe: a Mini-Library for Very Low Precision Computations with Correct Rounding[END_REF][START_REF] Hauser | SoftFloat release[END_REF][START_REF] Said | FPU Bit-Width Optimization for Approximate Computing: A Non-Intrusive Approach[END_REF][START_REF] Chatelain | Automatic exploration of reduced floating-point representations in iterative methods[END_REF]} are all good candidates. These combinations provide equivalent energy savings: 60.7%. Other criteria, such as the circuit area or the memory footprint, can be considered to guide the choice of final FPU configuration. According to our synthesis results, the best compromise for this application would be the [START_REF] Lefèvre | Sipe: a Mini-Library for Very Low Precision Computations with Correct Rounding[END_REF][START_REF] Moreau | A Taxonomy of General 770 Purpose Approximate Computing Techniques[END_REF][START_REF] Said | FPU Bit-Width Optimization for Approximate Computing: A Non-Intrusive Approach[END_REF][START_REF] Chatelain | Automatic exploration of reduced floating-point representations in iterative methods[END_REF] configuration, since it has the lowest circuit area overhead, which is estimated to be 27% higher than the standard RV64FD architecture.

Related Works

Our methodology takes advantage of ARP to improve energy efficiency. In this section we present some related 610 works from the literature.

FP Variable type optimization. Tools such as [START_REF] Graillat | Auto-tuning for floating-point precision with Discrete Stochastic Arithmetic[END_REF][START_REF] Lam | Automatically adapting programs for mixed-precision floating-point computation[END_REF][START_REF] Rubio-González | Precimonious: Tuning assistant for floating-point precision[END_REF] are designed to find, given an {application, input dataset, QoR constraint}, a variable type configuration of the original application that minimizes the number of high-precision variables and maximizes the number of low-precision variables. For some of these tools [START_REF] Lam | Automatically adapting programs for mixed-precision floating-point computation[END_REF][START_REF] Rubio-González | Precimonious: Tuning assistant for floating-point precision[END_REF], the objective is to optimize speed, whereas for others [START_REF] Graillat | Auto-tuning for floating-point precision with Discrete Stochastic Arithmetic[END_REF], the goal is to maximize the number of single-precision variables. As a consequence, results produced by [START_REF] Graillat | Auto-tuning for floating-point precision with Discrete Stochastic Arithmetic[END_REF] are more reproducible.

620

All these tools are based on delta-debugging heuristics [START_REF] Zeller | Simplifying and isolating failureinducing input[END_REF], which constitutes a scalability bottleneck when dealing with programs with numerous FP variables. Although [START_REF] Rubio-González | Precimonious: Tuning assistant for floating-point precision[END_REF] is relatively more scalable, the search space still depends on the number of variables in the application, and the transformations are specific to the x86 architecture. As a result, [START_REF] Lam | Automatically adapting programs for mixed-precision floating-point computation[END_REF][START_REF] Rubio-González | Precimonious: Tuning assistant for floating-point precision[END_REF] are architecture-and compiler-dependant.

Non-standard/arbitrary format support. The previous tools only support standard IEEE 754 [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF] formats, although [START_REF] Rubio-González | Precimonious: Tuning assistant for floating-point precision[END_REF] also supports extended precision formats (i.e. In-630 tel's 80-bit format implemented as long double in C). As explained in Section 2.2, our approach is complementary with these and other works [START_REF] Panchekha | Automatically improving accuracy for floating point expressions[END_REF][START_REF] Sanchez-Stern | Finding Root Causes of Floating Point Error with Herbgrind[END_REF][START_REF] Benz | A Dynamic Program Analysis 820 to find Floating-Point Accuracy Problems 10[END_REF], since it provides further computation power/energy/execution time optimization thanks to ARP and it can reuse their results to determine a first coarse-grained memory footprint and computation optimization.

Simulation of FP approximation impact. [START_REF] Ho | Efficient float-840 ing point precision tuning for approximate computing[END_REF] proposes an arbitrary precision impact simulation methodology and introduces an automatic source code transformation tool. 640 However, since it exploits the MPFR C library [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with 790 correct rounding[END_REF], it only allows variation in the mantissa bit-width. Other libraries such as FlexFloat [START_REF] Tagliavini | FlexFloat: A Software Library for Transprecision Computing[END_REF] (written in C) and FloatX [START_REF] Flegar | FloatX: A C++ Library for Customized Floating-Point Arithmetic[END_REF] (its C++ successor) were developed to make it possible to express variable precision behavior in software, to simulate the effects of arbitrary FP bit-width reduction on the QoR provided by applications. Using these libraries, developers are expected to modify all or some of their FP variable types using arbitrary-precision types as substitutes for the standard ones. The code transformation effort required can be very significant, especially for low-level specialized computational kernels in C/Assembly for example. Moreover, some of the classic languages do not support advanced features such as operator overloading. Adapting existing system stacks (hardware, compilers, standard math, scientific libraries, operating systems, etc.) to FP resizing consequently becomes a tremendous development effort, especially at the production stage.

Non-intrusive impact simulation. One of the most closely related works to the one described here is VPREC [START_REF] Chatelain | Automatic exploration of reduced floating-point representations in iterative methods[END_REF], a back-end system that enables non-intrusive variable precision simulation in the Verificarlo [START_REF] Denis | Verificarlo: Checking floating point accuracy through monte carlo arithmetic[END_REF] software toolchain. This tool was introduced to simulate in-time variable precision, specifically for iterative algorithms which benefit from adaptive gradually-increasing precision rather than fixed precision. Although the idea is very similar to ours, there are some key differences: 1) Verificarlo requires a specific instrumentation toolchain, making recompiling mandatory, and as a consequence the instrumented application is not "production-ready", 2) Instrumentation is performed on Valgrind, which means that the application was studied in a high-level operating system context. In contrast, here, we support virtually all contexts (including bare-metal). In addition, and most importantly 3) VPREC is only designed for simulation, and its conclusions cannot be readily transferred to hardware due to the lack of the necessary hardware support.

In contrast with all these approaches, ours operates on the final executable binary after all user-specific and compiler optimizations and can benefit from existing VTO tools. It can also be directly implemented in hardware to take advantage of the fast software design space exploration. To the best of our knowledge, our methodology is the first non-intrusive hardware-oriented approach to simulate the impact of FP approximation.

Discussion

Challenges. Approximating individual assembly instructions comes with a few drawbacks. Since approximation is applied locally, error propagation is not controlled globally because the algorithmic structure of the original program is practically lost. However, our approach considers the final production-ready executable binary as well as compiler optimizations. This feature is unique in that it reflects the expected production-stage behavior.

Limitations. The approach described in this paper is agnostic to the language (C, Fortran, etc.), to the context (bare-metal, OS/RTOS etc.), and to the FP analysis tools used upstream. It is therefore very powerful when compared to existing strategies. However, for a given FP program using complex functions, the results will depend on the underlying standard C library implementation. Indeed, the QoR variation behavior vis-a-vis the precision reduction might differ slightly between two implementations. For instance, transcendental functions are implemented in different ways across distinct C standard library implementations (e.g. Newlib, Glibc, MUSL, etc.). Their internal use of elementary FP operations (e.g., addition, multiplication) will differ, which explains the non-identical numerical results. The proposed methodology is thus C library-dependant when complex functions are invoked.

710

Another limitation of this approach is its datadependency. As for all other existing techniques, our strategy provides no guarantees that the application will produce an equivalent QoR for all possible input datasets. Hence, designers should carefully select representative data inputs to cover all intended application behaviors.

Future work. Future studies will focus on introducing variable precision in time in AxQEMU and propose an appropriate hardware architecture that benefits from operatorlevel precision tuning. A more complete hardware-level 720 evaluation taking system parameters into consideration will also be established.

Conclusion

This work introduces a hardware-friendly and code nonintrusive approach to optimize FPU bit-width sizing for a given embedded application, dataset, and quality of result level (QoR). A software implementation (AxQEMU) based on the well-known QEMU dynamic binary translator made it possible to perform a rapid evaluation of the impact of arbitrary reduced precision (ARP) on an 730 application's QoR, without modifying the software stack. Compared to other approaches, our methodology allows fine-grained bit-width tuning to obtain the best QoR vs. energy consumption trade-off. Our approach can be used either as a standalone process or in conjunction with variable type optimization tools for fine-grained tuning. The approach was tested on several benchmarks from the literature through case studies. Several hardware FPUs were designed on an ASIC 28-nm FD-SOI technology node to demonstrate the proposed methodology's efficiency. Ex-740 perimental results exhibited computational energy savings of between 19.4% and 60.7% when using ARP, compared to 7.7% when using the best current alternative tools.

Figure 1 :

 1 Figure 1: AxQEMU overview

#

  define PRECISE __attribute__((__section__(".precise"))) \ __attribute__((noinline)) double PRECISE compute_QoR(float result, double ref){ // Function body here } Listing 1: Example of a C macro for function annotation At run-time, the address of each assembly instruction instance is fetched from the Program Counter register (PC).

  350

4

 4

x 10 3 Figure 2 :

 32 Figure 2: QoR variation as a function of the exponent E d and mantissa M d bit-widths.

Figure 4 :

 4 Figure 4: Energy vs. QoR trade-offs

Figure 3 :

 3 Figure 3: QoR results for arclength, rectangle, and squareroot.

Table 1 :

 1 Benchmark summary

	Type	Configurations	Application	# Doubles	# Floats	QoR	Threshold	Candidate FPU Configurations
						10	(8, 3, 11, 44)
			arclength	8	1	8 6	(8, 3, 11, 36) (8, 3, 11, 32)
						4	(8, 3, 11, 28)
						10	(8, 3, 11, 36)
	V10	rectangle	4	3	8-6	(8, 3, 11, 28)
						4	(8, 3, 11, 16)
						10	(8, 3, 11, 32)
			squareroot	6	2	8 6	(8, 3, 11, 28) (8, 3, 11, 20)
						4	(8, 3, 11, 8)
			arclength	7	2	6 4	(8, 22, 11, 32) (8, 13, 11, 28)
	V6	rectangle	3	4	8-6 4	(8, 3, 11, 28) (8, 3, 11, 16)
			squareroot	0	8	6 4	(8, 15, 11, 4) (8, 7, 11, 4)
	V4	arclength rectangle	2 0	7 7	4 4	(8, 22, 11, 28) (8, 16, 11, 4), (8, 13, 11, 4)

The "d" in helper fadd d refers to double-precision, whereas the "s" in helper fadd s refers to single-precision.

http://promise.lip6.fr

RISC-V architecture with D and F extensions.

Acknowledgement

The authors would like to thank Stefan Mach (ETH Zurich) and Rodrigo Iga Jodue (TIMA Laboratory) for their valuable technical advice on the HW implementation.

CRediT authorship contribution statement N. Ait Said: Conceptualization, Methodology, Formal analysis, Software, Validation, Investigation, Data cu-750 ration, Visualization, Writing -Original Draft, Writing -

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could 760 have appeared to influence the work reported in this paper.


Paolo Massioni
email: paolo.massioni@insa-lyon.fr

Laurent Bako
email: laurent.bako@ec-lyon.fr

Gérard Scorletti
email: gerard.scorletti@ec-lyon.fr

Stability of uncertain piecewise-affine systems with parametric dependence

Keywords: Hybrid systems, Piecewise-affine systems, Linear matrix inequalities

INTRODUCTION

Finding a generic and systematic method for analyzing the properties of generic nonlinear system is a challenging problem in control theory. However, the seminal works in [START_REF] Johansson | Computation of piecewise quadratic Lyapunov functions for hybrid systems[END_REF]; [START_REF] Hassibi | Quadratic stabilization and control of piecewise-linear systems[END_REF] have shown that if we consider approximating the nonlinearities using a piecewise-affine (PWA) representation, then a systematic numerical approach can be developed. This approach is still an active topic of research in automatic control, as demonstrated by several recent works like [START_REF] Ameur | Lyapunov stability analysis of switching controllers in presence of sliding modes and parametric uncertainties with application to pneumatic systems[END_REF], [START_REF] Waitman | Incremental stability of Lur'e systems through piecewise-affine approximations[END_REF], [START_REF] Iervolino | Lyapunov stability for piecewise affine systems via conecopositivity[END_REF]. A piecewise-affine representation consists of a polyhedral partition of the state-space of the system into region with disjoint interiors, each of which is associated with an affine time-invariant dynamical model. Whereas there are systems that can naturally be described by PWA models (e.g. linear systems interconnected by static nonlinearities such as saturations, dead zones, etc.), in most of the other situations, PWA systems can be regarded merely as convenient approximations in view of easier control design or stability analysis. From this latter perspective, one must account for the uncertainty associated to the approximation. This can be done through, for example, the construction of piecewise differential inclusions that embed the dynamics of the nonlinear system as in [START_REF] Johansson | Piecewise linear control systems[END_REF] or more simply, by allowing the PWA model parameters to be parameter dependent.

We consider in this paper the problem of assessing the asymptotic stability for a class of PWA systems where the partition and the affine subsystem dynamical equations depend on a scalar time-varying parameter. The parameter in question, although unknown, is assumed to vary smoothly in a bounded interval. A comparable class of model has been defined in [START_REF] Lebel | Piecewise-affine parameter-varying control of wheeled mobile robots[END_REF] under the name of piecewise-affine parameter-varying model (PWA-PV), where the cells were limited to time-invariant. The class of systems that we propose in this work typically arises when approximating a nonlinear time-varying system with a PWA model. Also, assuming a variable partition and letting the matrices of the affine subsystems be parameter dependent can account for the uncertainty associated with the approximation of nonlinear systems. These situations are illustrated through some motivating examples in Section 3.

The proposed analysis method consists in searching for a quadratic Lyapunov function on the state-space with a polynomial dependency on the scheduling parameter. We show that this can be reduced to a feasibility problem in the decision variables under matrix sum of squares (MSOS) constraints. The latter can be made convex and solved efficiently using available numerical tools.

The paper is organised as follows. Section 2 introduces all the preliminary notions, including the notation, piecewiseaffine system, and sum of squares (SOS). Section 3 contains a discussion concerning the motivation for the work in this paper, whereas Section 4 explains the proposed approach, eventually leading to the main result of this paper. Section 5 shows the application of the main result to two cases, one derived from the motivating example and the second one derived from a famous example in the literature, and finally Section 6 draws the conclusions.

PRELIMINARIES

Notation

Let R denote the set of real numbers, R + the set of positive real numbers including zero, Z the set of the integers and N the set of positive integers. For a symmetric matrix A, the expression A 0 indicates that A is positive semidefinite. We also use the notation M (:) to indicate the entries of a matrix M taken one at a time, and for imposing component-wise constraints; for example M (:) 0 indicates that all the entries of M are positive or equal to zero. Definition 1. (Uncertain piecewise-affine system). In the context of this work, we call an autonomous dynamical system "uncertain piecewise-affine with parametric dependence" when it has a representation given by ẋ

(t) = A i (ρ)x(t) + a i (ρ) for x(t) ∈ X i (ρ)
(1) where ρ ∈ R := [ρ min , ρ max ] is a time-varying scalar parameter, continuous with respect to time and with bounded rate of change. The regions (or cells) X i (ρ), for i ∈ I := {1, . . . , L}, are closed convex polyhedral sets defined by

X i = {x ∈ R n | E i (ρ)x + e i (ρ)
0} with non-empty and pairwise disjoint interiors ∀ρ ∈ R, such that i∈I X i (ρ) = R n . Then, {X i (ρ)} i∈I constitutes a finite partition of R n . From the geometry of X i (ρ), the intersection X i (ρ) ∩ X j (ρ) between two different regions is always contained in a hyperplane. Let I 0 ⊆ I be the index set for cells that contain the origin for at least one value of ρ ∈ R, and let I 1 = I \ I 0 . We assume that a i (ρ) = 0 for i ∈ I 0 and for all ρ ∈ R (this is a necessary and sufficient condition for equilibrium at the origin).

The representation above can be simplified by extending the state vector x(t) as x(t) = [x(t) ⊤ , 1] ⊤ , which allows reformulating it as

ẋ(t) = A i (ρ)x(t) for x(t) ∈ X i (ρ) (2) and X i (ρ) = {x ∈ R n | E i (ρ)x 0}
(3) with opportune definitions of A i (ρ) and E i (ρ):

A i (ρ) = A i (ρ) a i (ρ) 0 0 , E i (ρ) = [ E i (ρ) e i (ρ) ] , (4)
where a i (ρ) = 0 and e i (ρ) = 0 for i ∈ I 0 (in the same way as in [START_REF] Johansson | Computation of piecewise quadratic Lyapunov functions for hybrid systems[END_REF]). In the remainder of the paper, we will assume that A i (ρ) are rational functions of ρ, i.e. their entries can be expressed as ratios of polynomials in ρ, and that E i (ρ) are polynomial functions of ρ. If A i (ρ) is rational, then we add the additional hypothesis that A i (ρ) is finite and defined for all ρ ∈ R (i.e., there is no denominator that is equal to zero for values of ρ ∈ R); notice also that there is no loss of generality in assuming E i (ρ) polynomial instead of rational: denominators in E i (ρ) that are always defined never assume a value of 0 for ρ ∈ R, so they always have the same sign and they can be simplified without changing the definitions of the cells. We also assume the absence of sliding modes, as it is the case for example if the right-hand side of (1) is continous across boundaries. In any case, sliding modes can be dealt with as explained in Section IIV of [START_REF] Johansson | Computation of piecewise quadratic Lyapunov functions for hybrid systems[END_REF].

Sum of squares

Definition 2. (Generic SOS problems, Parrilo (2003)). Let p(t) be a polynomial of degree up to 2d ∈ N in the variable t ∈ R. We call a sum of squares problem (SOS) the problem of finding whether there exists a finite number l of polynomials π i (t) such that

p(t) = l i=1 π i (t) 2 . ( 5
)
If such a decomposition into a sum of squares of polynomials exists, we say that p(t) is sum of squares (SOS), which implies that p(t) 0 for all t.

Let P(t) be a symmetric matrix of polynomials of degree up to 2d ∈ N in the variable t ∈ R. We call a matrix sum of squares problem (MSOS) the problem of finding whether there exists a finite number l of symmetric matrices of polynomials Π i (t) such that

P(t) = l i=1 Π i (t) ⊤ Π i (t). ( 6
)
If such a decomposition exists, we say that P (t) is matrix sum of squares (MSOS), which implies that P (t) 0 for all t.

The definition above is restricted to univariate polynomials (t ∈ R), as this is the case we consider in this paper. SOS and MSOS problems are convex problems. A derived class of problems is the one of feasibility problems under SOS or MSOS constraints, defined here. Definition 3. (SOS constraints feasibility, [START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF]). Let p i (t, θ), i = 1, . . . , q, be a set of q polynomials of degree up to 2d ∈ N in the variable t ∈ R, where θ ∈ R ρ is a vector of parameters or unknowns, with p i (t, θ) affine with respect to the entries of θ. A feasibility problem under SOS constraints consists in finding, if it exists, a value of θ = θ * for which p i (t, θ * ) is SOS, for i = 1, ..., q.

(7) If such a θ * exists, then the problem is feasible; otherwise it is unfeasible. Definition 4. (MSOS constraints feasibility, Chesi (2010)). Let P i (t, θ), i = 1, . . . , q, be a set of q symmetric matrices of polynomials of degree up to 2d ∈ N in the variable t ∈ R, where θ ∈ R ρ is a vector of parameters or unknowns, and the matrices P i (t, θ) are affine with respect to the entries of θ. A feasibility problem under SOS constraints consists in finding, if it exists, a value of θ = θ * for which

P i (t, θ * ) is MSOS, for i = 1, ..., q.
(8) If such a θ * exists, then the problem is feasible; otherwise it is unfeasible.

Feasibility problems under SOS and MSOS problems are convex optimisation problems, and they can be reformulated as linear matrix inequality (LMI) feasibility problems; this can either be done explicitly, or relying on automated procedures, like the one available in the Yalmip toolbox [START_REF] Löfberg | Pre-and post-processing sum-ofsquares programs in practice[END_REF] under Matlab.

The general S-procedure

By the name S-procedure, one commonly identifies a lemma which allows restricting a certain class of inequal-ities to a certain subset (see [START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF]). In this article, we rely on a very simple but general expression, whose proof is obvious, which will be specialised according to the cases. Lemma 1. (Generalised S-procedure). Let F (x), G(x) be (symmetric matrix) functions of the (vector) variable x. Let g(x) be a scalar function of x. The following implications hold.

F (x) 0 for G(x) 0 ⇐ F (x) -G(x)λ 0 for λ 0, ∀x (9) F (x) 0 for g(x) 0 ⇐ F (x) -g(x)Λ 0 for Λ 0, ∀x (10)
The terms Λ and λ are called multipliers, which can be chosen at one's convenience, i.e. they are decision variables subject to the positivity constraints above. When used in the context of polynomial problems as in the context of this paper, the lemma above is a direct consequence of a lemma known as Positivstellensatz (see Chesi ( 2010)) or p-satz, of which several versions exist in the literature. The multipliers can in this case have a polynomial dependence on x, which allows satisfying the positivity constraints by means of either SOS or MSOS constraints.

Lyapunov stability

Lemma 2. (Adapted from [START_REF] Kalman | Control system analysis and design via the "second method" of Lyapunov: I-continuous-time systems[END_REF]). Consider a parameter varying dynamical system ẋ = g(ρ, x), with g continuous depending on a continuously differentiable function ρ of time. If there exists a continuous and differentiable function V (x, ρ) = x ⊤ P (ρ, x)x such that α||x|| 2 V (x, ρ) β||x|| 2 (11) and V (x, ρ) -ε||x|| 2 (12) with α > 0, β > 0, ε > 0, then the system is asymptotically exponentially stable for any initial condition.

MOTIVATION

The motivation of this work, as stated in the introduction, is based on the need of an accurate description of nonlinearities that have some parametric dependence, for which in turn the parameter might be time-varying. Several examples of these kind of systems can be found in practical applications, we can for example mention varistors (see [START_REF] Eda | Zinc oxide varistors[END_REF]) and thermistors (see [START_REF] Scarr | Thermistors, their theory, manufacture and application[END_REF]); these electrical components have nonlinear current-versus-voltage curves depending on parameters (the exponent α for varistors and the temperature for the thermistors). Besides these, the other main motivation for a more accurate description of nonlinearities comes from friction and specifically dry friction (see for example [START_REF] Karnopp | Computer simulation of stick-slip friction in mechanical dynamic systems[END_REF]; [START_REF] Persson | Rubber friction: role of the flash temperature[END_REF]). Dry friction is a nonlinear phenomenon that is very difficult to model, and which depends on many factors and parameters. The dry friction coefficient µ, i.e. the ratio between friction force and contact force, depends not only on the relative sliding velocity but also on several parameters such as the temperature, the contact force itself, and even the acceleration of the motion. For example, from [START_REF] Persson | Rubber friction: role of the flash temperature[END_REF] one can derive a model of µ with a parametric dependence as in Fig. 1.

A first attempt to deal with time-varying and parametervarying friction coefficient has been done in [START_REF] Ameur | Lyapunov stability analysis of switching controllers in presence of sliding modes and parametric uncertainties with application to pneumatic systems[END_REF], with very gross and inaccurate differential inclusions. On the other hand, a curve like the one in Fig. 1 can be accurately described with a piecewise-affine model, if one only allows the boundary themselves to change as a function of the parameter. This model of friction, presented here as a motivation for the work, will be featured again in the first numerical example of this paper, in Section 5.1.

MAIN RESULT

In this section, we state the main result of the paper which is a stability theorem for the uncertain PWA system defined in (1). To assess the global asymptotic stability of that system we will be searching for a piecewise quadratic Lyapunov function V : R n × R → R + , which depends polynomially on the the same scalar parameter ρ as the dynamics of the system (1). More precisely, inspired by the work reported in Johansson and Rantzer (1998) we set V (x, ρ) to be of the form

V (x, ρ) = x⊤ P i (ρ)x if x ∈ X i (ρ)
for some i ∈ I (14) where x is defined from x as before and

P i (ρ) = P i (ρ) p i (ρ) p i (ρ) ⊤ r i (ρ) , (15)
of any chosen degree in ρ, with the constraints that p i (ρ) = 0 and r i (ρ) = 0 for i ∈ I 0 . Observe that V (x, ρ) is defined here on a partition of the state-space which coincides with that of the system to be analyzed, but in a more general setting one could consider a piecewise Lyapunov function candidate defined on a partition of the state-space which is completely different from that of the system to be analyzed.

Continuity of V : Similarly as in Johansson and Rantzer we parameterise the matrix functions P i (ρ) in the

for i ∈ I 0      P i (ρ) -E i (ρ) ⊤ W i (ρ)E i (ρ) -G i (ρ)η(ρ) -Iε is MSOS ν i (ρ) -A i (ρ) ⊤ P i (ρ) -P i (ρ)A i (ρ) -ρ d ∂P i (ρ) ∂ρ -Iε -E i (ρ) ⊤ U i (ρ)E i (ρ) -H i (ρ)η(ρ) is MSOS ν i (ρ) -A i (ρ) ⊤ P i (ρ) -P i (ρ)A i (ρ) + ρ d ∂P i (ρ) ∂ρ -Iε -E i (ρ) ⊤ U ′ i (ρ)E i (ρ) -H ′ i (ρ)η(ρ) is MSOS W i (ρ)(:), U i (ρ)(:), U ′ i (ρ)(:) are SOS, G i (ρ), H i (ρ), H ′ i (ρ) are MSOS for i ∈ I 1          P i (ρ) -E i (ρ) ⊤ W i (ρ)E i (ρ) -G i (ρ)η(ρ) -Iε is MSOS ν i (ρ) -A i (ρ) ⊤ P i (ρ) -P i (ρ)A i (ρ) -ρ d ∂P i (ρ) ∂ρ -Iε -E i (ρ) ⊤ U i (ρ)E i (ρ) -H i (ρ)η(ρ) is MSOS ν i (ρ) -A i (ρ) ⊤ P i (ρ) -P i (ρ)A i (ρ) -ρ d ∂P i (ρ) ∂ρ -Iε -E i (ρ) ⊤ U ′ i (ρ)E i (ρ) -H ′ i (ρ)η(ρ) is MSOS W i (ρ)(:), U i (ρ)(:), U ′ i (ρ)(:) are SOS, G i (ρ), H i (ρ), H ′ i (ρ) are MSOS (13) form P i (ρ) = F i (ρ) ⊤ T (ρ)F i (ρ) (16
) where T (ρ) is a parameterised polynomial matrix function of ρ of a chosen arbitrary degree, and the F i (ρ) are known matrix polynomial functions:

F i (ρ) = [ F i (ρ) f i (ρ) ] , (17)
of appropriate dimensions, satisfying F i (ρ)x = F j (ρ)x (18) for all pairs (i, j) ∈ I 2 of cells with non trivial intersection and all x such that x ∈ X i (ρ) ∩ X j (ρ) (with f i (ρ) = 0 for i ∈ I 0 , see again [START_REF] Johansson | Computation of piecewise quadratic Lyapunov functions for hybrid systems[END_REF]). With the parameterizations in ( 16) and the conditions (18), the Lyapunov function candidate V (x, ρ) expressed in ( 14) is automatically ensured to be continuous on R n when regarded as a function of x.

The Lyapunov functions candidate will need to satisfy three properties, of which a brief introductory account is given here, before stating the main theorem.

Positive definiteness of V (x, ρ): We must require V (x, ρ) to be positive-definite. For that it suffices that each piece of V is positive in the corresponding cell, i.e., x⊤ P i (ρ)x > 0 for all ρ ∈ R and all x ∈ X i (ρ), x = 0, i ∈ I. These conditions can be ensured by making use of the device of the S-procedure recalled in Lemma 1. For this purpose note that a cell X i (ρ) is contained in any set represented by the inequality x⊤ E i (ρ) ⊤ W i (ρ)E i (ρ)x 0 provided that W i (ρ) is a nonnegative matrix function, in the sense that all of its entries are nonnegative. Proceeding likewise the range R of ρ can be embedded in the set defined by the matrix inequality G i (ρ)η(ρ) 0 for some positive semidefinite matrix function G i (ρ) and with

η(ρ) = 1 4 (ρ max -ρ min ) 2 -ρ - ρ min + ρ max 2 2 (19) (η(ρ) is nonnegative when ρ ∈ R).
Wrapping it up, the existence for all i ∈ I of W i (ρ), G i (ρ) and ε > 0 such that

P i (ρ) -E i (ρ) ⊤ W i (ρ)E i (ρ) -G i (ρ)η(ρ) -εI is MSOS constitutes a sufficient condition for V (x, ρ) to be positive- definite.
Here

I = I n 0 0 0 , ( 20
)
with I n the identity matrix of order n (also written as I for a lighter notation).

Negativity of V (x, ρ): The negativity of the time derivative V (x, ρ) along the system trajectories can also be enforced per cell using the same technique of S-procedure. The details are in the proof on the theorem that we are now ready to state.

Theorem 1. Consider an uncertain PWA system according to Definition 1, with matrices defined according to ( 4), ( 16), ( 17), and assume | ρ| ρ d , with ρ d a given positive scalar.

The system is asymptotically exponentially stable for any initial conditions if for some ε > 0, the MSOS constraints in (13) are feasible (at the top of the page), for i = 1, ..., N . The decisions variables are the matrix coefficients of T (ρ) = T 0 +ρT 1 +ρ 2 T 2 ... and the polynomial multipliers

W i (ρ), G i (ρ), U i (ρ), U ′ i (ρ), H i (ρ), H ′ i (ρ) or W i (ρ), G i (ρ), U i (ρ), U ′ i (ρ), H i (ρ), H ′ i (ρ), of fixed arbitrary degree. The function η(ρ) is defined in (19), and ν i (ρ), with ν i (ρ) > 0 for ρ ∈ R, is a polynomial for which ν i (ρ)A i (ρ) (if i ∈ I 0 ) or ν i (ρ)A i (ρ) (if i ∈ I i ) are pure polynomials (in case there are denominators in A i (ρ) or A i (ρ), it simplifies them).
Proof. We show that V (x, ρ) is a valid Lyapunov function in each cell. By definition, V (x, ρ) is continuous with respect to x and ρ. Subsequently, we show that the conditions in (13) imply the conditions in Lemma 2; the proof is shown in details only for cells with i ∈ I 0 , but the extension to i ∈ I 1 follows exactly the same scheme. This proof is inspired by the approach to parameterdependent systems that can be found in [START_REF] Wu | SOS-based solution approach to polynomial LPV system analysis and synthesis problems[END_REF]. Consider a cell X i with i ∈ I 0 , and consider the first condition in (13). If the condition is satisfied, it means that the expression on the left hand side is always positive semidefinite; the terms E i (ρ) ⊤ W i (ρ)E i (ρ) and G i (ρ)η(ρ) are general S-procedure terms with W i (ρ) and G i (ρ) nonnegative multipliers (entry-wise for W i (ρ), positive semidefinite for G i (ρ), as implied by the fourth line of ( 13)); this implies that x ⊤ (P i (ρ) -Iε)x 0 when E i (ρ)x 0 and when η(ρ) 0, i.e. V (x, ρ) satisfies the left hand side inequality in (11) when x ∈ X i and ρ ∈ R. The right hand side of ( 11) is proven by the definition of V (x, ρ) as a quadratic function of x, with ρ bounded. Going now to the second and third inequalities in (13): by the same Sprocedure mechanisms explained above and subsequently simplifying the positive term ν i (ρ), they imply that

x ⊤ -A i (ρ) ⊤ P i (ρ) -P i (ρ)A i (ρ) -ρ d ∂P i (ρ) ∂ρ -Iε x 0 and x ⊤ -A i (ρ) ⊤ P i (ρ) -P i (ρ)A i (ρ) + ρ d ∂P i (ρ) ∂ρ -Iε x 0
for x ∈ X i and ρ ∈ R. Notice that any convex combination of the two inequalities above will still hold, so the two imply that

x ⊤ -A i (ρ) ⊤ P i (ρ) -P i (ρ)A i (ρ) - ρ ∂P i (ρ) ∂ρ -Iε x 0 for any ρ ∈ [-ρ d , ρ d ].
This last inequality implies V (x, ρ) -ε||x|| 2 , i.e. it implies (12), concluding the proof. Remark 1. (On moving boundaries). It might surprise some readers that no specific condition is set for the case when the border is moving, compared to the main result in [START_REF] Johansson | Computation of piecewise quadratic Lyapunov functions for hybrid systems[END_REF]. In fact, there is no difference whether the state crosses a boundary because it is moving itself, or because the boundary is moving: consider that crossing a boundary is an action that requires some time, whereas the conditions in Theorem 1 assure that instantly, at each time, the Lyapunov function is decreasing; this is true also across the boundary (see Figure 2, point x b ), regardless of how the state arrived at the boundary and regardless of how the boundary moves. The time derivative of the Lyapunov function includes a term due to the derivative of the state, and a term due to the time-change of the expression of the Lyapunov function itself ( Ṗi ); this last term includes also the effect of the boundaries, if any. The continuity of the Lyapunov function makes sure that there are no jumps of its value over time. Remark 2. (On necessity). One might wonder if, at least in some cases, the conditions of Theorem 1 can be necessary as well; in fact we are in the case of univariate polynomials, for which several necessary and sufficient results exist, like for example the powerful Markov-Lukacs's theorem [START_REF] Genin | Convex optimization over positive polynomials and filter design[END_REF]. Also in the matrix of square case the necessity would be obtained for a given (high) degree of the multipliers G i (ρ), [START_REF] Briat | Convex conditions for robust stabilization of uncertain switched systems with guaranteed minimum and mode-dependent dwell-time[END_REF][START_REF] Scherer | Matrix sum-of-squares relaxations for robust semi-definite programs[END_REF], but this degree cannot be known in advance and increasing it too much is often impossible once one implements the problem on a computer.

EXAMPLES

Motivating example (friction)

Let us consider a feedback positioning system with both viscous friction and dry friction, of equation:

m v = -ky -cv -N µ(v, ρ) (21
) where y is the position, v = ẏ is the velocity, k is the force feedback gain, c is a viscous friction coefficient, µ(v, ρ) is the static friction coefficient, N is the contact force and m is the mass. For µ(v, ρ), we consider a model as the one discussed in Section 3 (Fig. 1), which yields an uncertain PWA system depending on the time-varying parameter ρ. Taking the state vector as x = [y, v] ⊤ , the state-space is divided into 5 cells according to the value of the velocity v; the first cell is X

1 (ρ) = {x | -ρ v ρ} (1 ∈ I 0 ), where the dry friction coefficient is µ(v, ρ) = µmax ρ v; the second is X 2 (ρ) = {x | ρ v 2ρ}, where µ(v, ρ) = -µmax-µ∞ ρ v + (2µ max -µ ∞ ); and the third is X 3 (ρ) = {x | v 2ρ} (2, 3 ∈ I 1 ) where µ(v, ρ) = µ ∞ .
The fourth and fifth cells are just the symmetric of the second and third; considering that the dynamical system is invariant to a central symmetry, we do not need to worry about them: if a Lyapunov function exists for the first three cells, then it exists for all the state-space.

The cells are described by matrices

A 1 = 0 1 -k m -c m -µmaxN ρm , E 1 = 0, F 1 =      0 0 0 1 0 1 1 0 0 1      , (22)
A 2 =   0 1 0 -k m c m + (µmax-µ∞)N ρm -(2µmax-µ∞)N m 0 0 0   , E 2 = 0 1 -ρ 0 -1 2ρ 0 0 ρ , (23)
A 3 =   0 1 0 -k m -c m -µ∞N m 0 0 0   , E 3 = 0 1 -ρ 0 2 -4ρ 0 0 ρ , (24)
F 2 =      0 1 -ρ 0 -1 2ρ 0 0 ρ 1 0 0 0 1 0      , F 3 =      0 1 -ρ 0 2 -4ρ 0 0 ρ 1 0 0 0 1 0      . (25)
The numerical procedure in Theorem 1 is coded into Matlab with the help of Yalmip [START_REF] Löfberg | Pre-and post-processing sum-ofsquares programs in practice[END_REF][START_REF] Löfberg | YALMIP: a toolbox for modeling and optimization in MATLAB[END_REF], and Mosek (ApS, 2017) is used as solver. Setting

m = 1, N = 1, µ max = 2, µ ∞ = 1, k = 0.1, c = 0.1, | ρ| 1, ρ ∈ [0.01, 2]
, and ν i (ρ) = ρ, the procedure is successful even for T of zero degree in ρ, which yields a valid parameter-dependent Lyapunov function defined by the following matrices. The stability of the positioning servo-system is then assured, even in presence of high uncertainty in the friction model.

Example from the literature

We consider Example 1 on page 557 of [START_REF] Johansson | Computation of piecewise quadratic Lyapunov functions for hybrid systems[END_REF], and we add a time-varying uncertainty on the borders of the cells. The dynamics (second order) is the same as in the reference:

A 1 = A 3 = -ǫ ω -αω -ǫ , A 2 = A 4 = -ǫ αω -ω -ǫ , (29)
with α = 5, ω = 1 and ǫ = 3. The difference is introduced in the definition of the cells:

E 1 = -E 3 = -1 1 + ρ -1 -1 -ρ , (30)
E 2 = -E 4 = -1 1 + ρ 1 1 + ρ , (31)
and

F ⊤ i = [E ⊤ i , I].
We let ρ ∈ [-0.5, 0.5] and code the numerical procedure of Theorem 1 with a seconddegree matrix T (ρ) (whereas ν i (ρ) = 1). The procedure is successful in finding a Lyapunov function for small values of ρ d , i.e. it works roughly up to ρ d ≈ 0.25; for faster variations of the parameter, the system is not guaranteed to be stable, which was to be expected as a coordinated switching between the two different dynamics can lead the state to blow-up (i.e., there is no common Lyapunov function for both A 1 = A 3 and A 2 = A 4 ).

CONCLUSIONS

This work has demonstrated the possibility of systematically finding piecewise-quadratic function for piecewise-affine systems with parametric dependence, where not only the dynamics can change with respect to the parameter, but also the state partition itself. This work can be considered as preliminary and many possible extensions can be imagined, starting for example with piecewise polynomial Lyapunov function as in [START_REF] Ameur | Lyapunov stability analysis of switching controllers in presence of sliding modes and parametric uncertainties with application to pneumatic systems[END_REF] or [START_REF] Samadi | A unified dissipativity approach for stability analysis of piecewise smooth systems[END_REF]. In addition to this, although we have only considered here a scalar parameter dependency, the method can be generalized without many technical difficulties to dependency with respect to more than one parameter. Further research might also investigate the extension of the analysis from simple stability to several different indices of performance.

Fig. 1 .

 1 Fig. 1. Example of a curve of friction coefficient µ as a function of the sliding velocity v; the shape features a parametric dependence on ρ (simplified and adapted from Persson (2006)).

Fig. 2 .

 2 Fig. 2. Moving boundaries between two cells X 1 and X 2 .


S Gerbaud
email: lgerbaud@chu-clermontferrand.fr

Djezzar

Anne Batisse
email: anne.batisse@aphp.fr

Stéphanie Leger
email: stephanie.leger@math.univ-bpclermont.fr

Eric Vicaut
email: eric.vicaut@aphp.fr

Laurent Gerbaud

Samira Djezzar
email: samira.djezzar@aphp.fr

of psychoactive Substances among Youth Students (COSYS): a cross-sectional study in France Absrtract

Keywords:

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Drug use among young adults is a major concern and epidemiological monitoring is essential for public health policy makers. The beginning of high school is marked by a change in lifestyle and is often a time when individuals are first exposed to psychoactive substances. Particular attention must be focused on the student population, as their drug use is determined by their consumption during adolescence and also provides an insight to future use (i.e. transposition of consumption from the academic to the professional environment).

Increasingly, students must adapt themselves to new requirements and, for some, this is also in a very competitive environment. The use of psychoactive substances results in a form of neuroenhancement (NE). NE is the non-medical use (NMU) of psychoactive substances resulting in a subjective enhancement in psychological functioning and experience [START_REF] Brand | Drugs As Instruments: Describing and Testing a Behavioral Approach to the Study of Neuroenhancement[END_REF] . In COSYS (the French Study on COgnitive enhancement and consumption of psychoactive Substances among Youth Students), the behaviour of students is explored by analysing relationships between drug use (e.g. over-the-counter products, prescription drugs or illicit drugs) and motivations (e.g. improved cognitive performance, counteracting stress and fatigue, improved social interaction).

Methods

COSYS is a national cross-sectional survey of French students conducted by the Parisian

AddictoVigilance Centre. Further information on the research methodology can be found online. 9 Questionnaires were mailed to students between January and June 2017. Students were asked about their use of prescription and illicit drugs, in addition to their motivations for use. Participants reported their gender, age, physical activity habits, school, and disability status and level (from dyslexia to motor disorders). Firstly, descriptive statistics for all variables were examined. Subsequently, sociodemographic variables were included in a multivariable logistic regression model to determine conditional associations between each variable and different psychoactive drug use.

Results

This cross-sectional study included 46,203 respondents throughout France, mostly in universities (>60%) and 'grandes écoles' (20%). Respondents were mostly women (63.4%), with an average age of 21.4 years. Alcohol, tobacco and cannabis were the most frequently used substances, but the use of psychoactive medicines was also important: found in third position among women (prevalence = 22.48%) after alcohol/tobacco and in fourth position among men (prevalence = 15.14%) after alcohol/tobacco/cannabis (see table 1). In contrast, the prevalence of current use for 'illicit substances' is still much higher among men (28.12%) than women (17.58%) [p < 0.0001].

For all substances (with the exception of alcohol), physical activity appears to be a protective factor. Within psychoactive medicines, the most represented subgroup is opiates (15%) [e.g. paracetamolcodeine > codeine > morphine > dextromethorphan], followed by benzodiazepines (3.42%).

Two subgroups were identified to have a clear male predominance (psychostimulants and nicotine replacement therapy).

NMU of prescription drugs was reported, and these were obtained from the family medicine cabinet, a friend, a dealer or the Internet. NMU for stimulants and benzodiazepines was identified (60.8% modafinil, 32.9% methylphenidate, 27.5% bromazepam, 23.6% diazepam, 14.2% alprazolam and 14% Z-drugs). In total, 30% of students who used psychoactive medicines obtained them for free, through a friend or relative (mostly codeine). Moreover, nearly 8% of students using psychoactive medications stated that they had needed to increase their doses. Once again, opioid substances (e.g. codeine) are the most frequently quoted psychoactive medications. Opioids are most likely used in the preparation of 'Purple Drank' (a recreational drug beverage, prepared by combining prescriptiongrade cough syrup with a soft drink and hard candy): 15% of codeine users reported consuming half to a whole bottle. Similarly, more than one in six (1/6) students who took Z-drugs and benzodiazepines had to increase their prescribed doses.

The purpose for all illicit drug use is recreational ('to party') or experimental ('to experience new states'). 'Creativity' appears to be an effect of some products, such as LSD (lysergic acid diethylamide) [71%], hallucinogenic mushrooms (59%), ecstasy (26.5%), cannabinoids (32%) and cannabis (23%). 'When I'm bored' is a specific motivation for cannabinoids and cannabis (20.0% and 25.0%, respectively). Of course, the first most commonly cited reason for all psychoactive prescription drug use is 'prescribed treatment'. However, the same motivations (e.g. partying, identity and experimentation) are reported for psychostimulant medicine (up to 34%), or dextromethorphan (17%), codeine (up to 15%) and benzodiazepine/Z-drug (almost 4%) use.

Modafinil stands out in terms of reason for motivation with only 26% of cases stating 'treatment' as the primary reason.

In total, 18.6% of students (17.5% women vs 20.4% men) reported 'managing stress' as the motivation for consumption of psychoactive medications (14.5% with illicit drugs and 5.52% with prescription drugs). Among them, three-quarters used cannabis and half used benzodiazepines.

Likewise, 'managing sleep' was the intended effect for 14.1% of all students (8.9% with illicit drugs and 6.9% with prescription drugs) [note: there was no gender effect]. Among these individuals, twothirds used cannabis and one-quarter use benzodiazepines, including Z-drugs. Nearly 20% also reported managing sleep with codeine.

The results of this study indicated that NE in students is a problem, with 18.6% of students in the COSYS survey confirming use of psychoactive medications. Only 3.81% of students use psychoactive substances specifically during examination periods (4.66% women vs 2.36% men) and 95% of these are prescribed medications. Specifically, these are stimulants (piracetam 79%, modafinil 63% and methylphenidate 62%) or anxiolytic drugs (etifoxin in 56%, bromazepam 32%, zolpidem 29%). Nearly 4% of all students use psychoactive substances to improve concentration and alertness. Doping motivations with stimulants (both prescription drugs and illicit products) were cited by 3.1% of students. The search for both academic and social performance is emphasised by the use of medications and illicit drugs ('one has to succeed in all areas'). In fact, a very low prevalence for psychostimulant prescription drugs (0.57% of men), mostly NMU (nearly 67% of cases), was observed for participants aged >21 years or those having a disability status. Younger students do not appear to use stimulant medication. A significant trend of use corresponded with specific academic disciplines, including 'Health' and 'Law/Economy'. Multiple correspondence analysis yields three different profiles: doping candidate, experimenter and psychiatric profile (see table 2). The diversity of user profiles complicates the implementation of prevention programmes.

Discussion

NMU and conventional use of prescription psychoactive drugs is highly prevalent in French students, notably women. NMU prevalence is close to the prevalence of illicit drug use. NMU is robustly associated with problematic substance use, substance use disorders and psychopathology. [START_REF] Ty S Schepis 1 , Dalton L Klare 1 | Prescription drug misuse : taking a lifespan perspective[END_REF] NMU of opioids is associated with suicidality in adults aged ≥18 years [START_REF] Ashrafioun | The relationship between past 12-month suicidality and reasons for prescription opioid[END_REF] . Benzodiazepine NMU is also associated with substance use disorders, psychopathology and suicidality, [START_REF] Victoria R Votaw | The epidemiology of benzodiazepine misuse: A systematic review[END_REF] and stimulant NMU with neuropsychological impairment and psychopathology [START_REF] Wilens | Neuropsychological functioning in college students who misuse prescription stimulants[END_REF] . Thus, the implementation of intervention strategies are necessary.

Unlike the current situation in the US or Iceland, the NMU of stimulants is anecdotal in France.

Although improving academic performance is an often-cited motivation for NMU, stimulant NMU has not been shown to improve academic performance in those without a diagnosis of attention deficit hyperactivity disorder (ADHD) [START_REF] Peterkin | Cognitive performance enhancement: misuse or self-treatment?[END_REF] . Prescription stimulant NMU is believed to be a 'good study tool' by some students, which enables them to study. Moreover, because the use of stimulants as cognitive enhancers is, in part, a reaction to social pressures to seek high levels of achievement and to overcome personal shortcomings at all costs [START_REF] Cleveland | What's missing from the current smart drugs debate[END_REF] , changing this culture of high expectations could reduce the demand for performance-enhancing drugs. One study suggests that theory-guided and tailored anti-NMU campaigns have great potential in changing students' attitudes and behaviours towards NMU of medications [START_REF] Su | Study Natural" without Drugs: An Exploratory Study of Theory-Guided and Tailored Health Campaign Interventions to Prevent Nonmedical Use of Prescription Stimulants in College Students[END_REF] . Social norms and social media increase NMU of psychoactive substances, but also provide a potential platform for of anti-NMU campaigns.

Table 1 :

 1 Prevalence of psychoactive substances and motivations of use(2017)

	Psychoactive	Prevalence WOMEN	MEN	Nonmedical Use	Recreational	Improving	Managing
	substances		(%)	(%)	(NMU)	use	vigilance/concentration	stress/sleep (%/%)
							(%/%)
	Legal substances
	Prescribed
	medication
	Antidepressants		1.47%	0.79%			Fluoxetine (-/17%)	Fluoxetine (60%/34%)
	Benzodiazepines		4.04%	2.35%	Bromazepam	4%	Buspirone (-/19%)	Alprazolam (87%/75%)
					27.5%,			Bromazepam (77%/75%)
					Diazepam 23.6%,			Oxazepam (77% /75%)
		Women 22.48%			Alprazolam 14.2%, Z-drugs 14%			Clorazepam(75% /92%) Diazepam (65%/65%) Zopiclone (40%/89%)
	Other anxiolytics	Men 15.14%	1.76%	0.76%				Zolpidem (36%/94%) Etifoxine (90%/65%)
	Other hypnotics		1.89%	1.09%				Doxylamine (39%/94%)
	Stimulants		0.21%	0.57%	Modafinil 60.8%	34%	Modafinil (96%/88%)
					Methylphenidate		Piracetam (71%/86%)
					32.9%		Methylphenidate(71%/83%)
	Opioids		17.09%	11.53%	Codeine 30%	Codeine 15%		Codeine (16%/21%)
	Non-prescribed
	Alcohol		79.78%	84.74%
	Tobacco		35.52%	35.77%
	E-cigarette		6.55%	10.65%
	Poppers		6.22%	10.35%				(17%/-)
	Nitrous oxide		2.81%	5.65%
	Cannabis		16.86%	26.99%			(3%/ 6%)	(66%/42%)
	Illicit substances
	Ecstasy/MDMA	Women	2.98%	5.79%			(31%/ 11%)	(26%/-)
	Cocaine	8.12%	1.85%	3.57%			(50%/ 28%)	(16%/-)
	Hallucinogenic		0.96%	2.69%				(26%/2%)
	mushrooms

Table 2 :

 2 Psychostimulant medication users and Multiple Correspondence Analysis (MCA) a

	VARIABLES STUDIED	PROFILE 1	PROFILE 2	PROFILE 3
	Sex	p = 0.0003
	Women		48 (50%)	5 (13.51%)	4 (26.67%)
	Men		48 (50%)	32 (86.49%)	11 (73.33%)
	Disability p = 0.002
	No		71 (73.96%)	37 (100%)	13 (86.67%)
	Yes		25 (26.04%)	0 (0%)	2 (13.33%)
	Job	p = 0.0352
	No		75 (78.13%)	26 (70.27%)	7 (46.67%)
	Yes		21 (21.88%)	11 (29.73%)	8 (53.33%)
	Legally obtained
	prescription p < 0.0001
	No		33 (39.76%)	31 (88.57%)	8 (66.67%)
	Yes		50 (60.24%)	4 (11.43%)	4 (33.33%)
	Antidepressants
	No		91 (94.79%)	36 (97.3%)	4 (26.67%)
	Yes		5 (5.21%)	1 (2.7%)	11 (73.33%)
	Anxiolytics
	No		94 (97,92%)	35 (94.59%)	4 (26.67%)
	Yes		2 (2,08%)	2 (5.41%)	11 (73.33%)
	Hypnotics
	No		85 (88.54%)	33 (89.19%)	7 (46.67%)
	Yes		11 (11,46%)	4 (10.81%)	8 (53.33%)
	Benzodiazepines
	No		84 (87.5%)	26 (70.27%)	1 (6.67%)
	Yes		12 (12.5%)	11 (29.73%)	14 (93.33%)
	Legal substances (excluding
	alcohol, tobacco, E-
	cigarette)
	No		87 (90.63%)	8 (21.62%)	1 (6.67%)
	Yes		9 (9.38%)	29 (78.38%)	14 (93

.33%) Illicit substances (excluding cannabis) No Yes 80 (83.33%) 16 (16.67%) 2 (5.41%) 35 (94.59%) 2 (13.33%) 13 (86.67%) a

  Profile 1 = doping candidate; profile 2 = experimenter; profile 3 = psychiatric profile.

Funding

Authors declare sources of funding from MILDECA (INTERMINISTERIAL MISSION FOR COMBATING DRUGS AND ADDICTIVE BEHAVIOURS).

Author statements

Ethical approval

The COSYS project was approved by the 'Commission nationale de l'informatique et des libertés' (National Commission of Informatics and Liberties, Reference 2004-009). Students were informed online on the nature and purpose of the study. The study protocol was evaluated by the Local Committee for the Protection of Persons of Paris University which indicated that this research was outside the scope of the provisions governing biomedical research (reference 2016/28NI).

Competing interests

None declared.