Diego F Aranha
email: dfaranha@ic.unicamp.br

Pierre-Alain Fouque
email: fouque@irisa.fr

Benoît Gérard
email: benoit.gerard@irisa.fr

Jean-Gabriel Kammerer
email: jean-gabriel.kammerer@m4x.org

Mehdi Tibouchi
email: tibouchi.mehdi@lab.ntt.co.jp

Jean-Christophe Zapalowicz
email: jean-christophe.zapalowicz@inria.fr

GLV/GLS Decomposition, Power Analysis, and Attacks on ECDSA Signatures With Single-Bit Nonce Bias

Keywords: Elliptic Curve Cryptography, GLV/GLS Method, Bleichenbacher's ECDSA Attacks, Side-Channel Analysis

The fastest implementations of elliptic curve cryptography in recent years have been achieved on curves endowed with nontrivial efficient endomorphisms, using techniques due to Gallant-Lambert-Vanstone (GLV) and Galbraith-Lin-Scott (GLS). In such implementations, a scalar multiplication [k]P is computed as a double multiplication [k1]P + [k2]ψ(P ), for ψ an efficient endomorphism and k1, k2 appropriate half-size scalars. To compute a random scalar multiplication, one can either select the scalars k1, k2 at random, hoping that the resulting k = k1 + k2λ is close to uniform, or pick a uniform k instead and decompose it as k1 + k2λ afterwards. The main goal of this paper is to discuss security issues that may arise using either approach. When k1 and k2 are chosen uniformly at random in [0, √ n), n = ord(P ), we provide a security proofs under mild assumptions. However, if they are chosen as random integers of 1 2 log 2 n bits, the resulting k is slightly skewed, and hence not suitable for use in schemes like ECDSA. Indeed, for GLS curves, we show that this results in a bias of up to 1 bit on a suitable multiple of k mod n, and that this bias is practically exploitable: while lattice-based attacks cannot exploit a single bit of bias, we demonstrate that an earlier attack strategy by Bleichenbacher makes it possible. In doing so, we set a record by carrying out the first ECDSA full key recovery using a single bit of bias. On the other hand, computing k1 and k2 by decomposing a uniformly random k ∈ [0, n) avoids any statistical bias, but the decomposition algorithm may leak side-channel information. Early proposed algorithms relied on lattice reduction and exhibited a significant amount of timing channel leakage. More recently, constant-time approaches have also been proposed, but we show that they are amenable to power analysis: we describe a template attack that can be combined with classical lattice-based attacks on ECDSA to achieve full key recovery on physiscal devices.

Introduction

The GLV/GLS techniques. Many record implementations of elliptic curve cryptography in software, including, most recently, works such as [START_REF] Oliveira | Two is the fastest prime: lambda coordinates for binary elliptic curves[END_REF]5,[START_REF] Costello | Faster Compact Diffie-Hellman: Endomorphisms on the x-line[END_REF], rely on elliptic curves endowed with fast endomorphisms, as constructed by the methods due to Gallant-Lambert-Vanstone (GLV) [START_REF] Gallant | Faster Point Multiplication on Elliptic Curves with Efficient Endomorphisms[END_REF], Galbraith-Lin-Scott (GLS) [START_REF] Galbraith | Endomorphisms for faster elliptic curve cryptography on a large class of curves[END_REF], and generalizations thereof. In such implementations, the fast endomorphism ψ on the elliptic curve E/F q is used to speed up full size scalar multiplications [k]P by computing them as multi-exponentiation [k 1 ]P + [k 2 ]ψ(P ), where k 1 and k 2 are roughly half of the size of k. Indeed, on a prime order subgroup of E(F q ), ψ acts by multiplication by some constant λ, and thus, for a generator P of that subgroup, we have [k 1 ]P + [k 2 ]ψ(P ) = [k 1 + k 2 λ]P .

In order to compute random scalar multiplications with those techniques, two types of approaches have been considered, as far back as in the earliest presentations of the GLV method (such as Gallant's talk at ECC'99 [START_REF] Gallant | Efficient multiplication on curves having an endomorphism of norm 1[END_REF]).

On the one hand, k 1 and k 2 can simply be chosen uniformly at random in a suitable half-length interval. This appraoch, which we call the recomposition technique (since k is "recomposed" as k = k 1 + k 2 λ), results in a very simple implementation, and has been used in several implementation records including [START_REF] Oliveira | Two is the fastest prime: lambda coordinates for binary elliptic curves[END_REF], but Gallant expressed concerns about possible biases in the resulting scalar k. Such concerns have been partially vindicated by some numerical evidence provided by Brumley and Nyberg [START_REF] Brumley | On Modular Decomposition of Integers[END_REF], who also described a relatively general way to choose intervals for k 1 and k 2 so that the resulting choice of k is in fact secure (in the sense that it has high entropy). However, the Brumley-Nyberg method is a bit cumbersome, and no attack so far has been demonstrated against arbitrary half-length uniform choices of k 1 and k 2 , so that the security picture is somewhat unclear.

On the other hand, one can also pick k at random and subsequently deduce half-length values k 1 and k 2 , which eliminates concerns regarding possible biases in the distribution of k. This decomposition technique usually relies on lattice reduction in dimension 2 (or equivalently, continued fractions, a generalized Euclidean algorithm, etc.), as originally described in the GLV paper [START_REF] Gallant | Faster Point Multiplication on Elliptic Curves with Efficient Endomorphisms[END_REF], and is significantly more computationally demanding than recomposition. Simplifications of this method have later been proposed (particularly in [START_REF] Park | An Alternate Decomposition of an Integer for Faster Point Multiplication on Certain Elliptic Curves[END_REF]), as well as higher-dimensional generalizations [START_REF] Nguyen | Low-Dimensional Lattice Basis Reduction Revisited[END_REF] to tackle decompositions involving several endomorphisms (as recently used in [START_REF] Smith | Families of Fast Elliptic Curves from Q-curves[END_REF][START_REF] Guillevic | Four-Dimensional GLV via the Weil Restriction[END_REF] for instance).

ECDSA attacks. The success of GLV/GLS method in implementations lately makes it desirable to reconsider these decomposition and recomposition techniques from a security viewpoint. We do so in this paper in the context of ECDSA signatures, one of the most widely deployed elliptic curve cryptographic schemes, and an interesting target for the cryptanalyst (like other Schnorr-like signature schemes) due to its sensitivity to biases in the distribution of nonce values, as demonstrated by the powerful attack due to Howgrave-Graham and Smart [START_REF] Howgrave-Graham | Lattice Attacks on Digital Signature Schemes[END_REF] based on lattice reduction techniques, which breaks (EC)DSA when some of the most significant bits of the nonces are known. This attack was analyzed in further details by Nguyen and Shparlinski [START_REF] Nguyen | The Insecurity of the Digital Signature Algorithm with Partially Known Nonces[END_REF][START_REF] Nguyen | The Insecurity of the Elliptic Curve Digital Signature Algorithm with Partially Known Nonces[END_REF] and carried out in practice in many contexts, including against physical devices (see e.g. [START_REF] Naccache | Experimenting with Faults, Lattices and the DSA[END_REF][START_REF] Brumley | Cache-Timing Template Attacks[END_REF] for some examples). The basic idea is to express the key recovery problem as an instance of the Hidden Number problem (HNP), which reduces to the closest vector problem (CVP) in a suitable lattice. Since CVP is tractable in low-dimensional lattices, many practical instances of ECDSA can be broken depending on key size and the number of leaked nonce bits. The largest problem instance broken so far is the case of 2-bit nonce leaks on 160-bit curves, tackled by Liu and Nguyen [START_REF] Liu | Solving BDD by Enumeration: An Update[END_REF] using the most advanced known techniques for lattice reduction (BKZ 2.0 [START_REF] Chen | BKZ 2.0: Better Lattice Security Estimates[END_REF]). Breaking 2-bit leaks on 256-bit curves, or 4-bit leaks on 384-bit curves seems currently out of reach (see the discussions in [START_REF] Chen | BKZ 2.0: Better Lattice Security Estimates[END_REF][START_REF] Mulder | Using Bleichenbacher's solution to the hidden number problem to attack nonce leaks in 384-bit ECDSA: extended version[END_REF]).

In any case, there is a hard limit to what can be achieved using lattice reduction: due to the underlying structure of the HNP lattice, it is impossible to attack (EC)DSA using a single-bit nonce leak with lattice reduction. In that case, the "hidden lattice point" corresponding to the HNP solution will not be the closest vector even under the Gaussian heuristic (see [START_REF] Nguyen | Lattice-Based Fault Attacks on Signatures[END_REF]), so that lattice techniques cannot work. To break this "lattice barrier", the only known alternate attack is an algorithm due to Bleichenbacher [START_REF] Bleichenbacher | On the generation of one-time keys in DL signature schemes[END_REF] which predates the attack of Howgrave-Graham and Smart, but was generally considered of mostly theoretical interest until it was recently revisited by De Mulder et al. [START_REF] Mulder | Using Bleichenbacher's solution to the hidden number problem to attack nonce leaks in 384-bit ECDSA: extended version[END_REF] to attack 384-bit curves. Bleichenbacher devised his attack to demonstrate a vulnerability in DSS at the time, in which DSA nonces were generated by picking a random value of n bits, where n is the bit length of the group order n, and then to reduce it modulo n. Bleichenbacher showed that the resulting bias could be exploited in a very interesting way, obtaining a key recovery using about 2 41 signatures and about 2 47 time and 2 41 memory complexities. At that time, it was not possible to mount this attack and only simulations on reduced numbers were possible and the paper was never published.

In the first stage, Bleichenbacher's algorithm reduces the signatures from 160 bits to say 40 bits using linear combinations of the original signatures and then, during a second phase, a Discrete Fourier Transform is used to recover the most significant bits of the secret key. The bias of the reduced signatures is higher than the bias of the original signatures, that's the reason why Fourier technique is needed to extract this information. This algorithm is very similar to Blum, Kalai and Wasserman algorithms [START_REF] Blum | Noise-tolerant learning, the parity problem, and the statistical query model[END_REF][START_REF] Levieil | An Improved LPN Algorithm[END_REF] for solving LPN and LWE problems. For 384-bit order, the first stage of Bleichenbacher original attack is not sufficiently efficient to reduce the signatures and more advanced techniques based on LLL and BKZ are needed if the number of leaked bits is high enough [START_REF] Mulder | Using Bleichenbacher's solution to the hidden number problem to attack nonce leaks in 384-bit ECDSA: extended version[END_REF]. The modification of the first stage is not possible if less than one bit of nonces is available and we turn back to Bleichenbacher's original attack which requires a high number of signatures. out this attack on the standardized SECG P160 R1 elliptic curve. On this 160bit curve, we use 2 33 ECDSA signatures, and achieve a full key recovery in a few hours of wall-clock time on a 64-core workstation. The most time-consuming part of the attack is the first phase, in which a sorting algorithm is executed several times. This is the first key recovery from a single bit of bias, which paves the way to new applications. We stress again that this record cannot be achieved using lattice reduction techniques based on HNP problem, since even if the HNP lattice satisfies the Gaussian heuristic, a condition for finding the hidden lattice point is that the number of known bits of the nonce must be greater than log 2 ( πe/2) ≈ 1.0471 (hence at least 2) [START_REF] Nguyen | Lattice-Based Fault Attacks on Signatures[END_REF], irrespective of the underlying lattice reduction algorithm.

As a second contribution, we show a security proof for the recomposition method on curves obtained by the quadratic GLS method once the values k 1 and k 2 are uniformly distributed in the interval [0, √ n), where n is the prime group order. We prove that the statistical distance between this distribution and the uniform distribution in [0, n) is negligible. Furthermore, if k 1 and k 2 are taken at random in a small interval of the form [0, 2 m ), where m = 1 2 log 2 n , the bias on the distribution on k used in Bleichenbacher's attack is negligible. However, we show that the bias of the distribution on tk where t is the trace is sufficiently large and a Bleichenbacher's attack allows to recover the secret key. We also implement this attack and the complexities are similar to the previous part.

Finally, we study the decomposition technique proposed in GLV with the implementation described by Park et al. in [START_REF] Park | An Alternate Decomposition of an Integer for Faster Point Multiplication on Certain Elliptic Curves[END_REF]. To this end, we propose a very efficient side-channel attack that uses the leakage on the multiplication in order to recover some of the least significant bits of the nonces. Consequently, we can thus use lattice techniques to recover the secret key.

Preliminaries

Bias definition and properties

The measurement of the bias of random variables represents a significant part of our analyses. We thus recall the definition of the bias which was proposed by Bleichenbacher in [START_REF] Bleichenbacher | On the generation of one-time keys in DL signature schemes[END_REF].

Definition 1. Let X be a random variable over Z/nZ. The bias B n (X) is de- fined as B n (X) = E(e 2πiX/n ) = B n (X mod n),
where E(X) represents the mean. Similarly, the sampled bias of a set of points

V = (v 1 , • • • , v L ) in Z/nZ is defined by B n (V ) = 1 L L-1 j=0 e 2πivj /n .
The bias as defined above presents some useful properties we recall in Lemma 1.

Lemma 1. Let 0 < T ≤ n be a bound and X, Y random variables uniformly distributed on the interval

[0, T -1]. (a) If X is uniformly distributed on the interval [0, n -1], then B n (X) = 0. (b) If X and Y are independent, then B n (X + Y ) = B n (X)B n (Y ). (c) B n (-X) = B n (X)
where a denotes the conjugate of a.

(d) B n (X) = 1 T sin(πT /n) sin(π/n) and B n (X) is real-valued with 0 ≤ B n (X) ≤ 1.
(e) Let a be an integer with |a|T ≤ n and Y = aX, then B n (Y ) =

2: k $ ← [0, n -1] 3: (u, v) ← [k]P 4:
r ← u mod n; if r = 0 then goto step 2; 5:

s ← k -1 (H(m) + rx) mod n; if s = 0 then goto step 2; 6:

return (r, s) 7: end function

Bleichenbacher's Attack on single bit bias

In this part, we present our results on an ECDSA signature generation scheme where the nonce k is 1-bit biased. We demonstrate that an attack proposed some years ago by Bleichenbacher can succeed in retrieving the secret key in about 2 37 time and 2 33 memory complexities given 2 33 signatures, for 160-bit order. This attack was initially focusing on the DSA signature generation scheme but can be applied without any modification to ECDSA we consider in this paper.

The main idea consists in using the fact that the nonces k j are chosen from a biased random variable K, i.e. k are not randomly and uniformly generated on [0, n -1]. Because the values k j are biased and linked with the secret key x by the equations which are used for the signature computations, these signatures, correctly manipulated, also present a bias which will only be significant for the correct value of x. In other words the bias plays the role of the distinguisher in this attack.

Obviously, for cryptographic sizes, evaluating the bias for all values in [0, n-1] is impractical. However, Bleichenbacher observed that it is possible to "broaden the peak" of the bias in such a way that, with a value close the correct value of x, the bias will remain significant. Thus the bias computations can be performed on a more sparse set of candidates thanks to the Fast Fourier Transform. In return, it requires a non-negligible work on the signatures which reduces the bias, and the attack returns an approximation of the secret key, i.e. its most significant bits. The attack can be iterated to retrieve more bits of the secrets and as soon as sufficiently many bits of x are known, Pollard's lambda method [START_REF] Pollard | Kangaroos, monopoly and discrete logarithms[END_REF] can be used to derive the remaining bits. Algorithm 2 presents the main steps of the attack.

Algorithm 2 Bleichenbacher's attack given S ECDSA signatures. The parameters S, and ι have to be chosen accordingly to the bias.

Require: S biased ECDSA signatures (rj, sj) computed using a single secret key x. Ensure: The most significant bits of x.

1: Preprocessing

2: for j = 0 to S -1 do 3:
hj

← H(mj) • s -1 mod n 4: cj ← rj • s -1 j
mod n 5: end for 6: Reduction of the cj values (Sort-and-Difference Algorithm)

7: A ← [(cj, hj)] 0≤j≤S-1 8: for i = 1 to ι do 9:
Sort A by the cj values cj ≤ cj+1 10:

for j = 0 to S -ι do 11:

A[j] ← A[j + 1] -A[j] A[j] = (cj+1 -cj, hj+1 -hj) 12:
end for 13: end for 14: Only keep the pairs (cj, hj) such that cj < 2 15: Denote by L the number of such pairs 16: Bias computation using the inverse FFT 17: Z ← (0, • • • , 0) a vector of size 2 18: for j = 0 to L -1 do 19:

Zc j ← Zc j + e 2πih j /n 20: end for 21: W ← iFFT(Z)
Inverse FFT computation. The output is also a vector of complex numbers. 22: Find the value m such that |Zm| is maximal 23: return msb (mn/2 )

Attack analysis

We first explain why the bias can serve as a distinguisher and while doing so explain the goal of the preprocessing phase, as it was done in [START_REF] Mulder | Using Bleichenbacher's solution to the hidden number problem to attack nonce leaks in 384-bit ECDSA: extended version[END_REF], for the sake of completeness. For that purpose, consider S ECDSA signatures (r j , s j ) with biased nonces k j . We have the following relation due to step 5 of Algorithm 1:

k j = H(m j )s -1 j + r j s -1 j x mod n for 0 ≤ j ≤ S -1.
Now let h j = H(m j )s -1 j mod n and c j = r j s -1 j mod n. Then the set {h j + c j x} S-1 j=0 = {k j } S-1 j=0 will show a significant nonzero sampled bias. Moreover, for any w = x, the sampled bias from V w = {h j + c j w} S-1 j=0 will be relatively small. Since h j and c j are publicly computable, we thus have a way to determine the correct value of x by testing all the value w ∈ [0, n -1].

To have a practical test, we have to broaden the peak of the bias such that values of w close to the correct value x will also show a significant bias. The peak will be broad if the c j are relatively small. More precisely, by denoting 2 a bound such that 0 ≤ c j < 2 , then we can find an approximation of x by evaluating the sampled bias of 2 evenly-spaced values of w between 0 and n -1.

The reduction of the c j , second phase in Algorithm 2, can be done using a sort-and-difference algorithm. From S pairs (c j , h j ), we first sort them according to their first element. Then we subtract each c j from the next largest one and we take the differences of the corresponding h j as well. We thus obtain a list of S -1 pairs (c j , h j ) whose values c j are on average log(S) bits smaller. More details about the analysis of this reduction are given later. This reduction algorithm can be repeated in order to achieve the bound 2 . Now let w m = mn/2 , with m ∈ [0, 2 -1], be 2 evenly-spaced values between 0 and n -1. For sake of clarity, we keep the notation (c j , h j ) for the reduced pairs with c j < 2 and we consider having L such pairs. Then

B n (V wm ) = 1 L L-1 j=0 e 2πi(hj +cj mn/2 )/n = 2 -1 t=0 1 L {j|cj =t} e 2πihj /n e 2πitm/2 = 2 -1 t=0 Z t e 2πitm/2
with Z t = 1 L {j|cj =t} e 2πihj /n . B n (V wm ) can be viewed as the inverse Fast Fourier Transform of the vector Z = (Z 0 , • • • , Z 2 -1 ). Thus the multiple bias computations can be performed very efficiently using the FFT. From Step 17 to 20 in Algorithm 2, we compute this vector Z. Step 21 outputs a vector of the sampled bias for the 2 candidates, i.e.

iFFT(Z) = (B n (V w0 ), B n (V w1 ), • • • , B n (V w 2 -1 )).
Finally, the value of w m = mn/2 with the largest sampled bias should share its most significant bits with the secret key x.

Choosing the parameters. We first give some properties which will help to define the parameters for the attack. We can estimate the sampled bias for a wrong candidate w m , i.e. a value w m which do not share some most significant bits with the secret key x. More precisely, it can be shown that for w m either significantly larger or significantly smaller than x, we have B n (V wm ) ≈ 1 √ L , which corresponds to the average distance from the origin for a random walk on the complex plane.

The second property concerns the c j reduction phase and gives a relation between the number of signatures S and the number of reduced pairs L.

Proposition 1. Consider S ECDSA signatures of the form (c j , h j ) and γ ∈ Z. The percentage of signatures (c j , h j ) after the first application of the sort-anddifference algorithm such that c j < 2 log q-log S+γ can be approximated by 1-e -2 γ . Lemma 2. Let X 1 , . . . , X N be N independent uniformly distributed random variables over [0, 1], and for all i, denote by X (i) the i-th order statistic of the X j 's (namely, X (i) is the i-th smallest among the X j 's). Then, the random variables Y i = X (i+1) -X (i) for i = 1, . . . , N -1 are identically distributed, and all follow the beta distribution B(1, N ), of probability density function (hereafter pdf )

f (t) = N • (1 -t) N -1 . As a result, for any constant α > 0, we have Pr[Y i ≤ α/N ] = 1 -e -α + O(1/N ).
Proof. Indeed, a standard formula [11, 2.2.1] expresses the joint pdf of X (i) and X (i+1) as:

f i,i+1 (u, v) = N ! (i-1)!(N -i-1)! u i-1 (1 -v) N -i-1 for 0 ≤ u ≤ v ≤ 1, 0 otherwise.
Hence, the pdf f i of Y i is given by:

f i (t) = 1-t 0 f i,i+1 (u, u + t)dt for t ∈ [0, 1].
The change of variable u = (1 -t)w gives:

f i (t) = (1 -t) 1 0 f i,i+1 (1 -t)w, (1 -t)w + t dw = c(1 -t) 1 0 (1 -t) i-1 w i-1 (1 -w -t + wt) N -i-1 dw = c(1 -t) i 1 0 w i-1 (1 -t) N -i-1 (1 -w) N -i-1 dw = c(1 -t) N -1 1 0 w i-1 (1 -w) N -i-1 dw, where c = N ! (i-1)!(N -i-1)! .
In particular, we have f i (t) = c (1 -t) N -1 for some constant c and all t ∈ [0, 1], and since

1 0 f i = 1, we must have f i (t) = N (1 - t) N -1 = f (t) as required.
As a result, we obtain:

Pr Y i ≤ α N = α/N 0 N (1 -t) N -1 dt = 1 -1 - α N N = 1 -exp N • (-α/N + O(1/N 2 )) = 1 -e -α + O(1/N ).
This concludes the proof.

As an example consider a modulus n of size 160. Starting from 2 40 ECDSA signatures, after one iteration of the sort-and-difference algorithm, about 86.5% of them will have a value c j < 2 121 . The percentage drops to 22.1% if we consider only those ones with a value c j < 2 118 . Note that this proposition is only true for the first iteration of the algorithm where we really can consider variables as uniformly random and independently distributed. Clearly they are not after this: if after the first round variables were uniformly distributed, the ratio between γ = -2 and γ = 1 would be 0.125 = 1/2 3 where it is ≈ 0.255. Sadly, it appears that the ratio progress in our disfavor when we want to iterate, i.e. the ratio after ι iterations is less than (1 -e -2 γ ) ι . We thus do not have a lower bound. However the ratio can be experimentally determined and Table 1 gives an overview for different values of γ up to 6 iterations.

γ -2 -1 0 1 2 1 st iteration 0.22
0.39 0.63 0.86 0.98 2 nd iteration 0.031 0.12 0.36 0.75 0.94 3 rd iteration 3.2 10 -3 0.025 0.17 0.64 0.89 4 th iteration 3.0 10 -4 4.6 10 -3 0.069 0.53 0.84 5 th iteration 2.0 10 -5 6.7 10 -4 0.022 0.40 0.79 6 th iteration 2.8 10 -6 9.5 10 -5 6.5 10 -3 0.28 0.73 Table 1. Experimental ratio between the ECDSA signatures of the form (c j , h j ) such that c j < 2 log n-ι•(log S+γ) , and the S initial signatures, after ι iterations of the sortand-difference algorithm.

Given S signatures, we have to choose a pair (γ, ι) such that log n-ι•(log S + γ) = is sufficiently small to perform a FFT in 2 log time and 2 memory complexities. The algorithm complexity is O S log(S) + log( ) . Now a verification is necessary to be sure that this set of parameters will give a successful attack. Indeed denote by B n (K) the initial bias which is fully determined by the number of most (or least) significant bits of the k i which are known or set to zero (see Table 2 for some values). From properties (b) and (c) of the Lemma 1, each iteration of the sort-and-difference algorithm reduces the bias by raising it to the square of its norm (assuming that the variables are independant): indeed, let X, Y be uniformly distributed and independent random variables on [0

, n -1], then B n (X) = B n (Y ) and B n (X -Y ) = B n (X)B n (Y ) = |B n (X)| 2 . The final bias is then approximated by |B n (K)| 2 ι . Thus the following inequality holds since B n (V wm ) ≈ 1/ √ L: |B n (K)| 2 ι 1/ √ L,
where L represents as before the number of reduced pairs (c j , h j ) with c j < 2 .

Using Table 1 which gives the ratio L/S for different choices of pairs (γ, ι), we obtain a relation between S, ι, and n.

Note that contrary to previous reports in the literature [START_REF] Mulder | Using Bleichenbacher's solution to the hidden number problem to attack nonce leaks in 384-bit ECDSA: extended version[END_REF][START_REF] Bleichenbacher | On the generation of one-time keys in DL signature schemes[END_REF], we do not need to center the k j around 0. Indeed sort-and-difference algorithm performs only subtractions and does not mix subtractions and additions as is common with lattice reduction or generalized birthday algorithms. 2. Some values of bias for large n, when b most (or least) significant bits of k are known, using Property (d) of Lemma 1.

Implementation

We successfully implemented the attack. As our target, we chose the SECG P160 R1 curve, published in 2000 by the SECG consortium [START_REF]Standards for efficient cryptography, SEC 1: Elliptic curve cryptography[END_REF] and still considered secure. We fixed the most significant bit in the nonces and checked (with the help of the secret) that we indeed got the expected bias: ≈ 0.63662. Our C++ implementation was based on the RELIC toolkit (using its provided plain C integer arithmetic) [START_REF] Aranha | RELIC is an Efficient LIbrary for Cryptography[END_REF] and FFTW [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. We parallelized it in a straightforward manner (including (quick)sorting phases) and tested it on a multicore machine.

We generated 2 33 signatures and performed 4 sort-and-difference reduction phases. 450 millions (which is 52.5%) of our initial 2 33 signatures had their c j reduced down to 32 bits, as was expected from table 1. The bias after 4 reduction steps was 0.000743558 which is slightly greater than the expected 0.63662 2 4 ≈ 0.00072792. We then computed a FFT on 32 bits (we selected the reduced c j smaller than 2 32 ). The best candidate had a score approximately 35% greater than the second. Both corresponding MSB of the secret differed only by the 31 st and 32 nd most significant bits. The 3 rd and 4 th candidates were also very close to the two first ones, with score approximately 1/3 of the best candidate. Then, there was a number of random values with maximal score approximately 1/6 of the best one. We repeated the experiment several (5) times and got similar results, always finding at least the 30 MSB of the secret with the best candidate. We couldn't repeat it more because of the high computational resources involved.

The total memory used by the signatures and FFT tables was slightly more than 1 terabyte. To recover 32 bits of the private key, the attack took approximately 1150 CPU-hours, most of it being data exchange, which we can decompose as follows:

-70%: parallelised quicksort (the most memory-intensive phase) -18%: signature generation (approximately 250 to 430 kilocycles per signature depending on the CPU, excluding hash computations) -10%: candidate selection and FFT table preparation -≈ 1%: the FFT itself.

We did not use more parallelizable sorts like Batcher odd-even mergesort [START_REF] Batcher | Algorithm 5 Information propagation for one step of the multiplication inner loop[END_REF] but this would clearly be the next thing to do from a performance perspective.

Next steps of the attack to recover the following bits of the secret were done as in [START_REF] Mulder | Using Bleichenbacher's solution to the hidden number problem to attack nonce leaks in 384-bit ECDSA: extended version[END_REF]. Basically, it amounts to a replay of algorithm 2 on the initial signatures, putting the previously found MSB of the secret into h j . Write the private key x as x 0 2 m + x 1 where x 0 is the recovered m MSB at the first round. Then (h j ) + (c j )x = (h j + c j x 0 2 m ) + (c j x 1 ) and we want to recover the MSB of x 1 . We proceed as in the first round, except that we now keep the c j that are smaller than 2 +m instead of 2 (thus when = m we just have to stop the reduction one iteration earlier). Then we build the FFT table as Z[c j /2 m ] = Z[c j /2 m ] + e 2πih j /n . The FFT recovers the next most significant chunk of bits of the secret key. The computation restart makes it necessary to go back from the initial signatures, but there's no need to keep them in memory during the reduction. In practice we had barely enough memory to keep them, but in order to reduce memory usage they should either be stored on disk and retrieved to iterate the secret recovery, or tracked down through the reduction and rebuilt afterwards.

In practice, it is advisable to take a small security margin and reinject only 30 bits of the computed MSB of the secret to account for small variations of sampling around the peak. In any case, if we recurse with a wrong secret, the FFT will not detect any peak. Experiments indeed showed no peak in this case, with the highest score not being statistically different from the other ones. This paves the way for a time/memory tradeoff: suppose the hardware is limited in memory and can only work on (say) 2 31 signatures and 2 30 FFT size instead of the 2 33 needed for attacking 160 bits with 4 iterations with the previous algorithm. We first reduce the c j from 160 to 40 bits with 4 reductions as usual. We then simply guess the 10 MSB of the secret and build 2 30 -sized FFT tables accordingly. The guess will be correct on the only one FFT among the 2 10 which shows a significant peak. Since FFTs are particularly efficient, much more than sorting, this is of practical importance. Alternatively, if it's possible to compute 2 41 signatures, we can select only the expected 1/2 10 fraction of signatures whose corresponding c j have their 10 MSB already zeroes, that is to say that have 150 bits instead of 160 and can be reduced to 30 with 4 iterations. Finally, since the FFT table takes less memory than the signatures (a complex number occupies 16 bytes whereas a signature requires at least 40), we could improve the attack further by either carrying out several FFTs in parallel when guessing some bits of the secret, or by increasing the size of the FFT table slightly (with a corresponding increase of the selection bound on c j ). This would have two advantadges. Firstly, it would improve the sampling around the peak and reduce the uncertainty. Secondly, the bound increase implies that some signatures would be selected after the third round of reduction instead of the fourth, thus having a much better bias and hopefully revealing more precise information about the secret.

Our experiments targeted a 160-bit curve, but it should be pointed out that larger curves are susceptible to this attack as well. Roughly speaking, one can carry out the key recovery attack with 1-bit nonce bias on an N -bit curve in time ≈ 2 N/5 log 2 (N/5) and memory ≈ 2 N/5 . For example, a 256-bit curve can be attacked in time ≈ 2 58 and memory ≈ 2 52 : generate 2 52 signatures, perform 4 reduction steps (removing 4 • 51 = 204 bits on approximately 86% of the data), keep signatures with c j less than 2 52 and carry out the FFT on a table of size 2 52 . One signature is 64 = 2 6 bytes, so that the total memory needed for the attack is 2 18 terabytes of storage, which corresponds to 65536 of today's 4 TB disks. This does not appear to be out of reach of well-funded adversaries.

Security analysis of the recomposition technique

The results presented so far had no direct connection with GLV/GLS curves. We now turn to such curves, and first discuss in this section the security of what we called the "recomposition technique" for GLV/GLS coefficients (namely, choose k 1 and k 2 uniformly at random in some interval [0, K) to obtain k = k 1 + k 2 λ mod n), whereas the next section will focus on the "decomposition technique".

To fix ideas, we consider an elliptic curve E obtained by the quadratic GLS method over a prime field [13, §2.1]. In other words, there is an elliptic curve E 0 over the prime field F p such that E is the quadratic twist of E 0 over F p 2 . If we denote by p + 1 -t the order of E 0 (F p ) (where t is bounded as |t| ≤ 2 √ p by the Hasse-Weil theorem), the order n of E(F p 2 ) satisfies:

n = (p -1) 2 + t 2 . ( 1
)
We assume that this order n is prime, which is the main case of interest. Then, E is endowed with an efficient endomorphism ψ (obtained by conjugating the Frobenius map with the twisting isomorphism) which acts on the cyclic group E(F p 2 ) by multiplication by

λ ≡ t -1 (p -1) (mod n). (2)
In particular,

λ 2 ≡ (p -1) 2 /t 2 ≡ -t 2 /t 2 ≡ -1 (mod n).
In this setting, we first prove in §4.1 that if k 1 and k 2 are chosen uniformly at random in [0, √ n), then k = k 1 + k 2 λ is statistically close to uniform in Z/nZ, so that such a choice of (k 1 , k 2 ) can be used securely in any cryptographic protocol (and in particular ECDSA). On the other hand, we show in §4.2 that if k 1 and k 2 are chosen in [0, 2 m ) where m = 1 2 log 2 n instead, then k = k 1 + k 2 λ may not be close to uniform anymore, and we show that a variant of Bleichenbacher's attack can apply. In §4.3, we describe an implementation of that attack on a 160-bit GLS curve, similar to the attack of §3.

4.1 A secure choice of (k 1 , k 2 )
Let E be a curve of prime order n over F p 2 obtained by the quadratic GLS method as above. In view of (1), we have:

(p -1) 2 ≤ n ≤ (p -1) 2 + (2 √ p) 2 = (p + 1) 2 ,
and the inequalities are in fact strict, since n is prime. Thus, we have p -1 < √ n < p + 1, and it follows that the distribution of 2 . We will thus concentrate on the latter, and show that it is close to uniform in Z/nZ using the following lemma. Lemma 3. The following map is injective.

k = k 1 + k 2 λ for (k 1 , k 2 ) uniform in [0, √ n) 2 is statistically close to the distribution of the same k for (k 1 , k 2 ) uniform in [0, p -1)
F : [0, p -1) 2 -→ Z/nZ (k 1 , k 2 ) -→ k 1 + k 2 λ. Proof. Consider two distinct pairs (k 1 , k 2 ) = (k 1 , k 2 ) such that F (k 1 , k 2 ) = F (k 1 , k 2 ). We have: (x -x ) + (y -y )λ ≡ 0 (mod n) (x -x ) 2 ≡ λ 2 (y -y ) 2 (mod n) (x -x ) 2 + (y -y ) 2 ≡ 0 (mod n),
since λ 2 ≡ -1 (mod n). Thus, the positive integer (x-x ) 2 +(y -y ) 2 is divisible by n, and it is also smaller than 2(p -1) 2 < 2n, so we must have (x -x ) 2 + (y -y ) 2 = n. In other words, (x -x ) 2 + (y -y ) 2 is a decomposition of n as a sum of two squares. Now it is well-known that, as a prime number, n has at most one decomposition as a sum of two squares up to order and sign (see e.g. [20, §3.6]), and (p -1) 2 + t 2 is one such representation. As a result, we must have either x -x = ±(p -1) or y -y = ±(p -1), and neither is possible since those difference are bounded by p -2 in absolute value. Hence, F is injective as required. 2 is statistically close to the uniform distribution on Z/nZ. More precisely, the statistical distance:

Theorem 1. The distribution of the values k

= k 1 + k 2 λ for (k 1 , k 2 ) uniform in [0, p -1)
∆ 1 = k∈Z/nZ Pr k = k 1 + k 2 λ ; (k 1 , k 2 ) $ ← [0, p -1) 2 - 1 n
is given by ∆ 1 = 2t 2 /n, which is negligible.

Proof. Indeed, since the function F above is injective by Lemma 3, the probability Pr 2 for each of the (p -1) 2 points in the image of F , and 0 for each of the n -(p -1) 2 = t 2 points outside of that image. Therefore:

k = k 1 + k 2 λ ; (k 1 , k 2 ) $ ← [0, p -1) 2 is equal to 1/(p -1)
∆ 1 = (p -1) 2 • 1 (p -1) 2 - 1 n + t 2 • 0 - 1 n = 1 - (p -1) 2 n + t 2 n = 2t 2 n
as required. This is bounded above by 8p/(p-1) 2 , which is indeed negligible.

Remark 1. Theorem 1 means that it is secure, in any ECC protocol instantiated over the GLS curve E, to sample random scalars k by picking k 1 and k 2 uniformly in [0, p -1), or equivalently [0, √ n). As we can see, the proof relies on the particular arithmetic properties of the quadratic GLS method (mainly the fact that λ = √ -1 in Z/nZ), so that the result does not readily extend to different settings, like the GLV method on a curve of CM discriminant -3. And indeed, in that case, Brumley and Nyberg have provided evidence that choosing (k 1 , k 2 ) uniformly in [0, √ n) may not yield a close to uniform distribution for k [7, Example 3]. They suggest an alternate approach to select intervals to choose k 1 and k 2 from and still achieve high entropy in a more general setting, but since the quadratic GLS method is one of the most used variants of GLV/GLS, we believe Theorem 1 is of significant practical interest.

4.2 Breaking insecure choices of (k 1 , k 2 ) with Bleichenbacher's attack

In the quadratic GLS setting, we have just seen that choosing (k 1 , k 2 ) uniformly in [0, √ n) 2 yields a close-to-uniform distribution of k = k 1 + k 2 λ. However, we can reasonably suspect that if we choose k 1 and k 2 uniformly in [0, 2 m ), m = 1 2 log 2 n (i.e. uniform bitstrings of length just under half of the size of n), the distribution of k will no longer be uniform. This is not immediately visible on the bias, however.

Indeed, if we let T = 2 m and define K 1 , K 2 as independent uniform random variables over [0, T ) and K as the random variable in Z/nZ given by K = K 1 + K 2 λ, we have, by Lemma 1:

B n (K) = B n (K 1 ) • B n (λK 2 ) = 1 T sin(πT /n) sin(π/n) • 1 T sin(πλT /n) sin(πλ/n) .
The first factor is very close to 1, but the second factor is usually negligible. For example, on the 160-bit GLS curve (3) below, we have T = 2 79 and B n (λK 2 ) ≈ 1.52/T . As a result, Bleichenbacher's attack does not apply directly to this setting in general.

However, since λ ≡ t -1 (p -1) (mod n), we claim that there is a significant bias on the values t • k. Indeed, we have:

B n (tK) = B n (tK 1 ) • B n (p -1)K 2 = 1 T sin(πtT /n) sin(πt/n) • 1 T sin(π(p -1)T /n) sin(π(p -1)/n) = 1 T πtT /n + O((tT /n) 3 ) πt/n + O((t/n) 3 ) • 1 T | sin(π(p -1)T /n)| π(p -1)/n + O(((p -1)/n) 3 ) = 1 + O (tT /n) 2 + (p/n) 2 • sin(π(p -1)T /n) π(p -1)T /n .
The big-O in the first factor is negligible since tT /n = Θ(p 1/2 • p/p 2 ) = Θ(p -1/2 ) and p/n = Θ(p -1 ). On the other hand, (p -1)T /n ≈ T / √ n is roughly between 0.5 and 1 depending on how close n is to a power of two. Thus, the bias is significant in general, and is maximal when (p -1)T /n is smallest (close to 1/2), which happens when n is just under a power of two. The bias B n (tK) is then close to 1/(π • 1/2) = 2/π ≈ 0.637.

It is then straightforward to adapt Bleichenbacher's attack to this setting by targetting the values t•k instead of k. We can then break ECDSA signatures that use nonces of the form k = k 1 +k 2 λ above using that variant. An implementation of that attack is discussed in the next subsection.

Implementation of Bleichenbacher's attack in the GLS setting

We carry out the attack described above on the 160-bit GLS curve E defined as follows. Over the 80-bit prime field7 F p , p = 255 • 2 72 + 1, we define E 0 : y 2 = x 3 -3x/23 + 104. Then, the elliptic curve E is the quadratic twist of E 0 over F p 2 = F p ( √ 23), namely:

E : y 2 = x 3 -3x + 104 • √ 23 3 over F p 2 . ( 3
)
The order of E 0 (F p ) is p + 1 -t for t = 776009485427, and E(F p 2 ) is of prime order n = (p -1) 2 + t 2 . The theoretical value of the bias B n (tK), computed using the exact formula above, is then ≈ 0.634. We performed the recovery of 32 MSB of a private key as in section 3.2. We computed 2 33 signatures and unrolled the attack on (tc j mod n, th j mod n) instead of (c j , h j ). We checked the bias and obtained ≈ 0.634116 which is close to the theory. In practice the attack took about 2000 CPU-hours, with 56% for the signature generation, 37% for the four sort-and-difference reduction steps, 5% for the candidate selection and FFT table preparation and less than 0.5% for the FFT itself. In wall-clock time terms, except for the signature generation which took (much) longer, other phases were identical as 3.2. We attribute this unexpected increase in signing time to threshold effects: for example, representing elements on a prime field with ≈ 2 160 elements needs only 3 64-bit words, whereas a on F p 2 we needed 4 × 2 = 8 words.

In this section, we analyze the security of algorithms for computing the decomposition of the nonces used in the GLV method from a side-channel analytic perspective. Many techniques have been proposed, including [START_REF] Gallant | Faster Point Multiplication on Elliptic Curves with Efficient Endomorphisms[END_REF][START_REF] Park | An Alternate Decomposition of an Integer for Faster Point Multiplication on Certain Elliptic Curves[END_REF]. The original GLV method [START_REF] Gallant | Faster Point Multiplication on Elliptic Curves with Efficient Endomorphisms[END_REF] based on LLL reduction of a lattice that depends on the nonce k, and variants thereof, have an execution time that depends on k, and are therefore vulnerable to timing attacks.

Therefore, we examine the security of a potentially more secure approach, the Park et al. [START_REF] Park | An Alternate Decomposition of an Integer for Faster Point Multiplication on Certain Elliptic Curves[END_REF] decomposition technique, using more involved power analysis technique.

Decomposition Algorithm

Park et al. provide an alternative decomposition to the GLV paper [START_REF] Gallant | Faster Point Multiplication on Elliptic Curves with Efficient Endomorphisms[END_REF] which reduces the theoretical bound for the decomposition using the theory of µ-Euclidian algorithm and is a little bit faster. The algorithm requires two short and independent vectors v 1 and v 2 of the two-dimensional lattice L = {(x, y) :

x + yλ = 0 mod n}. We can find these vectors during a precomputation time using the Gauss reduction. The algorithm consists in finding a vector in the lattice L = Zv 1 + Zv 2 that is close to (k, 0) using linear algebra. Then, (k 1 , k 2 ) is determined by the equation:

(k 1 , k 2 ) = (k, 0) -( b 1 v 1 + b 2 v 2 ), where (k, 0) = b 1 v 1 + b 2 v 2 is an element of Q × Q.
Algorithm 3 Decomposition technique of Park et al. in [START_REF] Park | An Alternate Decomposition of an Integer for Faster Point Multiplication on Certain Elliptic Curves[END_REF] Require: k ≈ n, the shortest vectors v1 = (x1, y1), v2 = (x2, y2) Ensure: (k1, k2) such that k = k1 + k2λ(modn)

1: D = x1y2 -x2y1, a1 = y2k, a2 = -y1k 2: zi = ai/D for i = 1, 2 3: k1 = k -(z1x1 + z2x2), k2 = z1y1 + z2y2 return (k1, k2)
The decomposition technique depicted in Algorithm 3 makes many computations involving the sensitive nonce k. Particularly, the computation of a 1 (resp. a 2 ) is based on a multiplication of the nonce k by y 2 (resp. y 1 ) which is assumed to be known since it is a precomputed value obtained from public parameters using a deterministic algorithm.

Suppose now that we obtain the knowledge of the least significant byte of nonces k 1 , • • • , k . The best strategy for finding the secret key x consists in performing classical lattice attacks as proposed in [START_REF] Howgrave-Graham | Lattice Attacks on Digital Signature Schemes[END_REF][START_REF] Nguyen | The Insecurity of the Digital Signature Algorithm with Partially Known Nonces[END_REF][START_REF] Nguyen | The Insecurity of the Elliptic Curve Digital Signature Algorithm with Partially Known Nonces[END_REF]. For a 160-bit modulus, the lattice attack works consistently for 27. However the side-channel attack may sometimes fail, i.e. the returned byte of some k j can be a wrong value. Thus, by denoting 0 < c < 1 the confidence rate, the side-channel attack has to be performed on m > 27/c signatures. Then:

-Select 27 signatures at random among them.

-Perform the attack using these signatures.

-If the attack fails, goto the first step.

The probability of success at each iteration of the lattice attack is m•c 27 / m 27 . As an example, suppose we obtain m = 200 signatures, and can guess the least significant byte with 90% accuracy (c = 0.9). Then the probability of success of the lattice attack is about 4.7% and 21 lattice reductions have to be performed on average. Since LLL reductions are cheap, much lower success probabilities are tractable as well.

In the following, we discuss the side-channel attack that aims at recovering the first byte of the nonce targeting the two aforementioned multiplications. We present the attack in the particular case of a 8-bit implementation (that corresponds to the device we used in experiments). Note that this attack may also work for 16-bit implementation but in this case the computational cost will be larger and the success rate smaller.

Side-Channel Attack on this implementation

The details of the attack highly rely on the way the multiplication is implemented. Depending of the underlying algorithm, the attack may be more or less difficult. We present here the attack corresponding to the implementation we target but we will discuss adaptations to different algorithms. The multiplication we target is a schoolbook multiplication with the nonce being scanned in the outter loop. Algorithm 4 outlines the implementation of such multiplication for n -bit nonces and n/2 -bit b. The idea is to take profit of all operations involving the first nonce-byte in the inner loop to recover its value. This can be done by propagating a probability distibution from an operation to another and updating it with the corresponding leakages. Since the nonce bits have to be recovered using a single trace (the nonce is randomly generated for each signature) we place ourselves in the context of a profiled attack. The application of such an attack in a non-profiled setting is left as an open question.

Template Attack on One

Step. One step of the inner loop consists in a multiplication of the first byte of the nonce k 0 , a byte of the auxiliary input b i and the carry c i . This results in a value v i and a new carry c i+1 . We may obtain leakages for each of these variables. We denote by capital letters the output distributions of template exploitation corresponding to small letter variables. For instance, after processing the leakage corresponding to c i , the attacker gets a distribution C i = Pr(c i = 0), Pr(c i = 1), . . . , Pr(c i = 255) .

Since these variables may be manipulated more than once during the computation, different leakage points may be combined by multiplicating probabilities then normalizing the resulting distribution. More precisely, let l 1 , l 2 , . . . , l l be leakages corresponding to variable k 0 , then the distribution K 0 obtained from these leakages is computed as

Pr(k 0 = x) = 1 Z l j=1
Pr(k 0 = x|l j ), where the normalizing coefficient Z is given by x l j=1 Pr(k 0 = x|l j ).

Propagating and Updating Distribution. Let us now discuss how to take profit of all the leakages of the inner loop to gain information on the byte k 0 . The main idea is to gather all the information from all variables of a given step i into distribution K 0 and C i+1 then do the same at step i + 1 using the newly updated distributions. From a probabilistic point of view we should compute the joint distribution of variables of step i then compute marginalized distributions K 0 and C i+1 . The following algorithm updates the distributions K 0 and C i+1 according to the distributions of variables b i , v i and c i . The attacker starts with using Algorithm 5 for the first step. Then she uses the newly updated distributions K 0 and C 1 and the initial distributions B 1 , V 1 and C 2 as inputs of Algorithm 5 and so on . . . At the end, the attacker gets the final distribution K 0 from which she can derive the most likely value of the least significant bit (or more).

  Bn(K) 0.6366198 0.9003163 0.9744954 0.9935869 0.9983944 Table

Algorithm 4 7 : end for 8 :

 478 Multiplication v = kb of k = n /8 i=0 k i 2 8i times b = n/2 /8 i=0 b i 2 8i Require: n-bit k and n/2 -bit b two integers, v = 0 Ensure: v = k × b 1: v ← 0 2: for i = 0 to i < n/8 do 3: c0 ← 0 4: for j = 0 to j < n/2 /8 do 5: vi+j = (ki × bj + cj) & 0xFF 6: cj+1 = (ki × bj + cj)8 end for return v

This is an example of "optimal prime field" (OPF). See e.g.[START_REF] Wenger | An 8-bit AVR-based elliptic curve cryptographic RISC processor for the Internet of Things[END_REF].


Marc Dinh

Anton Korniienko
email: anton.korniienko@ec-lyon.fr

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

Gérard Scorletti Convex

Convex Hierarchical Analysis for the Performances of Uncertain Large-Scale Systems

Keywords: Large-scale system, hierarchical analysis, convex optimization. 1

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Since the beginning, robustness to uncertainties is a major issue. The sensitivity transfer function was introduced by Bode as a measure of the impact of uncertainties on the closed loop. Some design methods were developed with robustness at their core. It is the case for Quantitative Feedback Theory [START_REF] Houpis | Quantitative feedback theory: fundamentals and applications[END_REF] or H ∞ [28]. In the 80's-90's, µ-analysis [START_REF] Doyle | Analysis of feedback systems with structured uncertainties[END_REF][START_REF] Safonov | Stability margin of diagonaly perturbed multivariable feedback systems[END_REF] was developed to investigate the (H ∞ ) performances of Linear Time Invariant systems in the presence of structured uncertainties. This approach is based on the computation of the structured singular value µ of frequency dependent matrices, which was proved to be NPhard [START_REF] Braatz | Computational complexity of µ calculation[END_REF]. Fortunately, lower and upper bounds on µ can be efficiently computed; the µ upper bounds [START_REF] Fan | Robustness in the presence of mixed parametric uncertainty and unmodeled dynamics[END_REF] allow to guarantee a certain level of performances with some conservatism. By efficient, it is understood that the computation time is bounded by a power function of the problem size [START_REF] Garey | Computers and Intractability: A guide to the theory of NP-Completeness[END_REF]. A major interest of the use of µ upper bounds is to obtain a satisfying trade-off between the computation time and the conservatism of the obtained result.

Nevertheless, even if the computation of µ upper bound is efficient, its computation time can be important in the case of uncertain large-scale systems. The purpose of this paper is to propose a robustness analysis (of performances) method with a trade-off between computation time and conservatism adapted to large-scale systems. Our motivating example is the robustness analysis of Phase-Locked-Loop (PLL) networks [START_REF] Korniienko | Control law design for distributed multi-agent systems[END_REF], a challenging problem in Microelectronics.

In the case of large-scale systems with uncertainties defined by conic sector properties, that is unstructured uncertainties, Safonov approached the problem via an hierarchical analysis using recursively a propagation of conic sector properties result [START_REF] Safonov | Propagation of conic model uncertainty in hierarchical systems[END_REF]. To the authors' best knowledge, the work of [START_REF] Safonov | Propagation of conic model uncertainty in hierarchical systems[END_REF] is the only one devoted to the problem. Basically, a hierarchical system can be described as a tree with layers, each of them is composed of systems. Each system is modeled as the interconnection of sub-systems. These sub-systems are in fact the systems of the preceding layer. Assuming that the conic properties of the sub-systems are known, the propagation problem is to find 'propagated' conic properties for the system that are not trivially connected to the ones of the sub-systems. If one is able to find these 'propagated' properties, then the hierarchical analysis boils down to a recursive application of the propagation.

Unfortunately, its direct application to a (large-scale) system with structured uncertainties can lead to an overly conservative result. While keeping the same overall approach, we propose in this paper a new hierarchical analysis method which overcomes the disadvantages of [START_REF] Safonov | Propagation of conic model uncertainty in hierarchical systems[END_REF]: the method is efficient and adapted to structured uncertainties. Our method is based on a generalization of the propagation of conic properties to quadratic ones. Preliminary results of the present work can be found in [START_REF] Dinh | Embedding of uncertainty propagation: Application to hierarchical performance analysis[END_REF] where there was only one 'propagated' quadratic property (the conic sector one) with the uncertainties already described by quadratic properties. We present here other quadratic 'propagated' proper-ties. Combining them (their 'intersection'), a finer description can be made of the system. The user can then set the trade-off by choosing the 'propagated' properties used, the main contributor to the computation time being their number.

As in [START_REF] Dinh | Embedding of uncertainty propagation: Application to hierarchical performance analysis[END_REF], this propagation can be viewed as an embedding problem: find (simpler) sets that includes the uncertain system. Beyond the hierarchical analysis, this embedding has many interesting applications. For instance, in the Quantitative Feedback Theory [START_REF] Houpis | Quantitative feedback theory: fundamentals and applications[END_REF], the aim is to design a controller that achieves some level of performances in front of all the uncertainties in the system. In the method, it is assumed that it is possible to know the propagation of the uncertainties through the system. Another interesting application, is the use of the embedding to perform µ-synthesis in the integrated framework of [START_REF] Paw | Development and application of an integrated framework for small uav flight control development[END_REF]; the advantage being the computation time thanks to the embedding.

Our solution is based on a separation of graph theorem. First proposed in [START_REF] Safonov | Stability and Robustness of Multivariable Feedback Systems[END_REF] as a general approach to feedback system analysis, specialized forms were proposed in e.g. [START_REF] Iwasaki | Well-posedness of feedback systems: Insights into exact robustness analysis and approximate computations[END_REF][START_REF] Scorletti | A more praticle formulation for robustness analysis[END_REF][START_REF] Scherer | LPV control and full block multipliers[END_REF] for (uncertain) LTI system analysis. Since the µ upper bound proposed in [START_REF] Fan | Robustness in the presence of mixed parametric uncertainty and unmodeled dynamics[END_REF] can be interpreted as a particular application of the separation of graph theorem, this theorem was applied to extend µ-analysis to time-delay systems [START_REF] Scorletti | Approche Unifiée de l'analyse et la commande des systèmes par formulation LMI[END_REF] or time-varying/nonlinear systems [START_REF] Megretski | System analysis via integral quadratic constraints[END_REF], to reduce the conservatism of the µ upper bound [START_REF] Scorletti | Improved efficient analysis for systems with uncertain parameters[END_REF] to cite a few. We reveal here another interesting application of this powerful theorem.

Paper outline

Section 2 begins with definitions and fundamental properties of dissipative systems which are used afterward. It explicits then the uncertain large-scale system that is considered with the proposed approach. Section 3 proposed several dissipative properties that can be used practically. A numerical example on a PLL network is performed in Section 4. Section 5 concludes the paper.

Notations R (respectively C) denotes the set of real (resp. complex) numbers. R denotes R ∪ {-∞, +∞}. M R and M I stands for the real and imaginary parts of M .

R m×n (respectively C m×n ) denotes the set of real (resp. complex) matrices of dimension m × n. I r and 0 r denote the identity and the zero matrices of size r. M * (respectively M T ) stands for the transpose conjugate (resp. transpose) of M . For several matrices M i , i = 1, . . . , n, bdiag i (M i ) denotes the matrix

     M 1 . . . Mn      . RH ∞ (respectively RL ∞ )
denotes the set of matrices of stable (resp. non-causally stable) rational transfer function.

Throughout the paper, when dimensions are omitted, they are assumed to be clear from the context. Moreover, we consistently denote uncertainties by ∆ and interconnections by M . The set ∆ is referred to as the uncertainty set. We denote by ∆ M the set {∆ M, ∆ ∈ ∆}, with standing for the Redheffer star product. This set is also referred to as an uncertain system. For the uncertain system ∆ M , we further denotes the interconnection's partitioning of appropriate dimension by

M = A B C D .
Finally, we denote

L(M, Φ 11 , Φ 12 , Φ 22 , X, Y, Z) = M I *       -Φ 22 0 -Φ * 12 0 0 X 0 Y -Φ 12 0 -Φ 11 0 0 Y * 0 Z       M I .
2 Approach for Hierarchical Analysis of Performances

Definitions and preliminaries

An uncertain system is modeled as an interconnection ∆ M with ∆ ∈ ∆. Along with this definition, we assume: Introducing the internal signals and using the Fourier transform, we get for uncertain system ∆ M represented in the Fig. 1 and for all ω ∈ R:

p(jω) = ∆(jω) q(jω) q(jω) z(jω) = M (jω) p(jω) w(jω) = A(jω) B(jω) C(jω) D(jω) p(jω) w(jω) . (1)
The stability of an uncertain system is now defined.

Definition 2.1 An uncertain system ∆ M is said to be stable (respectively non-causally stable) if for any ∆ ∈ ∆, the system ∆ M is stable (respectively non-causally stable).

In this section, dissipative properties are used. They are defined below.

Definition 2.2 Let X(jω), Y (jω) and Z(jω) be 3 transfer functions of RL ∞ such that X(jω) = X(jω) * and Z(jω) = Z(jω) * . Then,

∆ ∈ ∆, the system ∆ M is {X(jω), Y (jω), Z(jω)} dissipative.
Note that when H is stable and

X(jω) Y (jω) Y (jω) * Z(jω) = -I 0 0 γ 2 I ,
H has an H ∞ norm strictly less than γ. More generally, dissipative properties can model performance indices. A dissipativity property can also be viewed as an inclusion of sets (of systems or signals). With the preceding example, the set composed of H is included in the ball of systems with an H ∞ norm strictly less than γ. This set of systems is convex which is proved by the following Lemma.

Lemma 2.1 Let X(jω), Y (jω) and Z(jω) be 3 transfer functions of RL ∞ such that X(jω) = X(jω) * and Z(jω) = Z(jω) * . Then, the set of systems

{H | H is {X(jω), Y (jω), Z(jω)} dissipative }
is convex, and thus connected.

Proof See Appendix A.

Fundamental properties are now given with a direct corollary. They states that new dissipativity properties can be generated from original ones. For ease of notation, they are stated for certain systems. The extension to uncertain systems is straightforward.

Lemma 2.2 Let H i be {X i (jω), Y i (jω), Z i (jω)} dissipative, i = 1, . . . , m. Then H = bdiag i (H i ) is {bdiag i (X i (jω)), bdiag i (Y i (jω)), bdiag i (Z i (jω))} dissipative. Lemma 2.3 Let H be {X k (jω), Y k (jω), Z k (jω)} dissipative, k = 1, . . . , n. Then for any τ k (jω) > 0, k = 1, . . . , n, H is { k τ k (jω)X k (jω), k τ k (jω)Y k (jω), k τ k (jω)Z k (jω)} dissipative. Corollary 2.1 Let H i be {X ik (jω), Y ik (jω), Z ik (jω)} dissipative, i = 1, . . . , m, k = 1, . . . , n. Then for any τ ik (jω) > 0, i = 1, . . . , m, k = 1, . . . , n, H = bdiag i (H i ) is {bdiag i ( k τ ik (jω)X ik (jω)), bdiag i ( k τ ik (jω)Y ik (jω)), bdiag k ( k τ ik (jω)Z ik (jω))} dissipative.
This corollary also defines a set of linearly parameterized dissipative properties. Let us denote it Φ(jω).

Hierarchical system description and proposed approach

From [START_REF] Safonov | Propagation of conic model uncertainty in hierarchical systems[END_REF], a large-scale system is described by a tree as illustrated in Figure 2 where a hierarchical structure arises naturally. Each branch of the tree is assigned an index. A branch, say i, is a two-way channel through which a signal w i (the input) ascends and another signal z i (the output) descends. The tree obtained by cutting branch i and retaining everything connected above is an uncertain system called T i with input w i and output z i . If a tree T i has other branches besides branch i then there is a single node denoted M i from which other branches ascend. If branch i is the only branch in the tree T i , then T i is called a leaf and is denoted by ∆ i . It is uncertain system belonging to the uncertain set ∆ i . Each M i and ∆ i is an LTI system. Furthermore, it is supposed that for each leaf ∆ i its dissipative properties are a priori known.

It should be noted from the practical point of view, that all uncertainties are represented by the leafs ∆ i while the certain part as well as the system interconnection structure are described by nodes M i . Assumption 2.2 Each ∆ i is a bounded and connected subset of RH ∞ and each M i belongs to RH ∞ . This is the counterpart of Assumption 2.1.

Assumption 2.3 For any i, ∆ i is an elementary uncertainty set: there exists a priori known X ik (jω), Y ik (jω) and Z ik (jω) such that ∆ i is {X ik (jω), Y ik (jω), Z ik (jω)} dissipative for k = 1, . . . , n and some n > 0. Examples of elementary uncertainty sets along with their dissipative properties can be found in [START_REF] Fan | Robustness in the presence of mixed parametric uncertainty and unmodeled dynamics[END_REF][START_REF] Scorletti | Robustness analysis with time delays[END_REF][START_REF] Scorletti | Improved efficient analysis for systems with uncertain parameters[END_REF]. The interest of this assumption will be made clearer in the discussion of the propagation theorem.

M 1 M 2 ∆ 5 M 3 ∆ 6 M 4 M 7 ∆ 8 ∆ 9 ∆
The proposed approach for the performance analysis of a hierarchical system as described in the introduction is based on a recursive application of the following propagation1 of dissipativity properties of uncertainties through an interconnection theorem.

Theorem 2.1 The uncertain system ∆ M is stable and {X(jω), Y (jω), Z(jω)} dissipative if and only if:

1. there exists ∆ 0 ∈ ∆ such that the system ∆ 0 M is stable; 2. there exists 3 transfer functions Φ 11 (jω), Φ 12 (jω) and

Φ 22 (jω) of RL ∞ , with Φ 11 (jω) = Φ 11 (jω) * and Φ 22 (jω) = Φ 22 (jω) * , such that: the uncertainty set ∆ is {Φ 11 (jω), Φ 12 (jω), Φ 22 (jω)} dissipative and ∀ω ∈ R, L(M (jω), Φ 11 (jω), Φ 12 (jω), Φ 22 (jω), X(jω), Y (jω), Z(jω)) > 0. ( 2
)
The uncertain system ∆ M is non-causally stable if and only if merely the condition 1 holds.

Proof See Appendix B.

Condition 1 is generally viewed as an assumption which is verified beforehand on the nominal system. This assumption is hopefully very mild and is extremely similar (even weaker in fact) to the one of µ-analysis [START_REF] Skogestad | Multivariable Feedback Control, Analysis and Design[END_REF].

The fact that the uncertainty set ∆ is {Φ 11 (jω), Φ 12 (jω), Φ 22 (jω)} dissipative in condition 2 can also be verified a priori using Assumption 2.3 and the set of linearly parameterized dissipativity properties Φ(jω) as defined by Corollary 2.1. Note however that this leads to sufficient conditions only. Note also that it is not compulsory as an assumption as it is not used in the proof. It is possible to find them directly as in [START_REF] Feron | Analysis and synthesis of robust control systems via parameter-dependent lyapunov functions[END_REF] for instance. This choice has been made to increase the efficiency of the approach.

As a consequence, Theorem 2.1 boils down in practise to verify condition (2).

Corollary 2.2 Let Φ(jω) be a set (defined by Corollary 2.1) such that for any

(Φ 11 (jω), Φ 12 (jω), Φ 22 (jω)) ∈ Φ(jω), the uncertainty set ∆ is {Φ 11 (jω), Φ 12 (jω), Φ 22 (jω)} dissipative.
Then the uncertain system ∆ M is non-causally stable and {X(jω), Y (jω), Z(jω)} dissipative if there exists

(Φ 11 (jω), Φ 12 (jω), Φ 22 (jω)) ∈ Φ(jω) such that ∀ω ∈ R, L(M (jω), Φ 11 (jω), Φ 12 (jω), Φ 22 (jω), X(jω), Y (jω), Z(jω)) > 0.
Moreover, ∆ needs to be a bounded and connected set only: it does not need to contain 0 as usually assumed. In our case it is an important fact from a practical point of view as the theorem is recursively applied so that the uncertainty set (the previous branches) does not necessarily contain 0.

When two uncertain systems ∆ 1 M 1 and ∆ 2 M 2 are homogenous, then they share the same dissipative properties: if

M 1 = M 2 and ∆ 1 = ∆ 2 , then ∆ 1 M 1 is {X(jω), Y (jω), Z(jω)} dissipative if and only if ∆ 2 M 2 is {X(jω), Y (jω), Z(jω)} dissipative.
This is the case of the PLL network example of Section 4.

For computational purposes, it may also be interesting to replace

M (jω) I by N (jω) D(jω) with M (jω) = N (jω)D(jω) -1 .
Finally, the theorem also states that the analysis can be performed by a frequency by frequency i.e. applying the frequency gridding approach. This fact will be used in the next section to drop the dependency on jω.

Let us now illustrate the use of Theorem 2.1 with the system presented in Figure 2. First, from the dissipative properties of ∆ 9 and ∆ 10 , find some dissipative properties of the branch T 7 using Theorem 2.1 with ∆ = bdiag(∆ 9 , ∆ 10 ) and M = M 7 . From these several dissipative properties and from the ones of ∆ 8 , use again Theorem 2.1 (and Corollary 2.1) with ∆ = bdiag(T 7 , ∆ 8 ) and M = M 4 to find dissipative properties of the branch T 4 . However, it is possible if T 7 is a bounded set (the connected part of the assumption is ensured by Lemma 2.1). This is the case if its dissipativity properties were well chosen: typically, a conic sector property is needed, see Section 3.1. And so on until branch T 1 where the dissipativity property is a performance index. The overall trade-off between conservatism and computation time then depends on the number of dissipative properties that are searched for at each step. The user can thus set this trade-off by setting the number of dissipative properties.

Practical Formulation of Dissipativity Propagation

In this section, we show how to find dissipative properties (referred to as 'propagated' in the introduction) for the uncertain system ∆ M from the ones of ∆: it is the propagation of dissipativity properties. For the sake of simplicity, we set the value of the frequency without loss of generality as the system is linear time-invariant, so that the dependency on jω is dropped. The problem can be stated as follows.

Problem 3.1 Let Φ be a set such that for any

(Φ 11 , Φ 12 , Φ 22 ) ∈ Φ, the uncertainty set ∆ is {Φ 11 , Φ 12 , Φ 22 } dissipative.
From the set Φ, find X, Y and Z such that the uncertain system ∆ M is {X, Y, Z} dissipative.

L(M, Φ 11 , Φ 12 , Φ 22 , X, Y, Z) > 0. Proof It is a direct consequence of Corollary 2.2.
Note that the optimization problem defined in Corollary 3.1 parameterizes all the possible propagated properties from the ones of ∆ in Φ: it is non conservative from a propagation perspective.

In the way the propagation is used, ∆ is either a leaf or a branch. In both case, either due to Assumption 2.3 or Corollary 2.1, the set Φ is of the form

i τ i Φ i 11 , Φ i 12 , Φ i 22 with a priori known (Φ i 11 , Φ i 12 , Φ i 22 )
. The optimization problem defined in Corollary 3.1 boils down to find τ i and X, Y and Z such that L(M, Φ 11 , Φ 12 , Φ 22 , X, Y, Z) > 0 with

Φ 11 = i τ i (Φ i 11 ), Φ 12 = i τ i (Φ i 12 ), Φ 22 = i τ i (Φ i 22 )
. It is thus an LMI optimization problem, is convex and can be solved efficiently. Corollary 3.1 defines a optimization problem with complex LMI constraints. For computational purpose, they can readily be converted as real LMI constraints [START_REF] Boyd | Method of centers for minimizing generalized eigenvalues[END_REF]: a complex matrix M is positive definite if and only if the real matrix

M R M I -M I M R is positive definite.
The three preceding remarks hold for all the optimization problems involved in this section.

To improve the overall conservatism of the hierarchical analysis, it is interesting to obtain the 'tightest' propagated dissipativity property. It is performed by interpreting the property in geometrical terms. For each geometrical interpretation, a notion of size is defined and one is interested in minimizing this size.

Conic sector: X < 0

In the case, when X < 0, let us define

H c = -X -1 Y and R * R = Z -Y * X -1 Y.
The dissipativity property of a system H writes then

(z -z c ) * (-X)(z -z c ) < w * R * Rw with z = Hw and z c = H c w. (3)
which defines the same set as a conic sector [START_REF] Safonov | Propagation of conic model uncertainty in hierarchical systems[END_REF] in which a system H is said to be in the conic sector (C, P, Q), with C the cone center, whenever

Q -1/2 (z -Cw) 2 < P 1/2 w 2 with z = Hw
The link is provided by

-X = Q -1 , H c = C and R * R = P.
For a SISO system, the inequality (3) defines a disk of center z c and radius w * (X/(R * R))w.

More generally, it is an ellipsoid. Indeed, the inequality can be rewritten as

z R -z cR -(z I -z cI ) T P z R -z cR -(z I -z cI ) < 1 with P = 1 w R -w I T (R * R) R (R * R) I -(R * R) I (R * R) R w R -w I -X R -X I X I -X R . (4)
Thus, for a given non null input w, the corresponding output signal z = z T R -z T I T , with zc = z T cR -z T cI T , belongs to the ellipsoid

P = {z | (z -zc ) T P(z -zc ) < 1}. ( 5
)
The volume is here evaluated as in [START_REF] Boyd | Linear Matrix Inequalities in Systems and Control Theory[END_REF].

Definition 3.1 The volume of the ellipsoid P defined by ( 5) and ( 4) is defined as

vol( P ) 2 = β det(P -1 )
where β is a positive scalar which depends on the size n z of the vector z -zc .

We are interested in finding the smallest one for all inputs such that w = 1.

-X R -X I X I -X R -1
and such that

1. L(M, Φ 11 , Φ 12 , Φ 22 , X, Y, Z) > 0 holds; 2.      I 0 0 0 0 0 0 0 0 0 I 0 0 0 0 0      ≥      Z R Y * R Z I Y * I Y R X R Y I X I -Z I -Y * I Z R Y * R -Y I -X I Y R X R      holds.
This optimization problem is a determinant maximization under linear matrix inequality constraints [START_REF] Vandenberghe | Determinant maximization with linear matrix inequality constraints[END_REF] and is convex.

Proof Problem 3.2 writes minimize maximize maximize vol( P ) 2 over X, Y, Z over ∆ ∈ ∆ over w = 1 subject to ∆ M is {X, Y, Z} dissipative
As the logarithm function is strictly increasing and as β is constant, the optimization problem is equivalent to minimize maximize maximize log(det(P

-1 )) over X, Y, Z over ∆ ∈ ∆ over w = 1 subject to ∆ M is {X, Y, Z} dissipative Now max w =1 log(det(P -1 )) is equal to max w =1 log det λ max -X R -X I X I -X R -1
with λ max the minimal value verifying

λ max I ≥ (R * R) R (R * R) I -(R * R) I (R * R) R .
As a dissipativity property is defined up to a strictly positive multiplicative coefficient and as {X, Y, Z} dissipativity defines the same ellipsoid as {τ X, τ Y, τ Z} dissipativity for any τ > 0 since

1 w R w I T τ (R * R) R τ (R * R) I -τ (R * R) I τ (R * R) R w R w I -τ X R -τ X I τ X I -τ X R = P,
one can search for X, Y and Z such that λ max = 1 without loss of generality. Thus, the optimization problem is equivalent to

minimize maximize log det -X R -X I X I -X R -1 over X, Y, Z over ∆ ∈ ∆ subject to ∆ M is {X, Y, Z} dissipative I ≥ (R * R) R (R * R) I -(R * R) I (R * R) R
Finally, condition 1 of Theorem 3.1 is obtained by applying Corollary 3.1 and condition 2 is obtained by applying Schur's lemma [START_REF] Boyd | Linear Matrix Inequalities in Systems and Control Theory[END_REF].

Half Planes: X = 0

Half plane A dissipativity property with X = 0 rewrites

ξ T z R -z I -η > 0 with ξ = 2 Y R Y I -Y I Y R w R -w I , η = - w R -w I T Z R Z I -Z I Z R w R -w I .
This last inequality express that the output signal belongs to a half plane defined by the hyperplane

z R z I ξ T z R -z I = η .
ξ is a vector normal to the hyperplane and η the 'signed distance' of the hyperplane to the origin (the dot product of any point of the hyperplane with ξ).

Remark: the passitivity kind of performance is a specific half plane with

X Y Y * Z = 0 I I 0 .
Band A band is the intersection of two half planes with the same normal direction but opposite way. As for the conic sector, we are interested in the smallest band for a given direction, that is for a given Y . The size of a band is defined by the distance between the two parallel hyperplanes.

Definition 3.2 Let ξ, η 1 and η 2 be of appropriate dimensions and define the two half planes

ξ T z R -z I -η 1 > 0 and -ξ T z R -z I -η 2 > 0.
The size of the corresponding band can be defined by

d Y = |η 1 + η 2 | = -η 1 -η 2 > 0.
We are interested in finding the smallest band for all inputs such that w = 1.

Problem
) ∈ Φ, Z 1 , λ 1 , (Φ 2 11 , Φ 2 12 , Φ 2 22 ) ∈ Φ, Z 2 , λ 2 and d that minimize d such that 1. L(M, Φ 1 11 , Φ 1 21 , Φ 1 22 , 0, Y, Z 1 ) > 0 holds; 2. L(M, Φ 2 11 , Φ 2 12 , Φ 2 22 , 0, -Y, Z 2 ) > 0 holds; 3. dI ≥ (Z 1 + Z 2 ) holds.
This optimization problem is a minimization of a linear cost under linear matrix inequality constraints [START_REF] Boyd | Linear Matrix Inequalities in Systems and Control Theory[END_REF] and is convex.

It is then possible to generate several bands with different normals as performed in the example of Section 4. It is also possible to search for the direction of the band by letting Y to be free.

Proof Problem 3.3 writes minimize maximize maximize |η 1 + η 2 | over X, Y, Z over ∆ ∈ ∆ over w = 1 subject to ∆ M is {0, Y, Z 1 } dissipative ∆ M is {0, -Y, Z 2 } dissipative with η 1 = - w R -w I T Z 1R Z 1I -Z 1I Z 1R w R -w I and η 2 = - w R -w I T Z 2R Z 2I -Z 2I Z 2R w R -w I .
Noting that, |η 1 + η 2 | = -η 1 -η 2 > 0 since the set ∆ M is not empty, the previous problem is equivalent to the optimization problem:

minimize maximize maximize (η 1 + η2 ) over Z 1 , Z 2 over ∆ ∈ ∆ over w = 1 subject to ∆ M is {0, Y, Z 1 } dissipative ∆ M is {0, -Y, Z 2 } dissipative with η1 = w R w I T Z 1R Z 1I -Z 1I Z 1R w R w I and η2 = w R w I T Z 2R Z 2I -Z 2I Z 2R w R w I . Now max w =1 (η 1 + η2 ) is equivalent to min d constrained by dI ≥ Z 1R Z 1I -Z 1I Z 1R + Z 2R Z 2I -Z 2I Z 2R .
After rewriting the LMI in complex form, the optimization is thus equivalent to

minimize maximize minimize d over Z 1 , Z 2 , d over ∆ ∈ ∆ over d subject to ∆ M is {0, Y, Z 1 } dissipative ∆ M is {0, -Y, Z 2 } dissipative dI ≥ (Z 1 + Z 2 )
Finally, conditions 1 and 2 of Theorem 3.2 are obtained by applying Corollary 3.1.

PLL network Example

Let us consider now a numerical example of hierarchical performance analysis of an uncertain large-scale system. One takes as an example the performance analysis of the active clock distribution network from [START_REF] Korniienko | Control law design for distributed multi-agent systems[END_REF] subject to technological dispersions. An active clock distribution network is composed of N = 16 mutually synchronized Phase-Locked-Loops (they constitute branches of the tree) delivering the clock signals to the chip (see the Fig. 3). To be able to synchronize the PLLs exchange the information on their relative phase through the interconnection network and the phase detectors. This example is particularly well adapted as the performance is measured in frequency domain with homogeneous PLLs.

PLL network description

Since the principal aim of the system is the synchronization, the PLLs are homogeneous, that is have a common interconnection and the same uncertainty set. Of course, during the manufacturing process, there are inevitable technological dispersions which can be represented in the form of parametric uncertainties belonging to the same set. We have thus ∀i ∈ {1, . . . , N }:

T i (jω) = k i (jω + a i ) -ω 2 + k i jω + k i a i (6)
where k i , a i are the real uncertain parameters defined as k i ∈ (0.76 • 10 4 , 6.84 • 10 4 ) and a i ∈ (91.1, 273.3). ω is the current frequency defined by gridding as suggested in the previous sections.

M net M P LL ∆ 1 M P LL ∆ i Branch/PLL i Branch/PLL 16 w z w net 1 w net i w net 16 q 1 q i z net 1 z net i z net 16 p 1 p i M P LL ∆ 16 q 16 p 16
Branch/PLL 1

Figure 3: PLL network topology

The exchange of information between the PLLs is represented in the Fig. 3 and is modeled by an interconnection matrix M net defined in [START_REF] Fan | Robustness in the presence of mixed parametric uncertainty and unmodeled dynamics[END_REF].

M net =                                0 1
3 0 0 1 3 0 0 0 0 0 0 0 0 0 0 0 1 3 1 3 0 1 3 0 0 1 3 0 0 0 0 0 0 0 0 0 0 0 0 1 3 0 1 3 0 0 1 3 0 0 0 0 0 0 0 0 0 0 0 0 1 2 0 0 0 0 1 2 0 0 0 0 0 0 0 0 0 1 3 0 0 0 0 1 3 0 0 1 3 0 0 0 0 0 0 0 0 0 1 4 0 0 1 4 0 1 4 0 0 1 4 0 0 0 0 0 0 0 0 0 1 4 0 0 1 4 0 1 4 0 0 1 4 0 0 0 0 0 0 0 0 0 1 3 0 0 1 3 0 0 0 0 1 3 0 0 0 0 0 0 0 0 0 1 3 0 0 0 0 1 3 0 0 1 3 0 0 0 0 0 0 0 0 0 1 4 0 0 1 4 0 1 4 0 0 1 4 0 0 0 0 0 0 0 0 0 1 4 0 0 1 4 0 1 4 0 0 1 4 0 0 0 0 0 0 0 0 0 1 3 0 0 1 3 0 0 0 0 1 3 0 0 0 0 0 0 0 0 0 1 2 0 0 0 0 1 2 0 0 0 0 0 0 0 0 0 0 0 0 1 3 0 0 1 3 0 1 3 0 0 0 0 0 0 0 0 0 0 0 0 1 3 0 0 1 3 0 1 3 0 0 0 0 0 0 0 0 0 0 0 0 1 2 0 0 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1

                               (7)
In this example, the transfer function between external signals w and z expresses the performance of the global PLL network and namely its ability to synchronize with periodic reference signal w. This reference signal is represented by its phase so that the PLL network has to track a ramp. More details can be found in [START_REF] Korniienko | Control law synthesis for distributed multi-agent systems: Application to active clock distribution networks[END_REF].

Hierarchical analysis set up

The proposed hierarchical analysis approach is applied in two steps for this PLL network:

1. obtain dissipativity properties of each individual PLL, each PLL being a branch. Here note that the PLL are homogeneous so that the dissipativity properties obtained for one PLL is valid for the others as well;

2. obtain the performance of the overall network through the interconnection of the 16 PLL branches and the matrix M net .

Individual PLL Each PLL can be readily written in the form of an interconnection, which leads after normalization of the uncertainties to:

T i (jω) = ∆ i (jω) M P LL (jω), ∆ i (jω) ∈ ∆
with ∆ of the form

δk i 0 0 δa i , with δ = δk i δa i ∈ R 2 such that δ ∞ ≤ 1 .
It is a standard elementary uncertainty set (the leaves) of the form

{ bdiag i (δ i I n i ) , with δ = [δ i ] ∈ R r such that ∞ ≤ 1 }
representing parametric uncertainties. The dissipativity property of the uncertainty set can then be chosen of the form

Φ 11 Φ 12 Φ * 12 Φ 22 = -D G G * D where D = bdiag i (D i ), with D i = D * i > 0, and G = bdiag i (G i ), with G i = -G * i .
This corresponds to the usual D-G scalings of the µ-analysis ( [START_REF] Fan | Robustness in the presence of mixed parametric uncertainty and unmodeled dynamics[END_REF]). The L scaling was introduced in [START_REF] Scorletti | Improved efficient analysis for systems with uncertain parameters[END_REF] to reduce the conservatism in the case of r ≥ 2. Indeed, it is possible to represent the branches T i with a non standard uncertainty set ∆ (and the appropriate interconnection matrix M ) of the form

{ δ ⊗ I , with δ = [δ i ] ∈ R r such that δ ∞ ≤ 1 } .
The dissipativity property of the uncertainty set can then be chosen of the form

Φ 11 Φ 12 Φ * 12 Φ 22 = -D + jL G G * D i with L =       0 V 1,2 . . . V 1,r -V 1,2 0 . . . . . . . . . . . . . . . V r-1,r -V 1,r . . . -V r-1,r 0       (8)
where

D = bdiag i (D i ), with D i = D * i > 0, G = [. . . G i . . .], with G i = -G * i , and V i,j = V T i,j
are real matrices. The elementary dissipativity properties of the uncertainty set are thus chosen of the form [START_REF] Feron | Analysis and synthesis of robust control systems via parameter-dependent lyapunov functions[END_REF].

As for the dissipativity properties of the PLL itself, we chose:

• a conic sector alone (for comparison with the result obtained in [START_REF] Dinh | Embedding of uncertainty propagation: Application to hierarchical performance analysis[END_REF]) or with;

• 4 bands (vertical, horizontal, and with a slope of +/-45 deg):

Y ∈ { 1, j, 1 + j, 1 -j }.
This choice has been made a priori, without particular knowledge on a PLL frequency response.

Network performance

The network performance is measured by its frequency response magnitude. The performance dissipativity property is thus chosen of the form

-I 0 0 γ 2 I

Results

Individual PLL For illustration purpose, Figure 4 displays the obtained dissipativity properties of a PLL viewed as system embeddings for different frequencies. The red circle (the red star is its center) and lines represent the embeddings where as the green stars and purple circles represents the systems for some values of the uncertainties.

Network performance

We are now interested in the performance of the PLLs network displayed in Figure 5 while Table 1 displays the characteristics of the different analysis results (the number in brackets for the hierarchical analysis columns is a comparison with the µ-analysis results).

All the analysis results reveal that the PLL network is able to track a ramp since the slope of the max magnitude at low frequencies is +40 dB/dec. Table 1 illustrates the trade-off between conservatism and computation time that can be set by the user with the hierarchical analysis approach: when using the conic sector alone, the result is conservative but is obtained really quickly; when using the conic sector with the bands, the result is much less conservative but is obtained in much more time. For this last hierarchical set up, the 0.998 difference in the maximal peak value with µ-analysis is +0.1 dB, that corresponds to 1.2 % of ratio, which is negligible; the result was obtained in 60 % of the time needed for µ-analysis.

1 1.002 1.004 1. 006 1.008 -8 -7 -6 -5 -4 -3 -2 -1 0 1 x 10 -3 Real
It should be noted that the computation time decrease of the proposed hierarchical analysis compared to the classical µ-analysis is expected even higher for the large scale system with bigger number of the subsystem. The difference (in the dB-scale) between the maximum gain given by µanalysis and by proposed hierarchical approaches is represented in the Fig. 6. It shows that the peak is not attained at the same frequency.

For most frequencies, µ-analysis performs better in terms of conservatism. Surprisingly, hierarchical analysis performs better in the frequency range (160,300) rad/s with a less demanding computation load at the same time.

Conclusion

In this paper, a hierarchical analysis approach has been proposed for the performance of uncertain large-scale systems. It relies on the propagation of dissipativity properties of subsystems through an interconnection; this propagation result is recursively applied leading to a multi steps analysis. The aim is to propose a trade-off adapted to these large-scale systems when a one-step approach as µ-analysis can lead to a large computation time. A numerical example on a PLL network illustrated the new achieved trade-off. • find other dissipative properties that can be used. We think to a cone as proposed in [START_REF] Tits | Robustness under bounded uncertainty with phase information[END_REF] for instance;

• further assess the achieved trade-off for other examples, especially MIMO ones;

• assess the evolution of the achieved trade-off as a function of the dissipativity properties used.

Another direction is to use differently the propagation result. It is used here in a multi steps approach; it could also be used to lead to a one step approach. This results in a less conservative approach but with more demanding computation. This trade-off can also be assessed.

[28] G. Zames. Feedback and optimal sensitivity: model reference transformations, multiplicative seminorms, and approximate inverses. IEEE Trans. Aut. Control, AC-26(2):301-320, April 1981.

A Proof of Lemma 2.1

Let H 1 and H 2 be {X(jω), Y (jω), Z(jω)} dissipative. We prove that, for any λ ∈ (0, 1), Then we have to prove that

λH 1 + (1 -λ)H 2 is {X(jω), Y (
f (λ) = (α + β -(γ + γ * )) λ 2 + (-2β + (γ + γ * )) λ + β with α = z 1 (jω) w(jω) * X(jω) Y (jω) Y (jω) * Z(jω) z 1 (jω) w(jω) , β = z 2 (jω) w(jω) * X(jω) Y (jω) Y (jω) * Z(jω) z 2 (jω) w(jω) , γ = z 1 (jω) w(jω) * X(jω) Y (jω) Y (jω) * Z(jω) z 2 (jω) w(jω) .
is strictly positive on (0, 1). f (λ) is a real-valued polynomial with real coefficients of degree 2. Note that, as w(jω) is non null otherwise z(jω) would also be null, we have f (0) = β > 0 and f (1) = α > 0. The only case where f (λ) could be negative for some value of λ is when (α + β -(γ + γ * )) > 0, df dλ (0) < 0 and df dλ (1) > 0. However

df dλ (λ) = 2 (α + β -(γ + γ * )) λ + (-2β + (γ + γ * )) .
The minimum of f is then attained for

λ 0 = 2β -(γ + γ * ) 2 (α + β -(γ + γ * )) . Then 4 (α + β -(γ + γ * )) f (λ 0 ) = 3 (-2β + (γ + γ * )) 2 + 4 (α + β -(γ + γ * )) β
is strictly positive as (α + β -(γ + γ * )) and β are strictly positive. Thus f (λ) is strictly positive on (0, 1).

B Proof of Theorem 2.1

It is a corollary of Lemma B.1 for stability and Lemma B.2 for dissipativity by noticing that condition [START_REF] Jakubovic | The S-procedure in nonlinear control theory[END_REF] implies condition [START_REF] Houpis | Quantitative feedback theory: fundamentals and applications[END_REF].

Lemma B.1 The uncertain system ∆ M is stable if and only if:

1. there exists ∆ 0 ∈ ∆ such that the system ∆ 0 M is stable;

and for any ω ∈ R M (jω)

I *       -Φ 22 (jω) 0 -Φ 12 (jω) * 0 0 X(jω) 0 Y (jω) -Φ 12 (jω) 0 -Φ 11 (jω) 0 0 Y (jω) * 0 Z(jω)       M (jω) I > 0. (12)
C Proof of Lemma B.1

Non-causal stability We begin by proving that the non-causal stability of ∆ M is equivalent to condition 2. As M is stable, we only need to prove the non-causal stability of the feedback ∆ A. By definition, we have: for any ∆ and ω det (I -A(jω)∆(jω)) = 0.

We prove it by contradiction. So let us assume that condition 2 is verified but there exists ∆ c and ω c such that det (I -A(jω c )∆ c (jω c )) = 0.

Equivalently, there exists a non null q c (jω c ) such that (I -A(jω c )∆ c (jω c )) q c (jω c ) = 0.

Let p c (jω c ) = ∆ c (jω c )q c (jω c ), then q c (jω c ) = A(jω c )p c (jω c ). By definition of dissipativity of ∆, we have:

p c (jω c ) q c (jω c ) * Φ 11 (jω c ) Φ 12 (jω c ) Φ 12 (jω c ) * Φ 22 (jω c ) p c (jω c ) q c (jω c ) > 0.
Post et pre multiplying of (10) by p c (jω c ) yields

q c (jω c ) p c (jω c ) * -Φ 22 (jω c ) -Φ 12 (jω c ) * -Φ 12 (jω c ) -Φ 11 (jω c ) q c (jω c ) p c (jω c ) > 0
which is a contradiction.

Necessity It is evident from the non-causal stability equivalence.

Sufficiency As M is stable, we only need to prove the stability of the feedback ∆ A.

We show it by contradiction. So let us assume that there exists ∆ u ∈ ∆ such that ∆ u A is unstable while conditions 1 and 2 are met. We show that this assumption with condition 1 leads to the existence of ∆ c and ω c such that det (I -A(jω c )∆ c (jω c )) = 0.

Which is a contradiction from the non-causal stability proof.

As ∆ u and A are stable, applying the Nyquist criterion, it shows that the Nyquist curve defined by det (I -A(jω)∆ u (jω)) , ω ∈ R encircles the origin in the complex plan. At the same time, as ∆ 0 A is stable, the curve det (I -A(jω)∆ 0 (jω)) , ω ∈ R does not encircle the origin. Now, as ∆ is a connected set, it is possible to find a continuous path inside the set that link ∆ 0 (jω) to ∆ u (jω), that is, a continuous function Ψ ω : satisfying constraints [START_REF] Iwasaki | Well-posedness of feedback systems: Insights into exact robustness analysis and approximate computations[END_REF] and [START_REF] Jakubovic | The S-procedure in nonlinear control theory[END_REF]. For convenience of writing, we drop the dependency to jω in this part of the proof. By definition: for any ∆ ∈ ∆, the uncertain system

[0, 1] -→ ∆(jω) such that Ψ ω (0) = ∆ 0 (
∆ M is {X(jω), Y (jω), Z(jω)} dissipative if p w * C D 0 I * X Y Y * Z C D 0 I p w > 0 (13) such that p = ∆ q q = A B p w
This last equality can be rewritten as

I -∆ I 0 A B p w = 0. (14)
From Finsler's lemma [START_REF] Boyd | Linear Matrix Inequalities in Systems and Control Theory[END_REF], condition (13) holds for p * w * * defined by ( 14) if and only if there exists τ such that2 :

C D 0 I * X Y Y * Z C D 0 I + τ I 0 A B * I -∆ * I -∆ I 0 A B > 0.
That is: there exists τ such that for any ∆ ∈ ∆,

     C D 0 I I 0 A B      *      X Y Y * Z 0 0 τ I -∆ * I -∆           C D 0 I I 0 A B      > 0.
Let µ 1 µ 2 be such that

µ 1 µ 2 ⊥ =      C D 0 I I 0 A B     
, then, by applying Finsler's lemma, we get the equivalent condition: there exists τ and η such that for any ∆ ∈ ∆,

     X Y Y * Z 0 0 τ I -∆ * I -∆      + η µ * 1 µ * 2 µ 1 µ 2 > 0
It is equivalent by Schur's lemma to: there exists τ and η such that for any ∆ ∈ ∆,

         X Y Y * Z + ηµ * 1 µ 1 > 0 τ I -∆ * I -∆ + ηµ * 2 µ 2 -ηµ * 2 µ 1 X Y Y * Z + ηµ * 1 µ 1 -1 ηµ * 1 µ 2 > 0
Thus there exists τ , η and > 0 such that for any ∆ ∈ ∆,

τ I -∆ * I -∆ + ηµ * 2 µ 2 -ηµ * 2 µ 1 X Y Y * Z + ηµ * 1 µ 1 -1 ηµ * 1 µ 2 -I > 0 Let us define Φ = ηµ * 2 µ 2 -ηµ * 2 µ 1 X Y Y * Z + ηµ * 1 µ 1 -1
ηµ * 1 µ 2 -I.

Then, using Finsler's lemma, the condition there exists τ such that for any ∆ ∈ ∆, τ I -∆ * I -∆ + Φ > 0 is equivalent to for any ∆ ∈ ∆,

∆ I * Φ ∆ I > 0.
That is ∆ is {Φ 11 , Φ 12 , Φ 22 } dissipative.

For the remaining part, let us notice that

-Φ + ηµ * 2 µ 2 -ηµ * 2 µ 1 X Y Y * Z + ηµ * 1 µ 1 -1 ηµ * 1 µ 2 > 0.
Then by Schur's lemma, it is equivalent to Now, the uncertain system ∆ M is {X(jω), Y (jω), Z(jω)} dissipative if and only if condition ( 13) is verified for the set Λ. Due to the inclusion, this is implied by condition [START_REF] Korniienko | Control law design for distributed multi-agent systems[END_REF] being verified for the set Γ. Using S-procedure [START_REF] Jakubovic | The S-procedure in nonlinear control theory[END_REF][START_REF] Boyd | Linear Matrix Inequalities in Systems and Control Theory[END_REF], it is the case if This is exactly condition (12) after factorization.

E Direct Additive Uncertainty Embedding

There exist:

• a matrix G nom ∈ C nz×nw ;

• invertible matrices W o ∈ C nz×nz and W i ∈ C nw×nw ; such that the uncertain system ∆ M is included in the set G u defined

G u = { G nom + W o ∆ u W i , σ(∆ u ) < 1 }
if and only if there exists a solution to one of the two following problems:

1. First: there exist

• a matrix P nom ∈ C nz×nw ;

• a positive definite matrix P o ∈ C nz×nz ;

• a positive definite matrix P i ∈ C nw×nw ;

• a matrix Φ ∈ C (no+n i )×(no+n i ) ; such that ( 15) and ( 16

Assumption 2 . 1 ∆Figure 1 :

 211 Figure 1: Uncertain system

Figure 2 :

 2 Figure 2: Uncertain linear large scale system

Corollary 3 . 1

 31 Problem 3.1 is solved by the following optimization problem: find (Φ 11 , Φ 12 , Φ 22 ) ∈ Φ and X, Y and Z such that

ImaginaryFigure 4 :

 4 Figure 4: Dissipativity properties of a PLL viewed as embeddings

Figure 5 :

 5 Figure 5: Performance of the PLLs network

  Difference of achieved performance between hierarchical analysis and µ-analysis

Figure 6 :

 6 Figure 6: Difference in achieved perfromance between hierarchical analysis and µ-analysis

  jω), Z(jω)} dissipative. Let z(jω) * w(jω) * be non null such that z(jω) = H(jω)w(jω) and define z 1 (jω) = H 1 (jω)w(jω) and z 2 (jω) = H 2 (jω)w(jω). Let us also denotef (λ) = z(jω) w(jω) * X(jω) Y (jω) Y (jω) * Z(jω)z(jω) w(jω) .

2 Necessity

 2 jω) and Ψ ω (1) = ∆ u (jω). As the determinant is a continuous function, the function which at ∆(jω) associates det (I -M (jω)∆(jω)) composed with Ψ ω is a continuous function of λ ∈ [0, 1]. Thus there exists ∆ c and ω c such that det (I -A(jω c )∆ c (jω c )) = 0. D Proof of Lemma B.The necessity is proved by construction of Φ(jω) = Φ 11 (jω) Φ 12 (jω) Φ 12 (jω) * Φ 22 (jω)

= 0

 0 constraint[START_REF] Jakubovic | The S-procedure in nonlinear control theory[END_REF] after rearrangement.Sufficiency It follows from the interpretation of a dissipativity property in terms of inclusion of sets. The {Φ 11 (jω), Φ 12 (jω), Φ 22 (jω)} dissipativity of ∆(jω) is equivalent to p(jω) q(jω) = 0 ∃∆(jω), p(jω) = ∆(jω)q(jω) ∃∆(jω), p(jω) = ∆(jω)q(jω), q(jω) = [ A(jω) B(jω) ]

2 .

 2 P o = W o W * o , P i = W *i W i and P nom = G nom . Second: there exist• a matrix P nom ∈ C nz×nw ;• a positive definite matrix P o ∈ C nz×nz ;• a positive definite matrix P i ∈ C nw×nw ;• a matrix Φ ∈ C (no+n i )×(no+n i ) ; such that (15) and (17) are satisfiedP o = W - * o W -1 o , P i = W * i W i and P nom = P o G nom .

  Problem 3.2 Let Φ be a set such that for any (Φ 11 , Φ 12 , Φ 22 ) ∈ Φ, the uncertainty set ∆ is {Φ 11 , Φ 12 , Φ 22 } dissipative. From the set Φ, find X, Y and Z such that: 1. the uncertain system ∆ M is {X, Y, Z} dissipative; 2. they minimize max ∆∈∆ max w =1 vol( P ) 2 . Problem 3.2 is solved by the following optimization problem: find (Φ 11 , Φ 12 , Φ 22 ) ∈ Φ and X, Y and Z that minimize log det

	Theorem 3.1

  3.3 Let Φ be a set such that for any (Φ 11 , Φ 12 , Φ 22 ) ∈ Φ, the uncertainty set ∆ is {Φ 11 , Φ 12 , Φ 22 } dissipative. Let Y be a matrix of appropriate dimension.From the set Φ and Y , find Z 1 and Z 2 such that:1. the uncertain system ∆ M is {0, Y, Z 1 } dissipative; 2. the uncertain system ∆ M is {0, -Y, Z 2 } dissipative;3. they minimize max ∆∈∆ max w =1 d Y .

	Theorem 3.2 Problem 3.3 is solved by the following optimization problem: find (Φ 1 11 , Φ 1 21 , Φ 1 22

Table 1 :

 1 Characteristics of the analysis results

		4 dB)	6.2 dB (+0.1 dB)	6.1 dB
	Computation time	72 sec (6 %)	767 sec (60 %)	1279 sec

  2. there exists 3 transfer functions Φ 11 (jω), Φ 12 (jω) and Φ 22 (jω) of RL ∞ , with Φ 11 (jω) = Φ 11 (jω) * and Φ 22 (jω) = Φ 22 (jω) * , such that: the uncertainty set ∆ is {Φ 11 (jω), Φ 12 (jω), Φ 22 (jω)} dissipative The uncertain system ∆ M is non-causally stable if and only if merely the condition 1 holds. The uncertain system ∆ M is {X(jω), Y (jω), Z(jω)} dissipative if and only if there exists 3 transfer functions Φ 11 (jω), Φ 12 (jω) and Φ 22 (jω) of RL ∞ , with Φ 11 (jω) = Φ 11 (jω) * and Φ 22 (jω) = Φ 22 (jω) * , such that: the uncertainty set ∆ is {Φ 11 (jω), Φ 12 (jω), Φ 22 (jω)} dissipative

						(9)
	and for any ω ∈ R
	A(jω) I	*	-Φ 22 (jω) -Φ 12 (jω) * -Φ 12 (jω) -Φ 11 (jω)	A(jω) I	> 0.	(10)
	Lemma B.2

The term propagation is kept in reference to[START_REF] Safonov | Propagation of conic model uncertainty in hierarchical systems[END_REF] even if the meaning is slightly different.

In fact, τ should depend on ∆, that is τ ∆ . As shown in[START_REF] Chou | Stability multipliers and µ upper bounds: connections and implications for numerical verification of frequency domain conditions[END_REF], it can be used a continuous function τ (∆) on the closure of ∆. But as ∆ is bounded, it can be selected independent of ∆ (take the maximum on the closure of ∆). This fact will be used several times.