Jean-Pierre Dussault
email: pierre.dussault@usherbrooke.ca.

Jean Charles Gilbert
email: jean-charles.gilbert@inria.fr

Exact computation of an error bound for the balanced linear complementarity problem with unique solution

Keywords: Balanced linear complementarity problem, complexity, data bitlength, error bound, extreme diagonal matrix, matrix inverse norm, P-matrix, rowwise convex combination of matrices, separable function, strong duality 15A09, 15A60, 49N15, 65F35, 65Y20, 90C33, 90C46, 90C60

This paper considers the balanced form of the standard linear complementarity problem with unique solution and provides a more precise expression of an upper error bound discovered by Chen and Xiang and published in 2006. This expression has at least two advantages. It makes possible the exact computation of the error bound factor and it provides a satisfactory upper estimate of that factor in terms of the data bitlength when the data is formed of rational numbers. Along the way, we show that, when any rowwise convex combination of two square matrices is nonsingular, the ℓ ∞ norm of the inverse of these rowwise convex combinations is maximized by an extreme diagonal matrix.

Introduction

Error bounds play a prominent role in the analysis of mathematical problems and the algorithms to solve them, in particular in numerical optimization [START_REF] Pang | Error bounds in mathematical programming[END_REF]. This paper focuses on error bounds discovered by Chen and Xiang [8; 2006] for the linear complementarity problem with a P-matrix and simplifies the expression of its upper factor. The paper also deduces some consequences of this new expression.

In its standard form [START_REF] Cottle | The Linear Complementarity Problem[END_REF], the linear complementarity problem (LCP) reads

0 x ⊥ (M x + q) 0, (1.1)
where the unknown is x ∈ R n (the set of real vectors with n components), while M ∈ R n×n (the set of real matrices of order n) and q ∈ R n are data. Inequalities on vectors must be understood componentwise (for example x 0 in (1.1) means x i 0 for all i ∈ [1 : n], the set of the first n integers). The compact writing of the problem in (1.1) means that one has to find a vector x in R n + := {x ∈ R n : x 0} such that M x + q 0 and x T (M x + q) = 0 (the superscript " T " is used to denote vector or matrix transposition).

A matrix M ∈ R n×n is said to be a P-matrix if all its principal minors are positive (i.e., the determinant det M II > 0, for all I ⊆ [1 : n]; by convention det M ∅∅ = 1). One denotes by P the class of P-matrices. It is known that problem (1.1) has a unique solution, whatever q is, if and only if M ∈ P [32; 1958]. There are many other characterizations of the P-matricity [START_REF] Cottle | The Linear Complementarity Problem[END_REF], including algorithmic ones [START_REF] Ben Gharbia | An algorithmic characterization of P-matricity[END_REF][START_REF] Ben Gharbia | An algorithmic characterization of P-matricity II: adjustments, refinements, and validation[END_REF].

For the sake of generality and for taking advantage of its symmetric formulation, which allows us to shorten some proofs, this paper considers an LCP in a slightly more general form than (1.1), namely 0 (Ax + a) ⊥ (Bx + b) 0, (

where A, B ∈ R n×n and a, b ∈ R n are the data. By comparison with problem (1.1), we name this problem the balanced LCP. This is a special case of the so-called (extended) vertical LCP, which uses more than two matrices and vectors in its formulation (see [START_REF] Cottle | A generalization of the linear complementarity problem[END_REF][START_REF] Sznajder | The generalized order linear complementarity problem[END_REF]). Throughout this work, we assume that problem (1.2) has a unique solution x.

An error bound associated with a set S is an estimate of the distance to S by quantities that are easier to evaluate than this distance, usually those that are used to define the set. The set considered in this paper is the solution set of the LCP (1.2), which has been said to be reduced to the singleton {x}, while the quantity used to estimate the distance to x is defined as follows.

Let • denote an arbitrary norm on R n . The natural residual [START_REF] Luo | Error bound and convergence analysis of matrix splitting algorithms for the affine variational inequality problem[END_REF][START_REF] Mangasarian | New improved error bounds for the linear complementarity problem[END_REF] associated with the linear complementarity problem (1.2) is the function r : R n → R n whose value at x ∈ R n is given by r(x)

:= min(Ax + a, Bx + b), (1.3)
where the minimum operator "min" acts componentwise (for two vectors

u, v ∈ R n and i ∈ [1 : n]: [min(u, v)] i = min(u i , v i )). It is clear that x solves (1.
2) if and only if r(x) = 0, since min(Ax + a, Bx + b) = 0 if and only if, for all i ∈ [1 : n], (Ax + a) i 0, (Bx + b) i 0 and either (Ax + a) i or (Bx + b) i vanishes. Therefore r(x) is a possible measure of the proximity of x to x. In this paper, we consider error bounds of the form

∀ x ∈ R n : β r(x) x -x β r(x) ,
where β and β are positive constants (independent of x), that we call the lower and upper error bound factors, respectively. Error bounds for the LCP have been the subject of many contributions, see [START_REF] Robinson | Some continuity properties of polyhedral multifunctions[END_REF][START_REF] Pang | Inexact Newton methods for the nonlinear complementarity problem[END_REF][START_REF] Mathias | Error bounds for the linear complementarity problem with a P-matrix[END_REF][START_REF] Luo | Error bound and convergence analysis of matrix splitting algorithms for the affine variational inequality problem[END_REF][START_REF] Mangasarian | New improved error bounds for the linear complementarity problem[END_REF][START_REF] Luo | New error bounds for the linear complementarity problem[END_REF][START_REF] Chen | Computation of error bounds for P-matrix linear complementarity problems[END_REF][START_REF] Chen | Perturbation bounds of P-matrix linear complementarity problems[END_REF][START_REF] Dai | Error bounds for linear complementarity problems of DB-matrices[END_REF][START_REF] Li | Some new error bounds for linear complementarity problems of H-matrices[END_REF][START_REF] Chen | Error bounds for linear complementarity problems of MB-matrices[END_REF] for entry points.

For P and Q ∈ R n×n , we define

[P, Q] := {X ∈ R n×n : P X Q},
where the inequalities act again componentwise (i.e., P X Q means P ij X ij Q ij for all i, j ∈ [1 : n]). Hence, for the identity matrix I, [0, I] is a compact notation for the set of diagonal matrices with diagonal elements in the interval [0, 1]. Note also that the set of extreme points of [0, I], denoted by ext[0, I], is the set of diagonal matrices with diagonal elements in {0, 1} (see [31; p. 162] for the definition of an extreme point of a convex set; one can use [14; proposition 2.12] for a meticulous proof of this claim).

For D ∈ [0, I], we denote by

C D := (I -D)A + DB (1.4)
the rowwise convex combination of the matrices A and B ∈ R n×n . The LCP (1.2) has a unique solution whatever the vectors a and b are, if and only if [START_REF] Aganagić | Newton's method for linear complementarity problems[END_REF][START_REF] Pang | Newton's method for B-differentiable equations[END_REF][START_REF] Sznajder | Generalizations of P 0 -and P -properties; extended vertical and horizontal linear complementarity problems[END_REF][START_REF] Zhang | Global error bounds for the extended vertical LCP[END_REF]]

∀ D ∈ [0, I] : C D is nonsingular. (1.5)
Then, the following lower and upper error bounds hold [START_REF] Chen | Computation of error bounds for P-matrix linear complementarity problems[END_REF][START_REF] Zhang | Global error bounds for the extended vertical LCP[END_REF]:

∀ x ∈ R n : max D∈[0,I] C D -1 r(x) x -x max D∈[0,I] C -1 D r(x) , (1.6)
where • denotes a norm on R n and the induced matrix norm.

In this paper, we are interested in giving more precision on the way the lower and upper error bound factors appearing in (1.6) can be computed, when the ℓ ∞ norm is used. If the lower bound factor is easy to evaluate (see section 3.1), the upper bound factor

β := max D∈[0,I] C -1 D ∞ (1.7)
raises more difficulty. This concern makes perfect sense because, as far as we know, this upper error bound factor is the best one obtained so far for the LCP (1.1) with M ∈ P; in particular, it is smaller, hence better, than the one of Mathias and Pang [START_REF] Mathias | Error bounds for the linear complementarity problem with a P-matrix[END_REF] (see [8; theorem 2.3]). We shall show that the evaluation of β can be simplified since one has

β = max D∈ext[0,I] C -1 D ∞ . (1.8)
This extends to higher dimensions the simple observation that, when n = 1, the map

D ∈ [0, 1] → C -1 D ∞ is monotone, so that it attains its maximum on [0, 1] at a point in {0, 1}. For n > 1, however, D kk ∈ [0, 1] → C -1
D ∞ may be nonmonotone, so that an analysis along this line is troublesome. Furthermore, this map can be neither convex nor concave (see [START_REF] Dussault | Exact computation of an error bound for a generalized linear complementarity problem with unique solution -The full report[END_REF]). For these reasons, we shall present a specific, rather long and indirect, proof of (1.8). The simplification (1.8) of (1.7) may look minor at first glance, but it may be interesting for reasons that are discussed in section 4: it provides a way of computing β exactly, it simplifies its computation for small n and it may be crucial for giving an upper estimate of β in terms of the data bitlength in some complexity analysis.

The paper is organized as follows. The next section presents two results that will play an important role in getting the expression (1.8) of β: the first one deals with the norm of a matrix inverse and the second deals with min-max duality in optimization. Section 3 is dedicated to the proof of (1.8). We conclude by some thoughts on complexity issues.

This paper is an abridged version of the more detailed report [START_REF] Dussault | Exact computation of an error bound for a generalized linear complementarity problem with unique solution -The full report[END_REF].

Notation

The unit closed ball associated with a norm • is denoted by B := {x : x 1} and the unit sphere by ∂B := {x : x = 1}.

Preliminaries

This section presents two results that will play a major part in our strategy to get the desired result in section 3. The first one (lemma 2.1) gives an expression of A -1 , for a nonsingular matrix A ∈ R n×n , in terms of an optimization problem. Consequences of this expression for the ℓ ∞ norm are given in corollary 2.2 and in the technical lemma 2.3. The second result (lemma 2.4) highlights conditions to have strong duality on a product space X × [1 : p] for a pairing function that has a separable property.

Norm of a matrix inverse

For a given nonsingular matrix function z ∈ R p → A(z) ∈ R n×n , analyzing the map [A(•)] -1 is often more difficult than analyzing [A(•) . It is possible, however, to toggle from one map to the other thanks to the following lemma (formula (2.1a) can be found in [18; problem 5.6.P46]; a proof of the lemma is given in [START_REF] Dussault | Exact computation of an error bound for a generalized linear complementarity problem with unique solution -The full report[END_REF]).

Lemma 2.1 (norm of a matrix inverse) If A ∈ R n×n is a nonsingular square matrix and if • denotes a vector norm and its induced matrix norm, then

min v =1 Av = A -1 -1 .
(2.1a)

In addition, v solves the problem in the left-hand side of (2.1a) if and only if w := A -1 Av solves the problem in the left-hand side of

max w =1 A -1 w = A -1 . (2.1b)
In the sequel, the infinity vector and its induced matrix norms, both denoted by • ∞ , are used. For this reason, we consider this case in corollary 2.2 below and bring some precision. We denote by e i the ith basis vector of R n and set e := i∈[1 : n] e i , which is the vector of all ones. By definition and computation [17; § 5.6.5] (see also (2.7a)-(2.7d) in the proof of corollary 2.2 below), for a matrix A ∈ R n×n , one has

A ∞ := max w ∞=1 Aw ∞ = max i∈[1 : n] A i : 1 , (2.2)
where A i : := (e i ) T A denotes the ith row of A and

v 1 := i∈[1 : n] |v i | denotes the ℓ 1 -norm of v ∈ R n .
We also denote by "sign" the maximal monotone multifunction R ⊸ R that is the subdifferential of the absolute value function: it associates with t ∈ R the following set of R:

sign t :=    {-1} if t < 0 [-1, 1] if t = 0 {1} if t > 0. (2.3)
One finds other definitions of sign(0), in particular to make the map "sign" a single-valued function, but our choice of definition is important for the sequel, like in the formulas (2.4b) below. Recall that • 1 is the dual norm of • ∞ with respect to the Euclidean scalar product, which means that

v 1 = max w ∞=1 v T w = max w ∞=1 |v T w|. (2.4a)
The solution sets of these maximum problems are Arg max

w ∞ =1 v T w = (sign v) ∩ ∂B ∞ and Arg max w ∞=1 |v T w| = (± sign v) ∩ ∂B ∞ , (2.4b)
where, for a vector

v ∈ R n , sign v := (sign v 1 )ו • •×(sign v n ) ⊆ R n (hence sign 0 R n = B ∞ ), ± sign v := (sign v) ∪ (-sign v)
and the boundary ∂B ∞ of B ∞ is present only to deal with the case where v = 0. For a nonsingular square matrix A, we adopt the following notation

W ∞ (A) := Arg max w ∞=1 A -1 w ∞ = {w ∈ ∂B ∞ : A -1 w ∞ = A -1 ∞ }, (2.5a) V ∞ (A) := Arg min v ∞ =1 Av ∞ = {v ∈ ∂B ∞ : Av ∞ = A -1 -1 ∞ }. (2.5b)
The second equality in (2.5a) comes from the definition of the induced matrix norm

• ∞ in (2.
2), while the second equality in (2.5b) is deduced from the identity (2.1a). The next corollary gives other expressions of these sets.

Corollary 2.2 (ℓ ∞ -norm of a matrix inverse) Suppose that A ∈ R n×n is a non- singular matrix. Set β := A -1 ∞ and α := 1/β. Then, W ∞ (A) = ± sign(A -T e i ) : (A -1 ) i : 1 = β , (2.6a) V ∞ (A) = α A -1 W ∞ (A) . (2.6b) Proof. [(2.6a)] Observe first that β = max w ∞=1 A -1 w ∞ [definition of the matrix norm • ∞ ] (2.7a) = max w ∞=1 max i∈[1 : n] |(e i ) T A -1 w| [definition of the vector norm • ∞ ] (2.7b) = max i∈[1 : n] max w ∞=1 |(e i ) T A -1 w| [the max's commute] (2.7c) = max i∈[1 : n] A -T e i 1 [(2.4a)]. (2.7d)
We can now establish the identity (2.6a).

[⊆] If w ∈ W ∞ (A), w solves the problem in (2.7a)-(2.7b), by definition. Let ī ∈ [1 : n] be a solution to the inner problem max{|(e i ) T A -1 w| : i ∈ [1 : n]} appearing in (2.7b).

Then, the pair ( w,ī) maximizes the map (w, i)

∈ ∂B ∞ × [1 : n] → |(e i ) T A -1 w|.
It follows that ī solves to the problems in (2.7c)-(2.7d) and w is a solution to the inner problem max{|(e ī) T A -1 w| : w ∞ = 1} appearing in (2.7c). Hence, by (2.4b), w ∈ ± sign(A -T e ī) and, by (2.7d),

β = A -T e ī 1 = (A -1 ) ī : 1 . [⊇] Suppose now that w ∈ ± sign(A -T e ī) for some ī ∈ [1 : n] satisfying A -T e ī 1 = β.
-By this last identity, ī solves the problem in (2.7d), hence the problem in (2.7c).

-By the nonsingularity of A -T , one component of A -T e ī does not vanish, so that w ∈ ± sign(A -T e ī) ∩ ∂B ∞ . By (2.4b), this implies that w solves the problem max{|(e ī) T A -1 w| : w ∞ = 1}. It results from these last two observations and (2.7a)-(2.7c), that w solves the problem in (2.7a)-(2.7b). We have shown that w ∈ W ∞ (A).

[(2.6b)] This is a consequence of the last claim in lemma 2.1, according to which

v ∈ V ∞ (A) if and only if v = α A -1 w with w ∈ W ∞ (A). ✷
We conclude this section by synthesizing in the following lemma a mechanism that, despite its innocuous appearance, plays a major part in the proof of proposition 3.1 below. As shown in the lemma's proof, this mechanism is only operational when some element of A -1 vanishes, but this fact is revealed indirectly, through a property of a vector v ∈ V ∞ (A).

Lemma 2.3 (technical) Suppose that A ∈ R n×n is nonsingular, that α := A -1 -1 ∞ and that v ∈ V ∞ (A) has the property that |(Av) k | < α for some k ∈ [1 : n]. Then, there exists a v ′ ∈ V ∞ (A) such that (Av ′ ) k = 0.
Proof. Let β := 1/α. Since v ∈ V ∞ (A), the vector defined by w := βAv is in W ∞ (A), by (2.6b). By assumption, (Av) k ∈ (-α, α), so that w k ∈ (-1, 1). These two facts on w and (2.6a) imply that there must be some index i such that w ∈ ± sign(A -T e i ), (A -1 ) i : 1 = β and (A -1 ) ik = 0.

Define the vector w ′ ∈ R n by vanishing the kth component of w:

w ′ i := w i if i = k 0 otherwise.
Then, we also have w ′ ∈ ± sign(A -T e i ), implying that w ′ ∈ W ∞ (A). The sought vector is v ′ := αA -1 w ′ . Indeed, on the one hand, v ′ ∈ V ∞ (A) by (2.6b). On the other hand, Av ′ = αw ′ implying that (Av ′ ) k = 0, as desired. ✷

Strong duality for separable functions

Let be given a set X and p functions ϕ i : X → R. Usually, equality does not hold in the weak duality inequality [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms[END_REF][START_REF] Golshtein | Modified Lagrangians and Monotone Maps in Optimization[END_REF][START_REF] Bonnans | Numerical Optimization -Theoretical and Practical Aspects[END_REF][START_REF] Ch | Fragments d'Optimisation Différentiable -Théorie et Algorithmes. Lecture Notes[END_REF] inf

x∈X max i∈[1 : p] ϕ i (x) max i∈[1 : p] inf x∈X ϕ i (x). (2.8)
Take for example, X = R, p = 2, ϕ 1 (x) = (x + 1) 2 and ϕ 2 (x) = (x -1) 2 , in which case the left-hand side value is 1, while the right-hand side value is 0 (see [15; lemma 4.5] for a way of modifying (2.8) that ensures equality). The situation is very different, more elementary and more favorable, when X is a Cartesian product X = X 1 × • • • × X p of sets X i and each function ϕ i only depends on the ith component x i ∈ X i of x = (x 1 , . . . , x p ) ∈ X; then equality holds above with some other interesting properties. This particular situation, which occurs below, is analyzed in the next lemma. In this one, the problems

inf x∈X max i∈[1 : p] ϕ i (x i ) and max i∈[1 : p] inf x i ∈X i ϕ i (x i ).
are called the primal and dual problems, respectively. A primal (resp. dual) solution is a solution to this primal (resp. dual) problem.

The next lemma not only shows the lack of duality gap for the separable case, but also describes the sets of primal and dual solutions. It also shows how to construct a primal solution from the dual solutions, as well as a dual solution from the primal solutions.

Lemma 2.4 (strong duality for separable functions) Let X := X 1 × • • • × X p be the Cartesian product of nonempty sets X i and let ϕ i : X i → R, i ∈ [1 : p], be arbitrary functions. An x ∈ X is written x = (x 1 , . . . , x p ), with x i ∈ X i for i ∈ [1 : p]. 1) (No duality gap) The following identity holds

inf x∈X max i∈[1 : p] ϕ i (x i ) = max i∈[1 : p] inf x i ∈X i ϕ i (x i ).
(2.9a)

Denote by v the common value of the two sides of this identity. 2) (Set of primal solutions) The set of primal solutions is the possibly empty set

X := X1 × • • • × Xp , where Xi := {x i ∈ X i : ϕ i (x i ) v}.
(2.9b)

3) (Set of dual solutions) The set of dual solutions is the nonempty set

Ī := {i ∈ [1 : p] : ϕ i (x i ) v for all x i ∈ X i }. (2.9c) 4) (Saddle-point property)
The following properties are equivalent: (2.9e) 6) (Deducing a dual solution from the primal solutions) Suppose that X = ∅. Then, ī ∈ Ī if and only if, for all x ∈ X, ī

(i) (x,ī) ∈ X × Ī, (ii) (x,ī) is a saddle-point of the map (x, i) ∈ X × [1 : p] → ϕ i (x i ), meaning that ∀ (x, i) ∈ X × [1 : p] : ϕ i (x i ) ϕ ī(x ī) ϕ ī(x ī), (2.9d
maximizes i ∈ [1 : p] → ϕ i (x i ),
Proof. We only prove points 1 and 2, which are those intervening below. For a proof of the other points, see [START_REF] Dussault | Exact computation of an error bound for a generalized linear complementarity problem with unique solution -The full report[END_REF]. 1) By the weak duality property (2.8) and the fact that ϕ i only depends on the ith component of x, the inequality " " holds in (2.9a). Let us prove the reverse inequality. Let ε > 0. For any i ∈ [1 : p], there is an

x ε i ∈ X i such that ϕ i (x ε i ) inf x i ∈X i ϕ i (x i ) + ε. Therefore, max i∈[1 : p] ϕ i (x ε i ) max i∈[1 : p] inf x i ∈X i ϕ i (x i ) + ε. (2.10)
It is here that the separability assumption intervenes. Since the left-hand side of (2.10) is the value at x ε := (x ε 1 , . . . , x ε p ) of the function x = (x 1 , . . . , x p ) ∈ X → max i∈[1 : p] ϕ i (x i ), the following inequality certainly holds

inf x∈X max i∈[1 : p] ϕ i (x i ) max i∈[1 : p] ϕ i (x ε i ).
Combining with (2.10), we get

inf x∈X max i∈[1 : p] ϕ i (x i ) max i∈[1 : p] inf x i ∈X i ϕ i (x i ) + ε.
Since ε > 0 is arbitrary, the inequality " " holds in (2.9a).

2) The point x = (x 1 , . . . , xp ) is a primal solution if and only if

max i∈[1 : p] ϕ i (x i ) inf x∈X max i∈[1 : p] ϕ i (x i ) = v or ∀ i ∈ [1 : p] : ϕ i (x i ) v.
This fact reads x ∈ X, for the given X. ✷

Finitely computable error bounds for the LCP

This section focusses on the exact computation of the lower and upper error bound factors in (1.6). More specifically, its main result, proposition 3.1, shows that a solution to (1.7) can always be found in ext[0, I], the set of extreme points of [0, I].

Computation of the lower error bound factor

Before focusing in section 3.2 on the main objective of this paper, which is the simplification of the upper error bound factor in (1.6), let us mention that the lower error bound factor in (1.6), with the ℓ ∞ norm, namely

max D∈[0,I] C D ∞ -1 , (3.1)
can be easily computed. Observe first that, for two vectors u and v and a vector norm • , one has

max t∈[0,1] (1 -t)u + tv ) = max( u , v ) = max t∈{0,1} (1 -t)u + tv ). (3.2) Next, max D∈[0,I] C D ∞ = max D∈[0,I] max i∈[1 : n] (1 -D ii )A i : + D ii B i : 1 [(1.4) and (2.2)] = max i∈[1 : n] max D ii ∈[0,1] (1 -D ii )A i : + D ii B i : 1 [the max's commute] = max i∈[1 : n] max( A i : 1 , B i : 1 ) [(3.2)] = max( A ∞ , B ∞ ) [the max's commute]. (3.3)
This shows that (3.1) can be easily computed. Note that, in view of (3.3), the maximum in (3.1) is obtained for D ∈ {0, I}.

Computation of the upper error bound factor

By the compactness of [0, I] and the continuity of D ∈ [0, I] → C -1 D ∞ , the maximization problem (1.7), recalled below β := max

D∈[0,I] C -1 D ∞ , (3.4)
has a solution, say D. Since C D is nonsingular, β given by (3.4) is finite and positive. Then, one can also define the positive number

α := β -1 . (3.5)
The goal of this section is to show that the value β can be obtained by restricting the feasible set of problem (3.4) to ext[0, I], the set of extreme diagonal matrices of [0, I]: Before starting the analysis, let us observe that the objective of problem (3.4) is made of the composition of the nonlinear smooth function D → C -1 D and the convex function • ∞ , but this objective is maximized, not minimized, so that the theory developed for the class of composite problems [START_REF] Burke | Optimality conditions for non-finite valued convex composite functions[END_REF]29,[START_REF] Bonnans | Second order optimality conditions based on parabolic second order tangent sets[END_REF] does not apply. For this reason, we provide a specific proof of proposition 3.1. This one is postponed to page 11, after the following preliminary considerations on problem (3.4).

β = max D∈ext[0,I] C -1 D ∞ . ( 3
Part of the analysis is based on the following rewriting of β, defined by (3.4) (some more justifications are given after (3.7e), C D is defined by (1.4)):

max D∈[0,I] C -1 D ∞ = max D∈[0,I] min v ∞ =1 C D v ∞ -1 [(2.1a)] (3.7a) = min D∈[0,I] min v ∞ =1 C D v ∞ -1 (3.7b) = min v ∞ =1 min D∈[0,I] C D v ∞ -1 [the min's commute] (3.7c) = min v ∞ =1 min D∈[0,I] max i∈[1 : n] |(C D v) i | -1 [definition of • ∞ ] (3.7d) = min v ∞ =1 max i∈[1 : n] min D ii ∈[0,1] |(1 -D ii )(Av) i + D ii (Bv) i | -1 , (3.7e)
where we have been able to switch min D and max i from (3.7d) to (3.7e), without duality gap, thanks to point 1 of lemma 2.4 and the fact that

[0, I] = [0, 1] × • • • × [0, 1] (n times) is a Cartesian product and that |(C D v) i | = |(1 -D ii )(Av) i + D ii (Bv) i | only depends on D ii .
In (3.7c), we have a minimum in v (i.e., the infimum is attained), since by (3.7b) the function (D, v) → C D v ∞ has a minimizer ( D, v) on [0, I] × ∂B ∞ , which implies that v solves the problem in (3.7c) (this property of nested optimization problems is discussed around [14; corollary 1.10]). Let us some consequences of the identities in (3.7).

According to (3.4), the value of the left-hand side in (3.7a) is β > 0 and, according to (3.5), the optimal values of the optimization problems inside the parentheses in (3.7b)-(3.7c) is α > 0, so that

α = min D∈[0,I] min v ∞ =1 C D v ∞ , (3.8a)
= min

v ∞ =1 min D∈[0,I] C D v ∞ . (3.8b)
Therefore, one can write

D solves (3.4) ⇐⇒ ∃ v such that ( D, v) solves problems (3.8). (3.9)
We also have Pursuing along the vein that exploits (3.7), we see that the optimal value of the optimization problems inside the parentheses in (3.7d)-(3.7e) is also α > 0, so that, for a v such that ( D, v) solves the problems in (3.8) for some D ∈ [0, I], one has

D solves (3.4) v ∈ V ∞ (C D ) ⇐⇒ ( D, v) solves problems (3.8). (3.10)
α = min D∈[0,I] max i∈[1 : n] |(C D v) i |, (3.11a)
= max i∈[1 : n] min D ii ∈[0,1] |(1 -D ii )(Av) i + D ii (Bv) i |. (3.11b)
We shall also use the following implication:

( D, v) solves problems (3.8) D′ solves (3.11a) =⇒ ( D′ , v) solves problems (3.8) D′ solve (3.4). (3.12)
Indeed, by the left-hand side of the implication,

( D, v) minimizes C D v ∞ on [0, I] × ∂B ∞ and D′ minimizes C D v ∞ on [0, I]. Then, ( D′ , v) minimizes C D v ∞ on [0, I] × ∂B ∞ or,
equivalently, ( D′ , v) solves problems (3.8). Next, D′ solves (3.4), by (3.9). We conclude this preliminary discussion with an elementary lemma [START_REF] Dussault | Exact computation of an error bound for a generalized linear complementarity problem with unique solution -The full report[END_REF].

Lemma 3.2 (elementary) Suppose that ν and µ ∈ R, that α > 0 and that

min δ∈[0,1] |(1 -δ)ν + δµ| = α. (3.13)
Then, the solution set of the optimization problem in (3.13) is {0}, {1} or [0, 1].

Proof of proposition 3.1. Suppose that A and B ∈ R n×n satisfy (1.5) and that D solves the optimization problem in (3.4). Since the last claim of the proposition is clear, we only focus on the first part of it, assuming that Dkk ∈ (0, 1) for some k ∈ [1 : n]. By (3.9)-(3.10), there is a v ∈ V ∞ (C D) such that ( D, v) solves the problems in (3.8). The goal of the proof is now to show that one can replace Dkk by any value in [0, 1], to form a diagonal matrix D′ that is still a solution to (3.4). Sometimes (case 1 below), this goal will be reached with the chosen initial v ; other times (case 2 below), it will be necessary to change the optimal D and v ∈ V ∞ (C D) several times (infinitely often is not excluded, with a limit argument) to reach the goal. Before introducing these cases, we highlight the principal argument that is used in the proof.

Principal argument. Recall that the optimal value of (3.4) is denoted by β := C -1 D ∞ , which is positive, and that α : Therefore, by the principal argument described above, for the considered vector v solving (3.8b), the kth element of the optimal D can be {0}, {1} or any value in [0, 1]. Since Dkk ∈ (0, 1), by asumption, one has [a k , b k ] = [0, 1], which concludes the proof in this case.

= 1/β = C D v ∞ for the chosen v ∈ V ∞ (C D).
2) Case where min

D kk ∈[0,1] |(1 -D kk )(Av) k + D kk (Bv) k | < α. (3.16)
In that case, the interval [a k , b k ] defined by (3.14) is not guaranteed to contain 0 or 1.

By modifying the vector v, however, we show that one can find intervals of substitutes for Dkk , maintaining the optimality of the diagonal matrix, that cover all the interval [0, 1]; this is the desired result. It suffices to extend the interval [a k , b k ] of optimal values for Dkk to the left so that it contains 0, because, by symmetry, the interval [a k , b k ] can then also be extended to the right so that it contains 1 (switch A and B and replace D by I -D).

One can assume that a k > 0, since otherwise there is nothing to prove. This implies that α < A ∞ (because, by optimality of D, one has α

= C D v ∞ C 0 v ∞ = Av ∞ A ∞ and α = A ∞ since otherwise |(Av) k | Av ∞ = α and a k = 0 by (3.14)).
We do this extension by an iterative procedure whose iterates, indexed by j ∈ N, are pairs ( Dj , vj ) verifying ( Dj , vj ) solves the problems in (3.8),

) Dj ii = Dii for i = k, (3.17b) (1 -Dj kk )(Av j ) k + Dj kk (Bv j ) k = 0, (3.17c) 0 < Dj+1 kk (1 -α/(2 A ∞ )) Dj kk . (3.17d) (3.17a
The iterative process is interrupted as soon as 0 is in the interval

[a j k , b j k ] := {D kk ∈ [0, 1] : |(1 -D kk )(Av j ) k + D kk (Bv j ) k | α}, (3.18)
that is, as soon as a j k = 0. It will be clear from the construction of these intervals that their union will be formed of solutions for Dkk . Actually, the reasoning below does not control directly a j k but it controls Dj kk ∈ [a j k , b j k ], which tends to zero by (3.17d).

r Let us determine ( D0 , v0 ) and verify (3.17a)-(3.17c) for j = 0 ((3.17d) for j = 0 will be verified when D1 kk will be determined, in the next point). When (3.16) holds, point 2 of lemma 2.4 ensures that changing Dkk in order to have |(C D v) k | < α will not change the optimality of D, so that we can actually assume that

|(C D v) k | < α. Then, lemma 2.3 with A = C D and v = v ∈ V ∞ (C D ) tells us that one can find a v0 ∈ V ∞ (C D ) such that (C D v0 ) k = 0, which reads (1 -Dkk )(Av 0 ) k + Dkk (Bv 0 ) k = 0.
Therefore, setting D0 := D, we see that (3.17b) and (3.17c) hold. Furthermore, (3.17a) also holds since, by (3.10), the fact that D0 solves (3.4) and that v0 ∈ V ∞ (C D0 ) implies that ( D0 , v0 ) solves the problems in (3.8).

r Let us now show how to construct ( Dj+1 , vj+1 ) from ( Dj , vj ), if this is necessary.

Discussion

The simplification (1.8) of the error bound factor β given by (1.7) allows us to compute it by evaluating the map D ∈ [0, I] → C -1 D ∞ at the 2 n extreme points of [0, I], which are the diagonal matrices D with diagonal entries in {0, 1}. This is an improvement. Nevertheless, for large n, this exponential number of evaluations can make this exact exhaustive computation approach very time consuming. Now, it is not unlikely that, for special classes of matrices, the simplification (1.8) can yield an efficient way of computing the error bound factor. Finally, we are also exploring the possibility to simplify this comprehensive evaluation by a specific algorithm based on the developments made in this paper.

Another interest of the simplified formula (1.8) of β deals with the complexity analysis of some algorithms for solving the balanced linear complementarity problem (1.2) with matrices A and B verifying (1.5) (equivalent to the P-matricity of M if (A, B) = (I, M )) and integer (or rational) data. When the complexity is expressed in terms of the data bitlength, which is the number of bits required to represent the problem data in the computer memory, and when the error bound (1.6) intervenes, the question may arise to know whether the upper error bound factor can be bounded above by a formula using the data bitlength or the bitlength of the matrix A and B, denoted L(A, B) say, since the data bitlength is certainly larger than L(A, B). It is known from [25; paragraph straddling pages 209-210] (probably also implicit in [START_REF] Khachiyan | A polynomial algorithm in linear programming[END_REF]), that, for an arbitrary nonsingular matrix Without (1.8), formula (4.1) could fail to hold and the upper bound of max{ C -1 D ∞ : D ∈ [0, I]} could be in terms of L(C D), which could be infinite since the optimal diagonal matrix D could then have irrational numbers in some entries. Therefore, thanks to (1.8), for the balanced linear complementarity problem (1.2), with matrices A and B verifying (1.5), one has the error bound

M ∈ R n×n , M -1 ∞ n 2 L(M )+1 .
∀ x ∈ R n : x -x ∞ n 2 L(A,B)+1 min(Ax + a, Bx + b) ∞ , (4.2)
where x is the unique solution to the balanced LCP.

)

  (iii) xī minimizes ϕ ī on X ī and ī maximizes ϕ i (x i ) on [1 : p]. 5) (Deducing a primal solution from the dual solutions) Suppose that, for any dual solution ī ∈ Ī, the problem inf{ϕ ī(x ī) : x ī ∈ X ī} has a solution xī , then the primal problem has a solution x ∈ X satisfying xĪ = xĪ and Ī = Arg max i∈[1 : p] ϕ i (x).

. 6 )

 6 Proposition 3.1 (validity of (3.6)) Suppose that A and B ∈ R n×n satisfy (1.5) and that D solves the optimization problem in(3.4). Then, if Dkk ∈ (0, 1) for some k ∈ [1 : n], D remains optimal if Dkk is changed to any value in [0, 1]. particular, the value of β defined by(3.4) is also given by(3.6).

  This is because, when D solves(3.4) and v ∈ V ∞ (C D ) (i.e., v minimizes C Dv ∞ on ∂B ∞ by (2.5b)), ( D, v) solves the problems in (3.8). Reciprocally, when ( D, v) solves the problems in (3.8), then D solves (3.4) by (3.9) and v minimizes C D v ∞ on ∂B ∞ , which also reads v ∈ V ∞ (C D) by (2.5b).

  Now, we want to determine the other values that the kth diagonal element Dkk of the optimal D can take, if any, while keeping the optimality of the resulting diagonal matrix. Here is the mechanism that allows us to change Dkk . By point 2 of lemma 2.4, for the current v and for any value D′ kk taken in the interval[a k , b k ] := {D kk ∈ [0, 1] : |(1 -D kk )(Av) k + D kk (Bv) k | α},(3.14)the diagonal matrix D′ defined byD′ ii := D′ kk if i = k Dii otherwise,is a solution to problem(3.11a). By (3.12), we get that D′ is a solution to problem(3.4). In conclusion, for any v ∈ V ∞ (C D) with D solving (3.4), the interval [a k , b k ] defined by (3.14) is a set of optimal values for Dkk . These intervals depend on v. Our objective is to show that the union of these intervals [a k , b k ] for some well chosen v ∈ V ∞ (C D ) and solutions D to (3.4) is [0, 1] (the reasoning only holds when Dkk ∈ (0, 1)). In case 1 below, one has [a k , b k ] = [0, 1], immediately. In case 2 below, the objective is realized by changing v and D, alternatively, possibly infinitely often.By optimality of D for (3.11a), we have|(C D v) k | α (recall that (C D v) k = (1 -Dkk ) (Av) k + Dkk (Bv) k only depends on Dkk ). Therefore, either min{|(C D v) k | : D kk ∈ [0, 1]} = α or min{|(C D v) k | : D kk ∈ [0, 1]} < α.We now examine these two complementary cases.1) Case where minD kk ∈[0,1] |(1 -D kk )(Av) k + D kk (Bv) k | = α. (3.15) By lemma 3.2, with ν = (Av) k and µ = (Bv) k , the solution set of problem (3.15) is either {0}, {1} or [0, 1]. From (3.14) and (3.15), this solution set is also the interval [a k , b k ].

Thanks to the formula ( 1 . 8 )

 18 of β, the upper error bound factor is equal to C -1 D ∞ , for some D ∈ ext[0, I]. Therefore, the rows of C D defined by (1.4) are those of A or B. As a result, one certainly has L(C D) L(A, B). (4.1) As a result, with D ∈ ext[0, I] solving the optimization problem in (1.8), one has max D∈[0,I]C -1 D ∞ = C -1 D ∞ n 2 L(C D )+1 n 2 L(A,B)+1 .

Acknowledgments

We would like to thank the referees for their helpful comments and suggestions.

Assume that a j k > 0 (otherwise, there is no reason to pursue the iterative process). Then, 0 = (Av j ) k = (Bv j ) k (otherwise a j k = 0) and by definition of a j k in (3.18):

so that (3.17b) is verified with j + 1 instead of j. Adding side by side (3.17c) and (3.19), and using the definition (3.20) of Dj+1 , we get

Subtracting side by side (3.17c) from (3.21), using (Av

which is (3.17d).

We still have to determine

(the last membership comes from the fact that ( Dj+1 , vj ) solves the problems in (3.8) and the implication "⇐" in (3.10)) tells us that one can find a

The first membership implies (3.17a) with j + 1 replacing j, by the implication "⇒" of (3.10) (note that Dj+1 solves (3.4) by the implication "⇐" of (3.10)). The second identity reads (3.17c) with j + 1 replacing j.

By the two previous points, the iterative procedure defining ( Dj , vj ), for j ∈ N, is well defined, unless it is interrupted by the fact that a j k = 0 for some j ∈ N, which is a desirable property since then D is solution to (3. ✷


Dennis Van Der Meer
email: dennis.vandermeer@minesparis.psl.eu

Simon Camal
email: simon.camal@minesparis.psl.eu

Georges Kariniotakis
email: georges.kariniotakis@minesparis.psl.eu

Seamless intra-day and day-ahead multivariate probabilistic forecasts at high temporal resolution

Keywords: Scenarios, operational forecasting, renewable energy sources, photovoltaic generation, virtual power plant

High temporal resolution intra-day and day-ahead photovoltaic (PV) power forecasts are important to maximize the value of PV systems because they enable stakeholders to participate in both the energy and ancillary services markets. Whereas most day-ahead electricity markets feature an hourly temporal resolution, intra-day markets may require forecasts at 5-minute resolution. In addition, battery integration can improve power system management in isolated grids with high PV power penetration, but battery control requires high temporal resolution forecasts. We propose an efficient method based on pattern matching to generate multivariate probabilistic forecasts, approximated by trajectories, at high temporal resolution and without the need to separately forecast the marginals and estimate the covariance matrix. We compare the proposed method against quantile regression forests in combination with copula theory and show that our method reduces the forecast time by approximately 98% and simplifies the modeling chain while incurring a minor performance penalty.

I. INTRODUCTION

Murphy divides the "nature of goodness in weather forecasting" into three aspects: (i) consistency, (ii) quality and (iii) value [START_REF] Murphy | What Is a Good Forecast? An Essay on the Nature of Goodness in Weather Forecasting[END_REF]. In itself, a forecast does not have value until it is used in a decision-making process [START_REF] Murphy | What Is a Good Forecast? An Essay on the Nature of Goodness in Weather Forecasting[END_REF]. In order to maximize the monetary value of a renewable energy source (RES) forecast, it is important that the system owner or system aggregator be able to participate in both the energy and ancillary service markets on the intra-day and day-ahead levels. Electricity markets feature various requirements for parties to place bids; of particular interest to this study are the requirements on temporal resolution and forecast horizon. For instance, the California Independent System Operator (CAISO) requires, every 5 minutes, forecasts with a lead time of 7.5 minutes, a temporal resolution of 5 minutes and a forecast horizon of 65 minutes for their real-time economic dispatch [START_REF] Makarov | Incorporating uncertainty of wind power generation forecast into power system operation, dispatch, and unit commitment procedures[END_REF]. For its short-term unit commitment, CAISO requires forecasts with a lead time of 75 minutes, a temporal resolution of 15 minutes and a forecast horizon of 5 hours, updated every hour [START_REF] Makarov | Incorporating uncertainty of wind power generation forecast into power system operation, dispatch, and unit commitment procedures[END_REF]. For its day-ahead market, CAISO requires hourly forecasts with a lead time of 14 hours and a forecast horizon of 24 hours [START_REF] Makarov | Incorporating uncertainty of wind power generation forecast into power system operation, dispatch, and unit commitment procedures[END_REF]. Typically, research implementing RES forecasts in power system applications focuses on one of the aforementioned Parts of this research were carried in the frame of the Smart4RES project (No. 864337), supported by the Horizon 2020 Framework Program. markets. For instance, [START_REF] Camal | Optimal offer of automatic frequency restoration reserve from a combined pv/wind virtual power plant[END_REF] estimate the optimal automatic frequency restoration reserve (aFRR) offer that RESs can provide on the day-ahead market. Similarly, [START_REF] Camal | Scenario generation of aggregated wind, photovoltaics and small hydro production for power systems applications[END_REF] focus on offers of aFRR and energy on the day-ahead market. In contrast, [START_REF] Yang | Operational solar forecasting for the real-time market[END_REF] concentrate on CAISO's short-term unit commitment market described above.

The value of a forecast also relates to the continuity of balance between supply and demand. Especially in isolated grids featuring high RES penetration, the lack of synchronous machines introduces challenges in power system management [START_REF] Vasconcelos | Advanced control solutions for operating isolated power systems: Examining the portuguese islands[END_REF]. In such a scenario, high temporal resolution RES forecasts are important to inform the grid operator on future variability, as well as to enable battery control, which in turn allows for simultaneous provision of numerous services (e.g., peak shaving or frequency containment reserve) [START_REF] Pinson | Identification of new and enhanced forecasting products[END_REF]. In the context of battery control, [START_REF] Van Der Meer | An alternative optimal strategy for stochastic model predictive control of a residential battery energy management system with solar photovoltaic[END_REF] generate 15-min resolution forecasts up to 24 h ahead that are updated every 15 mins and subsequently fed into a model predictive control algorithm.

Besides value, Murphy defines consistency as follows: "a forecast should always correspond to a forecaster's best judgment" [START_REF] Murphy | What Is a Good Forecast? An Essay on the Nature of Goodness in Weather Forecasting[END_REF]. Since the forecaster faces several uncertainties , his or her judgment is inherently uncertain and the forecast should reflect that [START_REF] Murphy | What Is a Good Forecast? An Essay on the Nature of Goodness in Weather Forecasting[END_REF]. Moreover, sequential optimization problems, e.g., model predictive control, require trajectory forecasts that account for dependencies in time and space1 [START_REF] Pinson | Wind energy: Forecasting challenges for its operational management[END_REF]. The combination of a multitude of trajectory forecasts represents an approximation of the multivariate probabilistic forecast and can readily be employed in Monte-Carlo-type decision-making processes [START_REF] Pinson | Wind energy: Forecasting challenges for its operational management[END_REF].

The de facto method to issue trajectory forecasts is to forecast the marginal cumulative distribution functions (CDFs) over a forecast horizon and compute the probability integral transform (PIT) variables to estimate a covariance matrix that describes the dependence structure, e.g., [START_REF] Camal | Optimal offer of automatic frequency restoration reserve from a combined pv/wind virtual power plant[END_REF], [START_REF] Camal | Scenario generation of aggregated wind, photovoltaics and small hydro production for power systems applications[END_REF], [START_REF] Van Der Meer | An alternative optimal strategy for stochastic model predictive control of a residential battery energy management system with solar photovoltaic[END_REF], [START_REF] Pinson | Wind energy: Forecasting challenges for its operational management[END_REF]. However, as both the temporal resolution and the forecast horizon increase, the number of required marginal CDFs increases rapidly, as does the computational burden. For instance, a 48 h ahead forecast at a 5 min resolution requires 576 marginal CDFs, after which the covariance matrix needs to be estimated. To alleviate the computational burden and to simplify the model chain, we propose to generate trajectory forecasts using pattern matching in order to downscale hourly numerical weather prediction (NWP) ensemble forecasts to power forecasts at a 5 min resolution, similar to [START_REF] Yang | Operational solar forecasting for the real-time market[END_REF]. However, unlike [START_REF] Yang | Operational solar forecasting for the real-time market[END_REF]-who use Mueen's Algorithm for Similarity Search [START_REF] Mueen | The fastest similarity search algorithm for time series subsequences under euclidean distance[END_REF]-we employ k-dimensional tree (kd-tree, [START_REF] Bentley | Multidimensional binary search trees used for associative searching[END_REF]) to allow for multivariate input features as recommended in [START_REF] Yang | Post-processing in solar forecasting: Ten overarching thinking tools[END_REF]. This approach constitutes a seamless method for the generation of RES trajectories at intraday and day-ahead horizons, in contrast with the seamless RES forecasting method of [START_REF] Carriere | A Novel Approach for Seamless Probabilistic Photovoltaic Power Forecasting Covering Multiple Time Frames[END_REF] that generates probabilistic density forecasts without temporal correlation.

The contributions of this paper can be summed up as follows:

• We use a pattern matching algorithm to downscale NWP ensemble forecasts to high temporal resolution multivariate probabilistic PV power forecasts, which has not been attempted before as far as we are aware. • The proposed method significantly reduces the computational burden and greatly simplifies the forecast model chain.

• The use of kd-tree means that the proposed method can be applied to any RES power forecast model chain using any number of input features. The remainder this paper is organized as follows: Section II provides an overview of the state-of-the-art and describes the data, methods and scoring rules. Section III presents and discusses the results while Section IV concludes our study.

II. BACKGROUND AND METHODOLOGY

A. Problem formulation

In this study, we concern ourselves with temporal trajectory forecasts only and what follows is a description of how such trajectories are often generated. We follow [START_REF] Pinson | From probabilistic forecasts to statistical scenarios of short-term wind power production[END_REF] to estimate the covariance matrix that describes the dependencies across time. This method requires that the marginal predictive CDFs be calibrated, which is to say that the observations Y are statistically similar to samples from the CDFs F , i.e., Y ∼ F .

Mathematically, a series of forecasts with horizon k is said to be calibrated when P k is uniformly distributed, i.e., P k ∼ U [0, 1], where p k,i can be computed using the probability integral transform (PIT)

p k,i = F i+k|i (y i+k|i ), ∀i (1)
where i is the index of the training set that runs from 1, . . . , N . Likewise, t is the index of the testing set that runs from 1, . . . , T . Subsequently, it is possible to transform P k to a standard normal random variable Z k ∼ N (0, 1) using the inverse Gaussian CDF

z k,i = Φ -1 (p k,i ) . ∀i (2)
In case of k = 1, . . . , K forecast horizons, the random vector Z = (Z 1 , . . . , Z K )

⊤ follows a multivariate normal distribution, i.e., Z ∼ N (µ 0 , Σ) where µ 0 is a vector of zeros and Σ a covariance matrix that describes the dependencies between forecast horizons with ones on its diagonal [START_REF] Pinson | From probabilistic forecasts to statistical scenarios of short-term wind power production[END_REF]. An unbiased estimate of Σ is [START_REF] Pinson | From probabilistic forecasts to statistical scenarios of short-term wind power production[END_REF]:

Σ = 1 N -1 N i=1 Z i Z ⊤ i . (3)
The trajectory forecasts can then be generated by sampling S K-dimensional vectors using a multivariate normal random number generator with µ 0 and Σ. Subsequently, the sampled vectors can be transformed back to the uniform variable P k using the standard normal CDF Φ:

p s,k = Φ (z s,k ) . ∀s, k (4)
Once the marginal inverse predictive CDFs F -1 t+k|t have been generated on the testing set, the autocorrelated standard uniform samples p s,k can be used to generate trajectory forecasts f s,t+k|t :

f s,t+k|t = F -1 t+k|t (p s,k ) . ∀s, k (5)
The above-described method becomes computationally very demanding when the temporal resolution and the forecast horizon increase. For instance, a forecast horizon of 48 h at a temporal resolution of 5 min requires 576 marginal predictive CDFs, whereas a forecast horizon of 24 h at a temporal resolution of 1 min requires 1,440 marginal predictive CDFs. Clearly, an alternative method that performs similarly but significantly faster is preferred, which is what we describe in the next section.

B. Pattern matching

The pattern matching model (PMM) proposed here is conceptually straightforward and is based on the analog ensemble [START_REF] Alessandrini | An analog ensemble for short-term probabilistic solar power forecast[END_REF]. However, rather than searching for analogs for a single forecast horizon to approximate F t+k|t , we aim to search for S analog trajectories to approximate the multivariate predictive CDF F t . In addition, we intend to leverage the computational efficiency of kd-tree, which can improve the computational performance by three orders of magnitude compared to bruteforce search [START_REF] Yang | Ultra-fast analog ensemble using kd-tree[END_REF]. To that end, we modify the original similarity metric in [START_REF] Alessandrini | An analog ensemble for short-term probabilistic solar power forecast[END_REF] such that X t contains the query NWP forecast issued at testing time t organized as one vector. Similarly, A i is a single vector containing the analog NWP forecast issued at training time i. The similarity metric accepts any number of variables and is defined as [START_REF] Yang | Post-processing in solar forecasting: Ten overarching thinking tools[END_REF]

d(X t , A i ) = J j=1 w j x (j) t -x (j) i 2 , (6)
where for simplicity J is the total dimension and w j = 1 ∀ j.

The result is the approximate predictive distribution F t comprising detrended 5-min PV power measurements observed from time i+1 up to time i+K, such that F t ∈ R S×K and the corresponding observations are y t ∈ R K . In addition, we compute the clear-sky global horizontal irradiance (GHI) vector g t ∈ R K with which F t is multiplied to attain PV power measurements, see Section II-D. In order to value each feature similarly, A i is centered and scaled and these factors are used to center and scale X i . C. Benchmarks 1) Quantile regression forests (QRF): In random forests, the prediction is the weighted average of the observed response variables. In contrast, in QRF, output is the weighted distribution of the observed response variables [START_REF] Meinshausen | Quantile regression forests[END_REF]. During training, each tree is grown on a random sample of the training data, thus reducing the correlation between the trees. To further decrease this correlation, a random subset of the features is selected at each candidate split [START_REF] Meinshausen | Quantile regression forests[END_REF]. Note that QRF is combined with the Gaussian copula described above to generate state-of-the-art trajectory forecasts.

2) Multivariate probabilistic ensemble: The second benchmark is the naive multivariate probabilistic ensemble (MuPEn) of which the marginal predictive CDFs are identical to those of the complete-history persistence ensemble [START_REF] Van Der Meer | A benchmark for multivariate probabilistic solar irradiance forecasts[END_REF]. It can be constructed by gathering all N K-length vectors from the historical observations that start at the same time ("HH:MM") as the current forecast issue time t. The result is an N ×K matrix from which we randomly sample-without replacement-S trajectories such that the result is a multivariate predictive distribution F t ∈ R S×K .

D. Data

In this study we use two data sources, the first of which is the aggregated PV power measurements at a 1-min resolution from Rhodes, Greece, which we average to a 5-min resolution. The only available metadata is the installed capacity, which is 18,164 kW. The PV power measurements are subsequently detrended using clear-sky GHI generated by the McClear model [START_REF] Lefèvre | McClear: a new model estimating downwelling solar radiation at ground level in clear-sky conditions[END_REF] for a tilt of 25°and azimuth of 180°due south. These values constitute the optimal orientation on Rhodes and we assume that a large share of the PV systems are installed with this orientation [START_REF] Kambezidis | Estimation of the optimum energy received by solar energy flat-plate convertors in greece using typical meteorological years. part i: South-oriented tilt angles[END_REF]. However, we notice that the tails of the PV power measurements at sunrise and sunset do not coincide with the clear-sky GHI, which is likely caused by a number of systems orientated towards the east and west, respectively. The detrended PV power measurements are therefore set to 0 when the zenith angle is larger than 90°or when the clear-sky GHI is lower than 1 W/m 2 as these values would otherwise be significantly higher than can be reasonably assumed. The PV power is normalized with the installed capacity so that the errors can be presented as dimensionless numbers or as percentages.

The second data source is the hourly NWP forecasts issued by the European Centre for Medium-Range Weather Forecasts (ECMWF), issued daily at 00:00 UTC. Figure 1 presents the grid points in red at which we have collected NWP forecasts, i.e., 66 grid points in total to capture the large scale weather pattern. The variables used in this study are surface solar downward radiation (SSRD) converted to GHI and then clearsky index, total cloud cover (TCC), 10 meter u-and v-wind components (U10 and V10) and 2-meter temperature (T2M). In order to have a sense of the variability present in the forecasts, we use the ensemble prediction system featuring 50 ensemble members [START_REF] Leutbecher | Ensemble forecasting[END_REF]. Given the 48 h forecast horizon of interest, query X t has a length of 48×66×50×5 = 792, 000. To reduce the vector length, we summarize the grid points and ensemble members in two ways: (i) the mean and standard deviation of the ensemble members at each grid point, resulting in 31,680 features (referred to as case study MS); and (ii) the 0.01, 0.02, . . . , 0.99 quantiles over all ensemble members and grid points, resulting in 23,760 features (referred to as case study QS).

Given the relatively short length of the data set (1 year), we test on a monthly basis while the remaining months serve as training data. For instance, when forecasting January, we use February up to and including December as training data. In addition, it is important to note that forecasts are issued on a rolling basis with hourly intervals. This means there are approximately 700 forecasts per month, although this can be fewer in the case where the PMM is unable to find sufficient quality analogs. Finally, we train 48 QRF models instead of 576 and use each model to forecast the 12 horizons within the hour to ease the computational burden.

E. Forecast verification

Since consistency is one of the three aspects of a "good" forecast, it is encouraged in expectation by proper scoring rules [START_REF] Murphy | What Is a Good Forecast? An Essay on the Nature of Goodness in Weather Forecasting[END_REF]. Such scoring rules can be used to evaluate the quality of forecasts, which is Murphy's second aspect [START_REF] Murphy | What Is a Good Forecast? An Essay on the Nature of Goodness in Weather Forecasting[END_REF]. In this study, we use the continuous ranked probability score (CRPS) to evaluate the sharpness and calibration of the marginal predictive CDFs F t+k|t , which is defined as [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF] CRPS(F t+k|t , y t+k|t ) =

α 0 F t+k|t (x) -1{y t+k|t ≤ x} 2 dx, ( 7
)
where 1 is the indicator function and α is the maximum observed power. Subsequently, we average the CRPS over all forecast valid times where the zenith angle is smaller than 85°.

To evaluate the multivariate predictive CDFs F t , we use the energy score (ES) and variogram score (VS). ES is a generalization of CRPS and can be defined as [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF] ES

(F t , y t ) = E F ∥X t -y t ∥ - 1 2 E F ∥X t -X ′ t ∥, (8)
where X t and X ′ t are independent random vectors sampled from F t and ∥ ∥ represents the Euclidean norm. A disadvantage of ES is that it is unable to discriminate between poorly or correctly specified dependence structures [START_REF] Pinson | Evaluating the quality of scenarios of short-term wind power generation[END_REF]. Therefore, we include VS to compare the quality of the dependence structure specification, which is defined as [START_REF] Scheuerer | Variogram-based proper scoring rules for probabilistic forecasts of multivariate quantities[END_REF]:

VS p (F t , y t ) = K i,j=1 w ij (|y i -y j | p -E F |x i -x j | p ) 2 , ( 9
)
where x i and x j are components i and j of random vector X t , which is distributed according to F t . In this study, F t is approximated by S K-dimensional trajectories x (1) x (2) • • • x (S) ⊤ and E F |x i -x j | p can be approximated by [START_REF] Scheuerer | Variogram-based proper scoring rules for probabilistic forecasts of multivariate quantities[END_REF]:

E F |x i -x j | p ≈ 1 S S s=1 |x (s) i -x (s) j | p , i, j = 1, . . . , K, (10) where x
(s) i and x

(s) i are elements i and j of the s th trajectory forecast. Although weights w ij can be used to add or reduce importance between certain forecast horizons, we use identity weights and set p = 0.5 as recommended in [START_REF] Scheuerer | Variogram-based proper scoring rules for probabilistic forecasts of multivariate quantities[END_REF]. Finally, CRPS, ES and VS are averaged over the testing samples after which skill scores relative to MuPEn are computed.

III. RESULTS AND DISCUSSION

A. Experimental setup

Preliminary experiments reveal that the zenith angle along the entire forecast horizon provides valuable information to accurately determine sunrise and sunset timing. Similarly, U10, V10 and T2M do not contribute positively to the forecast accuracy and are therefore omitted for the PMM. In addition, we find that a larger number of trajectories S positively affects VS, likely because it is more challenging to approximate the prevailing dependencies with fewer trajectories. We therefore set S = 50 in our experiments. Furthermore, recall that the forecast horizon is 48 h and that the temporal resolution is 5 min, which means that K = 576. Finally, recall that the two ways of representing the spatio-temporal information from the NWP forecasts are abbreviated by MS and QS, resulting in four case studies (e.g., PMM+MS or QRF+QS).

B. Univariate forecasts

Figure 2 presents histograms of the PIT variables computed over the testing set using (1). In addition, Fig. 3 presents the CRPS in percent of the nominal capacity as a function of the forecast horizon. The mean and standard deviations are computed across the 12 testing months, as described above. Several interesting observations can be made from these figures. First, Fig. 2 shows that the PMM+MS and QRF+QS forecasts tend to be underdispersed. In contrast, forecasts generated by PMM+QS and QRF+MS are calibrated better in the main distribution although the lowest quantile deviates in case of the former and the most extreme quantiles deviate in case of the latter. Second, Fig. 3 shows that PMM+MS does not perform well and that is mainly caused by poor forecasts in December (not shown here). The poor performance during December is likely due to the inability of the PMM to find high quality analogs, although PMM+QS did not suffer from this problem. The poor performance in December could therefore be caused by the way the NWP ensemble information is summarized.

Third, CRPS of PMM+QS is constant over the forecast horizon, whereas CRPS of QRF+MS and QRF+QS increases sharply over the first hours and then stabilizes. In control problems, stability in the forecast error variance over the entire forecast horizon is preferable because it reduces the bullwhip effect [START_REF] Yang | Operational solar forecasting for the real-time market[END_REF]. For instance, underdispersive forecasts could result in an optimistic strategy by the control algorithm that could lead to the overutilization of available storage to correct for the forecast errors.

Finally, QRF+MS performs best in terms of CRPS. However, it is likely that more data would favor the PMM as it requires a long history to find quality analogs. In contrast, QRF has been proven to be robust against relatively small data sets. In terms of computation time, it is important to note that the PMM is approximately 98% faster than QRF. In hypothesis testing it is common to set the null hypothesis (H 0 ) such that there is no difference between the forecasts whereas the alternative hypothesis (H a ) is that there is a difference. The Diebold-Mariano (DM) test is commonly used to test H 0 . The test requires computation of the forecast error loss differential d t = ℓ(F 1,t , y t ) -ℓ(F 2,t , y t ). However, the DM test is designed specifically for a forecast horizon and consequently, we are uncertain it applies here. Furthermore, the paired t-test is not recommended when temporal dependence and contemporaneous correlation are present [START_REF] Gilleland | Testing the tests: What are the impacts of incorrect assumptions when applying confidence intervals or hypothesis tests to compare competing forecasts?[END_REF]. Instead, [START_REF] Gilleland | Testing the tests: What are the impacts of incorrect assumptions when applying confidence intervals or hypothesis tests to compare competing forecasts?[END_REF] recommend the Hering-Genton (HG) test and circular block bootstrapping. We choose to use block bootstrapping because the HG test, like the DM test, assumes that the series d is covariance stationary, whereas circular block bootstrapping is relevant for small testing sets [START_REF] Gilleland | Testing the tests: What are the impacts of incorrect assumptions when applying confidence intervals or hypothesis tests to compare competing forecasts?[END_REF]. The bootstrap is repeated 10,000 times with a block length of √ T and the results are presented in Table I as the mean plusminus the standard deviation and the confidence interval (CI, significance level α = 5%) in parentheses. As the table shows, 0 always lies within the CI and therefore we cannot reject H 0 .

C. Multivariate forecasts

Finally, we report the skill scores relative to MuPEn in Table II, computed as 1 -L model /L MuPEn and where L is the loss averaged over the testing set. Evidently, QRF+MS performs best on all scores. However, given that the marginal predictive CDFs of PMM+QS are calibrated slightly better, that no learning step but only proper data organization is required, and that computation time is reduced by approximately 98%, we argue that the PMM is at the very least a highly efficient and interpretable asset in a forecaster's toolbox.

IV. CONCLUSIONS

We proposed a forecast model based on a pattern matching model (PMM) and used it to generate multivariate probabilistic forecasts of aggregated photovoltaic (PV) power production, approximated by trajectories. We compared the PMM against the de facto method of generating trajectories, namely by separately forecasting the marginal predictive distributions and estimating the covariance matrix that describes the dependencies between forecast horizons. Quantile regression forests (QRF) were used to forecast the marginals and mostly outperformed the PMM, although block bootstrapping indicated that we could not reject the null hypothesis of an accuracy difference between the forecasts with α = 5%. We argue that because of its simplicity and efficiency, PMM can be a valuable asset in a forecaster's toolbox. Lastly, it is important to highlight that the limited amount of available data may have had a negative effect on the quality of the analogs found by the PMM, whereas QRF is known to be robust against small sample sizes. Furthermore, extensive feature selection could improve the PMM's performance further. These aspects will be investigated in future work.

Fig. 1 .

 1 Fig. 1. Map of the area of interest where the red points represent the grid points at which we have collected NWP forecasts.

Fig. 2 .

 2 Fig.2. Histograms of the marginal PIT variables combined over all forecast horizons and testing instances where the zenith angle is smaller than 85°.

Figure 4 presentsFig. 3 .Fig. 4 .

 434 Figure 4 presents ES and VS as monthly averages in addition to the total averaged scores. The figure clearly shows that PMM+MS performs poorly in December. Besides December, the figure shows that QRF+MS performs substantially better in February than the other models. Overall, QRF+MS outperforms the others but given the limited number of forecastverification pairs-recall that there are about 700 forecastverification pairs per month-it is worthwhile to test the significance of the difference between the forecasts.

TABLE I BLOCK

 I BOOTSTRAPPED LOSS DIFFERENTIAL PRESENTED AS µ ± σ (2.5% -97.5%).

				ES
		PMM+QS	QRF+QS	PMM+MS	QRF+MS
	PMM+QS	0±0 (0-0)	0.06±0.31 (-0.5-0.83)	-0.08±0.21 (-0.67-0.31)	0.13±0.32 (-0.4-0.82)
	QRF+QS	(-)	0±0 (0-0)	-0.11±0.29 (-0.78-0.46)	0.11±0.22 (-0.2-0.68)
	PMM+MS	(-)	(-)	0±0 (0-0)	0.15±0.27 (-0.31-0.74)
	QRF+MS	(-)	(-)	(-)	0±0 (0-0)
				VS
		PMM+QS	QRF+QS	PMM+MS	QRF+MS
	PMM+QS	0±0 (0-0) 397.02±1116.4 (-1557.52-3187.29)	-246.11±794.92 (-1996.46-1632.91)	722.43±1199.35 (-977.84-3611.67)
	QRF+QS	(-)	0±0 (0-0)	-533.09±1108.93 (-3073.65-1598.25) 534.45±1119.49 (-657.42-4011.29)
	PMM+MS	(-)	(-)	0±0 (0-0)	776.92±1091.41 (-770.52-3274.18)
	QRF+MS	(-)	(-)	(-)	0±0 (0-0)

TABLE II CRPS

 II , ES AND VS SKILL SCORES, RELATIVE TO MUPEN. NOTE THAT WE COMPUTE THE MEAN AND STANDARD DEVIATION (µ ± σ) OVER ALL FORECAST HORIZONS FOR CRPS.

	Model	CRPS	ES	VS
	PMM+MS 0.188 ± 0.250 0.237 0.384
	PMM+QS	0.241 ± 0.234 0.272 0.398
	QRF+QS	0.289 ± 0.221 0.308 0.452
	QRF+MS	0.348 ± 0.203 0.361 0.536

Note that the spatial dependency is disregarded in the present study because it focuses on aggregated power production.

ACKNOWLEDGMENT

The authors wish to thank ECMWF for the access to Numerical Weather Predictions and HEDNO for the provision of renewable production data.


Loreine Makki

Antoine Laspeyres

Corentin Darbas

Anne-Sophie Descamps

Marc-Anthony Mannah

Christophe Batard

Nicolas Ginot

Equivalent Circuit Model of a Pulse Planar Transformer and Endurance to Abrupt Dv/dt

Keywords: SiC MOSFET, driver circuits, pulse transformer, planar transformer, Altium Designer, Ansys Q3D, ElectroMagnetic Compatibility (EMC)

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

I. INTRODUCTION

I N high power applications (>100 kW), dedicated gate drivers are often used to operate power semiconductor devices. They must provide an optimal command for power semiconductors using an electrically insulated circuit. Gate drivers design can be divided into 4 functions [START_REF] Weckbrodt | Short pulse transmission for sic communicating gate driver under high dv/dt[END_REF]: the transmission of command orders through the galvanic isolation, the isolated power supply, the output stage and the protection functions such as short circuit protection.

A. Gate driver functions

The transmission of command orders through the galvanic isolation is often realized by digital isolators or pulse transformers. Fig. 1 illustrates a typical structure used in recent gate drivers such as the CMT-TIT8244 from Cissoid [START_REF]CMT-TIT8244 1700v half-bridge sic mosfet gate driver (optimized for 62mm power modules)[END_REF] or the CGD1700HB3P-HM3 from Wolfspeed [START_REF]CGD1700HB3P-HM3 dual channel differential isolated gate driver hm3 cpm3 sic half-bridge module companion tool[END_REF].

To drive power semiconductors, an open loop DC-DC converter is used for the isolated power supply. Transformer based converters using push-pull topologies are widely used. As the power consumption is low in gate driver applications (<5 W), these topologies' efficiency and good EMI while being simple to design make them attractive [START_REF] Kamath | Push-pull converter simplifies isolated power supply design in HEV/EV systems[END_REF].

This work was supported by the "Additive Manufactured high voltage high current diode less SiC smart Power Module" (AM-PM) RAPID project from the French ministry of defense, SAFRAN Tech and DEEP concept

The output stage circuit is placed as close to the power semiconductor as possible to drive it. Even if the mean power consumption is low, a gate current inrush of several amperes is needed to turn-on and turn-off power semiconductors. The output stage allows signal buffering in order to drive the gate voltage of the power transistor. To drive SiC MOSFETs, gate voltage applied by the output stage ranges from +20/-5 V for second generation SiC devices to +15/-4 V for third generation ones.

Fig. 1. Basic structure for gate drivers using pulse transformers [START_REF] Weckbrodt | Short pulse transmission for sic communicating gate driver under high dv/dt[END_REF] In recent gate drivers, various protection functions are implemented. This allows the gate driver to detect or prevent abnormal operations of the power component. In high power applications, three main features are often available. An overcurrent detection method detects short circuits and turns off the power component as fast as possible to avoid short circuit failure and the component destruction. Anti-overlap protection system prevents both power components to turn on at the same time to avoid short-circuit situations. Under voltage lockout (UVLO) systems are also used to detect low power supply in the gate driver which can lead to fail operation of the power component.

B. Isolation technology and EMC issues

Wide bandgap (WBG) semiconductors such as SiC MOS-FETs allow higher frequency switching power electronics for better efficiency [START_REF] Morya | Applications of wide bandgap (wbg) devices in ac electric drives: A technology status review[END_REF]. Nevertheless, as switching frequency increases, WBG poses new EMI issues such as switching oscillations [START_REF] Zhang | An overview of wide bandgap power semiconductor device packaging techniques for emi reduction[END_REF], [START_REF] Geramirad | Experimental emi study of a 3-phase 100kw 1200v dual active bridge converter using sic mosfets[END_REF], [START_REF] El-Sharkh | Investigation and reduction of near electric field emitted from a power helical inductor[END_REF]. Given the close proximity between gate drivers and power component, this EMI sources create common currents though the gate driver as illustrated in Fig. 2. To ensure common mode current does not trigger the output stage, close attention must be paid to reduce the parasitic capacitance of the transmission line through the galvanic isolation. CMTI defines the maximum tolerable rate of rise or fall of the common mode voltage applied between both isolated circuits which will not trigger the detection circuits. This ensures both transmitter side and receiver side of the isolated transmission line to function without error within the specified application. Fig. 2. Common mode current path in electrically isolated gate drivers in half bridge (HB) configuration [START_REF] Makki | Investigating the Shielding Effect of Pulse Transformer Operation in Isolated Gate Drivers for SiC MOSFETs[END_REF] In high power applications, opto-isolator technologies are rarely used given their parasitic capacitance and the difficulties for integration of optical fibers commands. Recent gate drivers mainly use two technologies: Digital isolators and pulse transformers.

Digital capacitive isolators allow unidirectional isolated data transmission. A single-ended input signal is split into two differential signal components to pass through the isolation barrier. Each signal component is differentiated to be compared to one another. These transient inputs are then converted into short pulses thanks to Schmitt triggers. As long as the positive input of a comparator has a higher potential than its negative input, the comparator output will present a logical high, thus converting an input transient into a short output pulse. The input signal is then reconstructed with a NOR-gate flip-flop (Fig. 3). Fig. 3. Block diagram of a digital capacitive isolator [START_REF] Kugelstadt | Designing with digital isolators[END_REF] In pulse transformer-based drivers, the transient signal is applied to the pulse transformer. The transient signal is then decoded on the other side of the pulse transformer. A serial capacitance is connected with the primary to prevent DC magnetic field in the pulse transformer which can lead to the core saturation. Bidirectional communication is possible by applying a longer pulse to the secondary winding of the pulse transformer. Error message is usually sent using this method in order to not increase the parasitic capacitance by adding a dedicated transmission line (Fig. 4). Fig. 4. Pulse transformer block diagram with error pulse circuit [START_REF] Weckbrodt | Short pulse transmission for sic communicating gate driver under high dv/dt[END_REF], [START_REF] Ginot | MOSFET et IGBT : circuits de commande, sécurisation et protection du composant à semi-conducteur[END_REF] In order to reduce the circuit size, new gate driver integrated circuits (ICs) using coreless transformer technology have been developed [START_REF] Varajao | Isolated gate driving solutions[END_REF]. This design allows smaller parasitic capacitance and a greater operating limit because of the core saturation. This study focuses on planar transformers technology in order to improve their CMTI while reducing the parasitic capacitance as low as possible by adding electrostatic shielding. The targeted application is the transmission line of high-power gate driver for a 3300 V, 500 A MOSFET SiC module. The operating voltage chosen is 2.5 kV. First, a design method for pulse planar transformers with Altium Designer is proposed. Then a simulation methodology of said pulse transformers using Ansys is exposed. Finally, experimental results are compared to Ansys simulation using 125 kV µs -1 dv/dt to simulate the switching of a power transistor connected to the gate driver.

II. PULSE TRANSFORMER DESIGN USING ALTIUM DESIGNER

Pulse transformer is one of the mainstream solutions used to drive power components in Half Bridge configuration. This transformer must provide galvanic isolation between the controller connected to the primary side and the power component connected to the secondary side. Moreover, pulse transformers are often used as a DC supply to drive the power components [START_REF] Varajao | Isolated gate driving solutions[END_REF]. Recent gate drivers (CMT-TIT8244, CGD1700HB3P-HM3) often use a dedicated DC power supply to power new features such as UVLO or short circuit detection. In such applications, pulse transformers do not need to supply DC energy reducing the risk of magnetic saturation. In aeronautical applications, planar design present two advantages : a reduced footprint and a good reliability due to the isolation properties of FR-4 epoxy [START_REF] Weckbrodt | Short pulse transmission for sic communicating gate driver under high dv/dt[END_REF].

Due to EMC issues and common mode current issues, the pulse transformer integrates electrostatic shielding between the primary and the secondary windings. Electrostatic shielding act as a Faraday cage to block the effect of an electric field. They can reduce the electrostatic field created by common mode currents while being transparent to pulses as they are ineffective against magnetic field. These screens are commonly linked to their respective ground in order to drain common mode currents and reduce EMC noises. Planar topology allows easier screen design and integration. The pulse transformer structure is detailed in Fig. 5.

A. Standard requirement

Pulse transformers must provide galvanic isolation between the primary side and the secondary side. As the secondary side is connected to the power component, its potential can rise from 0 V to the power component maximum potential (2.5 kV in this study). To ensure people safety who are connected to the primary side, the pulse transformer must respect electrical standard isolation. In planar applications, electrical standards can be divided into three measures : clearance, creepage distance and solid isolation. Fig. 6 highlights the three standard parameters.

Clearance is the shortest distance in the air between two conductors. This standard value ensures conductors will not arc through the air. In planar transformer applications, one must know ferrite cores are considered as conductors. Thus, the distance between two conductors must be clearance distance plus ferrite cores' width. Creepage distance is the shortest distance to another conductor along the surface of the insulating material of the PCB. This distance ensures conductors will not arc if conductive dust particles were to settle on the PCB. Finally, solid isolation gives the peak voltage solid isolators such as FR4-Epoxy must hold to ensure conductors will not arc. When dealing with solid isolation, planar transformers design usually considers the direct path between both conductors. Red arrow in Fig. 7 illustrates this path. As stated previously, ferrite cores offer a zero length path for electricity as their are electrical conductors. Thus, one must take into consideration the electrical path going through the ferrite core. One example of such a path is shown in yellow in Fig. 7. The pulse transformer design focuses on a 3300 V, 500 A SiC MOSFET module. The operating voltage chosen is 2.5 kV. The European standard EN 50178 [START_REF]Electronic equipment for use in power installations[END_REF] is used as a requirement for aeronautical purposes. The reinforced isolation was used. The European standard gives a 22.9 mm clearance as well as a 25 mm creepage distance. According to the standard, solid isolation must withstand a 18.4 kV impulse voltage. The MCL-E-679F (J) FR4 epoxy multi-layer material is used for solid isolation because of its high dielectric strength (≥ 40 kV mm -1 ). Therefore, the solid isolation between each conductor must be higher than 460 µm.

2) Magnetic core requirements: In pulse transformer applications, a rectangular shape pulse is applied across the transformer winding. If a voltage V is applied to a coil of N turns during ∆t time, a magnetic flux density ∆B will build up. To avoid saturation of the magnetic core, the core area A must be chosen to respect the basic pulse transformer (1) [START_REF] Lord | Pulse transformers[END_REF].

V • ∆t = N • A • ∆B (1)
In this paper, 25 ns long pulses are used under a 15 V voltage constant (see Fig. 10). To reduce the transformer size as much as possible, the smallest planar core from ferroxcube was used (E14/3.5/Core 5/R). Assuming a maximum magnetic flux density ∆B = 50 mT to reduce core loss, (1) gives a minimum of turns N = 0.7.

3) Coil turn requirements: To simply pulse transformers transient behavior analysis for operation with rectangular pulse, IEEE standardized the equivalent circuit of pulse transformers [START_REF]Ieee standard for pulse transformers[END_REF]. In this simplified model 2 resulting criteria are mainly discussed : the overshoot and the rising time of the rectangular pulse. In their study, Bortis and al [START_REF] Bortis | Design procedure for compact pulse transformers with rectangular pulse shape and fast rise times[END_REF] discussed about the transformer leakage inductance and capacitance impact regarding these criteria. To have the fastest rising time possible, leakage inductance must be kept as low as possible. Nevertheless, if the transformer leakage inductance is too low, high voltage overshoot are to be expected.

As, leakage inductance is proportional to N 2 [START_REF] Bortis | Design procedure for compact pulse transformers with rectangular pulse shape and fast rise times[END_REF], a balance must be achieved to keep a fast rising time while minimizing the resulting overshoot. Several designs with N ∈ [1;5] were simulated using Ansys. A number of turns N=3 was finally chosen.

4) Resulting design:

The pulse transformer is composed of three turns with an electrostatic screen 210 µm away from the primary side. The secondary side is a mirror of the primary side to get a winding ratio of 1 between primary and secondary circuit. As solid isolation is critical in planar transformer application, a safety factor of 2 was applied to the standard impulse voltage resulting. Thus, each side is separated from the other by a 1 mm FR-4 epoxy distance to respect standard. A minimum distance of 1mm was also applied between conductors and the ferrite according to Fig. 7. To design such transformer a glued double E core is used as it is not a typical PCB thickness of 1.6 mm.

III. PULSE TRANSFORMER SIMULATION METHODOLOGY

The employed simulation process in ANSYS included alternative software tools. The model being primarily implemented in Altium designer was imported into ANSYS SIWAVE using specific data exchange files. After being verified in terms of layer stack up, thickness, and material assignment, it was exported into ANSYS Q3D Extractor which can be dynamically linked with ANSYS Circuit Design for further analysis.

Being the leading extraction parasitic tool, ANSYS Q3D Extractor was employed to simulate and analyze the planar pulse transformer model. The electromagnetic field simulations conducted by this software are vital for the extraction of resistance, conductance, partial inductance and capacitance, known as RLGC parameters [START_REF] Kolobov | Research of planar transformer properties using ansys software[END_REF]. The central objective is to ensure compliant simulation and experimental results which will establish a reliable foundation for imminent experimentations. Q3D Extractor exploits the Finite Element Method (FEM) and method of moments to provoke an electromagnetic field solution. It relies on the simplification of Maxwell's equations, termed quasistatic approximation by anticipating that the size of the analyzed design is diminutive compared to the wavelength of maximum sought frequency where the coupling between magnetic and electric fields can be ignored [START_REF] Makki | Investigating the shielding effect of pulse transformer operation in isolated gate drivers for sic mosfets[END_REF]. As depicted in Fig. 9, the modelled design in Altium designer was imported into ANSYS SIWAVE to verify layer stack up and material assignment before being exported into ANSYS Q3D Extractor. The complete analysis performed in the latter allows the determination of RLGC parameters. The analyzed model can be dynamically linked with ANSYS Circuit design as it grants coordinated simulation of the subcircuit with the project original simulator. This will permit the determination of the transformer Scattering (S) parameters and effectively analyze the transformer equivalent circuit model which will be demonstrated in the sections thereafter. The model reveals the copper windings and shielding layers which are integrated in an FR4 Epoxy dielectric material, all secured by an EE 3F36 ferrite core.

Import

IV. TRANSFORMER SCATTERING PARAMETERS

A. Transformer Eligible Frequency Band

As the pulse transformer was exported into ANSYS, a spectral analysis is need to estimate the frequency band used in short pulse transmission. The method used to pass data through the galvanic isolation can be described as follow: a positive pulse on the primary winding means a rising edge of the PWM input and a negative pulse means a falling edge. With WBG semiconductors and the rising in the switching frequency, short pulses were chosen to maximize the gate driver switching frequency. Based on the pulse generator (UCC27517DBV) used, 25 ns long short pulse was targeted. Pulse parameters can be seen in Fig. 10. The signal of Fig. 10 was implemented on Matlab and a frequency bandwidth analysis was performed using the obw function. Fig. 11 shows the function result. A 50 kHz PWM was chosen with a duty cycle of 50 %. Sampling time is 0.7 ns. The signal processing analysis gave back a 99 % occupied bandwidth of 84.3 MHz ranging from 149.7 kHz to 85.75 MHz (Fig. 11).

B. S Parameter Simulation

Intending to study the electrical behavior of the planar transformer and demonstrate a reliable simulation criterion, S parameter simulation analysis was implemented and compared to experimental results. In order to achieve this evaluation, the pulse planar transformer was simulated in ANSYS Q3D Extractor by assigning a frequency sweep in accordance to the eligible frequency band ranging from 100 kHz till 85 MHz. Furthermore, the model was analyzed based on designating all conductors as nets which allows the computation of GC matrices, whereas AC/DC RL analysis requires the assignment of source and sink excitations which signify where the current enters and exits the windings respectively. After a successful analysis, the model was dynamically linked to ANSYS Circuit design as a four port sub-circuit demonstrated with the source and sink excitations of the primary and secondary winding as shown in Fig. 12. In order to reproduce an experimental examination, 50 Ω impedance interface ports were connected to the source pin of each winding, and distinctive ground ports to the sink pins.

ANSYS Q3D ANSYS Circuit

A frequency domain analysis is realized by setting up a Linear Network Analysis (LNA) over the previously specified frequency band. This final step will enable the computation of the frequency-dependent scattering parameters which will be compared to the experimental measurements revealed in section IV-C.

C. S Parameter Experimental Validation

To authenticate the simulation results, the pulse planar transformer was fabricated in accordance to the design modelled in Altium designer. S parameter measurements were conducted using a two-port, fifty-ohm Vector Network Analyzer (VNA). Fig. 13 correlates the agreeable S parameter results over a frequency band from 100 kHz to 100 MHz. The dashed lines represent the simulated results, whereas solid lines depict the experimental ones. S11 and S22 refer to the reflection coefficients of ports 1 and 2, respectively, yet S21 and S12 identify the transmission coefficients from ports 1 to 2 and contrarily, respectively [START_REF] Jargon | The relationship between switch-term-corrected scattering-parameters and wave-parameters measured with a two-port vector network analyzer[END_REF]. The symmetrical structure of the planar transformer explains the parity of S11 and S22, and S12 and S21. It is imperative to proclaim that this experimental verification was performed to reveal result conformity which grants reliable reliance on the simulation analysis implemented. Thus, fortifying the ability to propose an equivalent circuit model of the planar transformer relying on ANSYS Q3D Extractor RLGC parameters conferred thereafter.

V. ELECTROMAGNETIC COMPATIBILITY ADHERENCE

Medium-Voltage (MV) Silicon Carbide (SiC) Metal Oxide Semiconductor Field Effect Transistors (MOSFETs) acquire significant intrinsic properties, making it viable to intensify the switching frequency of converters from a few hundred kHz to the MHz region [START_REF] Rothmund | Highly compact isolated gate driver with ultrafast overcurrent protection for 10 kv sic mosfets[END_REF], [START_REF] Østergaard | Calculation of planar transformer capacitance based on the applied terminal voltages[END_REF]. Nonetheless, with every upgrade, researchers encounter challenges and restrictions to achieve enhanced power converter designs. Advantageous rapid switching frequency realization becomes detrimental when resulting in high dv/dt amounting to 100 kV µs -1 , generating electromagnetic disturbances in gate driver circuits due to the passage of Common Mode (CM) currents through parasitic capacitances, known as interwinding capacitances [START_REF] Saket | Low parasitics planar transformer for llc resonant battery chargers[END_REF]. Hence, being responsible for invoking semiconductor switching orders, critical attention is required when designing pulse planar transformers. This section will reveal the endurance of the pulse planar transformer to a steep dv/dt occurrence by performing an innovative simulation test of an equivalent circuit model and complying the results to experimental examination, whilst ensuring Electromagnetic Compatibility (EMC) standards.

A. Pulse Transformer Equivalent Circuit Model

A generic and symmetrical equivalent high-frequency circuit model of the pulse planar transformer is presented in Fig. 14. The electrical model comprises an ideal transformer with a turn's ratio m, magnetizing inductance L m , core resistance R c , leakage inductances L P and L S , winding resistances R P and R S and stray capacitances C w , C mc1 and C mc2 . This model mainly focuses on the values of the interwinding capacitances C mc1 and C mc2 as they are the common mode current passage route connecting the primary and secondary winding which will be quantifies in the section thereafter.

B. Pulse Transformer Parameter Extraction

In order to extract the parasitic parameters of the simulated pulse planar transformer and establish an equivalent circuit model, ANSYS Q3D Extractor tool was utilized [START_REF]ANSYS, Inc. ANSYS Electromagnetics Suite[END_REF]. As previously illustrated in Fig. 12, the model is set up by assigning material characteristics, net identifications and conductor excitation allocations. RLGC parameters were quantified over a sufficient frequency band extending from 100 kHz till 85 MHz. The reduce matrix feature in ANSYS Q3D was used to evaluate the resistance, inductance and capacitance values of the transformer model. The main objective is to investigate the susceptibility of the equivalent circuit model of the planar transformer to an application of high dv/dt switching transition. span which holds remarkable significance. The variation of R P and R S is common due to skin effect phenomenon, however, negligible over the tested frequency range. Several simulations were conducted to create a definitive model of the transformer independent of the applied frequency. The examination was implemented based on two approaches. The first approach was by averaging the value of each element over the frequency span. The second approach was testing the model at a frequency of 10 MHz and 50 MHz separately. Nevertheless, the simulation results of all the aforementioned propositions yielded to an agreeable outcome.

Hence, the equivalent circuit model illustrated in Fig. 14 at a peculiar frequency of 50 MHz is considered with m=1, R C =500 Ω and C w =21.37 pF. The aim is to perform a high dv/dt application on this model and compare the results to an actual experimented observation which will be subsequently discussed.

C. dv/dt Simulated Immunity Test Setup

ANSYS Circuit design software was used to simulate a dv/dt application on the modelled pulse transformer [START_REF]Ansys Circuit Online Help[END_REF]. To replicate a severe switching condition of power semiconductors, a steep dv/dt of 125 kV µs -1 was applied between the primary and secondary windings' disparate grounds as revealed in Fig. 15. Firstly, on the primary side short 25 ns pulses with a 15 V amplitude are connected to transmit turn on and turn off signals intended to switch a 1200V SiC MOSFET. At the secondary side, pulse selection is essential to segregate positive and negative pulses, and a filter stage necessary to curtail the signal noise for a 0-3.3 V logic at the output. After setting up and verifying the design, a transient analysis from 0 till 2000 ns was applied to evaluate the results.

D. dv/dt Experimental Immunity Test Setup

The pulse generator is composed of 2 NAND gates which receive respectively the input signal (IN) and its inverted The pulse planar transformer was tested at IETR lab using the circuit described in Fig. 15. A 50 kHz PWM signal is generated with a 0.5 duty cycle by the function generator. A Haefely generator is connected to both ground to perform a steep dv/dt immunity test between the primary side and the secondary side of the pulse transformer. A series of bursts are generated in sync with the PWM signal to simulate the power semi-conductor switching. A 125 kV µs -1 voltage slope is used which is the theoretical maximum switching speed of the power module multiplied by a safety factor of 4. The installation of the experimental setup of the planar pulse transformer CMTI test is revealed in Fig. 16.

E. dv/dt Immunity Test results

Static CMTI monitors the output stage of the pulse planar transformer when a CMT strike happens. Testing conditions cover the input logic condition and the CMT waveform. In high power gate driver, parasitic turn on is a major issue as it can lead to short circuit situations. In this application, the input is tied in logic low (i.e power component turned off). At 0.8 µs, a transient signal is sent in the input to simulate a switch on. Then, at 1 µs, a common mode transient voltage (CMTV) of 125 kV µs -1 is generated between the primary and the secondary side. Noise sustainability results are taken from the SET input of the NOR-gate Flip-Flop in Fig. 15 to check if the CMTV triggers the output stage. Fig. 17 In order to compare the shielding effect on the pulse transformer sustainability, another transformer was produced without shielding. Ansys simulation gives an increase of the C mc1 and C mc2 values by a factor 2. CMTI results show again a good correlation between the experimental results and the simulated ones, with near 1 V CMT parasitic voltage and 0.9 V, respectively in the Flip-Flop input. In the end, both transformers give minor voltage spikes under a 125 kV µs -1 applications with respect to the threshold voltage of the Schmitt-trigger logic gate circuitry. In the shielded pulse transformer case, a higher CMTI is achieved.

VI. CONCLUSION

A planar pulse transformer using electrostatic shielding for 3.3 kV applications with respect to European standard EN 50178 was proposed in this paper. The simulation analysis used alternative software tools, Altium Designer to design, and ANSYS Q3D with dynamic links to further analyze the transformer.

After performing a spectral analysis on the applied pulse signal, an eligible frequency band was identified to extract RLGC parasitic parameters in ANSYS Q3D Extractor that revealed coherent results. The study contributed to construct an equivalent circuit model of the transformer at a specific frequency and investigate the latter under abrupt dv/dt conditions which simulate the rapid switching of power semiconductors. Convenient CM noise results reveal EMC conservancy and simulation analysis reliability.

In this paper, 2 pulse transformers were designed : one with shielding and one without. In both cases, a CMTI test was performed using a 125 kV µs -1 dvdt. The shielded transformer show a better dvdt immunity than the non shielded one by a factor 2.

Fig. 5 .

 5 Fig. 5. Pulse transformer overview with Altium Designer

  cl e a ra n ce C re e p a g e d is ta n ce S o li d is o la ti o n

Fig. 6 .

 6 Fig. 6. Pulse transformer standard considerations

Fig. 7 .

 7 Fig. 7. Pulse transformer solid isolation considerations

Fig. 8 .

 8 Fig. 8. Pulse transformer design and overview with Altium Designer

  Fig 8 summarises the planar transformers characteristics.

Fig. 9 .

 9 Fig. 9. Simulation methodology of the analysed pulse transformer

Fig. 10 .

 10 Fig. 10. Pulse parameters used in this study

Fig. 11 .

 11 Fig. 11. Power spectral density of the pulse signal

Fig. 12 .

 12 Fig. 12. Simulation arrangement of S parameter evaluation

Fig. 13 .

 13 Fig. 13. Experimental vs. simulated S parameter results of a typically shielded transformer

Fig. 14 .

 14 Fig. 14. Common-mode current modeling of the pulse planar transformer

Fig. 15 .

 15 Fig. 15. Circuit to test pulse planar transformer's susceptibility to an application of 125 kV µs -1 dv/dt between both primary and secondary ground

Fig. 16 .

 16 Fig. 16. Experimental setup of pulse transformer sensitivity to a high dv/dt

  displays the experimental CMTV generated by the Haefely generator. To match the experimental setup, a similar CMTV was used in the simulated test.

Fig. 17 .

 17 Fig. 17. Experimented vs. Simulated steep dv/dt application Fig. 18 demonstrates convenient experimental and simulation results after applying the input signal. CMTI results show a good correlation between the experimental results and the simulated ones, with 0.55 V CMT parasitic voltage and 0.45 V, respectively in the Flip-Flop input. The values are minor with respect to the threshold voltage (VT+ = 1.88 V) of the Schmitttrigger logic gate circuitry [25].

Fig. 18 .

 18 Fig. 18. Pulse planar transformer sustainability to abrupt 125 kV µs -1 dv/dt with shielding

Fig. 19 .

 19 Fig. 19. Pulse planar transformer sustainability to abrupt 125 kV µs -1 dv/dt without shielding

  Table. I provides additional information of each solutions.

  TableIIdemonstrates the consistency of the parameters L P , L S , L m , C mc1 and C mc2 over the desired frequency

		Pulse generator	15V							3V3	Ch1	NOR-gate Flip-Flop
	IN	XOR	15V	Transformer	1kΩ	33Ω 33Ω 120pF	12Ω 12Ω	Set Reset	NOR	OUT
					1kΩ		120pF			NOR
					Pulse selection Filter			3V3
		COM 1	COM 1		COM 2			COM 2
				Passive dv/ dt controller	dv/dt choice (125 kV/µs) HaefelyP90.1Burst Tester (HighdV/dt Generator )				Oscilloscope on battery
		Functiongenerator	Peak Voltagechoice						Ch1
								TABLE II
					ANSYS Q3D EXTRACTOR RLC PARAMETERS
					Frequency		ANSYS Q3D Parameter Extraction
					(MHz)	C mc1	C mc2	R P	R S	L P	L S	Lm
						(pF)	(pF)	(Ω)	(Ω)	(nH)	(nH)	(µH)
					0.1	0.47	0.61	0.47	0.47	249	249	10.1
					1	0.47	0.61	0.52	0.52	246	246	10.1
					10	0.46	0.58	0.99	0.99	230	230	10.1
					20	0.45	0.58	1.3	1.3	226	226	10.1
					30	0.45	0.58	1.55	1.55	225	225	10.1
					40	0.45	0.57	1.77	1.77	224	224	10.1
					50	0.45	0.57	1.96	1.96	223	223	10.1
					60	0.45	0.57	2.12	2.12	223	223	10.1
					70	0.45	0.57	2.28	2.28	222	222	10.1
					80	0.45	0.57	2.42	2.42	222	222	10.1

ACCEPTED MANUSCRIPT -CLEAN COPY

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TPEL.2022.3160278, IEEE Transactions on Power Electronics IEEE TRANSACTIONS ON POWER ELECTRONICS