Denis Laxalde

Mathias Legrand

Christophe Pierre

Nonlinear Modal Analysis of Mechanical Systems with Frictionless Contact Interfaces

This paper investigates the nonlinear analysis, in the form of nonlinear modes, of mechanical systems undergoing unilateral frictionless contact. The nonlinear eigenproblem is introduced by a Rayleigh quotient minimization with inequality constraints formulated in the frequency domain. An augmented Lagrangian procedure is used for the calculation of the nonlinear contact forces. The efficiency of the proposed method for large-scale mechanical systems involving non-smooth nonlinear terms is shown. An industrial application consisting of a compressor blade in contact with a rigid casing is proposed. Sensitivity of the nonlinear modal parameters to contact is illustrated.

Introduction

The concept of nonlinear modes is now commonly accepted as a relevant approach for the analysis of nonlinear dynamical systems [START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF]. While the variety of analytical methods developed during the past decades have contributed to strengthen the theoretical background on nonlinear modes, new challenges appear when it comes to applying this tool to large scale systems [START_REF] Kerschen | Nonlinear normal modes, part I: A useful framework for the structural dynamicist[END_REF]. An encouraging perspective is however brought by numerical methods such as invariant manifold approach [START_REF] Pesheck | A new galerkin-based approach for accurate non-linear normal modes through invariant manifolds[END_REF], asymptotic-numerical methods [START_REF] Arquier | Two methods for the computation of nonlinear modes of vibrating systems at large amplitudes[END_REF], shooting and continuation [START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF] or Fourier methods [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF].

Dealing with non-smooth nonlinearities yields other challenges both on numerical analysis and phenomenological aspects. Concerning contact nonlinearities, difficulties come from their essential nonsmooth nature. As contact is defined by a multi-valued law which takes the form of inequality constraints, solution methods are often related to optimization [START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF][START_REF] Mijar | Review of formulations for elastostatic frictional contact problems[END_REF]. Dual formulations then provide physical sense to Lagrange multipliers. When it comes to numerical formulations, a popular and robust approach to deal with contact constraints is to introduce augmented Lagrangian functionals [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF][START_REF] Laursen | A continuum-based finite element formulation for the implicit solution of multibody, large deformation-frictional contact problems[END_REF][START_REF] Pietrzak | Large deformation frictional contact mechanics: continuum formulation and augmented lagrangian treatment[END_REF].

In contact dynamics, most integrated approaches are based on time integration and alternatives are unfortunately seldom explored [START_REF] Woo | Application of the harmonic balance method to ground moling machines operating in periodic regimes[END_REF][START_REF] Leib | Experimental and numerical analysis of two vibro-impacting beams[END_REF]. Thus, considering analysis of contact dynamics in the framework of nonlinear modes seems an interesting yet challenging approach. In this paper, we investigate the dynamics of linear elastic and autonomous flexible structures with contact interfaces. The target application concerns rotating components in turbomachinery which can be subject to intermittent contacts (or impacts) with surrounding stationary structures. Complex phenomena often occurs in such applications [START_REF] Azeez | Numerical and experimental analysis of a continuous overhung rotor undergoing vibro-impacts[END_REF][START_REF] Lesaffre | Stability analysis of rotating beams rubbing on an elastic circular structure[END_REF][START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF] and a modal approach appears interesting in order to describe their essential nature in particular.

The proposed approach is first described in a continuum mechanics formalism in which the problem of finding eigensolutions can be regarded as a mixed frequency-time boundary value problem. The solution is sought in the form of a Fourier series and the eigenvalue problem reduces to finding critical values of a generalized Rayleigh quotient. An inequality constrained optimization problem then arises which is solved using an augmented Lagrangian approach. Then, approximations for numerical treatments and algorithm are described.

A numerical application is proposed in which the first nonlinear mode of a turbomachinery compressor blade with contact occurring at its tip edge is explored. Effects of contact on modal parameters are investigated.

Theoretical formulation

A flexible structure in intermittent contact with a rigid foundation is investigated. Since attention is given to eigensolutions of the nonlinear problem [START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF], no external forces are applied such as the final system of equations is autonomous. Usual linear damping is neglected as well as friction.

Two main difficulties arise. First, the non-differentiability of the non-smooth contact conditions with respect to the displacement field, makes any linearization impossible. Second, as in every modal formulation, displacements and frequency of the motion are to be determined.

Unilateral frictionless contact formulation

A flexible component in contact with a rigid foundation is considered as illustrated in Fig. 1. Contact

u (x) Ω Γ d Γ c
g(u) = u(x)•n -g 0 (x) (1)
where g 0 (x) represents the initial positive gap, n, the outward normal to Γ c , and u = (u 1 , u 2 , u 3 ) the displacement field of a material point x = (x 1 , x 2 , x 3 ). Contact conditions [START_REF] Wriggers | Application of augmented lagrangian techniques for non-linear constitutive laws in contact interfaces[END_REF], referred to as Kuhn-Tucker optimality conditions are such that for all x ∈ Γ c :

τ N ≤ 0, g (u) ≤ 0, g (u) • τ N = 0 (2)
where τ N = σ • n stands for the contact pressure acting on Γ c .

Boundary value problem

The linearized strain tensor mapping displacements to deformations is defined as ǫ kl = 1 2 (∂u k /∂x l + ∂u l /∂x k ) and the constitutive linear Hooke law, as σ ij = λδ ij ǫ kk + 2µǫ ij . The strong formulation of the problem with Dirichlet and Signorini boundary conditions takes the following form:

ρü -div σ(u) = 0 on Ω × R + * (3a) u = 0 on Γ d × R + * (3b) g (u) ≤ 0, τ N ≤ 0 and g (u) • τ N = 0 on Γ c × R + * (3c)
Since the subsequent developments involve variational forms, we introduce the space of admissible displacement fields:

V = u = (u 1 , u 2 , u 3 ) , u i ∈H 1 (Ω) , u = 0 on Γ d (4)
equipped with a classical inner product (denoted ., . in the following), and the usual linear forms on V ×V :

m (u, v) = Ω ρuv dx (5a) k (u, v) = Ω σ (u) : ǫ (v) dx (5b)

Mixed frequency-time eigenvalue problem

We are interested in finding periodic and non-trivial solutions of problem (3). As periodic but not necessarily , it seems natural to expand the displacement field in the frequency-domain using Fourier series:

u(t) = n∈Z ûn e jnωt with ûn = 1 T T u(t)e -jnωt dt (6)
in which ω = 2π/T is the fundamental eigenfrequency of the motion to be determined.

The eigenvalue problem consists in finding {ω, û}, with û = {û n , n ∈ Z} such as,

-div σ(û n ) = (nω) 2 ρû n on Ω × Z (7a) ûn = 0 on Γ d × Z (7b) g (û) ≤ 0, τ N ≤ 0 and g (û) • τ N = 0 on Γ c × [0, T ] (7c)
in which the gap function is now expressed using Eq. ( 6) as:

g (û) = n∈Z ûn • n e jnωt -g 0 (8)
Note that while Eqs. (7a) and (7b) are formulated using frequency-domain fields, Eq. (7c) involves time-domain fields only. Thus Problem [START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF] can be regarded as a mixed time-frequency domains boundary value problem.

Rayleigh quotient and optimization

Based on the previous eigenvalue problem and using the linear forms defined in Eqs. ( 5), we introduce the following generalized Rayleigh quotient:

r (û) = n∈Z k (û n , ûn ) m (û n , ûn ) (9)
An eigenvector û is a critical point of this Rayleigh quotient and the associated eigenvalue is such that ω 2 = r (û). Here, due to the Signorini conditions, critical points of the Rayleigh quotient are solutions of a constrained minimization problem:

min ûn∈V r (û) s.t. g (û) ≤ 0 on Γ c ∀t ∈ [0, T ] (10)
Optimality conditions for this constrained minimization problem imply that there exists a Lagrange multiplier time-function λ(t) > 0 such that:

∇ ûn r (û) + T λ(t), ∇ ûn g (û) dt = 0 ∀n ∈ Z (11a) λ(t) > 0 and λ(t), g (û) = 0 ∀t ∈ [0, T ] (11b)
The constraint g (u) ≤0 has been transferred to the Lagrange multiplier λ(t) and the integral in Eq. (11a) can be regarded as a continuous equivalent to traditional weighted sum of individual constraints by the Lagrange multipliers. Expanding Eqs. [START_REF] Pietrzak | Large deformation frictional contact mechanics: continuum formulation and augmented lagrangian treatment[END_REF] leads to a variational formulation:

Find {û n , n ∈ Z} such as ∀v ∈ V 2 m (û n , ûn ) -1 k (û n , v) -n 2 ω 2 (û) m (û n , v) + T λ, v e jnωt dt = 0 ∀n ∈ Z (12a) λ > 0 and λ, g (û) = 0 ∀t ∈ [0, T ] (12b)

Augmented Lagrangian

The technique of augmented Lagrange multipliers converts an initially constrained minimization problem into a saddle point problem for which the addition of penalty terms improves convergence and benefits computationally due to subsequent convex properties. The respective Uzawa-like solution algorithm performs descent on the primal displacement variable and ascent on the dual Lagrange multipliers (coinciding with contact forces when convergence is reached).

For problem [START_REF] Laursen | A continuum-based finite element formulation for the implicit solution of multibody, large deformation-frictional contact problems[END_REF], an augmented Lagrangian can be defined as [START_REF] Powell | Algorithms for nonlinear constraints that use lagrangian functions[END_REF]:

L κ (û, θ) = r (û) + 1 2 T √ κ (g (û) + θ) 2 + dt (13)
in which κθ represents a Lagrange multiplier and κ, a diagonal matrix of positive penalty parameters.

During an iteration along the primal variable, the solution algorithm holds estimated values of dual parameters θ and κ and searches for û (θ, κ) that minimizes the augmented Lagrange functional (13) such as:

∇ ûn L κ (û, θ) = 0 ∀n ∈ Z (14)
or, equivalently:

Find {û n , n ∈ Z} such as ∀v ∈ V 2 m (û n , ûn ) -1 k (û n , v) -n 2 ω 2 (û) m (û n , v) + T κ (g (û) + θ) + , v e jnωt dt = 0 ( 15
)
Parameter θ is then adjusted as follows:

θ ← θ + max x,t (g (û) , -θ) (16)
while the penalty coefficient κ can be updated if condition [START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF] did not improve the constraints.

Mode normalization

Since the frequency of the periodic motion is an unknown of the eigenvalue problem, Eq. ( 15) has more unknowns than equations. This under-determination can be overcome by conducting a mode normalization through one (since the system is conservative) master coordinate which amplitude defines the modal amplitude and which governs all remaining dependent variables. This amplitude is denoted by q and modal parameters ûn (q) and ω(q) are now fully defined by its value.

Numerical approximations and algorithm

In the framework of finite-element method, we consider a discretized displacement field U together with consistent mass M and stiffness K structural matrices based on [START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF]. Contact constraints (1) can be re-written as:

g (U) = AU -g 0 ( 17
)
where rectangular matrix A restricts the displacement vector U to contact degrees-of-freedom. Furthermore, a uniform time discretization is introduced:

t = {t k = kT /m, k = 1, . . ., m} (18)
such as ūk = u (t k ) are discrete time values of u and the Fourier series of Eq. ( 6) is truncated up to N -th order. Hence, frequency-domain and time-domain variables are now related to each other in a discrete form: In practical applications, it is usually required to have N significantly smaller than m.

ûn = 1 T m k=1 ūk e -j 2πkn
The final eigenvalue problem is obtained from Eq. ( 15) in which space, time and frequency approximations are introduced and then reduces to 2N + 1 coupled sub-problems:

2M -1 K -(nω) 2 M Ûn + A T m k=1 κ g Û + θ k + e j 2πkn m = 0 (20)
The proposed algorithm 1 based on the one detailed in [START_REF] Arora | Multiplier methods for engineering optimization[END_REF] also features Powell's measure of convergence:

G = max dof max k (g (ū k ) , -θ k ) (21)
and an adaptive procedure for updating penalty parameters and Lagrange multipliers.

Algorithm 1 Frequency-time domain augmented Lagrangian

Set i = 0, i max , G (0) = ∞, θ = θ (1) , κ = κ (1) Choose α > 1, β > 1, ǫ > 0 while i < i max do Find {û n } n=-N,...,N solution of Eq. (20) Evaluate {g(ū k )} k=1,...,m Set L = {l : | max (g(ū k ), -θ k ) |≥G (i) /α} Set Ḡ = max dof |max k (g (ū k ) , -θ k )| if Ḡ ≥ G (i) then ∀l ∈ L, κ k,l ←βκ k,l and θ k,l ←θ k,l /β else i←i + 1 θ (i) k ←θ k , κ (i) k ←κ k and G (i) ← Ḡ if G (i) ≤ ǫ then break end if θ k ←θ (i) k + max (g (ū k ) , -θ k ) if G (i) ≥ G (i-1) /α then ∀l ∈ L, κ k,l ← βκ k,l end if end if end while

Analysis of contact in compressor blades

As an illustration of the proposed technique, we consider a turbomachinery compressor blade with a rigid contact interface at its tip edge. The finite element model of this blade is shown in Fig. 2; nodes in red define the contact interface. This model is reduced by means of a Craig-Bampton procedure composed of 72 constraints modes associated with displacements of 24 retained nodes belonging to contact interface and 40 component normal modes. Contact occurs in the radial direction only and a uniform initial gap is assumed. In the sequel, all displacements and velocities are normalized with respect to this initial gap and nonlinear frequencies are normalized with respect to the first linear frequency of the blade.

Figs. 3 display modal parameters of the first nonlinear mode (associated with first bending mode of the blade). In Fig. 3a, which represents the evolution of eigenfrequency as function of modal amplitude, two particular regions with distinct dynamical behaviors can be distinguished. The first refers to modal amplitude values such as the gap is always open and is characterized by a linear dynamics and constant frequency. The second, when the modal amplitude is such that the contact can occur at the interface, is highly nonlinear. The frequency-dependent amplitude phenomenon with stiffening effects is retrieved. Fig. 3b shows a manifold for this mode, that is the displacement of one nonlinear coordinate plotted versus the displacement and velocity of a linear one (the master coordinate here). The effect of contact, which limits the amplitude of the nonlinear coordinate, is noticeable. In order to investigate in more details the effects of contact on modal parameters, Figs. 4 display the evolution of radial displacements of all edge nodes as function of time (one period) and for three selected values of the modal amplitude (with respect to Fig. 3a). Note that contact occurs when the radial displacements reach a unit value due to the chosen normalization. In Fig. 4a, the modal amplitude has not reached the limit contact amplitude yet and no contact occurs. For the amplitude considered in Fig. 4b, contact occurs on the two last two nodes of the edge and during about 15% of time in one period. In the largest modal amplitude depicted Fig. 4c, contact occurs on the last 4 nodes and also on the first 5 nodes of the edge and these two contact region are out of phase. These observations bring further explanations to the stiffening effects of contact which was noticed on Fig. 3a. Indeed, stiffening is related to both the size of the actual contact interface (i.e. the number of nodes for which contact occurs) and the duration of contact.

Finally, Figs. 5 display the mode shapes on the finite element mesh for the three previous modal amplitude values. The changes in shape due to contact are substantial in particular during initial times of the period and can result in stress concentrations.

Conclusions

A methodology for modal analysis of elastic structures with contact interfaces was proposed. Based on the definition of nonlinear modes as periodic solutions of the autonomous dynamical system, the eigenvalue problem is defined in a continuum mechanics framework with Signorini boundary conditions. A generalized Rayleigh quotient functional is introduced, which minimization under constraints is formulated using Lagrange multipliers and solved using an augmented Lagrangian approach. Numerical approximations in space, time and frequency are then introduced for the sake of numerical treatments.

As an application, a turbomachinery compressor blade in contact with a rigid casing is considered. The first nonlinear mode is studied using the proposed approach and the influence of contact (which depends on the amplitude of the motion) is investigated on modal parameters such as frequency and shapes.

Based on these encouraging developments, future works will focus on further investigations on nonlinear dynamics to include stability and bifurcations analysis and also more complex systems which include several deformable bodies in respective contacts. This tool can be used to predict dangerous interactions zones between blades and surrounding casings. In the future, friction effects and as well as wear of contact interfaces can also be studied.

Figure 1 :

 1 Figure 1: Notations for contact formulation forces acting on interface Γ c are defined according to a continuous gap function separating the possibly contacting bodies: g(u) = u(x)•ng 0 (x)(1)

Figure 2 :Figure 3 :

 23 Figure 2: Blade model and details of the contact interface

Figure 4 :

 4 Figure 4: Radial displacements of edge nodes for three modal amplitudes: (a) 5.7, (b) 9.3 and (c) 31.3

Figure 5 :

 5 Figure 5: Mode shapes for three modal amplitudes: (a) 5.7, (b) 9.3 and (c) 31.3 and at t = 0 and t = T /2.

Acknowledgment

Thanks go to Snecma for its technical and financial support. This work takes place in the framework of the MAIA mechanical research and technology program sponsored by CNRS, ONERA and SAFRAN Group.


Espen Gaarder Haug
email: espenhaug@mac.com

Spavieri Gianfranco

Espen Gaarder

Gianfranco Spavieri
email: gspavieri@gmail.com

Mass-charge metric in curved spacetime

Keywords: 04.20.-q, 04.20.Jb, 04.70.-s Mass-charge, general relativity, strong gravitational fields, black holes, Reisser-Nordström metric

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

Spacetime curvature arises as an e↵ect due to mass, or equivalently energy. Because of the equivalence mass-energy, the energy produced by fields of the same form but di↵erent nature, is expected to provide the same e↵ect on the corresponding curvature regardless of the fields di↵erent origins.

In pre-relativistic classical mechanics and Maxwell's electromagnetism we have the relations,

F New m = G M r 2 = M r 2 F Coul q = k Q r 2 = Q r 2 , ( 1 )
where the first equation represents Newton's law of gravity with G the constant of gravity, m and M the interacting masses and F New the interacting force, while the second equation represents Coulomb's law with k the constant, q and Q the interacting electric charges and F Coul the interacting force. The last terms in (1) correspond to the choice of units where G = k = 1. If a point charge q is brought from x = 1 to point x in a region where there is an electric scalar potential (x), the work done on the charge (representing its potential energy) is W = q (x). For a system composed of point charges of total charge Q localized on the mass M , the electric potential energy and corresponding energy density are,

W Q = 1 8⇡ ˆ|r | 2 d 3 x = 1 8⇡ ˆ|E| 2 d 3 x w Q = 1 8⇡ |E| 2 , (2)
where for spherical symmetry and at the distance r from the localized mass M , the electric field is E r = Q/r 2 . Considering that analogous results are obtained for the gravitational potential energy, 1 the total energy density of the system M + Q can be conveniently expressed as,

w Q M = 1 8⇡ |E Q M | 2 = 1 8⇡ Q 2 M r 4 = 1 8⇡ M 2 r 4 + 1 8⇡ Q 2 r 4 , (3)
where

Q 2 M = M 2 + Q 2 Q M = p M 2 + Q 2 . ( 4
)
The contribution to the total energy densities in (3) that can never disappear is the gravitational, because present even for a perfectly neutral body. Since electromagnetic is much stronger than gravitational interaction, the term with M 2 in (3) and ( 4) is generally neglected in most physical situations where electric fields are present. However, in the context of black holes, astrophysical objects with a large mass M are generally almost neutral and in proximity of the mass M where r is small the energy density due to M may exceed the one due to Q.

The stress-energy tensor T µ⌫

With the total energy density given by ( 3) we can form a stress-energy tensor T µ⌫ that takes into account the field energy density due to M and Q. Stress-energy tensors that take into account di↵erent forms of energy are common in physics in several contexts outside or within general relativity. As an example, we may consider the tensor T µ⌫ = ✓ µ⌫ + S µ⌫ + 0 U µ U ⌫ , used in the literature [START_REF] Møller | The Theory of Relativity[END_REF][START_REF] Spavieri | A non-standard expression for the force and torque on a magnetic dipole[END_REF][START_REF] Aharonov | Comment on "proposed Aharonov-Casher e↵ect: Another example of an Aharonov-Bohm e↵ect arising from a classical lag[END_REF][START_REF] Spavieri | E↵ective interaction force between an electric charge and a magnetic dipole and locality (or nonlocality) in quantum e↵ects of the aharonov-bohm type[END_REF] in the context of electrodynamics of charged bodies with internal stresses. This tensor is complemented by the continuity equation, @ µ T µ⌫ = 0 while the term ✓ µ⌫ represents the usual electromagnetic tensor, S µ⌫ the stress tensor, and 0 the proper density of the proper mass.

In our case, we do not consider the stresses inside the mass M as we are interested in the solution of Einstein's field equations outside M . For our purpose we can simply exploit the fact that, for spherical symmetry, the form of the gravitational and electric fields is the same, as shown in [START_REF] Aharonov | Comment on "proposed Aharonov-Casher e↵ect: Another example of an Aharonov-Bohm e↵ect arising from a classical lag[END_REF]. Hence, we can formally adopt the stress-energy tensor T µ⌫ = ✓ µ⌫ (Q M ) where Q M is given by (4). Thus, our approach reflects the one by Reisser-Nordström [START_REF] Reissner | Über die eigengravitation des elektrischen feldes nach der Einsteinschen theorie[END_REF][START_REF] Nordström | On the energy of the gravitation field in Einstein's theory[END_REF] and we use here the formalism adopted by Zee [START_REF] Zee | Einstein gravity in a nutshell[END_REF] using geometric units c = G = k = 1.

The formal properties of T µ⌫ and the related field tensor F µ are unchanged, being,

T µ⌫ = 1 µ 0 (F µ F ⌫↵ g ↵ 1 4 g µ F F ).
The corresponding Einstein's [START_REF] Einstein | Näherungsweise integration der feldgleichungen der gravitation[END_REF] field equations are ,

R µ⌫ 1 2 g µ⌫ R + 8⇡T µ⌫ = 0.
In the context of black holes and for spherically symmetric spacetime the metric has the form,

ds 2 = A(t, r)dt 2 + B(t, r)dr 2 + r 2 (d✓ 2 + sin 2 (✓)d' 2 ) ( 5 )
In the solution of the Reisser-Nordström form, the relevant component of T µ⌫ is T 22 = (1/2µ 0 )r 2 F 01 F 01 . As we have only a radial component of the electric field, we may write:

E r = F 01 = F 10 = Q M /r 2 .
Then, T 22 = (1/2µ 0 )r 2 (Q 2 M /r 4 ) and we find that A(t, r) and B(t, r) obey the equations (see [START_REF] Zee | Einstein gravity in a nutshell[END_REF]),

A(r, t) = f (t) + C(t) r + Q 2 M r 2 AB = f (t).
So far, the solution for ( 5) is expressed in terms of mere algebraic quantities and does not represent physical reality. In order to be used to reflect physical reality and provide a falsifiable theory that can be tested experimentally, we need to relate the algebraic quantities to observable physical quantities. The quantities f and C can be determined by taking the limit to the Schwarzschild metric Q M ! 0 and,

g 00 = 1 2M r f (t) = 1 C(t) = 2M ⌘ r S
Finally, by means of ( 4) we find the horizon function = A(r),

= A(r) = 1 2M r + M 2 r 2 + Q 2 r 2 = 1 2M r (1 M 2r ) + Q 2 r 2
and the resulting metric is,

ds 2 = dt 2 + 1 dr 2 + r 2 d⌦ 2 , (6)
which, for almost neutral bodies where the charge Q ' 0 can be neglected, can be expressed as,

ds 2 = [1 2M r (1 M 2r )]dt 2 + [1 2M r (1 M 2r )] 1 dr 2 + r 2 d⌦ 2 (7)
Expressing the energy of gravitational fields as mass contributions to M . A possible interpretation of the additional mass term in the new metric ( 7) is obtained by comparing the Schwarzschild mass M S with the mass M (1 M 2r ) = M M 2 2r of our metric ( 6) or [START_REF] Zee | Einstein gravity in a nutshell[END_REF]. Integrating the gravitational energy density given by (3) starting from the observer position r up to r = 1 we find,

W M = 4⇡ ˆr=1 r ( 1 8⇡ M 2 r 04 )r 02 dr = M 2 2r , (8)
which is the mass, or equivalent energy, stored in the gravitational field M/r 2 outside r. The stored field energy in the space from r to r = 1 must come from the original mass corresponding to absence of field energy. Hence, the e↵ect of the field energy from r to r = 1 is to reduce the Schwarzschild mass M S (and corresponding space curvature) by M 2 2r . Consequently, the resulting metric [START_REF] Nordström | On the energy of the gravitation field in Einstein's theory[END_REF], or [START_REF] Zee | Einstein gravity in a nutshell[END_REF], can be obtained directly from the Schwarzschild (or Reisser-Nordström) metric by inserting into it the contribution to M S given by [START_REF] Einstein | Näherungsweise integration der feldgleichungen der gravitation[END_REF]. Our metric indicates that the energy of gravitational fields plays a role analogous to that of electromagnetic fields. Nevertheless, the additional mass term in our metric has to be considered in a wider context than that of simple analogy between Newton's gravity law and Coulomb's law. Of course, we have to go beyond Newton's gravity theory but we know today that the analogy between features of gravitational fields and those of electromagnetic fields is holding and is supported by developments in gravitational waves. A relevant confirmation is provided by the observed orbital decay of the Hulse-Taylor binary pulsar that matched the decay foreseen by general relativity as gravitation is lost to gravitational radiation, see [START_REF] Taylor | A new test of general relativity: Gravitational radiation and the binary pulsar psr 1913+16[END_REF]. Bondi et al. [START_REF] Bondi | Gravitational waves in general relativity[END_REF][START_REF] Bondi | Gravitational waves in general relativity, vii. waves from axi-symmetric isolated system[END_REF] in 1960 were likely the first to suggest that gravitational objects could lose mass due to the emission of gravitational waves, see also [START_REF] Penrose | On cosmological mass with positive ⇤[END_REF][START_REF] Hollands | Bondi mass cannot become negative in higher dimensions[END_REF]. Since the detection of gravitational waves from the merger of a binary black hole system by LIGO in 2016, the relationship between mass loss and gravitational waves has garnered increased interest, see for example [START_REF] Cao | Bound on the rate of bondi mass loss[END_REF][START_REF] Saw | Mass-loss of an isolated gravitating system due to energy carried away by gravitational waves with a cosmological constant[END_REF][START_REF] Saw | Mass loss due to gravitational waves with > 0[END_REF] for more information. In conclusion, the added mass term in the exact solution to Einstein's field equation given by our metric [START_REF] Zee | Einstein gravity in a nutshell[END_REF] is not surprising because based essentially on the equivalence between mass and field energy, regardless of whether the field source is the charge Q or the mass M. Therefore, our metric could potentially provide new insights in this area and, as shown below, in the interpretation of observational data in the context of strong fields.

Funding

No funding was received for this research.

Availability of data and materials

The study contains no new data.

.

Horizons of the metric

The metric gives the following horizons:

In the case of a neutral charge (Q = 0), the metric only yields one horizon, which is equal to half the Schwarzschild radius, i.e., r = M , instead of r = r s = 2M . For comparison, the horizons in the Reissner-Nordström metric are given by r ± = M ± p M 2 Q 2 , as shown in [START_REF] Hobson | General Relativity, An Introduction for Physicists[END_REF][START_REF] Bambi | General Relativity[END_REF][START_REF] Aretakis | Dynamics of Extremal Black Holes[END_REF] for example. The well-known extremal solution of the Reissner-Nordström metric (see [START_REF] Zee | Einstein gravity in a nutshell[END_REF]) occurs when Q = M . In this case, it predicts the same outcome as our metric when the charge is neutral.

Predictions of the metric in strong fields

The metric ( Eq. 7) predicts that the event horizon for black holes without charge is only half of what is predicted by the Schwarzschild metric. This would result in a stronger predicted gravitational redshift for light emitted just outside the horizon (from the accretion disk). The accretion disk will be able to approach closer to the center of the black hole compared to the Schwarzschild metric. This finding could potentially help explain the absence of observed velocity time dilation in high-z quasars (see [START_REF] Hawkins | Time dilation and quasar variability[END_REF][START_REF] Hawkins | On time dilation in quasar light curves[END_REF]), as more of the red-shift could be simply due to gravitational redshift. Although a detailed study would be necessary before drawing any definite conclusions, it seems that the new metric can be of help for explaining these observed phenomena.

Declarations Ethical Approval

Not applicable as a pure theoretical study.

Competing interests

The authors have no competing interests to report.

Authors' contributions

Both authors contributed to the research writing and review of the manuscript.


C Feutry

P Piantanida

P Duhamel

LEARNING SEMI-SUPERVISED ANONYMIZED REPRESENTATIONS BY MUTUAL INFORMATION

This paper addresses the problem of removing from a set of data (here images) a given private information, while still allowing other utilities on the processed data. This is obtained by training concurrently a GAN-like discriminator and an autoencoder. The optimization of the resulting structure involves a novel surrogate of the misclassification probability of the information to remove. Several examples are given, demonstrating that a good level of privacy can be obtained on images at the cost of the introduction of very small artifacts.

INTRODUCTION

The problem we address, i.e. removing private information from a set of images, while still allowing any other search in the set is a problem that is supervised with respect to the private label, but unsupervised with respect to the original image. Hence the name of semi-supervised anonymization.

Related work. Some previous works addressed similar problems, very often in a more constrain setup. Reference [START_REF] Chen | Vgan-based image representation learning for privacypreserving facial expression recognition[END_REF] proposes a tool that allows to perform some emotion preservation while transferring identities using a variational generative adversarial network. This work is fully supervised contrary to our work, that only use private labels during training. Several other works, such as [START_REF] Edwards | Censoring representations with an adversary[END_REF] or [START_REF] Raval | Protecting visual secrets using adversarial nets[END_REF] address the problem of hiding tags (text label, QR code,...) in images. This problem is somewhat different from ours, because their target is to hide some regions of the image while we focus on hiding inherent private data. Another very related approach is image to image translation, such as [START_REF] Zhu | Unpaired image-to-image translation using cycle-consistent adversarial networks[END_REF][START_REF] Liu | Unsupervised image-to-image translation networks[END_REF]. Instead of transferring the images to another domain, we provide a single "identity domain" representation where it is as hard as possible to track the initial identity (assuming that the private information is the identity).

De-identification is similar to anonymization but consists in the process of preventing someone personal ID from being revealed. The main difference between this concept and data anonymization is that some identifying information can be preserved in order to be relinked only by a trusted party or by the original data operator, whereas in the case of anonymization, no re-identification should be possible (e.g., [START_REF] Blaz Meden | Face deidentification with generative deep neural networks[END_REF][START_REF] Li | De-identification without losing faces[END_REF]).

Contributions. This paper presents a semi-supervised framework for removing a private information from image datasets, while keeping the distortion on the resulting images as small as possible. This model was used on two image datasets: a dataset of handwritten digits and a dataset of facial expression. It provided interesting results in the context of anonymization without a target task. Two cases should be distinguished: the case where the targeted private features are strongly connected to the other constitutive features of the image (denoted as intertwined) or not. Obviously, the later case is less difficult. Both cases are addressed below.

PRESENTATION OF THE PROBLEM

Since we wish to be able to perform the same classifications (except the private one) on the processed images and on the original images, the metric used to evaluate the output must be a mix between a fidelity metric (to ensure as much information as possible remains) and a private metric. However, the similarity metric cannot be a fixed metric (e.g., a traditional distortion). Indeed, a metric that performs well on a kind of data may not be as efficient on another data. Preliminary tests using fixed metrics yielded poor results. This directed us to use a GAN-like discriminator as the flexible metric that is trained concurrently with an autoencoder. This has the advantage of keeping a wide range of applications and yet to obtain accurate and neat results. This flexibility is also useful to train an anonymizing autoencoder where images will need to change significantly to remove private features. The resulting architecture is displayed on Fig. 1. Random variables are in capital letters; b • means computed over empirical distribution or estimate when over respectively a function or a random variable. Let T n be a training set of size n, where each element (x i , z i ) is composed of x i 2 X , a real vector and of z i 2 Z, the private label of this sample.

Our goal is to learn the parameters of the three networks: (i) q X|X the deep autoencoder; (ii) q Ŷ |V the discriminator; and (iii) q Ẑ| X the private classifiers. The selector S is before the discriminator and it allows a selection between x and x. The input v of the discriminator is the output of S which is driven by a variable

Y . If y = 1 then v = x otherwise v = x.
Learning with anonymization constraints. Consider the following constrained classification problem:

min (q X|X ,q Ŷ |V )2F
n P e q X|X , q Ŷ |V : min q Ẑ| X : X !P(Z)

P e q X|X , q Ẑ| X 1 " o , (1)
for a prescribed probability 1/|Z|  " < 1, where the minimization is over the set of restricted encoders and classifiers (q X|X , q Ŷ |V ) 2 F according to a model class F.

The above expression requires representations with (1 ")-approximate guarantees (over all possible classifiers) w.r.t. the misclassification probability of the private labels. Obviously, " can be replaced by a suitable positive multiplier ⌘ (") yielding a relaxed version of the objective.

min n P e q X|X , q Ŷ |V • P e q X|X , q ? Ẑ| X o , (2)
where q ? Ẑ| X is the minimizer of P e q X|X , q Ẑ| X . Expression (2) does not lead to a tractable objective for training (q X|X , q Ŷ |V ). Yet, it suggests a competitive game between two players: an adversary trying to infer the private labels Z from our representations x, by minimizing P e q X|X , q Ẑ| X over all possible q Ẑ| X , and a fidelity learner predicting the labels Y (i.e. the representations' fidelity), by optimizing a classifier q Ŷ |V over a prescribed model class F. We can trade-off these two quantities via the autoencoder model q X|X . This idea will be further developed in the next section.

Cross-entropy loss: Given two distributions q X|X : X ! P(X ) and q Ŷ |V : V ! P(Y), define the average (over representations) cross-entropy loss as:

` q X|X (•|x), q Ŷ |V (y|•) := ⌦ q X|X (•|x), log q Ŷ |V (y|•) ↵
(3) As usual, we shall measure the expected performance of (q X|X , q Ŷ |V ) via the risk:

L(q Ŷ |V , q X|X ) := E pXY ⇥ ` q X|X (•|X), q Ŷ |V (Y |•) ⇤ . (4)
It is not difficult to show from Fano's inequality that the lower bound on the misclassification task of the private labels is a monotonically decreasing function of the mutual information I(Z; X). This implies that any limitation of the mutual information between private labels Z and representations X will bound from below the probability of misclassification of private labels, whichever the chosen classifier q Ẑ| X . Besides, it is well-known that the cross-entropy provides a surrogate to optimize the misclassification probability of Ŷ , which motivates the cross-entropy loss. These information-theoretic bounds provide a mathematical objective in order to browse the trade-off (1) between all feasible misclassification probabilities P e q X|X , q Ŷ |V as a function of the prescribed (1 ") probability. Therefore, the learner's goal is to select an autoencoder q X|X and a classifier q Ŷ |V by minimizing jointly the risk and the mutual information. Yet, since p XY Z is unknown one cannot directly measure neither the risk nor the mutual information. It is usual to measure the agreement of a pair of candidates using a training dataset whose empirical distribution pXY Z is known. This yields an informationtheoretic objective, being a surrogate of expression of eq. ( 2):

min L emp (q Ŷ |V , q X|X ) + • b I(Z; X) , (5)
for a suitable multiplier 0, where L emp (q Ŷ |V , q X|X ) denotes the empirical risk of eq (3) (i.e. averaged w.r.t. pXY ). The mutual information must be empirically evaluated using qZ| X as being the posterior according to q X|X pXZ . Eq. ( 5) may be independently motivated by a different problem studying distortion-equivocation trade-offs [START_REF] Villard | Secure multiterminal source coding with side information at the eavesdropper[END_REF].

Reconstruction learning with anonymization. The initial experiments performed with a similar training objective as the one introduced by [START_REF] Ganin | Unsupervised domain adaptation by backpropagation[END_REF] led to an unstable training and a poor trade-off between the degree of anonymity and the quality of the representation. This guides us to use instead a new adversarial training objective given below.

An examination of expression [START_REF] Liu | Unsupervised image-to-image translation networks[END_REF] shows that it cannot be optimized since the posterior distribution qZ| X is not computable in high dimensions. We will looser this surrogate by upper bounding the empirical mutual information b

I(Z; X) = b H(Z) b H(Z| X).
The empirical entropy of Z can be upper bounded as follows:

b H(Z)  E pZ ⇥ log q Ẑ (Z) ⇤ (6)  E pZ E q X ⇥ log q Ẑ| X (Z| X) ⇤ (7)
⌘ E pZ E pX ⇥ ` q X|X (•|X), q Ẑ| X (Z|•) ⇤ (8)
:= L obj emp (q Ẑ| X , q X|X ), (9)
where ( 6) follows since the relative entropy is non-negative; [START_REF] Li | De-identification without losing faces[END_REF] follows by the convexity of t 7 ! log(t) and ( 8) follows from the definition of the cross-entropy loss. We will also resort to an approximation of the conditional entropy b H(Z| X) by learning an adequate empirical cross-entropy risk:

b H(Z| X) ⇡ E pXZ ⇥ ` q X|X (•|X), q Ẑ| X (Z|•) ⇤ ⌘ L emp (q Ẑ| X , q X|X ), (10)
which assumes a well-selected classifier q Ẑ| X , i.e., the resulting approximation error given by D qZ| X kq Ẑ| X |q X w.r.t. the exact q Ẑ| X is small enough. By combining expressions [START_REF] Ganin | Unsupervised domain adaptation by backpropagation[END_REF] and [START_REF] Alimoglu | Methods of combining multiple classifiers based on different representations for pen-based handwriting recognition[END_REF], and taking the absolute value, we obtain : b I(Z; X) -L obj emp (q Ẑ| X , q X|X ) L emp (q Ẑ| X , q X|X ) , [START_REF] Aneja | Modeling stylized character expressions via deep learning[END_REF] that leads to our tractable learning objective, which is an approximation of expression ( 5), being the surrogate of (2):

L (q Ŷ |V , q Ẑ| X , q X|X ) := L emp (q Ŷ |V , q X|X ) + • L obj emp (q Ẑ| X , q X|X ) L emp (q Ẑ| X , q X|X ) , (12)
for a suitable classifier q Ẑ| X and multiplier 0 which tunes the trade-off between the discriminator and the private task.

The data representations provided by the autoencoder q X|X must blur the private labels Z features from the raw data X while preserving original data's relevant features. Note that (9) corresponds to the loss of a 'random guessing' classifier in which the representations X are independent of private labels Z. Thus, training encoders q X|X to minimize (12) enforces the best classifier q Ẑ| X to get closer -in terms of loss -to the random guessing classifier.

Training objectives. The terms of loss function ( 12) are successively the binary cross-entropy of discriminator, the cross-entropy of the private branch and the objective crossentropy which expressions follow:

L emp (q Ŷ |V , q X|X ) = 1 n n X i=1 h y i log q Ŷ |V (•|v i ) +(1 y i ) log(1 q Ŷ |V (•|v i )) i , (13)
L emp (q Ẑ| X , q X|X ) =

1 n n X i=1 ⌦ e(z i ), log q Ẑ| X (•|x i ) ↵ , (14)
L obj emp (q Ẑ| X , q X|X ) = 1 n n X i=1 ⌦ pZ , log q Ẑ| X (•|x i ) ↵ . (15)
e(z i ) are "one-hot" vectors (z i component is 1 and the others 0) of the true labels of sample i = [1 : n]. pZ is the empirical distribution of the private labels. Each branch, i.e., q Ŷ |V and q Ẑ| X , is trained to minimize the associated cross-entropy loss (respectively (13) and ( 14)), whereas the autoencoder q X|X is trained to minimize eq. ( 12) (minimizing the cross-entropy loss with respect to Ŷ predictor while maximizing the new adversarial loss defined with respect to the b Z predictor). Training procedure. First the autoencoder is pre-trained using a pixel to pixel metric loss for a few epochs. Then the discriminator is pre-trained by alternating between: (i) updates on the discriminator parameters while autoencoder's parameters are frozen, with inputs going through the selector being a combination of autoencoded and original samples; (ii) updates on the autoencoder while discriminator's parameters are frozen. This allow the discriminator to learn the difference between both distributions. Finally, updates between the autoencoder on one side and the discriminator and the private branch on the other side are done alternatively. It is during this step that the privacy/fidelity trade-off occurs in the autoencoder.

SIMULATION RESULTS

The datasets used are doubly labelled: one is the private task classification and the other is used for the regular task classification in the performances assessment.

The Pen-digits dataset, from [START_REF] Alimoglu | Methods of combining multiple classifiers based on different representations for pen-based handwriting recognition[END_REF] provides the coordinates of digitally acquired pen movements of 44 persons (30 are in the training set and 14 in the test-set) writing digits from 0 to 9. We only used the training set (only 30 identities) which was randomly split into training, validation and test sets (size 5494, 1000 and 1000, respectively), sharing images of the same 30 persons. After being drawn, each image was down-sampled into a 20x20 grey-scale image.

The FERG (Facial Expression Research Group) dataset [START_REF] Aneja | Modeling stylized character expressions via deep learning[END_REF] contains 55767 annotated face synthetic images of six stylized characters modeled with the MAYA software. These images depict the seven following facial expressions: 'neutral', 'anger', 'fear', 'surprise', 'sadness', 'joy' and 'disgust'. Original 256x256 colour images have been pre-processed into re-sized 8-bit grey-scale 50x50 images.

Pen-digits results The private task labels in this dataset correspond to the writer identity. Hence, we want to transform handwritten digit images into other handwritten digit images while removing identity features and keeping relevant features about the original image digit.

Once the model is trained, original images are autoencoded into transformed images. Assessing the performance in a fair and automated manner is hard for it implies comparing and judging pictures. We processed images of the whole dataset and associated them with their original labels. Two as- sessments networks are then trained on the processed images, trying to recover both labels (ID and digit). The evaluation, plotted on Fig. 2, gives the accuracy on both tasks. This assessment method judges a semi-supervised transformation in a fully supervised manner, yet it is a good compromise to quantify the effect on the ID features and on some other features that we wanted to keep. This testing is realistic since we used all the available ID information to train a ID features extractor from the transformed images. It is the best case scenario for an attacker wanting to recover the ID from the processed dataset. Note that even if testing on only one non-private feature is not ideal, the fidelity assessment is performed on a feature not used during the training of the model.

FERG dataset. The private features to remove from the representations are the users identity labels.

Intertwined task method. The standard procedure explained before failed on the FERG dataset. Despite producing images that fool human eyes, they fail to fool a neural network. This is caused by the intertwined nature of the data. Identity traits and emotions are quite mixed. Our previous method led to a network that can produce visually anonymized results, but not a computer validated anonymization. To improve the anonymization on this kind of datasets, we derived a second method which specificity is restricting to one chosen identity the discriminator's train samples. This method although providing a small gain in anonymization was not sufficient. To perform a good anonymization we provided a protocol of substitution. Once all images are transformed into the chosen identity, each is compared to the chosen identity's real samples. Using pixel to pixel euclidean distance images are substitute by the closest chosen identity's training set image. The substitution is made really easy because the encoder shifts the representation of the original images very close to the selected identity while maintaining the underlying structure of the emotion features.

FERG Results. As stated above, two methods are associated with this particular framework. The first method restricts the discriminator to only one chosen ID and produces qualitative results shown on the Figs. 3 and4. Whereas the substitution method produces measurable results shown in Table 1.

Remark. One could argue that even if the dataset is anonymized, the chosen identity is still disclose because every sample wears the face of one preexisting identity. This is not a problem at all considering one can add an artificially identity in the databse to be the template of the anonymized (3,4), [START_REF] Raval | Protecting visual secrets using adversarial nets[END_REF][START_REF] Edwards | Censoring representations with an adversary[END_REF]. Other shift emotions [START_REF] Chen | Vgan-based image representation learning for privacypreserving facial expression recognition[END_REF][START_REF] Zhu | Unpaired image-to-image translation using cycle-consistent adversarial networks[END_REF]. faces, while still using our framework. Indeed the fact that the intertwined method allows to choose the final identity from the training set make it possible.

CONCLUSION

We succeed in providing an anonymization framework where only the private labels are used to train the model. This is possible via the introduction of a new loss (expression (12)), an adapted architecture and training protocol. To train this loss, we took advantage of and reused widely known and efficient methods such as GAN discriminator and autoencoder, and we managed to perform a smart combination of both to design a working architecture that was suited to our loss.

The results of this paper are made even more impressive because of the visual representation. The network we designed manages to affect the hidden representation of the data in the autoencoder to shift private bit of information toward a dataset universal template with minimum cost to the other parts of the information. The generalized information that we can deduce is that this method allows a structure adapted spatial shift in the learnt features space: shifting from identity to another without major loss of sub-cluster structure (here emotions or digits). The usage of the learnt discriminator allows a neat output image, pixel perfect in most cases.

Fig. 1 :

 1 Fig. 1: Architecture (S is the selector block).

Fig. 2 :

 2 Fig. 2: Accuracy as a function of 2 [0, 4.0] on Pen-digits dataset. Horizontal black dashed line is the random guessing classifier over the user-ID (3.33%).

Fig. 3 :

 3 Fig. 3: Evolution of the reconstructed images for a training and a validation samples (respectively top and bottom line). The resulting identity is fixed by the operator. The first image is the end of the pre-training phase: the autoencoder is trained to reproduce its entry (due to the pre-training's performance, this image is identical to the input sample). Going right follows the increase of the training epochs one by one. The most impressive changes occur over these first epochs. The following epochs (not show here) only refine the details of the encoded image to make it look genuine, resulting in the far right images that are the end training results.

Fig. 4 :

 4 Fig. 4: 24 couples of images composed of the sample and corresponding model output. All are randomly taken from the test set. Some representations are showing artifacts (missing an eye), e.g., representation at positions (line,column): (3,4),(3,2). Other shift emotions (1,4).

Table 1 :

 1 Accuracies for anonymization (supervised during training) with different targeted IDs and with the accuracies of emotions recognition (unsupervised during training) to measure conservation of unlabeled intertwined features. These are obtained by training a classifier on the processed dataset.

	targeted identity	0	2	3	4	5
	identity accuracy (%)	42.47 38.25 36.21 32.77 31.62
	emotions accuracy (%)	72.45 77.46 72.32 66.79 72.95

Acknowledgement

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 792464.


Mohammadhadi Shateri
email: mohammadhadi.shateri@mail.mcgill.ca

Francisco Messina
email: francisco.messina@mail.mcgill.ca

Pablo Piantanida

Fabrice Labeau

Deep Directed Information-Based Learning for Privacy-Preserving Smart Meter Data Release

Keywords: Privacy-preserving mechanism, Deep learning, Adversarial training, Recurrent Neural Networks, Directed information, Smart meters data

The explosion of data collection has raised serious privacy concerns in users due to the possibility that sharing data may also reveal sensitive information. The main goal of a privacypreserving mechanism is to prevent a malicious third party from inferring sensitive information while keeping the shared data useful. In this paper, we study this problem in the context of time series data and smart meters (SMs) power consumption measurements in particular. Although Mutual Information (MI) between private and released variables has been used as a common information-theoretic privacy measure, it fails to capture the causal time dependencies present in the power consumption time series data. To overcome this limitation, we introduce the Directed Information (DI) as a more meaningful measure of privacy in the considered setting and propose a novel loss function. The optimization is then performed using an adversarial framework where two Recurrent Neural Networks (RNNs), referred to as the releaser and the adversary, are trained with opposite goals. Our empirical studies on real-world data sets from SMs measurements in the worst-case scenario where an attacker has access to all the training data set used by the releaser, validate the proposed method and show the existing trade-offs between privacy and utility.

I. INTRODUCTION

In recent decades, there has been an explosive progress in data collection tools that can measure, analyze and disseminate users data. These advances have led to a large impact in several different fields such as electrical power systems, health care, digital banking, etc. However, users are generally unwilling to share their data due to the possibility that a third party could infer their personal private information, e.g., living habits, economic status, health history. Therefore, guaranteeing users privacy while preserving the benefits of data collection is an important challenge in modern data science [START_REF] Nelson | Security and privacy for big data: A systematic literature review[END_REF]. In this paper, we focus on this problem when the data has a time series structure and, in particular, we consider different privacy scenarios motivated by the deployment of Smart Meters (SMs) in electrical distribution networks [START_REF] Giaconi | Privacy-aware smart metering: Progress and challenges[END_REF]. The SMs are devices that can register a fine-grained electricity consumption of users and communicate this information to the utility provider in almost real-time. The utility of the SMs data is diverse [START_REF] Wang | Review of smart meter data analytics: Applications, methodologies, and challenges[END_REF], [START_REF] Depuru | Smart meters for power grid -challenges, issues, advantages and status[END_REF]. It can be used for power quality monitoring, timely fault detection, demand response, energy theft prevention, etc. However, widespread usage of SMs data can lead to serious leakages of consumers' private information, e.g., a malicious third party could use the data to detect the presence of residents at home as well as their personal habits [START_REF] Giaconi | Privacy-aware smart metering: Progress and challenges[END_REF]. This problem can have a serious impact on the deployment pace of SMs and, more broadly, in the development of smart electrical grids. Thus, it is critical to ensure that SMs data are sanitized before being released.

A very simple strategy that has been proposed in the context of SMs is to use pseudonyms rather than the real identities of users for data publishing purposes [START_REF] Efthymiou | Smart grid privacy via anonymization of smart metering data[END_REF]. However, this approach implicitly assumes that a trusted anonymizer is available. Another simple technique suggested in the literature is downsampling of the data, where the sampling rate is reduced to a level that does not pose any privacy threat [START_REF] Cárdenas | Privacy-aware sampling for residential demand response programs[END_REF], [START_REF] Mashima | Authenticated down-sampling for privacy-preserving energy usage data sharing[END_REF]. Although this approach may be effective from the privacy point of view, it could also limit seriously the utility of the SMs data for some applications requiring a timely response. More sophisticated and recent approaches exploit the presence of renewable energy sources and rechargeable batteries in homes to modify the actual energy consumption of users in order to hide the sensitive information [START_REF] Giaconi | Privacy-aware smart metering: Progress and challenges[END_REF], [START_REF] Backes | Differentially private smart metering with battery recharging[END_REF]- [START_REF] Erdemir | Privacy-cost trade-off in a smart meter system with a renewable energy source and a rechargeable battery[END_REF]. Some of these works use ideas from the well-known principle of differential privacy. However, recent articles suggest that the utility loss of differential privacy may be significant in practice [START_REF] Mendes | Privacy-preserving data mining: methods, metrics, and applications[END_REF]. It should be noted that our approach to the problem, which works only with the power measurement data, does not preclude the use of methods that change the energy consumption patterns by using physical resources. In fact, it should be viewed as a complementary approach that could even be used on top of the above mentioned methods.

In an information-theoretic context, privacy is generally measured by the Mutual Information (MI) between the sensitive and release variables [START_REF] Giaconi | Privacy-aware smart metering: Progress and challenges[END_REF], [START_REF] Li | Information-theoretic privacy for smart metering systems with a rechargeable battery[END_REF], [START_REF] Erdemir | Privacy-cost trade-off in a smart meter system with a renewable energy source and a rechargeable battery[END_REF], [START_REF] Sankar | Smart meter privacy: A theoretical framework[END_REF]. Some of these studies aim to find a privacy-utility trade-off using ideas from rate-distortion theory [START_REF] Sankar | Smart meter privacy: A theoretical framework[END_REF], [START_REF] Tripathy | Privacy-preserving adversarial networks[END_REF]. More specifically, the theoretical framework of the privacy-utility problem was proposed in [START_REF] Sankar | Smart meter privacy: A theoretical framework[END_REF], where a hidden Markov model for the power measurements of SMs is considered in which the distribution is assumed to be controlled by only the state of the home appliances. The privacy-utility trade-off is then found for a stationary Gaussian model of electricity load with MI between 1 arXiv:2011.11421v3 [cs.LG] 24 Nov 2021 release and private sequence of variables as a privacy measure. Besides the limitation of the Gaussian model, it is noted that the MI is not well-suited to capture the causal structure of the time series data.

In this paper, an information-theoretic cost function for privacy-preserving data release of time series is proposed. In order to take into account the time series structure and causality of the data in the privacy measure, we use the Directed Information (DI) [START_REF] Massey | Causality, feedback and directed information[END_REF] between the sensitive time series and an estimation of it. Then, a cost function is derived for the releaser mechanism based on an upper bound of the DI. To optimize and validate our cost function without imposing constraints on the data distribution, two recurrent neural networks, named as releaser and adversary networks, are employed. This approach is based on the framework of Generative Adversarial Networks (GANs) [START_REF] Goodfellow | Generative adversarial nets[END_REF], [START_REF] Huang | Contextaware generative adversarial privacy[END_REF], where two neural networks are trained simultaneously with opposite goals. We will show that by controlling the relative weight between a distortion measure and the DI privacy measure, we can control the utility-privacy trade-off of SMs power measurements. A similar approach for the privacy problem, but for different applications, was considered in [START_REF] Tripathy | Privacy-preserving adversarial networks[END_REF], [START_REF] Feutry | Learning Anonymized Representations with Adversarial Neural Networks[END_REF] considering the standard MI and independent and identically distributed (i.i.d.) data, where the authors use two deep feed-forward neural networks for the releaser and adversary. However, to the best of our knowledge, this is the first work to consider a DI privacy measure for time series data in the general privacy-preserving context and in SMs applications in particular.

This paper is organized as follows. In Section II, we present the theoretical formulation of the problem. Then, in Section III, a privacy-preserving data release method based on Long-Short Term Memory (LSTM) Recurrent Neural Networks (RNNs) is introduced along with the training algorithm. Results for two different applications based on SMs data are presented in Section IV. Finally, some concluding remarks and a discussion about future work are given in Section V.

Notation and conventions

A sequence of random variables (X 1 , . . . , X T ) of length T is denoted by X T , while x T = (x 1 , x 2 , . . . , x T ) denotes a realization of X T and x (i)T = (x

(i) 1 , x (i) 2 , . . . , x (i)
T ) denotes the i th sample in a minibatch used for training. Mutual information [START_REF] Cover | Elements of information theory[END_REF] between variables X and Y is represented as I(X; Y ) and the entropy as H(X). We use X ⌦ Y ⌦ Z to indicate that X, Y and Z form a Markov chain. The expectation of a random variable X is denoted as E[X].

II. PROBLEM FORMULATION AND TRAINING OBJECTIVE

A. Main definitions

Consider the private variables X T (such as occupancy label, household identity, or acorn family type), useful variables Y T (such as actual electricity consumption of household), and observed variables W T (which could be a combination of private and useful variables). We assume that X t takes values on a discrete alphabet X , for t 2 {1, . . . , T }. A releaser R ✓ (this notation is used to denote that the releaser is controlled by its parameters ✓) produces the release variables as Z t based on the observation W t , for each time t 2 {1, . . . , T }, while an adversary A attempts to infer X t based on Z t by finding an approximation of p X T |Z T which we shall denote by p XT |Z T . Thus, the Markov chain (X t , Y t ) ⌦ W t ⌦ Z t ⌦ Xt holds for all t 2 {1, . . . , T }. In addition, due to causality, the distribution p Z T XT |W T can be decomposed as follows:

p Z T XT |W T (z T , xT |w T ) = T Y t=1 p Zt|W t (z t |w t )p Xt|Z t (x t |z t ).
(1) The goal of the releaser R ✓ is to minimize the flow of information from the sensitive variables X T to their estimation XT while simultaneously keeping the distortion between the release variables Z T and the useful variables Y T below some given value. On the other hand, the goal of the adversary A (this notation is used to denote that the adversary is controlled by its parameters ) is to estimate X T as accurately as possible.

To take into account the causal relation between X T and XT , the flow of information is quantified by the DI [START_REF] Massey | Causality, feedback and directed information[END_REF]:

I(X T ! XT ) = T X t=1 I(X t ; Xt | Xt 1 ), (2)
where I(X t ; Xt | Xt 1 ) is the conditional mutual information between X t and Xt conditioned on Xt 1 [START_REF] Cover | Elements of information theory[END_REF]. The expected distortion between Z T and Y T is defined as:

D(Z T , Y T ) , E[d(Z T , Y T )], (3)
where d : R T ⇥ R T ! R is any distortion function (i.e., a metric on R T ). In order to ensure the quality of the release we shall impose the following constraint: D(Z T , Y T )  " for some given " 0. In this work, we will consider the normalized squared error as in [START_REF] Sankar | Smart meter privacy: A theoretical framework[END_REF], i.e.,

d(z T , y T ) , 1 T T X t=1 (z t y t ) 2 . (4)
Nevertheless, it should be noted that other distortion measures can also be relevant for the SMs data. For instance, demand response programs usually require an accurate knowledge of peak power consumption, so a distortion function closer to the infinity norm would be more meaningful for this particular application. This brief discussion simply illustrates that the distortion function should be properly matched to the intended application of the release variables Z T in order to preserve the characteristics of the useful variables Y T that are considered essential. Since the goal of this paper is mainly to introduce a new privacy measure and privacy-preserving data release framework, we will not further investigate different fidelity measures.

Therefore, the problem of finding an optimal releaser subject to the aforementioned adversary and distortion constraint can be formally written as follows:

min ✓ I ⇣ X T ! XT ⌘ , s.t. D(Z T , Y T )  ". (5)
Note that the solution of this optimization problem is a function of p XT |Z T , the conditional distributions that represent the adversary A .

B. Novel training objective

The optimization problem ( 5) can be directly used to define a loss function for R ✓ . However, note that the cost of computing the DI term is O(|X | T ), where |X | is the size of X . Thus, for the sake of tractability, DI will be replaced with the following surrogate bound:

I ⇣ X T ! XT ⌘ = T X t=1 H( Xt | Xt 1 ) H( Xt | Xt 1 , X t ) (i)  T X t=1 H( Xt | Xt 1 ) H( Xt | Xt 1 , X t , Z t ) (ii) = T X t=1 H( Xt | Xt 1 ) H( Xt |Z t ) (iii)  T log|X | T X t=1 H( Xt |Z t ), (6)
where (i) is due to the fact that conditioning reduces entropy; equality (ii) is due to the Markov chains X t ⌦ Z t ⌦ Xt and b X t 1 ⌦ Z t ⌦ b X t ; and (iii) is due to the trivial bound H( Xt | Xt 1 )  H( Xt )  log|X |. Therefore, the loss function for R ✓ can be written as

L R (✓, , ) = D(Z T , Y T ) T T X t=1 H( Xt |Z t ), (7)
where 0 controls the privacy-utility trade-off and the factor 1/T has been introduced for normalization purposes. It should be noted that the value of in [START_REF] Mashima | Authenticated down-sampling for privacy-preserving energy usage data sharing[END_REF] indirectly controls the achievable " in [START_REF] Efthymiou | Smart grid privacy via anonymization of smart metering data[END_REF], which means that we can control the privacy-utility trade-off by varying . For = 0, the loss function L R (✓, , ) reduces to the expected distortion, being independent from the adversary A . In such scenario, R ✓ offers no privacy guarantees. Conversely, for very large values of , the loss function L R (✓, , ) is dominated by the upper bound on the DI, so that privacy is the only goal of R ✓ . In this regime, we expect the adversary A to completely fail in the task of estimating X T , i.e., to approach to random guessing performance.

On the other hand, the adversary A is a classifier which optimizes the following cross-entropy loss: where the expectation should be understood w.r.t.

L A ( ) = 1 T T X t=1 E h log p Xt|Z t (X t |Z t ) i , (8)
p XtZ t . Notice that 1 T T X t=1 H(X t |Z t )  L A ( ). (9)
Therefore, if the adversary is ideal (i.e., p Xt|Z t = p Xt|Z t for all t), the releaser network, by maximizing

1 T P T t=1 H( Xt |Z t )
, prevents the adversary to infer private data.

III. PRIVACY-PRESERVING MECHANISM

Based on the previous theoretical formulation, an adversarial modeling framework consisting of two RNNs, a releaser R ✓ and an adversary A , is considered (see Fig. 1). Note that independent noise U T is appended to W T in order to randomize the released variables Z T , which is a popular approach in privacy-preserving methods. In addition, the available theoretical results show that, for Gaussian distributions, the optimal release contains such a noise component [START_REF] Sankar | Smart meter privacy: A theoretical framework[END_REF], [START_REF] Tripathy | Privacy-preserving adversarial networks[END_REF]. For both networks, a LSTM architecture is selected (see Fig. 2), which was shown to be successful in several problems dealing with sequences of data (see [START_REF] Goodfellow | Deep Learning[END_REF] and references therein for more details). The training of the suggested framework is performed using Algorithm 1 which uses k gradient steps to train A followed by one gradient step to train R ✓ . Note that k should be large enough to ensure that A is a strong adversary during training. This is in fact the common practice for effectively training two networks in an adversarial framework [START_REF] Goodfellow | Generative adversarial nets[END_REF] and, in particular, in privacy scenarios [START_REF] Tripathy | Privacy-preserving adversarial networks[END_REF]. It should be recalled that, after the training of both networks is completed, an attacker network is trained in order to test the privacy achieved by the releaser network. It should be clarified that this attacker network is distinct from the adversary network used during training and illustrated in Fig 1: the attacker used in testing mimics a real-world attacker that would try to deduce the private date from the release data.

IV. RESULTS AND DISCUSSION

A. Description of datasets

In this study, the Electricity Consumption & Occupancy (ECO) and Pecan Street data sets are used. The ECO data set, collected and published by [START_REF] Beckel | The eco data set and the performance of non-intrusive load monitoring algorithms[END_REF], includes 1 Hz power consumption measurements and occupancy information of five Algorithm 1 Algorithm for training privacy-preserving data releaser neural network. Input: Data set (which includes sample sequences of useful data y T , sensitive data x t ); seed noise samples u T ; seed noise dimension m; batch size B; number of steps to apply to the adversary k; gradient clipping value C; L 2 recurrent regularization parameter .

Output: Releaser network R ✓ .

1: for number of training iterations do 2:

for k steps do Compute the gradient of L A ( ), approximated with the minibatch B, w.r.t. to .

5:

Update the adversary by applying the RMSprop optimizer with clipping value C. Compute the gradient of L R (✓, , ), approximated with the minibatch B, w.r.t. to ✓.

9:

Use Ridge(L 2 ) recurrent regularization with value and update the releaser by applying RMSprop optimizer with clipping value C. 10: end for Fig. 2. LSTM recurrent network cell diagram. The cell includes four gating units to control the flow of information. All the gating units have a sigmoid activation function ( ) except for the input unit (that uses an hyperbolic tangent activation function (tanh) by default). The parameters b, V, W are respectively biases, input weights, and recurrent weights. In the LSTM architecture, the forget gate ft = (b f + K f h t 1 + V f wt) uses the output of the previous cell (which is called hidden state h t 1 ) to control the cell state Ct to remove irrelevant information. On the other hand, the input gate gt = (b g + K g h t 1 + V g wt) and input unit adds new information to Ct from the current input. Finally, the output gate ot = (b o +K o h t 1 +V o wt) generates the output of the cell from the current input and cell state.

houses in Swiss over a period of 8 months. In this study we resampled the data to have hourly samples. On the other hand, the Pecan Street data set contains hourly SMs data of houses in Texas, Austin and was collected by Pecan Street Inc. [START_REF]Dataport: the world's largest energy data resource[END_REF]. Pecan Street project is a smart grid demonstration research program which provides electricity, water, natural gas, and solar energy generation measurements for over 1000 houses in Texas, Austin. In order to model time dependency over each day (with a data rate of 1 sample per hour), the data was reshaped to sample sequences of length 24. For the ECO and Pecan Street data set, a total number of 11225 and 9120 sample sequences are used, respectively. The data is splitted into train and test sets with a ratio of roughly 85 : 15 while 10% of the training data is used as the validation set. It should be noted that in this study, we assume that the attacker has access to all the training data used by the releaser, which can be considered as a worst-case scenario study.

B. Inference of households occupancy

The first practical case of study regarding privacy-preserving in time series data is the concern of inferring presence/absence of residents at home from the total power consumption collected by SMs [START_REF] Kleiminger | Household occupancy monitoring using electricity meters[END_REF], [START_REF] Jia | Human-factor-aware privacy-preserving aggregation in smart grid[END_REF]. For this application, the electricity consumption measurements from the ECO data set are considered as the useful data, while occupancy labels are considered as the private data. Therefore, our privacy-preserving data release method aims to minimize a trade-off between the distortion of the total electricity consumption incurred and the probability of inferring the presence of an individual at home from the release signal. The releaser and adversary networks used for the training consist of 4 LSTM layers with 64 cells and 2 LSTM layers with 32 cells, respectively where a tanh activation function used. In addition, recurrent regularizer with parameter = 1.5 was used in each layer of the release network. The values of the other hyperparameters (B, k,m) were set to (128, 4, 8), respectively. Finally, after training, a strong attacker is used, consisting of 3 LSTM layers.

Based on the target data Y T and the release data Z T , the normalized root mean-square-error (NRMSE) is defined by

NRMSE , v u u u t E h Y T Z T 2 i E h Y T 2 i . (10)
Fig. 3 shows the empirically found privacy-utility trade-off for this application. It can be seen that by adding more distortion on the released data, the attacker is pushed toward a random guessing classifier.

In order to provide more insights about the release mechanism, the Power Spectrum Density (PSD) of the input signal and the PSD of the error signal (defined as the difference Fig. 3. Privacy-utility trade-off for house occupancy inference application. Since in this application the attacker is a binary classifier, the random guessing (balanced) accuracy is 50%. The fitted curve is based on an exponential function and is included only for illustration purposes.

between the actual power consumption and the released signal) for four different cases along the privacy-utility trade-off curve of Fig. 3 are estimated using Welch's method [START_REF] Stoica | Spectral Analysis of Signals[END_REF]. It should be mentioned that two stationary tests including the Augmented Dickey-Fuller test [START_REF] Dickey | Distribution of the estimators for autoregressive time series with a unit root[END_REF] and the Kwiatkowski, Phillips, Schmidt, and Shin (KPSS) test [START_REF] Kwiatkowski | Testing the null hypothesis of stationarity against the alternative of a unit root: How sure are we that economic time series have a unit root?[END_REF] applied to our data set indicates that there is enough evidence to suggest the data is stationary, supporting our PSD analysis.

C. Inference of house identity

The second practical case of study regarding the privacypreserving in SMs measurements is identity recognition from total power consumption of households [START_REF] Efthymiou | Smart grid privacy via anonymization of smart metering data[END_REF]. It is assumed that the attacker has access to total power consumption of different households in a region (training data) and then attempts to determine identities of the households using the new released data (test data). Thus, our model aims at generating release data of total power consumption of households in a way that prevents the adversary to perform the identity recognition while keeping distortion on the total power minimized. For this task, total power consumption of five houses is used. For this application, the releaser consists of 6 LSTM layers each includes 128 cells and adversary has 4 LSTM layers with 32 cells. Similarly, a tanh activation function was applied and = 2 was used in each layer of the release network. The values of the other hyperparameters (B, k,m) are set to (128, 5, 3), respectively. Finally, after training, an attacker, consisting of 4 LSTM layers was used. The empirical privacyutility trade-off curve obtained for this application is presented in Fig. 5. Comparing Fig. 5 with Fig. 3, we see that a high level of privacy is expensive. For instance, in order to obtain an attacker accuracy of 30 %, the NRMSE should be approximately equal to 0.30. This is attributed to the fact that this task is harder from the learning point of view than the one considered in Section IV-B.

PSD analysis was also performed for this application, yielding the results of Fig. 6. Once again, we see that the release network provides privacy-utility trade-off by mainly distorting the harmonics on the actual electricity consumption signal.

V. CONCLUSION

We have presented a new method to train privacy-preserving mechanisms controlling the privacy-utility trade-off in time series data. This lead us to define the directed information between sensitive variables and their estimation as a more suitable privacy measure than previous proposals in the literature. A tractable upper bound was then derived and a deep learning adversarial framework between two recurrent neural networks was introduced to optimize the new loss function. Our method was validated with two well-known privacy problems in smart meters data using two different open data sets. For both privacy problems we considered the worst-case where an attacker has access to all the training data used by the releaser. In future work, we will consider alternative formulations of the problem such as different distortion measures and a more general loss function in order to attempt to provide universal privacy guarantees (i.e., independent of the attacker structure and computational power).

Fig. 1 .

 1 Fig. 1. Privacy-Preserving framework. The seed noise U T is generated from i.i.d. samples according to a uniform distribution: Ut ⇠ U [0, 1].

3 :⌘ ; b = 1 , 2 ,

 312 Sample minibatch of B examples: B = n w (b)T = ⇣ x (b)T , y (b)T , u (b)T

7 :

 7 Sample minibatch of B examples: B = n w (b)T = ⇣ x (b)T , y (b)T , u (b)T ⌘ ; b = 1, 2, .., B o .8:

  For each case, we use 10 release signals and average the PSD estimates. Results are shown in Fig. 4. Looking at the PSD of the input signal (useful data) some harmonics are visible. The PSD of the error signals show that the model controls the trade-off in privacy-utility by mainly modifying the distortion on these harmonics.

Fig. 4 .

 4 Fig. 4. PSD of the actual electricity consumption and error signals for the house occupancy inference application.

Fig. 5 .

 5 Fig. 5. Privacy-utility trade-off for house identity inference application. Since in this application the attacker is a five-class classifier, the random guessing (balanced) accuracy is 20%. The fitted curve is based on an exponential function and is included only for illustration purposes.

Fig. 6 .

 6 Fig. 6. PSD of the actual electricity consumption and error signals for the house identity inference application.

ACKNOWLEDGMENT This work was supported by Hydro-Quebec, the Natural Sciences and Engineering Research Council of Canada, and McGill University in the framework of the NSERC/Hydro-Quebec Industrial Research Chair in Interactive Information Infrastructure for the Power Grid (IRCPJ406021-14). This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 792464.