Felix Koeth
email: felix.koeth@g2elab.grenoble-inp.fr

Nicolas Retiere
email: nicolas.retiere@univ-grenoble-alpes.fr

Spectral properties of dynamical power systems

The integration of new renewable energies sources in the power system will stress the power system even more. To assure save operation, the stability of the power system has to be studied in detail. A well known feature in the power system are the inter-and intra-area oscillations of generators. This work investigates the origin and properties of this oscillations from a fundamental point of view, applied to a realistic power system test case.

I. INTRODUCTION

The power system is probably the worlds most complicated human made system and of vital importance for modern society. Stable operation of this system is required to ensure safety and well being a large group of people. But this system faces challenges due to climate change, as the default means of producing energy have to change to renewable energies with lower carbon footprints. This changes further complicate the problems of stable operation. So, more than ever, the stability of the power system has to be investigated and understood. The power system consist of a large number of different, coupled devices. The many different scales and complicated, time dependent dynamics result in a very complex system. To assess the interplay between the dissimilarity and the coupling of the system, simplified power system models are often used. The stability of this systems is usually investigated by finding a synchronized solution [START_REF] Dörfler | Synchronization in complex networks of phase oscillators: A survey[END_REF] and investigating the transient stability of this solutions after large perturbations [START_REF] Ji | Basin stability of the Kuramoto-like model in small networks[END_REF], [START_REF] Dörfler | On the Critical Coupling for Kuramoto Oscillators[END_REF], [START_REF] Rohden | Impact of network topology on synchrony of oscillatory power grids[END_REF]. Also, the small signal stability from small perturbations around the operating point has to be taken into account [START_REF] Motter | Spontaneous synchrony in power-grid networks[END_REF].

An interesting property found in stability investigations are the power system oscillations [START_REF] Obaid | Power system oscillations and control: Classifications and PSSs' design methods: A review[END_REF]. The fast intra-area oscillations describe the dynamics of neighboring generators, which swing against each other. The intra-area oscillators are slower and deal with large scale interactions of coherent regions of generators. For large-scale control of the power system, these clusters are often bundled to a aggregated generator and the intra-area modes describe the oscillations between bundled oscillators. Understanding the formation and properties of these clusters and the coherent dynamics are thus vital for the control and design of the power system. The objective of this work is to investigate how the dynamical parameters and topology of a generator network influences the spectral properties and the eigenvectors, which are the basis for the formation of coherent clusters.

*Institute of Engineering Univ. Grenoble Alpes

In this work, the spectral properties of a specific IEEE test case are investigated. At first, modeling background and theoretical results are introduced. Then, numerical experiments on the test case are carried out. The main question is the influence on the spectral properties of the parameters from the system, mainly the generator properties and the topology of the network. Also, the robustness of the spectrum under parameter changes is investigated. In a last step, the existence of coherency in the eigenvectors is investigated.

II. MODELING AND THEORY

In this section, the simplified power system model is introduced. The chosen model is the structure preserving model (SP-Model). The model is linearized and the mathematical problem statement is explained. Also, a review of the known analytical properties of the mathematical problem is given.

A. SP Model and quadratic eigenvalue problem

A common way to investigate fundamental properties of power systems and the interplay between the topology, parameters and the system dynamics is to use simplified power system models. Many aspects of complex power system models, as the reactive power flow, load dynamics and active control effects are ignored. The differential equations describing the dynamics of the phases θ i at the generator (V G ) and load (V L ) nodes i are given by [START_REF] Dorfler | Synchronization in complex oscillator networks and smart grids[END_REF], [START_REF] Filatrella | Analysis of a power grid using a Kuramoto-like model[END_REF], [START_REF] Nishikawa | Comparative analysis of existing models for power-grid synchronization[END_REF]:

M i θi + D i θi = ω i - j a ij sin (θ i -θ j ) i ∈ V 1 (1)
D i θi = ω i - j a ij sin (θ i -θ j ) i ∈ V 2 (2)
With the inertia M , the damping D, the power input/output ω and the adjacency matrix A with entries a ij of the weighted graph G, which is given in power system nomenclature as a ij = U i U j Im Y ij with the bus Voltages U and the nodal admittance matrix of the network Y . For the remainder of this work, the damping of the load nodes is assumed to be zero. The resulting equations for the load nodes are thus purely algebraic. The SP-Model was extensively used in the studies linked in the introduction, especially for transient stability studies. It is said to capture the dynamics of the first swing of generators reasonably well. Usually, no external control is used. The damping term D θ incorporates not only mechanical damping, but also the rotor-stator flux effects and machine specific controls.

For small perturbations around a fixed point, the small angle approximations leads to the following linear system:

M θ + D θ + Lθ = ω (3)
With L being the Laplacian matrix of the Graph G. Depending on the model, the diagonal matrices M and D have have zero entries for the nodes i ∈ V G . The main interest is the dynamical behavior of the generators, so the responsefree system with ω = 0 is considered. The solutions to this equation can be found by solving the following quadratic eigenvalue problem (QEP):

λ 2 M + λD + L x = 0 (4)
In accordance to the eigenvalue problem, we call λ k the eigenvalues and x k the (right) eigenvector of the QEP. Due to the quadratic nature of the problem, we have 2n solutions to the quadratic eigenvalue problem. Using γ = [γ 1 . . . γ 2n ] as a vector of arbitrary constants related to the initial values of the problem, the solution θ can then be calculated as [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]:

θ(t) = 2n k γ k x k exp (λ k t) (5)
From equation [START_REF] Motter | Spontaneous synchrony in power-grid networks[END_REF], coherency in the phases of the power system can be defined. If, for a given mode λ k , some components of the eigenvector x k are (roughly) identical, this modes will show similar response to a excitation of this mode. As the eigenvectors of the quadratic eigenvalue problem are generally complex, coherent behavior can be observed if the complex phase and magnitude of eigenvector components are similar. The quadratic (and higher order) eigenvalue problems are usually calculated by using the companion form, where the QEP with n × n matrices is transformed to a generalized eigenvalue problem with 2n × 2n matrices [START_REF] Walsh | Quadratic eigenvalue problems[END_REF].

B. Kron reduction

As previously mentioned, the assumption of undamped loads results in algebraic equations for the load nodes. This load nodes will thus not exhibit a dynamic response, but change infinity fast. This can be used to simplify the systems equations, by writing the Laplacian matrix L and the inertia M in block form and solving the eigenvalue problem for the M, L matrix pencil1 :

L G L GL L T GL L L v G v L = λ 2 M G 0 0 0 v G v L
The eigenvector for the load nodes can be removed by inserting the second row of the matrix equation in the first row, resulting in:

L G v G + L GL v L = L G v G -L GL L -1 L L T GL v G =Lredv G = λ 2 M G v G
Which is identical to the Shur complement of the Laplacian matrix with respect to the sub-matrix of load nodes. The reduced Laplacian L red is of size n G ×n G . The reduced matrix is still a Laplacian matrix [START_REF] Dorfler | Kron Reduction of Graphs With Applications to Electrical Networks[END_REF].

This process is identical to the well known Kron reduction, where all load nodes are removed and the admittances between the generators are augmented. This process results in a complete graph. The underlying topological properties are encoded in the changed edge-weights of the complete graph. The other matrices of the QEP are now the submatrices M G and D G , which are full diagonal sub matrices of the generator machine properties.

C. Theoretical background

The main known properties of the quadratic eigenvalue problem are discussed in [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]. For the reduced eigenvalue problem stated here, the following properties of the QEP are respected:

1) M G is non-singular -2n G finite eigenvalues 2) M G , D G , L red are real -eigenvalues are real or complex conjugate pairs 3) M G hermitian and positive definite, D G and L red positive semi-definite -The real parts of the eigenvalues are non-positive A fourth condition is known as the overdamped condition, named after the related phenomenon in mechanical systems. For this condition, M G and D G have to be symmetric and positive definite and L red is positive semi-definite. Additionally, γ (M, D, L) > 0, with γ given as:

γ (M, D, L) = min x 2=1 (x * Dx) 2 -4 (x * M x) (x * Lx) (6)
This requirement can be studied by a simple approximation to the optimization problem. An upper γ ≥ γ up bond for the condition can be verified to be:

γ up (M, D, L) = λ 2 D,min -4λ M,max λ L,max > 0 (7)
Here, λ X,min/max is the smallest/largest eigenvalue of the matrix X. Overdamped behavior can be archived by increasing the damping coefficient of the generators or by lowering the inertia or connectivity of the network.

III. NUMERICAL RESULTS

A. The IEEE 145 system

The objective of this work is to investigate how the parameters of the network and the generators influence the spectral and coherent properties of the quadratic eigenvalue problem. For this, a suitable test case and meaningful parameters have to be found. While the nodal admittance matrix Y , the power in-and outputs ω and the voltages U i are known or can be calculated by the default load-flow equations, the dynamic parameters are not always well known. In this study, the methods to calculate parameters found in [START_REF] Nishikawa | Comparative analysis of existing models for power-grid synchronization[END_REF], [START_REF] Motter | Spontaneous synchrony in power-grid networks[END_REF] are used. The main test system is dynamic 50 generator test case (sometimes called IEEE 145 system) [START_REF]Power Systems Test Case Archive -UWEE[END_REF]. It features 50 generators and has 145 buses. The system is mainly used for dynamical studies and some generator data is known in literature.

B. The Laplacian spectrum

In figure 1, the spectrum λ(L) of the Laplacian matrix L is shown, ranked by the absolute value. Not included is the zero eigenvalue of the Laplacian matrix. It can be observed that the first few eigenvalues are relatively small and increase slowly, which is typical for Laplacian matrices of power systems.

Also included in Graph 1 is the spectrum of the Kron reduced Laplacian L red , as described in the theory. The resulting reduced system is less complex (as it has less nodes).

The full system has a larger maximum ranked eigenvalue, but for small ranked eigenvalues, the reduced spectrum is larger, corresponding to a higher connectivity of the network.

C. Inertia

At first, the influence of the inertia in the spectrum is investigated. For that, the damping matrix is set to zero at first. To assess how the parameter of the inertia (e.g. the dissimilarities of the diagonal entries of M ) influence the spectrum, the inertia M is perturbed to M by α = (α 1 , . . . , α N G ), with M = diag (M 1 α 1 , . . . , M N G α N G ) and uniform distributed α with α ∈ [0.5, 2.0]. Figure 2 shows the original spectrum and spectra of 200 realizations with different M i . It can be noted that for smaller ranks, the differences in spectra are relatively small. At the same time, the range of eigenvalues grows larger for higher ranked eigenvalues.

The highest ranked eigenvalues are linked to intra-area oscillations, where individual generators oscillate against each other. This localized dynamics is stronger influenced by the inertia, which is a local effect on a single node.

D. Damping

One major simplification up to this point is the fact that the generators were undamped. Realistic generators are damped, to reduce the instabilities caused by perturbations from their operating point. Here, the interest is how this intrinsic damping influences the spectra of power systems.

Similarly to the previous section, figure 3 shows the spectra for different damping realizations around the original damping ratio D 0 . Additionally, spectra without any damping is included. Here, the distribution of α is chosen larger as α ∈ [0.25, 4.0]. Two observations are directly clear: even though the random range is larger, the difference seems, for the most part, to be smaller for different damping perturbations. Also, the difference in the spectrum without damping and the default damping parameter is very small. As before, one can notice the difference for larger and smaller damped eigenvalues. Larger damped eigenvalues seem more affected by the localized parameters of D, while the lower ranked eigenvalues are almost unchanged and more depended on the structure of the underlying network.

The observations indicate that the spectrum is not at all overdamped. From equation [START_REF] Dorfler | Synchronization in complex oscillator networks and smart grids[END_REF], one can estimate that overdamped behavior would be reached if the eigenvector of D is multiplied by a factor of α OD ≈ 1130. This is very far from the assumed conditions in the power system and one can assume that power systems are not overdamped. The spectra for under-and overdamped behavior is shown in figure 4, where the damping matrix D is calculated by D 0 α, with different values for α. Again, it is clear that low to moderate values of α does not change the spectrum significantly. For very large α, the observations from the theory are observable. The eigenvalues turn purely real and the gap between the first n = 50 and last n eigenvalues is clearly present.

E. Sensitivity

The previous sections show that the introduction of inertia and damping changes the results of the quadratic eigenvalue problem, but the general shape of the spectra seems to be dominated by the topology of the network. To answer which parameter, damping or inertia, influences the system stronger, the following experiment was carried out:

1) The spectra λ 0 (M 0 , D 0 , L red ) was calculated.

2) The diagonal matrices M 0 and D 0 were perturbed by random, uniform distributed vectors r M and r D to

M i = M 0 r M,i and D i = D 0 r D,i . 3) The spectra λ M i (M i , D 0 , L red ), λ D i (M 0 , D i , L red ) and λ M,D i (M i , D i , L red )
were calculated for 200 different realisations. 4) The average differences from (λ i -λ 0 ) and variances were calculated for all cases of λ i . The resulting graphs are shown in figure 5. All perturbations r i are uniform distributed 0.5 < r i < 1.5. From the figure, it appears that changing both M and D has a very similar effect to only changing the inertia. Only changing the D has a much smaller effect, as the variance and the difference to the unperturbed values is much lower. While the behavior of the individual ranks are not directly clear, especially the highest ranked eigenvalue has the largest difference and variance, for all considered eigenvalues.

Another tool to investigate the global sensitivity of matrices, and in more general matrix pencils, are the pseudospectra [START_REF] Tisseur | Structured Pseudospectra for Polynomial Eigenvalue Problems, with Applications[END_REF], which gives insights on the region in which eigenvalues change under matrix perturbations with a certain strength. For the power system used here, the pseudo-spectrum is given in figure 6. In this figure, the markers show the real-and imaginary-part of the quadratic eigenvalue problem. The colors correspond to a perturbation strength of the matrices M , D and L and the perturbed eigenvalues will lie in a region corresponding to the color.

From the pseudo-spectrum, it is clear that the eigenvalues with a larger negative real part can be perturbed with weaker perturbations. At the same time, changing this value will not fundamentally alter the dynamics of the full system. The more interesting eigenvalues are those with a small real part. To perturb this eigenvalues, much larger perturbations are necessary.

From this observations, the spectral properties seem to depend mostly on the structure of the underlying network, with the exception of the highest ranked eigenvalues. In general, the inertia seems to have a bigger effect on the spectrum then the damping and for small damping values, the damping can be neglected. The next question is how do the parameters influence the coherency, or the eigenvector structure of the system.

IV. COHERENCY OBSERVATIONS

The general structure of the eigenvalues for all modes of the given problem is shown in figure 7. The color indicates the magnitude of the eigenvector components for each mode.

For smaller modes, most nodes are active. Similar values can be found for some bands in the first modes. At higher modes (> 60), the dynamics are more isolated. Only a few nodes oscillate at that mode at all. On the graph, the phases and magnitudes of the eigenvector components of the second smallest eigenvalue (called Fielder vector) are shown in figure 8, where the color of the nodes corresponds to the complex phase while the size of the nodes is the complex magnitudes. To improve observability, only edge-nodes with a edge weight2 over 0.5 are shown in this figure. The isolated nodes without any connections are only weakly connected to the rest of the system. Nodes with a similar color and size will be coherent. Two observations are interesting: The complex phases are all very close to ±π/2. In the time response, this corresponds to a phase shift of π, so both clusters oscillate anti-phasic. Also, most nodes are in one, large cluster, with different magnitudes. Only a few nodes, included in a small, interconnected region are antiphasic. Those have smaller magnitude, but the magnitudes are not identical over all oscillators.

In contrast, the 38h eigenvector components are displayed in figure 9. The magnitudes of a large number of nodes is almost zero, meaning that the node does not react to excitations of this mode. Some single nodes, especially larger blue circle, dominate this mode. Coherency is almost not found, except for the fact that many modes not oscillate at all.

The eigenvector properties can be investigated in a polar projection, as shown in figure 10. As seen earlier, most eigenvector components for the higher ranked eigenvalue are centered around zero, indicating no oscillation in this modes. Single nodes have a large, individual share of the dynamics in this mode. Fourth mode, the clustering of the most modes can be observed at a small magnitude around 260 degrees. The 18th eigenvector has most modes active at small magnitudes, centered around zero. The 38th mode has, as seen before, only a few active nodes, while most modes are almost zero. Also, noteworthy is the alignment of most eigenvalues around the ±π/2 axis (90 • to 270 • ). This corresponds to a low damping value in the system, as the undamped system is a hermitian generalized eigenvalue problem, where the eigenvectors will be purely imaginary.

V. CONCLUSION

In this work, the influence of the inertia and damping on the spectral properties of a linearized power system model was investigated. It was shown that the damping has, for physical relevant values, only small effect on the structure of the eigenvalues. The influence of the local parameters is mostly noted at the higher ranked eigenvalues, while the lower ranked eigenvalues are dominated by the global properties of the Laplacian matrix. Also, general questions of the robustness of the structure of the spectrum was investigated, with the result that the spectral properties are relatively robust for moderate changes in the relevant parameters.

A theoretical explanation of slow coherency is given in [START_REF] Ayazifar | Graph spectra and modal dynamics of oscillatory networks[END_REF]. There, the case of D = 0 was investigated analytically. It was shown that a network which can be partitioned in q clusters, which are weakly connected to each other, will show slow coherency in the q slowest modes. In the reduced realistic system as used in this work, the conditions for this theorem are not directly observed. As this condition is only sufficient, coherency can be observed in different systems as well.

In general, the structure of eigenvectors are a very difficult mathematical topic and only little analytical results are known to the authors. The most recent advance is proving an upper bound for the number of nodal domains (regions where the sign of the eigenvectors does not change) in homogeneous graphs [START_REF] Urschel | Nodal decompositions of graphs[END_REF]. The extension to the quadratic eigenvalue problem (and their complex nature) and to weighted, non-homogeneous graphs is not directly possible.

Even though the system does not respect the theoretical properties, the system shows coherent behavior, as shown in the last section. Future work will focus on the structure of the eigenvectors. While general theoretical results might be impossible to get, the idea is to investigate the undamped system (which simplifies the system from a QEP to a generalized eigenvalue problem) of a homogeneous complete graph with a general weight matrix M . The results will be perturbed by including weights in the complete graph, which should allow for insights of the structure of the eigenvectors of the undamped system.

Fig. 1 .

 1 Fig. 1. Spectrum of Laplacian Matrix L and reduced matrix L red .

Fig. 2 .

 2 Fig. 2. Spectrum of eigenvalues of matrix pair (M, L) for IEEE 145 system with randomly perturbed inertia M . The grey lines indicate perturbed spectra.

Fig. 3 .

 3 Fig. 3. Spectrum of eigenvalues of matrix pencil (M, D, L) for IEEE 145 system with randomly perturbed damping D. The grey lines indicate perturbed spectra.

Fig. 4 .

 4 Fig. 4. Normalized spectra for different damping factors α.

Fig. 5 .

 5 Fig. 5. Average spectral difference and variance of perturbed spectra for perturbed M , D or both.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Pseudo-spectrum of quadratic eigenvalue problem. Marker indicate the original eigenvalues, the color indicate the region in which the eigenvalue can fall after a perturbation of strength corresponding to the color.

Fig. 8 .

 8 Fig. 8. Phase and magnitude of fourth ranked eigenvalue.

Fig. 9 .

 9 Fig. 9. Phase and magnitude of 38th ranked eigenvalue.

Fig. 10 .

 10 Fig. 10. Polar projection of eigenvector components for different modes.

For simplicity, the damping is ignored here, but the calculations are identical for damped generators.

The edge weight between two nodes is proportional to the inverse of the admittance between the two nodes.

VI. ACKNOWLEDGMENTS

This project has received funding from the European union's Horizon 2020 research and innovation programme under Marie Sklodowska-Curie grant agreement No 675318


Ibrahima Diop

Mathieu Carbone

Sebastien Ordas

Yanis Linge

Pierre Yvan Liardet

Philippe Maurine

Collision for estimating SCA Measurement Quality and Related Applications

If the Signal to Noise Ratio (SNR) is a figure of merit commonly used in many areas to gauge the quality of analogue measurements, its use in the context of Side-Channel Attacks (SCA) remains very difficult because the nature and characteristics of the signals are not a priori known. Consequently, the SNR is rarely used in this latter area to gauge the quality of measurements or of experimental protocols followed to acquire them. It is however used to quantify the amount of leakage in a set of traces regardless of the quality of the measures. This is a surprisingly lack despite the key role of measurements and experiments in this field. In this context, this paper introduces a fast and accurate method for estimating the SNR. Then, simple and accurate techniques are derived. They allow to process some daily tasks the evaluators have to perform in a pragmatic and efficient manner. Among them one can find the analysis of the electrical activity of Integrated Circuit (IC) or the identification of the frequencies carrying some information or leakages.

Introduction

In recent years, emphasis has been put on the localization of leaking content in space [START_REF] Réal | Enhancing correlation electromagnetic attack using planar near-field cartography[END_REF], [START_REF] Dehbaoui | Enhancing electromagnetic attacks using spectral coherence based cartography[END_REF], [START_REF] Sauvage | Blind cartography for side channel attacks: Cross-correlation cartography[END_REF] and time [START_REF] Bhasin | NICV: normalized inter-class variance for detection of side-channel leakage[END_REF], but also to increase our understanding of the latest or again to enhance our extraction capability [START_REF] Kocher | Differential power analysis[END_REF], [START_REF] Brier | Correlation power analysis with a leakage model[END_REF], [START_REF] Gierlichs | Mutual information analysis[END_REF]. Efforts have also been invested to assess the efficiency of attacks resulting in notions like the success rate and the guessing entropy.

Based on [START_REF] Standaert | A Unified Framework for the Analysis of Side-Channel Key Recovery Attacks[END_REF], one generally considers two types of evaluation metrics for assessing the leakage emitted by cryptographic devices. First, information theoretic metrics aim to capture the amount of information available in a sidechannel leakage, independently of the adversary exploitation. Second, security metrics aim to quantify how this information can be exploited by adversary using the notions of Success Rate (SR) and Guessing Entropy (GE). In SCA context, these two types of metrics are clearly related but they are also surrounded by the quality of measurements for which less attention has been paid to define criteria allowing to compare experimental practices that are central.

The Signal-to-Noise Ratio (SNR) is a simple, intuitive and relevant metric to characterize any analog phenomena, such as side channel measurements in order to gauge the quality of acquisition campaign. However the first definition of SNR [START_REF] Mangard | Hardware countermeasures against dpa ? a statistical analysis of their effectiveness[END_REF] in the SCA context has some drawbacks making its practical use difficult given that the signal and its characteristics are a priori not known. Some workarounds have been proposed in the literature but most of them rely in converting the problem into leakage quantification or Points of Interest (PoI) detection rather than focusing the quality of the acquisition campaigns. Among these works, one may include the approach proposed in 2011 [START_REF] Guilley | Quantifying the quality of side channel acquisitions[END_REF] that offers a method for estimating the SNR which is based on a linear model of the leakage but also on empirical considerations related to the use of Digital Sampling Oscilloscopes (DSO). However, the SNR definition considered in this work is relatively far from the standard electrical engineering definition and finally closer to that of a Leakage to Noise Ratio (LNR) concept.

This latter concept appears explicitly in [START_REF] Tiran | A model of the leakage in the frequency domain and its application to CPA and DPA[END_REF] with an estimation method based on a modeling of the leakage in the frequency domain. Alongside the introduction of the LNR, authors in [START_REF] Bhasin | NICV: normalized inter-class variance for detection of side-channel leakage[END_REF] propose the Normalized Inter-Class Variance (NICV). This approach based on a clever and effective use of the statistical tool of ANOVA (ANalysis Of VAriance) where the F-test allows to identify the PoI thanks to the knowledge of public information (plaintext or ciphertext) (i.e. allows to detect interesting time samples carrying potential leakage based on the manipulation of the inputs). Besides this important application, it is also suggested that the NICV provides an estimate of the SNR. However, again the SNR definition considered in this paper is more an estimate of the relative amount of leakage in the samples of traces that an estimate of the SNR of a measurement set.

In this context, the main contribution of this paper to the State of the Art is a method, based on signal collisions, to estimate the SNR. This method requires a small amount of measures to be applied: two in theory and up to few tens for convenient practice. The second contribution is to show that this SNR estimate is also an effective way to gauge the quality of traces in view of the application of SCA. Additionally, it is then shown how to exploit the above technique to:

analyze the temporal behavior of a circuit from a small amount of measures, distinguish the frequencies carrying information related to the signal from the frequencies essentially made up of noise.

The distinction of frequencies carrying information is crucial. Indeed, it enables, without any a priori knowledge about the characteristics of the signal, to rationally and adaptively guide the filtering procedures. This paper is organized as follows: in section 2, the meaning of the term signal considered in this paper, and related SNR is defined. Section 3 recalls the standard electrical engineering definition of the SNR and the corollary definition thereof in the case of zero mean signals. These definitions recalled, section 4 demonstrates theoretically how to get an estimate of the SNR from signal collisions and more particularly from a bounded collision detection criterion introduced in [START_REF] Diop | Collision based attacks in practice[END_REF]. It is then experimentally demonstrated in section 5 that the obtained SNR estimate is actually a figure of merit for gauging the quality of traces in a SCA context. Similarly, section 6 shows how this fast estimate method of the SNR allows analyzing SCA traces in a simple, effective and complementary way to the approaches introduced in [START_REF] Bhasin | NICV: normalized inter-class variance for detection of side-channel leakage[END_REF] and [START_REF] Tiran | A model of the leakage in the frequency domain and its application to CPA and DPA[END_REF]. Section 7 shows how to accurately distinguish with a small amount of traces the harmonics of the signal from harmonics mainly related to noise. Finally, a conclusion is drawn in section 8.

Preamble

In this paper, the term signal is devoted to the deterministic evolution (repeated and measured n times with perfect equipments in the absence of any noise source, the signal is unique) of a physical quantity over time, evolving in response to an unique stimulus. For example the evolution of the magnetic field at the coordinate (X, Y, Z) above a given IC computing the cipher of a given text with a given key is a signal. The changing of one of these parameters by solely one bit gives a new signal.

In the classical literature such as [START_REF] Mangard | Power analysis attacks: Revealing the secrets of smart cards[END_REF], [START_REF] Bhasin | NICV: normalized inter-class variance for detection of side-channel leakage[END_REF], [START_REF] Mangard | Hardware countermeasures against dpa ? a statistical analysis of their effectiveness[END_REF], [START_REF] Tiran | A model of the leakage in the frequency domain and its application to CPA and DPA[END_REF], the power consumption of a measure (trace) is defined as the sum of the exploitable power consumption from an SCA point of view (P exp ), the power of the noise ( considered in [START_REF] Mangard | Power analysis attacks: Revealing the secrets of smart cards[END_REF] as the sum of the electronic noise P el.noise and of the switching noise P sw.noise ) and a constant power consumption (P const ). In the present paper, we consider the power consumption of a measure as the sum of the power consumptions of the signal (P S ) to be measured and of the noise (P N ). P S is considered regardless the power consumption of signal is exploitable or not. As a result, P S includes but is not reduced to P exp . This choice is done in order to evaluate the quality of the measurement protocol independently from the fact that traces contain exploitable information.

None of these meanings is preferable. They are adopted for different purposes. In our case, a paradigm shift is investigated to assess the quality of a measurement. This also allows to assess the quality of the related experimental protocol, as an alternative of quantifying the amount of leakage present in a set of SCA traces. However we make a parallel with our approach and this point later in the paper.

SNR definition and the related problem

With the definition of the signal adopted in this paper, the SNR is an objective figure of merit for gauging the quality of the measurement M (t) = [m 1 , ..., m q ] of a signal S(t) = [s 1 , ..., s q ], which depends on the surrounding noise, but also of the quality of the equipments used to collect it. The standard electrical engineering definition of the SNR is the ratio between the power P S of the signal to be measured and the power P N of noise

N (t) = [η 1 , ..., η q ]: SN R = P S P N = 1 q • q i=1 (s i ) 2 1 q • q i=1 (η i ) 2 (1)
If the signal and the noise are zero mean, the numerator and the denominator of Equation 1 are the signal and noise variances, respectively. This leads, considering S(t) and N (t) as random variables over time (horizontal random variables or random processes), to the following compact expression of the SNR:

SN R = σ 2 S σ 2 N = 1 q • q i=1 (s i -0) 2 1 q • q i=1 (η i -0) 2 (2)
with σ S and σ N denoting the standard deviations (measured horizontally and not vertically) of the signal S(t) and noise N (t), respectively. It should be noted that even if the SNR takes the form of a variance ratio, this is really a power ratio and must be perceived as it. This SNR definition is particularly interesting in the context of SCA exploiting the Electro-Magnetic (EM) channel. Indeed, the means of the EM signals are always very close to zero, or can even be fixed to zero by forcing the usual voltage amplifiers to work in small signals or by compensating their input offset voltage. This definition is therefore considered in the rest of this paper. However, given that in the context of SCA, the signal characteristics are not known, this definition is extremely difficult to apply in this form. A solution must therefore be found to extract from measures the SNR values ( Equation 2).

From Collisions to SNR estimations

If in the context of SCA the signal S(t) is not known, other degrees of freedom are available to adversaries. One of these degrees of freedom is the capability to stimulate the circuit n times in the same manner to collect n measures of a same execution or calculus. However, this degree of freedom is not unlimited. Indeed, some security protocols limit to few units the possibility to repeat a same execution. In addition, if this is not the case, repeating the same execution is theoretically possible but often tedious because of misalignment of traces which implies the use of some more or less costly realignment techniques. However, it seems possible to generate few colliding traces and exploit them to estimate the SNR as explained below.

From Bounded Collision Detection Criterion to SNR

The use of signal collisions is a recognized technique to attack the implementations of symmetric [START_REF] Schramm | A new class of collision attacks and its application to DES[END_REF] or asymmetric [START_REF] Fouque | The doubling attack -Why Upwards Is Better than Downwards[END_REF], [START_REF] Yen | Power analysis by exploiting chosen message and internal collisions -vulnerability of checking mechanism for rsadecryption[END_REF] algorithms. For this purpose, a choice of plaintexts is done in order to induce a collision between intermediate values according to the targeted secret. The adversary finally derives the secret from the occurrence or the absence of the collision. However, because in practice SCA traces are noisy, detection is not straightforward and some signal processing techniques are necessary to exploit collisions. Recently in [START_REF] Diop | Collision based attacks in practice[END_REF], a criterion that automates the detection of collisions was proposed. This criterion, called Bounded Collision Detection Criterion (BCDC), takes values in [0,1] and is defined by:

BCDC(M 1 , M 2 ) = 1 √ 2 × σ (M1-M2) σ (M1) (3)
where

M i = [m i 1 , • • • , m i q ]
is a measurement, collected with a Digital Sampling Oscilloscope (DSO), representing a measure of the signal S(t) i characterized by a (horizontal) standard deviation σ (Mi) . The idea behind this criterion is that in case of collision the numerator is close to zero, while in the absence of collision the BCDC is close to one.

Assuming that all samples m i j of these measurements (with j ∈ [1, q]) can be expressed by m i j = s i j + η i j , i.e. is the sum of: • a deterministic value (s i j ) which is a sample of S(t) emitted by the Device Under Test (DUT), at time t, during its operation,

• and of the realization of the noise (η i j ), which is a random variable, drawn in a normal distribution with zero mean and an unknown standard deviation, it is then possible to express the difference between M 1 and M 2 :

∆M = M 1 -M 2 = [δ m1 , • • • , δ mq ] (4)
with

δ mj = m 1 j -m 2 j = s 1 j + η 1 j -(s 2 j + η 2 j ). ( 5
)
Considering that all terms of Equation 5 are independent, it appears that the denominator of the BCDC, Equation 3, can be expressed as:

σ (M1-M2) = σ 2 ([s 1 1 -s 2 1 , • • • , s 1 q -s 2 q ]) + σ 2 ([η 1 1 -η 2 1 , • • • , η 1 q -η 2 q ]), (6)
expression that can be re-written:

σ (M1-M2) = 2 • σ 2 S + 2 • σ 2 N = √ 2 • σ 2 S + σ 2 N (7)
Starting from this expression, one can calculate the asymptotic BCDC value when M 1 and M 2 are two measures of the same signal, i.e. in case of a collision. Indeed, in that case, because s 1 j = s 2 j and thus because:

σ (M1-M2) = √ 2 • σ 2 N , (8)
equation 3 becomes:

BCDC(M 1 , M 2 ) = 1 √ 2 √ 2 • σ N σ 2 S + σ 2 N = 1 1 + σ 2 s σ 2 N = 1 1 + SN R(M 1 , M 2 ) . (9)
As a result, the quick and pragmatic estimate of the SNR becomes possible in the context of SCA. To do this, simply collect n (> 1) traces corresponding to the execution of the same calculus by the DUT; then calculate the 1 2 • n • (n -1) BCDC values and finally deduce as much values SN R(M i , M j ):

SN R(M i , M j ) = 1 BCDC 2 (M i , M j ) -1 (10)
Finally, the SNR of a set of traces is estimated by averaging these 11). In the remainder of the paper this estimation is denoted SN R.

1 2 • n • (n -1) SN R(M i , M j ) values (see Equation
SN R = 2 n • (n -1) n i=1 n j=i+1 ( 1 BCDC 2 (M i , M j ) -1) (11)
This mean value as well as the shape of the SN R(M i , M j ) distribution constitute figures of merit for gauging the quality of a SCA measures and that of the experimental protocol that has been followed to collect them: the higher SN R is, the better the experimental protocol is. As an illustration, typical SN R(M i , M j ) distributions are given Figure 2c and Figure 2d in section 5. Next paragraphs give an illustration in a school book case.

Example

To provide practical and convincing elements about the accuracy and the simplicity of this fast SNR estimation method, which is applicable when the signal and the data are not known, simulated traces were generated. They relate to the measure of a sinusoidal signal of amplitude equal to 1 and with an horizontal variance over a period of 2 (σ S = √ 2 on a period) in an environment producing a normal measurement noise with zero mean and standard deviation σ N .

Table 2 gives the exact SNR value wrt σ N as well as SN R wrt the number of measurements (n = {2, ..., 50}) used to estimate it with the proposed method. One may observe that with only five colliding measures (i.e. 10 BCDC and SN R(M i , M j ) estimates) SN R n=5 is really close to the exact SNR value. Indeed, the relative absolute value is lower than 5%. This result appears even more accurate when σ SN R is low. For n = 50 measures (i.e. 1275 BCDC values) and for σ N = 1.9, SN R n=50 is equal to 0.1385 and σ SN R is equal to 0.0048 i.e. 3.45% of the exact SNR value. These results demonstrate the practicability as well the accuracy of the proposed method.

Discussion

At this stage, a fast SNR estimation method to be applied when the signal is not known has been introduced and its accuracy demonstrated. It is based on the SNR definition commonly used to assess the quality of analog measures or 0.1 0.3 0.5 0.7 0.9 1.1 1.3 1.5 1.7 1.9 SNR 50 5.56 2.00 1.02 0.62 0.41 0.30 0.22 0.17 0.14 SN Rn=2 49.38 5.58 2.03 1.00 0.64 0.41 0.28 0.24 0.17 0.16 SN Rn=5 49.76 5.51 1.99 1.03 0.60 0.41 0.29 0.21 0.17 0.14 SN Rn=10 50.41 5.56 1.99 1.02 0.62 0.41 0.29 0.21 0.18 0.14 SN Rn=20 49.91 5.52 2.00 1.03 0.62 0.41 0.30 0.22 0.17 0.14 SN Rn=50 50.08 5.56 2.01 1.01 0.62 0.41 0.30 0.22 0.17 0.14 the quality of a communicating channel. This method is the main contribution of this paper.

If the definition considered in this paper is different to that adopted in [START_REF] Mangard | Hardware countermeasures against dpa ? a statistical analysis of their effectiveness[END_REF], [START_REF] Mangard | Power analysis attacks: Revealing the secrets of smart cards[END_REF] and [START_REF] Bhasin | NICV: normalized inter-class variance for detection of side-channel leakage[END_REF], in which the authors aim at quantifying the amount of exploitable (with CPA) leakage in set of traces, nothing forbids to assess whether the definition considered in this paper is also an indirect measure of the amount of leakage present in traces. After all, measures of high quality are they not susceptible to convey more leakage than bad quality measurements?

This can be envisaged since the SNR definition adopted in this paper is the ratio between the signal and noise powers, and thus because the signal power is necessarily an upper bound of the exploitable leakage power Q L , i.e. the power consumed by the cryptographic operation targeted by the adversary:

σ 2 S σ 2 N ≥ σ 2 Q L σ 2 N ≥ 0 ( 12
)
If the SNR definition considered in this paper is indeed a figure of merit of SCA trace quality, the idea we propose to quantify the quality of a set of SCA traces is simple. It consists in acquiring a small amount of traces (2 to 50 for instance) of a same computational activity (or pairs of the same activity) in order to evaluate SN R and if necessary to draw the SN R(M i , M j ) distribution. This value (or the estimated distribution law) will be considered as a quality measure of SCA traces but also of the experimental protocol applied to collect the traces.

Estimated SNR wrt CPA efficiency

To assess if the SNR definition considered in this paper, dedicated to the evaluation of the quality of measures and of the experimental protocol used to obtain it, is also an indicator of the amount of leakage in a set traces of SCA and more precisely an upper bound, various experiments were conducted. One challenge was to define an experiment allowing to vary the SNR.

Experiments

In order to vary the SNR, different experimental campaigns were conducted, and different sets of EM traces, characterized by different SNR values were collected. These campaigns have consisted in measuring the EM radiations of an AES mapped into an FPGA. The EM sensor was placed at different distances Z from the IC surface and this without changing the settings of the DSO (the vertical caliber and the time base) so that to reduce the amplitude (the power) of the signal while increasing the noise. For each considered Z value, 50 EM traces of 20000 samples corresponding to ciphering of the same plaintext were acquired to compute SN R; 5000 others traces corresponding each to the ciphering of a random plaintext were also acquired with the intent to perform CPA. This experimental procedure was defined to vary the SNR but also because measuring the EM radiations so far from the IC (when it is possible at Z = 0) is a really bad experimental protocol to collect EM traces. Figure 1 illustrates the principle of these experiments and reports 25 EM measurements obtained at Z = 0 and Z = 500µm while the AES ciphers the same plaintext. The black thick traces are the mean of the 25 measures. One can observe that the amplitudes of the measurements performed at Z = 500µm are lower than the ones at Z = 0, but above all, that some large amplitude noises alter some traces. These spurious noises render the reading of traces less intelligible. Indeed, contrarily to measures at Z = 0, the rounds of the AES are less visible, and completely disappear in the noise at Z = 4000µm.

traces @ Z=0 traces @ Z=500µm

Fig. 1. 25 EM Measures EM collected at Z = 0 and Z = 500µm above an AES mapped into an FPGA

Analysis of the experimental results

A CPA was applied to each set of 5000 traces collected at Z values ranging between 0 and 5mm. Figure 2 gives the evolution with Z of σ S and SN R. As expected the SN R decreases rapidly when Z increases; the distributions of the estimated SN R values (with Equation 2) for Z = 0 and Z = 2000µm are given Figures 2c andFigures 2d, respectively. For measures done at Z = 0, SN R is equal to 15.11 and σ √ SN R is equal to 4.80 while at Z = 2000µm, SN R is equal to 7.31 and σ √ SN R to 1.64. This demonstrates that the SNR is, as expected, a measure of the quality of the experimental protocol.

Figure 2b gives the evolution of the global Guessing Entropy (gGE) with respect to SN R. One can observe that the lower the SN R value is, the less efficient the CPA is. In addition, as soon as SN R is lower than 5, the CPA does not succeed in disclosing the correct key with 5000 traces and more over the gGE remains close to 128. Nevertheless, Figure 2b highlights the existence of a link between gGE and the SNR. The standard electrical engineering definition of the SNR (Equation 2) is therefore a figure of merit of the quality of EM traces from an adversary point of view, as the SNR definition considered in [START_REF] Mangard | Hardware countermeasures against dpa ? a statistical analysis of their effectiveness[END_REF] and [START_REF] Bhasin | NICV: normalized inter-class variance for detection of side-channel leakage[END_REF]. The higher the SN R associated to a set of traces is, the leakier it is expected to be. However, it should be noticed that this figure of merit does not indicate possible problems that can encounter adversaries to exploit this information. These difficulties are for instance:

the presence or the absence of countermeasures such as masking [START_REF] Chari | Towards sound approaches to counteract power-analysis attacks[END_REF], the choice of the adequate distinguisher [START_REF] Gierlichs | Mutual information analysis[END_REF], [START_REF] Brier | Correlation power analysis with a leakage model[END_REF], the identification of the accurate leakage model [START_REF] Mangard | Power analysis attacks: Revealing the secrets of smart cards[END_REF], which are problems independent from the measurement quality.

Collision Criterion for analyzing leakage traces

At that stage, the BCDC has been used for estimating the SNR (Equation 2).

Then it has been experimentally shown that this criterion is also a figure of merit to gauge the quality of a set of EM traces with the intent to apply an SCA on it.

The BCDC can also be used to analyze the computational activity of a DUT over time. Indeed, it is quite possible to calculate BCDC values for the different time windows ([t i , • • • , t k ] with i ≥ 1 and k ≤ q) constituting the traces:

SN R(t i , • • • , t k ) = 1 BCDC 2 (t i , • • • , t k ) -1 (13)
and to deduce SN R values so that to distinguish windows characterized by different SNR values. This leads to draw the estimated SN R evolution over time, considering that each type of electrical activity (or instruction) consumes a specific amount of power and is thus characterized by a specific SNR value. Figure 3 gives 25 EM traces related to the ciphering of 25 random plaintexts (traces with no collision), the standard deviations of all the samples constituting the traces (thus computed with only 25 values), and finally the evolution of SN R estimated from 25 colliding traces using windows of 250 consecutive This result may come from two main reasons. The first reason could be the fact that the computation of BCDC values involves the calculus of standard deviations which are performed on a set of 250 consecutive samples rather than 25 values in case of the use of the vertical variance as analysis criterion. This limits the impact of noise. The second could also be related to the reduction of the impact of noise, and more particularly of the effect of outliers. Indeed, 25 traces allows computing the mean of the SNR distribution, SN R, with 300 different estimates of the SNR.

Figure 4 shows the same results than Figure 3 for 25 traces collected at Z = 3000µm i.e. for traces characterized by a lower SNR value. In this condition the AES is no more visible on the mean trace (in black Figure 4a). In addition, the transitions of a synchronization signal on an IO pad are much more visible; the EM probe is now above the related bonding.

If the results presented above, obtained on an hardware AES embedded on an FPGA, are interesting, the full effectiveness of the technique appears during the analysis of the electrical activity of a micro-controller which executes a software implementation of the AES. To illustrate this, Figure 5 gives 25 traces EM (and their mean) corresponding to the computation of the two last rounds of the AES by a 32-bits processor (cortex M4) designed in 90nm process technology and clocked at 40MHz. As shown, neither the mean nor the variance of traces allow to distinguish any pattern related to the two rounds. The vertical mean and variance look like a really noisy trace. This is not the case with the evolution of SN R on which two quite similar (but not exactly the same) patterns are visible. One can even distinguish similar shapes in the two rounds.

Collision Criterion for adaptive filtering

In the preceding paragraphs, the BCDC has been exploited to analyze the evolution of the electrical activity of two DUT, an FPGA and a micro-controller, using a sliding windowing solution. Such an approach can be adopted to distinguish in the frequency domain, the frequencies carrying information related to the signal from those that convey essentially noise. At that end, a frequency bandwidth, δf , is chosen. The colliding traces are then translated in the frequency domain using for example the FFT. This done, for each frequency f of the spectrum, the harmonics falling out of the frequencies bandwidth of interest (f ± δf ) are canceled before to translate back the traces in the time domain and compute SN R(f ± δ) at each frequency f of the spectrum, as in the preceding sections.

It should be noticed that δf could be set to the spectral resolution of the FFT which depends on the number of points constituting the traces. This approach with the minimal possible δf value has a larger computational cost and returns in computing SN R(f ) i.e. in keeping only one frequency in the bandwidth. Contrarily to σ N (f ± δf ), SN R(f ± δf ) seems to decrease following a 1 f function as indicated by the leakage model in the frequency domain introduced in [START_REF] Tiran | A model of the leakage in the frequency domain and its application to CPA and DPA[END_REF] and as formerly observed in [START_REF] Mateos | A new correlation frequency analysis of the side channel[END_REF]. However, it should be noticed that peaks of SN R(f ± δf ) appear at several multiples of the frequencies (50MHz). This is probably due to the influence of the clock tree circuit, one of the largest energy consumers at these frequencies in IC.

More importantly, it should be observed that SN R(f ± δf ) remains high (> 5) for f < 275MHz and particularly high for f < 52MHz. It therefore seems wise to only keep harmonics below 275MHz or even 52MHz during SCA.

In order to validate that SN R(f ± δf ) is an interesting figure of merit allowing to adapt at best the filtering steps usually applied prior any SCA, CPA were first launched on the rough EM traces, then on these same traces after canceling of harmonics higher than 52MHz and 275MHz, respectively. Figure 7 gives the evolutions of the gGE for these three CPA in the case of traces collected at Z = 0 and at Z = 4200µm.

Figure 7 shows that the filtering of traces characterized by a high SNR value (traces collected at Z = 0) has a really limited interest as it was expected. By contrast, for traces collected at Z = 4200µm, characterized by a SN R value lower than 3, keeping harmonics with the highest SN R(f ± δf ) values allows to enhance significantly the efficiency of the CPA. Indeed, keeping only the f values lower than 52MHz allows reaching a gGE lower than 20 after the processing of 3000 traces, value to be compared to 90 without application of any filtering procedure.

Conclusion

The SNR, as defined in engineering, is a key figure for gauging the quality of analogue measurements or communication channels. However, its use requires the knowledge of the signal to be measured. In the context of SCA, the nature and characteristics of the measured signals are not known. The use of the SNR is therefore particularly difficult. In this context the paper proposes a method to estimate the SNR, when the signal is not known. This method is based on a collision detection criterion. Its advantages are its easiness of use, its accuracy, and its low computation time. It has also been shown in this paper that the SNR is an interesting figure of merit to gauge the quality of SCA traces, to analyze the behavior of IC, or to guide rationally (and therefore in an efficient way) the filtering step commonly used as pre-processing of traces. There are also other applications of this fast SNR estimation method. Among them it use for the spatial analysis IC activity has given interesting results not shown in this paper by lack of place. However, we believe that the main interest of the proposed SNR estimation technique is that it is a sufficiently simple, low cost and convenient technique to establish the electrical engineering SNR as a standard measure for comparing the various experimental practices of SCA of both academic and industrial laboratories.

Fig. 2 .Fig. 3 .

 23 Fig. 2. (a) Evolutions with Z of σS and SN R, (b) Evolution of the gGE wrt the SN R (c) histogram of the √ SN R values estimated with traces collected at Z = 0, (d) histogram of the √ SN R values estimated with traces collected at Z = 2000µm

Fig. 4 .

 4 Fig. 4. (a) 25 EM traces collected @ Z = 3000µm over an AES (b) evolution of SN R

Fig. 5 .

 5 Fig. 5. (a) 25 EM traces (and their mean in black) collected above a 32bit microcontroller processing an AES (b) evolution σ ( mi) and (c) evolution of SN R

Fig. 6 .

 6 Fig. 6. SN R(f ) and σN (f ) computed with 8 EM traces

Figure 6

 6 Figure 6 gives the evolution with f of SN R(f ± δf ) estimated with only 8 EM traces collected above the AES mapped into the FPGA, considering a frequencies bandwidth δf = 2MHz. The quite flat evolution of σ N (f ± δf ) is also drawn in this figure.

Fig. 7 .

 7 Fig. 7. Evolution of gGE with the number of processed traces (collected either at Z=0 or at Z = 4200µm when keeping respectively all harmonics (in blue), harmonics below 275MHz (in red) and harmonics below 52MHz (in black).

Table 1 .

 1 Evolution of the mean SN R value wrt σN and the number n of measurements used to estimate it σ N


Denis Teissandier

Yann Ledoux

Santiago Arroyave-Tobon

Vincent Delos

Jean Marc Linares

Effect of form errors on the positioning precision of over-constrained systems

Keywords: Precision, Positioning, Geometric modelling

This paper analyses the influence of form deviations on over-constrained mechanical systems. A mechanical clamp was equipped with LVDT sensors to measure relative displacements between its two parts to derive their relative position. Each part was measured using CMM filtering its form deviations. Part feature variations (position, orientation and dimension) and assembly clearances were aggregated on a set of 6d linear constraints. From the measured local displacements another polytope is derived. These two polytopes are compared to quantify the influence of form error on the geometrical behaviour of the assembly. To do this evaluation new metrics are presented and discussed.

Introduction

The objective of a tolerance analysis is to aggregate the geometric error variations of the constitutive parts of a mechanical system in order to check the functional requirements. Such an analysis process is based on: a model of geometric error variations, a set of relations between them derived from the topological structure of the mechanical system and numerical tools to simulate its geometrical behavior. This tolerance analysis can be done according to worst case or statistical approaches [START_REF] Cao | A comprehensive review of tolerance analysis models[END_REF]. The major question is how to select the suited values of geometrical product specifications in order to control the geometric error variations (form, orientation, position and dimension) to guarantee the expected product quality. To simplify the tolerance analysis, some assumptions are usually made, among which form errors are neglected and mechanism parts are considered as infinitely rigid bodies. These assumptions could lead to inaccuracies in the computation results as shown by [START_REF] Grandjean | On the role of form defects in assemblies subject to local deformations and mechanical loads[END_REF]. Only a few research works, summarized in [START_REF] Xingyu | Review and Comparison of Form Error Simulation Methods for Computer-Aided Tolerancing[END_REF], address the effect of form error variations on the geometrical behavior of an assembly. This is a key point in mechanism analysis that should be quantified. Some works are based on experimental devices to compare the simulated and the real behaviors of the assembly. However, these works were dedicated to analyze, experimentally, isostatic assemblies.

In this work, a practical case study is used to demonstrate the influence of form errors in the positioning precision on over-constrained mechanical systems. Error metrics are proposed to quantify that. In the context of Industry 4.0, these metrics are computed using a digital twin of the mechanical system. Their uses in production control are discussed in the manufacturing field. In a first part of this paper, a geometrical model for over constrained mechanisms is set up assuming no error from variations. An example of a mechanical clamp with 2 parts and their several parallel joints is treated. In the second part, based on rigid body restrictions, the influence of the form error variations on the relative position between two parts is evaluated. This evaluation is based on the comparison of the calculated and the measured displacements due to geometrical errors. To do this evaluation four new metrics are presented and discussed.

Materials and methods

The parts were modeled as infinitely rigid bodies. The displacement restrictions induced by the non-interpenetration conditions between two surfaces potentially in contact are formalized by a set of linear constraints [START_REF] Fleming | Geometric relationships between toleranced features[END_REF][START_REF] Giordano | Clearance space in volumic dimensioning[END_REF][START_REF] Davidson | A new mathematical model for geometric tolerances as applied to round faces[END_REF]. Each set of contact constraints is represented by a polyhedron shown in the next section [START_REF] Arroyave Tobon | Tolerance analysis with polytopes in HV-description[END_REF].

Geometrical model for over-constrained systems

Let us consider a contact restriction between surfaces 1,j (surface j of part 1) and 2,j (surface j of part 2), see fig. 1. The general definition of the contact restriction is given by (1). This definition is applied to a set of Nk points (with 1≤k≤kmax and kmax is the number of facets for the discretization of each contact surface). At each point, the contact restriction is applied along the normal vector nk where dk is the distance between 1,j and 2,j. The intersection of this set of restrictions defines the polyhedron Pj. P 𝑗 = ⋂(𝐭 𝑁𝑘_𝟏,𝒋/𝟐,𝒋 ). 𝐧 𝒌 ≤ 𝑑 𝑘 ⟺ 𝑘 P 𝑗 = ⋂ (𝐭 𝑶_𝟏,𝒋/𝟐,𝒋 + 𝐍 𝐤 𝐎 × 𝐫 𝟏,𝒋/𝟐,𝒋 ). 𝐧 𝒌 ≤ 𝑑 𝑘 𝑘 [START_REF] Cao | A comprehensive review of tolerance analysis models[END_REF] Pj is a polyhedron of dimension 6 since each local contact restriction is a half-space k H -in the deviation space of dimension 6 [START_REF] Arroyave Tobon | Tolerance analysis with polytopes in HV-description[END_REF]. The 6 parameters (Rx,Ry,Rz,Tx,Ty,Tz) are coming from the components of r1,j/2,j (rotation vector) and tO˗1,j/2,1 (translation vector expressed at point O). The point O is any point assumed to be rigidly linked to the set of points Nk. Let us consider a mechanical assembly composed of two parts: part 1 and part 2, in multiple relations by several contacts in parallel between surface 1,j and 2,j (1≤ j ≤ jmax): see fig. 2. A datum R1 is associated to the surfaces 1,j of part 1 by measurements. Following the same method, a datum R2 is associated to the surfaces 2,j of part 2. The determination of the local distances dk from (1) takes into account the location deviations of 1,j with respect to R1 (d1,j/R1) and 2,j with respect to R2 (d2,j/R1). The relative position between R1 and R2 is defined by Pc. Pc is the intersection of the n polyhedron Pj (1≤ j ≤ jmax; where jmax is the number of contact surfaces). If the contacts between parts 1 and 2 suppress the 6 degrees of freedom, Pc is a polytope (i.e. a bounded polyhedron). Then the vertices of Pc define the extremal values of the parameters (Rx,Ry,Rz,Tx,Ty,Tz) in a deviation space. Pc is the result of the intersection between (jmax×kmax) local restrictions (i.e. half-spaces of dimension 6). Finally, we will use the general expression (2), commonly named Hdescription, to define the polytope Pc. According to the Minkowski-Weyl theorem [START_REF] Ziegler | Lectures on polytopes[END_REF], an equivalent definition (3) can be used for Pc where vci are the vertices of Pc. The definition ( 3) is commonly named the V-description of Pc, where Pc is defined by the convex hull of its vertices vci. We use this theorem in order to directly measure the relative position between parts 1 and 2 by a measured polytope Pm. A mechanical system was instrumented by sensors which measures a finite number of relative positions between parts. Each position corresponds to a set of parameters defined by a point devu in the deviation space. Finally, Pm can be estimated by the convex hull of points devu according to (4) where umax is the number of measures to generate the points devu.

P c = ⋂ P 𝑗 𝑗 = ⋂ H ̅ 𝑢 -with: 1 ≤ 𝑢 ≤ 𝑗max × 𝑘max 𝑢 (2)
P c = Conv(𝑣 𝑐𝑖 ) (3)
P 𝑚 = Conv(𝑑𝑒𝑣 𝑢 ) 1 ≤ 𝑢 ≤ 𝑢max (4)
Figure 2. Contact restrictions at the assembly level.

Experimental protocol

Our aim is to evaluate the influence of the form error variations on the relative position between two parts in contact in the case of an over constrained assembly. This evaluation is based on the comparison between the calculated polytope Pc and the measured polytope Pm. The surfaces 1,j and 2,j are modeled by substitute surfaces through the measurement process, in order to eliminate the form deviations by least square best-fit method (i.e. gaussian filters). So, the polytope Pc is the aggregation of orientation, position and dimension variations. The convex hull of points (devu) derives from direct measurements of the relative position between parts. The polytope Pm is then the aggregation of all the possible sources of geometrical variations of surfaces 1,j and 2,j (i.e. orientation, position, dimension and then form errors).

The study is carried out on a mechanical system inspired by a clamp (fig. 3a andb). The constitutive parts 1 and 2 are in multiple contacts by three couples of pin 1,j -hole 2,j (1≤ j ≤ 3) with a nominal diameter equal to Ø10. These couples are evenly distributed on a nominal circle (Ø80) (fig. 3c). Two additional datums 1,4 and 2,4 respectively on part 1 and part 2, are nominally located on the center of this circle. These respective center holes O1,4 and O2,4 with a fixed square on sensors part 1 are used to measure the relative position between datum R1 and datum R2 respectively linked to parts 1 and 2. We define contact restrictions assuming the contact zones between pins and holes lie in the same plane (i.e. the common plane between parts). As a consequence, the relative position between R1 and R2 is defined by a small rotation along a normal to the common plane and two orthogonal translations in this plane. The polytopes Pc and Pm are both 3d polytopes. The origin of datum R1 is defined by O1,1 and one axis by the square. The origin of datum R2 is defined by O2,1 and one axis by a line (O2,1, O2,4). The five steps of the analysis protocol are depicted in fig. 4. First, part 1 is measured to estimate the location of surfaces 1,j in R1 and their diameters ØD1,j (step 1). By analogy, surfaces 2,j in R2 and their diameters ØD2,j are estimated on the part 2. Each surface is measured by 8 points evenly distributed with a CMM machine (Renishaw® TP20, measurement uncertainty: 7 µm in volume). Although they are not used in the computation of polytope Pc, the form deviations of each contact surface are also evaluated. In step 2, part 1 is assembled on part 2 and fixed in an arbitrary location such as their common planes are in contact. Then, in step 3, the square and the hole 1,4 are measured in R2 and the LVDT sensors are initialized in this configuration. This is a major step in the best-fit of the origins of polytopes of Pc and Pm. In step 4, an operator manually moved part 1 with respect to part 2 reaching as many positions as possible. The motions between parts are measured by the three LVDT sensors (DP5 probes of Solartron Metrology® with accuracy of 0.5 µm). Finally, all the data are saved to be used for the computation of polytopes Pc and Pm (step 5).

Computed polytope Pc

The contact restrictions between pin 1,j and hole 2,j derive from definition (1), and Pj can be expressed as [START_REF] Giordano | Clearance space in volumic dimensioning[END_REF].

𝐏 𝒋 = ⋂ (𝐭 𝑶_𝟏,𝒋/𝟐,𝒋 + 𝐍 𝐤 𝐎 × 𝐫 𝟏,𝒋/𝟐,𝒋 ). 𝐧 𝒌 ≤ 𝑑 𝑘 𝒌 [START_REF] Giordano | Clearance space in volumic dimensioning[END_REF] with 𝑑 𝑘 = ( 𝐷 2,𝑗 -𝐷 1,𝑗 2 ) + 𝐎 1,𝑗 𝐎 2,𝑗 • 𝐧 𝑘

The couples Nk , nk are deduced from a mesh built on the CAD model of the clamp. This couple and the point O control the normal of each half-space k H -. An operand Pj is generated with kmax=36 such that the maximal deviation between the mesh and the surfaces of pin 1,j and hole 2,j is less than 1 µm [START_REF] Arroyave Tobon | Tolerance analysis with polytopes in HV-description[END_REF]. The other parameters to compute a local distance dk are coming from the diameters and the location deviations of pin 1,j and hole 2,j. Their values derive from the measurements performed in steps 1 and 2 of the synoptic (see fig. 4). Least square filtering is used to remove the form deviations to others for pins 1,j and holes 2,j. From (2) the polytope Pc can be computed as shown in fig. 5.

Measured polytope Pm

From local measures of the sensors, we obtain a set of relations using equation ( 6). These relations can be expressed at a common point O assumed to be rigidly linked with the points Mi. Then the parameters Rz, Tx and Ty can be determined from ( 6) with ( 7) and ( 8) by [START_REF] Ushakov | Optimization of the Hausdorff distance between convex polyhedrons in R3[END_REF].

{𝒕 𝑴𝒊_𝑹𝟏/𝑹𝟐 . 𝒔 𝒊 = 𝑚 𝑖 } 1 ≤ 𝑖 ≤ 3 (6)
𝒕 𝑴𝒊_𝑹𝟏/𝑹𝟐 = 𝒕 𝑶_𝑹𝟏/𝑹𝟐 + 𝒓 𝑹𝟏/𝑹𝟐 × 𝑴 𝒊 𝑶 (7) 𝒔 𝒊 = (𝑠 𝑖𝑥 , 𝑠 𝑖𝑦 , 0) 𝜪𝜧 𝒊 = (𝑂𝑀 𝑖𝑥 , 𝑂𝑀 𝑖𝑦 , 0) 𝒓 𝑹𝟏/𝑹𝟐 = (0,0, 𝑅 𝑧 ) 𝒕 𝑶_𝑹𝟏/𝑹𝟐 = (𝑇 𝑥 , 𝑇 𝑦 , 0)

{(-𝑠 𝑖𝑥 . 𝑂𝑀 𝑖𝑦 + 𝑠 𝑖𝑦 . 𝑂𝑀 𝑖𝑥 )𝑅 𝑧 + 𝑠 𝑖𝑥 . 𝑇 𝑥 + 𝑠 𝑖𝑦 . 𝑇 𝑦 = 𝑚 𝑖 } 1 ≤ 𝑖 ≤ 3(9)

Finally, from umax records of triplets mi we can deduce a cloud of points devu in a deviation space (Rz,Tx,Ty).

Measurement of R 1 such that part 1 is fixed on part 2 Step 1

Measurements of part 1 and part 2 individually

O 1,j in R 1 , ØD 1,j O 2,j in R 2 , ØD 2,j O 1,4 in R 2 square in
Step 2

Step 3

Step 4

Step 5

R z T x T y P 1 P 2 P 3

Results and discussions

Measured polytope Pm

Computed polytope Pc

The results of the CMM measurement of parts 1 and 2 are displayed in table 1. The vertices of Pc result from the intersections of its half-spaces initially generated on the nodes of a regular mesh, see fig. 5 and7.

Comparison between polytopes Pm and Pc

The difference between the shape of polytopes Pc and Pm is due to the form deviations. Based on equation ( 5), the distance dk* between a pin 1,j and a hole 2,j is redefined by (10).

𝑑 𝑘 * = 𝑑 𝑘 + 𝑓𝑣 𝑘 (10)
The parameter fvk is the accumulation of form error deviation of a pin 1,j with respect to a hole 2,j at a node Nk along the normal nk, see fig. 1. It was filtered by the best-fit process after the CMM measurement of parts 1 and 2 and it is not taken into account in [START_REF] Giordano | Clearance space in volumic dimensioning[END_REF]. In consequence, form error variations generate a translation fvk of the half-space k H -along its normal. In general, the parameters fvk have an influence on the boundary of the intersection of these half-spaces. Each vertex of a polytope defines an extremal position between the parts 1 and 2 and is created by a minimum of 3 concurrent hyperplanes in dimension 3. These hyperplanes are the boundaries of the half-spaces deriving from some points Nk which define the contact zone between the parts 1 and 2. The relation between the half-spaces and the vertices ensures the complete traceability between nodes Nk, nk and the extremal positions extracted from its faces [START_REF] Arroyave Tobon | Tolerance analysis with polytopes in HV-description[END_REF]. The comparison between polytopes Pm and Pc cannot be done based on properties of their lattices. Therefore, the direct estimation of the values fvk is not possible. In response to this, new geometrical metrics are proposed to compare Pm and Pc such as the bounding boxes, the mass centers, the volumes and the Hausdorff distances [START_REF] Ushakov | Optimization of the Hausdorff distance between convex polyhedrons in R3[END_REF]. The bounding box of Pm is included inside the bounding box of Pc. According to the parameters (Rz,Tx,Ty), the ratios of the three extends of the bounding boxes are 88%, 87% and 67%, see table 2. If we want to limit the global motions of one part compare to the other one, this bounding box computation is a useful quantification. It is also, a well indication of the positioning precision of the parts (small volume meaning a small motion). These results do not integrate the correlations between the parameters (Rz,Tx,Ty). Assuming that the relative position of parts is randomly distributed in the polytope, the vector GcGm then gives the average of the relative location between datum R1 and R2 in terms of (Rz,Tx,Ty) in the 3d deviation space, see table 2. The magnitude and the direction of this vector are strongly correlated to the distribution of fvk on the contact surfaces. Furthermore, the polytope Pm is not included inside the polytope Pc but the volume of Pc ∩ Pm is very closed to the volume of Pm, see table 4. These volumes indicate the inclusion rate and traduce if the polytope Pc increase or not the extremal displacements of Pm. These two-last metrics (GcGm and volumes) can help a designer to validate a design choice in a preliminary design stage.

Finally, the deviations from the tessellated boundary of Pm with respect to the boundary of Pc can be deduced by a set of distances, see fig.

Conclusion

A complete analysis protocol to evaluate the influence of the form error variations on the relative position between parts of an over constrained assembly was presented. This study proposes new metrics to detect and localize assembly interferences induced by form deviations. Moreover, it brings to the fore the influence of form errors in the positioning precision. The four proposed metrics were: Bounding boxes, Mass centers, Volumes of polytopes and Hausdorff distances.

Figure 1 .

 1 Figure 1. Local (a) and global (b) contact restrictions.

Figure 3 .

 3 Figure 3. Experimental device (a); CAD part (b); geometric setting of the clamp (c).

Figure 4 .

 4 Figure 4. Experimental protocol.

Figure 5 .

 5 Figure 5. Polytope Pc = P1 ∩ P2 ∩ P3.

Figure 6 .

 6 Figure 6. Cloud of 2000 points (a) and its convex hull Pm (b).

Fig. 6a gives

  Fig. 6a gives an example with 2000 records of triplets mi. Fig 6b shows the polytope Pm resulting from the convex hull of 2000 points. The vertices of Pm are the extremal displacements between parts 1 and 2 reached by the manual operation. The point density on the boundary of Pm is very heterogeneous.

Figure 7 .

 7 Figure 7. Projections of polytopes Pc and Pm in canonical views.

8 .

 8 Each distance is the Hausdorff distance between a point from the boundary of Pm and Pc. Each signed distance indicates if the form deviation must be added or subtracted along a given direction in the neighborhood of an extremal position between the parts 1 and 2. These distances and their directions can lead to detect what half-space k H -are on the boundary of Pm. Thus, it is possible to investigate about the evolution of the contact zones with and without form error variations in a geometrical model. This diagnosis can lead to detect precisely where are the influent zones on the relative position between parts 1 and 2. Furthermore, this diagnosis can be useful to analyze the load transfer between the parts in order to predict the mechanical behavior of the assembly. These new metrics were integrated in the open source software (i2m.u-bordeaux.fr/politopix).

Figure 8 .

 8 Figure 8. Hausdorff distances between tessellation of Pm with respect to the boundary of Pc.

Table 1 .

 1 Measurement results of individual parts 1 and 2 in R2 in mm

	Surface	1,1	1,2	1,3	2,1	2,2	2,3
	Form dev	0.014	0.011	0.012	0.007 0.009	0.007
	Diameter 10.004 10.008 10.006 10.554 10.558 10.557
	Loc x	-0.273	34.904	-34.422	0.000	34.663	-34.619
	Loc y	39.973	20.070	20.056	40.000	19.995	19.995

Table 2 .

 2 Bounding

	boxes of polytopes
	Bounding box	Rz (10 -3 rad)	Tx (mm)	Ty (mm)
	Pc	[-7.321; 6,384]	[-0.510; 0.023]	[-0.302; 0.207]
	Pm	[-6.190; 5.860]	[-0.456; 0.007]	[-0.217; 0.133]
	Ratio Pm / Pc	88%	87%	67%

Table 3 .

 3 Mass

	centres of polytopes
	Mass center	Rz (10 -3 rad)	Tx (mm)	Ty (mm)
	Gc (Pc)	-0.468	-0.247	-0.056
	Gm (Pm)	-0.165	-0.211	-0.065
	GcGm	0.303	0.026	-0.009

Table 4 .

 4 Volumes of polytopes in 10 -3 mm 2 .rad

		Pc	Pm	Pc ∩ Pm
	Volume	46.04	25.58	25.46

In future, the metric uncertainties will be evaluated in order to take into account the effect of assumptions (rigid body, small screw displacement, surface discretization,…) and measure uncertainties. These metrics can help designers to take decisions in geometrical product specification. Additionally, they could be used in Industry 4.0, for example in smart part pairing, for choosing the constitutive parts that warrantee an assembly with no interference. This could be done, in real time in the production chain, by digital twin assembly simulation. The knowledge of the assembly interference localization opens the way to automatic processes to repair high added value parts.