Thibaut Barbier
email: thibaut.barbier@mines-paristech.fr

Robin Girard

Nicolas Kong

François-Pascal Neirac

Georges Kariniotakis

Elena Magliaro

AN OPTIMISATION ALGORITHM FOR MATCHING LARGE SCALE DATABASES ON CUSTOMERS FOR IMPROVED CHARACTERISATION OF ELECTRICITY CONSUMPTION

Keywords: database matching, electricity consumption, electricity customers, housings, optimization algorithm

This paper presents a method that permits to match customer information from the French DSO Enedis and housing information from the French population census institute INSEE. Our method allows having a list of housings linked to each customer in order to add household and building information to customers. We show with our method improvements in predictions of aggregated load curve indicators compared to the traditional method that averages socio demographic indicators from housing information of the zone covered by measurements. Our results indicate that the proposed algorithm is able to capture efficiently the information of housings in some feeders. This permits to combine the databases of the DSO with external databases that exist from census or other processes. Enriching the information at the level of clients through the proposed automated way is a cost effective approach given the number of customers served by a DSO. This enhanced information can be then the basis to model, analyse and simulate demand in a bottom up approach which can be useful for planning purposes of the distribution networks.

Introduction

Most of scientists agree that human being has influenced the climate system, especially due to anthropogenic emissions of greenhouse gases. IPCC [START_REF]Intergovernmental Panel on Climate Change[END_REF] showed that in recent decades, climate change has damaged natural and human ecosystems in all continents, oceans and atmosphere. These impacts will continue for centuries, even if anthropogenic emissions are stopped. Climate change will exacerbate current issues in many aspects (with very high confidence), such as health problems, migrations, weather extreme events, and this phenomenon of amplification will increase quickly during the 21st century. Adaptation is a strategy for managing these issues. Also mitigating the magnitude of global warming by lowering those emissions is crucial, in order to delay the risk of severe, widespread and irreversible impacts globally [START_REF]Intergovernmental Panel on Climate Change[END_REF].

Electricity represents 18% of the total final energy consumption in the world in 2013 [START_REF]Key word Energy Statistics[END_REF] and is expected to rise up to 25% by 2040 [3]. Electricity and heat was 42% of world CO2 emissions in 2013, i.e. 13.8 Giga tons of CO2 [START_REF]Key trends in CO2 emissions[END_REF]. Such awareness has led to consider larger and larger shares of the electricity demand as controllable. This offers many opportunities [START_REF] Siano | Demand response and smart grids. A survey[END_REF] such as an increased potential for energy saving, the possibility to adapt electricity demand to intermittent renewable energy, or a lowered of peak demand to prevent costly investments on the electricity grid. In order to be able to exploit these opportunities and integrate them in the management and especially the planning processes of the power system, it is necessary to have an in depth knowledge of the characteristics of the demand not only in terms of energy but also in terms of power. However, electricity consumption, especially at the city scale, is the sum of a considerable number of loads, the characteristics of which are often unknown.

Characterizing each contributor of the electricity consumption is not a recent topic. Swan and Ugursal [6] present examples of data required to develop residential energy model. It includes information on the physical characteristics of the dwellings, occupants and their appliances, historical electricity consumption, climatic conditions, and macroeconomic indicators.

DSOs (Distribution System Operators) dispose large amounts of data on energy consumption and aggregated load curves. However, when it goes down to the layer of buildings, information is lacking [6]. Given the number of clients it is very costly to implement such data collection procedure. This type of data is often publicly available by external organizations. It is a challenge to match the information in such external databases to the data in the DSO databases.

In this paper, we present a method that permits to match customer information from the French DSO Enedis and housing information from the French population census institute (INSEE).

We define housing as houses or dwelling places thought of as a group. A housing occupied as a principal residence has a household.

Electricity customer (or customer) materializes a unique contract with the DSO who is paid in order to provide an electricity supply service at a point of common coupling where the customer belongs.

Our method allows having a list of housings linked to each customer in order to add household and building information to customers. We call the information about building and household "socio demographic information". When this information is gathered at district scale or at a zone covered by an electric line we call it "socio demographic indicator".

We show with our method improvements in predictions of aggregated load curve indicators compared to method chosen as reference that averages socio demographic indicators of the zone covered by measurements. For each district we can have from external databases, like the ones resulting from population census, the list of the housings. Electricity customers can be either residential or not. If a customer is residential, it implies that it is linked to one or multiple housings. But also a single housing can have different electricity contracts and so can be linked to different customers.

Illustration of the problem

Our method proposes to match customers and housings, in order to have a list of housing linked to residential customers. At least, as we know from the DSO which customer is linked to which feeder we can have a list of housing per electric feeder 3 Problem solving method

3.1

Estimation of common parameters Electricity customers and housings have different characteristics listed in Table 1. This table explains the difference between customer and housing in terms of information. We can notice that housings and electricity customers have no common parameters except that we both know the district where they are.

The idea is to create common parameters between the two databases in order to measure a distance and be able to match them.

We choose to estimate the yearly energy consumption of housings as a linear model of surface, presence of an electric heating, age of building, number of occupants and number of rooms. The model is fitted with aggregated information on each district, which was presented in a previous work [START_REF] Barbier | A novel approach for electric load curve holistic modelling and simulation[END_REF].

3.2

Matching data with the common parameters The problem set is as following: we want to match the housings and the residential electricity customers. Housings' annual electricity consumption was estimated with the method described in the previous section. We define below the list of parameters to be matched from the two considered databases: 𝑥 𝑖∈{1,2,…,𝑛} represents the estimated electricity annual consumption of the housings, n is the number of housings in the considered district. 𝑥 𝑖 is the i th housing. 𝑦 𝑗∈{1,2,…,𝑝} represents the residential annual consumption of the electricity customers. p is the number of customers in the considered district. 𝑦 𝑗 is the j th electricity customer.

The problem of matching can be formulated as a distance (here chosen to be quadratic) minimization problem between the matched entities:

𝛼 𝑖𝑗 ∈ {0,1}, ∀ (i, j) ∈ [1: n]x[1: p] (2) ∑ 𝛼 𝑖𝑗 𝑛 𝑗 ≥ 1 ∀i ∈ [1: n] (3) ∑ 𝛼 𝑖𝑗 𝑛 𝑖 =1 ∀j ∈ [1: p] (4)
where: [1:p] is the matching unknown:

 (α ij ) (i,j)∈[1:n]x
α ij = 1 if x i is matched with y j , 0 else (2); 
(3) is the constraint that each housing is linked to at least an electricity customer;  (4) is the constraint that each electricity customer is linked an unique housing. This problem is a classical quadratic function to minimize (1). However, ( 2) is a non-continuous constraint. Solving method by relaxing the constraint (2) in an interval and putting additional constraints to force solution to be closed to boundaries 0 and 1 where tested, but without any success: solutions where not closed to the boundaries.

Testing each possibilities of (α ij ) (i,j)∈[1:n]x[1:p] leads to too long calculations as the number of possibilities is 2 𝑛.𝑝 . The approximate method of resolution chosen is as following: we first set the initial state by ordering data to match and make a random gathering in order to have the same number of data to match 𝑦 𝑗∈{1,2,…,𝑝} and 𝑥 ̃𝑗∈{1,2,…,𝑝} , as illustrated figure 2. Then we optimize the global distance by moving the 𝑥 𝑖 that minimizes the global distance ∑ (𝑥 𝑗 ̃-𝑦 𝑗 )² 𝑝 𝑗=1

. We finally stop when the global distance remains the same. As shown figure 3, housings are now linked to customers. Data and reference of comparison Data were collected every 10 minutes by Enedis (main French DSO) for the whole year 2012. We used 59 feeder measurements from the region of Lyon, in France. For each feeder we have the list of customers supplied, whose characteristics are listed in Table I. The characteristics of the feeders are as following: For each district supplied by the feeders, we have also a list of housings with some characteristics listed in Table I. This data are provided by the French institute of statistics (INSEE).

The application of the matching algorithm between customers and housings for each district permits to derive a list of housings per feeder. For example in figure 3, the matching procedure for Feeder1 permits to identify that this feeder supplies housings H1, H3, H4, H5, and H6. With this list of housings it is possible to characterize socio-demographically each feeder. This is done through averaging the values of the parameters (as per Table I) of the associated feeder housings. At the scale of the feeder we also calculate reference sociodemographic indicators, which are the sum of average values of socio demographic indicators for each district weighted by the ratio of residential customers in each district.

Model to explain thermo sensitivity using socio-demographic indicators

Thanks to the previous section we have associated to each feeder additional useful information resulting from the reference calculation and the matching. Now we use this information to model the thermo-sensitivity per feeder. The thermo-sensitivity chosen here is the linear fit of consumption points in ordinate and temperature points on the abscissa, by only selecting points whose temperature are below a certain threshold of nonheating. Figure 4 shows the thermo-sensitivity plot of a feeder. In this figure we see the link between electricity and temperature, which represents well the French behaviour: because of electric heaters the consumption rises with temperature decrease. This behaviour is similar to other intermediate temperature European countries (averages temperature between 9°C and 14°C) [START_REF] Bessec | The non-linear link between electricity consumption and temperature in Europe: A threshold panel approach[END_REF]. This information can be useful for a DSO in order to size its grid as in this case peak of consumption is due to cold temperatures.

Quantitative interest of matching to explain thermo-sensitivity

We display prediction errors made by the linear model fitted on the data available.

We compare the errors of thermo-sensitivity predicted by using socio-demographic indicators taken as reference (average values in each district covered by a given feeder) and by using socio-demographic indicators from matching algorithm. As shown figure 6, the error of prediction is 0.8% less

Conclusion

In this work we propose a method to match housing information with residential customers. We test the interest of this method: we compare errors of prediction of thermo-sensitivity of electric feeders with sociodemographic data from average values of the zone and from data resulting from the matching algorithm. We show that the matching algorithm improves prediction, meaning that the matched information is closer to real customers' housing information than the averaged values. Figure 7 shows the distribution of the relative improvement for each feeder. Our results indicate that the proposed algorithm is able to capture efficiently the information of housings in some feeders. This permits to combine the databases of the DSO with external databases that exist from census or other processes. Enriching the information at the level of clients through the proposed automated way is a cost effective approach given the number of customers served by a DSO. This enhanced information can be then the basis to model, analyse and simulate demand in a bottom up approach which can be useful for planning purposes of the distribution networks. Further work will focus on investigating the local characteristics of the feeders for which thermo-sensitivity is well predicted, on testing the matching algorithm on higher amount of important with matched data compared to reference average data. data and multiple years and regions to validate the results. This method allows also estimating thermosensitivity for scales from districts to regions.

Figure 1

 1 Figure 1 illustrates an area of 5 districts where the electricity customers are supplied by the different feeders of a same substation. This is typical information we can have from DSO's databases. Some customers have the same address (e.g. C3-C4 and C11-C12): it represents different customers in the same building.

Fig. 1 :

 1 Fig.1: available information on the left (DSO) and on the right (population census database)

  min (𝜶 𝒊𝒋 ) (i,j)∈{1,2,…,n}x{1,2,…,p} [ ∑ (∑ 𝛼 𝑖𝑗 𝑥 𝑖 -𝑦 𝑗 𝑛

Fig. 2

 2 Fig. 2 Illustration of the different steps of the matching between customers and housings

Fig. 4

 4 Fig. 4 example of thermo-sensitivity plot of a feeder In order to quantify the impact of socio-demographic indicators on the consumption thermo-sensitivity, we use a linear model as shown figure 5.

Fig. 5

 5 Fig. 5 Illustration of the different steps of the matching between customers and housings

Fig. 7

 7 Fig. 7 Histogram of the relative amelioration of errors with data from matching method compared to reference average data.

Table 1

 1 List of the different characteristics of customers and housings

	Database	Useful characteristics to characterize
	entity	electricity consumption
	Electricity	Subscribed	power,	annual
	customer	consumption, text address, district,
		type of customer (e.g. residential,
		tertiary artisan or shopkeeper).
	Housing	Number of inhabitants, type (flat,
		house), surface, number of rooms, age
		of the building, category (vacant,
		secondary	residence,	principal
		residence), type of heating (electric,
		fuel, gas, wood, other), district.

Table 1

 1 Characteristics of the feeders data set used

	Indicator	Min	Mean	Max
		value	value	value
	Number of customers	330	3400	7000
	supplied
	Fraction of residential	42%	67%	88%
	consumption
	Number of districts	2	5	13
	supplied

Acknowledgements

We thank Enedis for providing the data for this work.


Simon Jones
email: jones5@rose-hulman.edu

Mathias Legrand
email: mathias.legrand@mcgill.ca

VIBRATIONS OF AN AXIAL BAR EXPERIENCING PERIODIC UNILATERAL CONTACT USING THE WAVELET BALANCE METHOD

Efficiently predicting the vibratory responses of flexible structures which experience unilateral contact is becoming of high engineering importance. An example of such a system is a rotor blade within a turbine engine; small operating clearances and varying loading conditions often result in contact between the blade and the casing.

The method of weighted residuals is a effective approach to simulating such behaviour as it can efficiently enforce timeperiodic solutions of lightly damped, flexible structures experiencing unilateral contact. The Harmonic Balance Method (HBM) based on Fourier expansion of the sought solution is a common formulation, though it is hypothesized wavelet bases that can sparsely define nonsmooth solutions may be superior. This is investigated herein using an axially vibrating rod with unilateral contact conditions. A distributional formulation in time is introduced allowing periodic, square-integrable trial functions to approximate the second-order equations. The mixed wavelet Petrov-Galerkin solutions are found to yield consistent or better results than HBM, with similar convergence rates and seemingly more accurate contact force prediction.

INTRODUCTION

Efficiently predicting the vibratory responses of flexible structures which experience unilateral contact is becoming of high engineering importance. For example, consider aircraft engines where slender, twisted blades rotate at high rotational velocity within stationary casings where minimal clearance is desirable for turbine energy efficiency. Simulating contact between the blades and casing is not a trivial exercise since unilateral contact is usually described by inequalities and complementary conditions [START_REF] Wriggers | Computational Contact Mechanics[END_REF]. In the time domain, structural displacements and velocities which satisfy these contact conditions are known to respectively feature absolute continuity and bounded variation only [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous galerkin methods: Definition and outlook[END_REF]. This implies displacements are not necessarily differentiable everywhere in the domain and velocities may exhibit jumps; these types of problems are generally referred to as nonsmooth [START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF][START_REF] Leine | Dynamics and bifurcations of non-smooth mechanical systems[END_REF].

This class of unilateral problems can also be approached using periodic vibration theory. This allows the original initial-value formulations to be transformed into partial differential equations periodic in time [START_REF] Laxalde | Nonlinear modal analysis of mechanical systems with frictionless contact interfaces[END_REF]. In structural mechanics, two families of numerical techniques can efficiently predict periodic solutions in time. The first is the shooting method which consists of finding the initial conditions that result in periodic motion. The overall approach relies on time integration of the governing equations over one period using a nonlinear solver which iterates on the initial conditions, and possibly the period. The second technique is based on weighted residual formulations, which are of interest in the current investigation. This method involves approximating the solution using a set of time-dependent basis functions, called trial functions, and enforcing the respective residual error to be orthogonal to an independent set of weighting functions [START_REF] Lapidus | Numerical Solution of Ordinary Differential Equations[END_REF][START_REF] Reddien | Projection methods for two-point boundary value problems[END_REF]. Unlike the shooting method which can become numerically sensitive to possible jumps in the velocity field, weighted residual techniques directly enforce the periodicity conditions while the remaining unilateral contact constraints and governing local equations of motion are satisfied in a weak integral sense [START_REF] Wang | Dynamics of unsymmetric piecewiselinear/non-linear systems using finite elements in time[END_REF]. It is worth noting that the well-known Harmonic Balance Method (HBM) is a weighted residual approach where both the trial and weighting functions are identical Fourier series. The main goal of the current work is to explore relevant basis functions whose order of smoothness can be adapted to a particular system to attain accurate approximations and rapid convergence.

The paper is broken into five sections: Section 2 introduces the system of interest; Section 3 reviews the weighted residual formulation in time; Section 4 provides an overview of the different basis functions under investigation; Section 5 specifies the unilateral contact model; and Section 6 presents and compares results for a number of basis function combinations.

System of interest

A schematic of a simple unilateral contact system showing a one-dimensional bar clamped on the left is depicted in Fig. 1; such a system is known to be a well-posed unilateral contact problem [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF]. A gap g exists between the tip and a rigid foundation. When the tip displacement is sufficiently large, the bar comes into contact with the rigid foundation and unilateral contact conditions must be satisfied. The existence of periodic solutions of period T are assumed where T is the period of the external forcing f .x; t / acting on the bar. Accordingly, the unknown displacement has to satisfy the following complementary boundary value problem:

Here signifies density, S cross-sectional area, and E elastic modulus. The dot superscript represents a temporal derivative, whereas the subscript x represents a spatial derivative. Structural damping is incorporated later by assuming ˇ-damping of the elastic modulus, as detailed in Section 2.

The periodicity conditions result in a problem which can be formulated on a circle in time. Without loss of generality, the basis functions are taken from an L 2 .S 1 / N Hilbert space [START_REF] Holschneider | Wavelet analysis on the circle[END_REF], where the period of the steady-state solution T has been normalized to 1 [START_REF] Berberian | Introduction to Hilbert space[END_REF]. The quantity .t / in Eq. ( 4) is a contact force, which stems from the enforced non-penetration condition g u.`; t / 0 and is necessarily positive (by convention). The complementarity condition .t / g u.`; t / D 0 states that the contact force .t / and the distance g u.`; t / separating the rod's end-tip from the rigid foundation may not be zero at the same time. These three conditions are such that the mathematical object pairing the contact force to the displacement is not a function in the usual sense. It is also known that the displacement field may be non-differentiable at given instants within the period of motion and the velocity field may be discontinuous. This motivates the derivation of numerical techniques capable of efficiently handling nonsmoothness. As a first approach, the unilateral contact inequalities from Eq. ( 4) are simplified and replaced by a penalty function; an exponential spring K p is used to approximate the contact forces. The gap g in Fig. 1 separates the tip of the rod and the exponential spring. The penalty function is of the form f c .u.`; t // D max a c .e ˛.u.`;t / g 0 / 1/; 0 :

(

Specific values for the penalty function are provided with the other model parameters in Section 4.

Weighted residual formulations in time

The method of weighted residuals is a classic method of obtaining numerical solutions to boundary value problems by expanding the sought solution as a finite sequence of basis functions (i.e. trial functions) in a proper functional space [START_REF] Lapidus | Numerical Solution of Ordinary Differential Equations[END_REF]. The subsequent residual is rendered orthogonal to a set a linearly independent functions of the same space (i.e. weightnig functions) through an inner product. Trial and weighting functions may contain the same basis function but not necessarily: The Galerkin method is a special case where the weighting functions are taken from the same function basis as that of the trial functions; the Petrov-Galerkin method involves the selection of weighting function comprising a basis which is independent of the trial function [START_REF] Durran | Numerical methods for wave equations in geophysical fluid dynamics[END_REF]. The functional space into which the solution is sought is such that the boundary conditions in space and time are satisfied. In this study, the displacement of the investigated rod is numerically sought in various Hilbert spaces of square integrable functions defined on the circle S 1 .

To solve Eqs. ( 1), ( 2), (3), and (4), the unknown displacement is expanded into a truncated series of N functions separated in space and time:

u.x; t / D N X i D1 #.x/u i .t / (6)
The standard Finite Element Method is implemented for the spatial variable using two-node linear rod elements [START_REF] Bathe | Finite element procedures in engineering analysis[END_REF][START_REF] Liu | Finite element method: a practical course[END_REF]. This yields the following vector ordinary differential equation of size N :

M R u.t / C C P u.t / C Ku.t / D f ext .u.t /; t / (7)
together with the remaining periodicity conditions in time and the unilateral contact condition. It should be noted that Eq. ( 7) is an approximate equation of a vibrating structure which may be of one, two or three dimensions. The system is constrained so that condition (3) is now satisfied. Here M and K are the standard mass and stiffness matrices for a rod, and ˇ-damping is enforced such that C D ˇK to account for structural damping. The displacement vector u.t / in Eq. ( 7) stores the temporal unknowns u i .t /, i D 1; : : : ; N , where N is the number of spatial degrees-of-freedom. Similarly, f ext .u.t /; t / stores the external forcing functions for each degree-of-freedom as well as the contact forces stemming from Eq. ( 5). The remainder of the derivation involves a weighted residual formulation in the time dimension. Three forms of the weighted residual method are discussed: the strong integral form, the weak form, and the distributional formulation. For this discussion it is understood that u and v are time-dependent vectors, i.e. u.t / and v.t /, but for the sake of clarity the .t / is omitted below.

Strong integral form

The standard weighted residual formulation of a differential equation is commonly termed the strong form. Taking the inner product of Eq. ( 7) with a weighting function stored column-wise in v results in the strong integral form of the equation: find u 2

H 2 .S 1 / N such that 8v 2 L 2 .S 1 / N : Z S 1 v > M R u C v > C P u C v > Ku v > f ext .u; t / dt D 0 (8)
and the superscript T denotes a transpose. This strong form of the equation is not necessarily the best framework for obtaining a solution [START_REF] Liu | Finite element method: a practical course[END_REF]; for this example involving a vector ordinary differential equation of order 2, the solution must be at least H 2 , limiting the permissible basis of trial functions.

Weak form

The respective weak form of the weighted residual statement can be obtained by performing one integration by parts over the domain S 1 for all terms containing a double time derivative in Eq. ( 8). This results in: find u 2 H 1 .S 1 / N such that 8v 2

H 1 .S 1 / N , Z S 1 P v > M P u C v > C P u C v > Ku v > f ext .u; t / dt D 0 (9)
The integral form of the weak formulation offers the advantage of shifting a portion of the functional smoothness requirement from the trial functions onto the weighting functions. More precisely, both the trial and weighting functions must now be H 1 . This allows the trial functions to be chosen from a wider permissible space [START_REF] Liu | Finite element method: a practical course[END_REF].

Formulation in a distributional sense

The above procedure of obtaining the solution from the weak formulation by performing integration by parts on the weighted residual statement can be extended one step further. Theoretically, this extension is not necessary since it is known that displacements should be absolutely continuous in time, though it may assist in the numerical derivations and allow very simple basis functions to be considered [START_REF] Studer | Numerics of unilateral contacts and friction[END_REF]. A weaker formulation is proposed by integrating again the terms involving time derivatives of the trial functions. This formulation can be understood in the sense of distributions, also known as generalized functions, ie:

find u 2 L 2 .S 1 / N such that 8v 2 H 2 .S 1 / N : Z S 1 R v > Mu P v > Cu C v > Ku v > f ext .u; t / dt D 0 (10)
Here the double time differential on the field variable is transferred to the weighting function and the continuity requirement on the trial function is reduced. As discussed later, the desired displacement functions u can now be described using a series of constant piecewise functions for instance.

Time discretization

Let each nodal displacement u i .t / be approximated by a combination of M linearly independent trial functions q k .t /, k D 1; : : : ; M defined on t 2 OE0; 1/. Similarly, let each v i .t / be approximated by a combination of linearly independent weighting functions p k .t /, k D 1; : : : ; M also defined on t 2 S 1 . Collectively, this leads to

u i .t / D M X kD1 a i k q k .t /; i D 1; : : : ; N (11)
and

v i .t / D M X kD1 b i k p k .t /; i D 1; : : : ; N (12)
or, equivalently in a vector form:

u.t / D '.t /a and v.t / D .t /b (13)
In Eq. ( 13), '.t / and .t / are rectangular matrices of dimension N NM ; a and b are vectors of size NM 1.

Each of the strong, weak, and distributional formulations can now be discretized with proper basis functions. The overall goal of this paper is to discuss the numerical properties of certain wavelet bases to efficiently solve the problem of interest. Depending on the selected formulation, the order of derivation acting on u.t / and v.t / (i.e. '.t / and .t / respectively) will affect the admissible basis functions. The corresponding discretized versions are:

Strong integral form:

Z S 1 > M R ' C > C P ' C > K' a > f ext .'a; t / Á dt D 0 (14)
Weak form:

Z S 1 P > M P ' C > C P 'C > K' a > f ext .'a; t / Á dt D 0 (15)
Distributional integral form:

Z S 1 R > M' P > C' C > K' a > f ext .'a; t / Á dt D 0 (16)
The resulting nonlinear equations can generically be recast in the following form

Ga f.a/ D 0 (17)
where Ga and f.a/ respectively stand as the linear internal and the nonlinear external contributions in Eq. ( 14), [START_REF] Studer | Numerics of unilateral contacts and friction[END_REF], and ( 16).

Trial and weighting function bases

The selection of functional bases to be used in the above approaches is an important factor in approximation accuracy and computational efficiency [START_REF] Bathe | Finite element procedures in engineering analysis[END_REF][START_REF] Liu | Finite element method: a practical course[END_REF]. A priori selection of the optimal bases for unilateral contact problems is not a trivial matter due to potential nonsmoothness of the response. In the current investigation a number of bases are investigated to compare the quality of approximation, including Fourier functions, B-spline wavelets, Daubechies wavelets, and Haar wavelets.

Harmonic balance and Fourier series

The harmonic balance method (HBM) is a special case of the weighted residual method where the Fourier basis is used for both the trial and weighting functions [START_REF] Mickens | Comments on the method of harmonic balance[END_REF]. This technique is particularly effective when dealing with smooth nonlinear systems; convergence is often reached with very few terms. It has been used to study steady state response of turbine engine blades with friction dampers using a multiterm approximation [START_REF] Wang | Investigation of the vibration of a blade with friction damper by HBM[END_REF]. This approach was extended to unilateral contact and friction conditions [START_REF] Legrand | n-dimensional Harmonic Balance Method extended to non-explicit nonlinearities[END_REF] through an Alternating Frequency/Time domain strategy proposed by Cameron et al. [START_REF] Cameron | An alternating frequency/time domain method for calculating the steadystate response of nonlinear dynamic systems[END_REF] and Pilipchuk [START_REF] Pilipchuk | Application of special nonsmooth temporal transformations to linear and nonlinear systems under discontinuous and impulsive excitation[END_REF]. It is worthy to note that the HBM formulations of Eqs. ( 14), [START_REF] Studer | Numerics of unilateral contacts and friction[END_REF], and ( 16) are identical.

The Fourier basis for L 2 .S 1 / is defined as f1g [ fcos.2m t /; sin.2m t /jm 2 N g where m signifies the harmonic number of the function. The first six functions are shown in Fig. 2. Fourier basis functions feature an infinite degree of

0 1 -1 0 1 (a) constant term 0 1 -1 0 1 (b) cos.2 t / 0 1 -1 0 1 (c) sin.2 t / 0 1 -1 0 1 (d) cos.4 t / 0 1 -1 0 1 (e) sin.4 t / 0 1 -1 0 1 (f) cos.6 t /

FIGURE 2. FIRST SIX FUNCTIONS OF THE FOURIER BASIS

smoothness. While this property can be beneficial in some cases, it is unclear whether Fourier functions are optimal when simulating potentially nonsmooth problems. It is possible that a large number of harmonics would be required to accurately capture the nonsmooth response, or the approximation may exhibit Gibbs phenomenon at localized discontinuities in the sought displacement and velocity fields. Accordingly, other basis functions featuring a lesser degree of functional smoothness are explored.

Brief review of discrete orthogonal wavelets

Discrete orthogonal wavelet families are composed of highly localized, oscillatory functions which provide a basis of L 2 .R/ and can be adapted to the periodic domain L 2 .S 1 / [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. These localized characteristics, or compact support, allow sparse repre-sentation of piecewise signals including transients and singularities [START_REF] Burrus | Introduction to wavelets and wavelet transforms: a primer[END_REF]. This makes them useful functions for use in the Galerkin approach when nonsmooth solutions are predicted [START_REF] Strang | Wavelets and Filter Banks[END_REF]. There are a large number of wavelet families and definitions thereof; both Mallat [START_REF] Mallat | A wavelet tour of signal processing[END_REF] and Strang [START_REF] Strang | Wavelets and Filter Banks[END_REF] provide excellent introductions to wavelet theory and history. Galerkin methods using appropriate discrete wavelet families as the trial functions have been shown to accurately approximate the solutions to both ordinary and partial differential equations [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF][START_REF] Beylkin | On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases[END_REF][START_REF] Pernot | A wavelet-Galerkin procedure to investigate time-periodic systems: transient vibration and stability analysis[END_REF][START_REF] Williams | Introduction to wavelets in engineering[END_REF][START_REF] Christon | The numerical performance of wavelets for PDEs: the multi-scale finite element[END_REF][START_REF] Burgos | Direct solution of differential equations using the wavelet-Galerkin method[END_REF].

The discrete wavelet family is built from scaling functions .t / and wavelet functions .t /. These functions are analogous to the low-pass and high-pass filters of a filter bank [START_REF] Strang | Wavelets and Filter Banks[END_REF]; decomposition using scaling functions will give a "smoothed" approximation of the original signal, while decomposition using wavelet functions provides the details of the signal, or high-frequency content.

The exact decomposition of a continuous time signal y.t / can be written

y.t / D X k g k J;k .t / C m X j DJ X `hj;` j;`. t / (18)
and

J;k .t / D 2 J =2 .2 J t k/ (19)
J;k .t / D 2 J =2 .2 J t k/ (20)
where J; k 2 Z; J is the dilation parameter (i.e. level), k is the translation parameter, and m is the maximum resolution given by the sampling rate of the function y.t /.

The span of the scaling functions at level J is commonly denoted V J , while the wavelet span is denoted W J . For orthogonal wavelet families, W J is the orthogonal complement to

V J in V J C1 V J C1 D V J ˚WJ : (21)
Provided the wavelet family is orthogonal [START_REF] Strang | Wavelets and Filter Banks[END_REF], the space of square-integrable functions on the real line L 2 .R/ can be decomposed using multiresolution analysis as a nested sequence of closed subspaces [START_REF] Mallat | A wavelet tour of signal processing[END_REF]

V 2 V 1 V 0 V 1 V 2 L 2 .R/ (22) such that lim J !1 V J D L 2 .R/: (23)
This implies that V J , the subspace composed of the set of scaling functions at level J , can approximate any function in L 2 .R/. A reduced orthonormal basis of L 2 .R/ is constructed by truncating the wavelet terms, resulting in

u.t / D X k u k J;k .t / (24)
The accuracy of this approximation increases as the level J is increased. This property of the scaling functions makes them excellent trial and weighting functions in weighted residual methods because they can be adapted to the accuracy level required. The reduced orthonormal basis of scaling functions is used in the current investigation.

Periodization of wavelet families

Standard wavelet definitions (i.e. scaling and wavelet functions) are commonly built on the real line. The functions can be adapted to periodic functions of L 2 .S 1 / by utilizing a standard periodization technique [START_REF] Mallat | A wavelet tour of signal processing[END_REF][START_REF] Erlebacher | Wavelets theory and application[END_REF][START_REF] Meyer | Ondelettes[END_REF]. Let .p/ .t / be the periodized form of the scaling function .t / defined on R

.p/ J;k .t / D X `2Z J;k .t `/ 0 Ä k Ä 2 J 1: (25)
This is equivalent to "wrapping around" the support R on S 1 through summation. The finite size of the interval results in the condition J 0

V 0 V 1 V 2 V 3 L 2 .S 1 / (26)
such that lim

J !1 V J D L 2 .S 1 /: (27)
A number of periodic discrete wavelet families exist [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. The investigation considers three families to determine how they perform in unilateral, nonsmooth contact problems: B-spline, Daubechies, and Haar.

Orthogonal cubic B-spline scaling functions

The scaling function for the orthogonal cubic B-spline wavelet family built on R is given as [START_REF] Blatter | Wavelets a primer[END_REF] .

t / D X k c k B 3 .t k/ (28)
where cubic B-spline B 3 .t / can be written using the following formulas

B 3 .t / D 8 < : t 3 =6 .0 Ä t Ä 1/ t 3 =2 C 2t 2 2t C 2=3 .1 Ä t Ä 2/ t 3 =2 4t 2 C 10t 22=3 .2 Ä t Ä 3/ t 3 =6 C 2t 2 8t C 32=3 .3 Ä t Ä 4/ 0 .otherwise/: (29)
The coefficients c k can be determined using [START_REF] Blatter | Wavelets a primer[END_REF]

c k D c k D 1 Z 0 cos.k / p p 3 .cos / d .k 0/ ( 30
)
where k is an integer. The cubic polynomial p 3 is given as

p 3 .y/ D 1 630 .y 3 C 60y 2 C 297y C 272/: (31)
Numerical simulations showed that truncating the summation in Eq. ( 28) at 50 Ä k Ä 50 is sufficient; larger k terms add negligibly to the summation. After periodization to S 1 and normalization by R 1 0 J;k .t /dt D 1 [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF], sample orthogonal cubic B-spline scaling functions for J D 0; 1; : : : ; 5 are shown in Fig. 3

FIGURE 3. EXAMPLES OF THE FIRST SIX LEVELS OF PERI-ODIZED ORTHOGONAL CUBIC B-SPLINE SCALING FUNCTIONS

Daubechies Wavelets

The Daubechies wavelet family is defined by a set of L filter coefficients fp `W `D 0; 1; : : : ; L 1g, where L is an even integer. The scaling function is defined by the fundamental twoscale equation [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF][START_REF] Daubechies | Ten lectures on wavelets[END_REF] .

t / D L 1 X `D0 p ` .2t `/ (32)
which has fundamental support over the finite intervals OE0; L 1. This equation can be used to determine the value of the scaling function at dyadic points t D n=2 J , n D 0; 1; : : : using the algorithm provided by Chen et al. [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF].

The corresponding scaling functions are highly nonsmooth and fractal in nature: as one increases the resolution of the functions, the shape does not converge but rather continues to increase in complexity. This makes accurately estimating the inner products of such scaling functions with each other prone to error when numerical integration is used [START_REF] Restrepo | Inner product computations using periodized Daubechies wavelets[END_REF]. When Daubechies scaling functions are used in a Galerkin approach, it is necessary to derive the inner products of the scaling function with itself or derivatives of itself. The exact solution to these inner products can be found by using the recursive nature of the fundamental equation on L 2 .R/ [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF][START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF]; the solution to these inner products are commonly referred to as connection coefficients. When Daubechies scaling functions are periodized on S 1 , the wrapping procedure results in functions which are no longer scale-invariant at low J values (scale-invariance requires wavelets at any scale to be a pure dilation of the mother-wavelet). Fig. 4 provides examples of the periodized scaling functions for L D 6. Interestingly, the lack of scale-invariance for small J values does not invalidate the connection coefficient algorithms derived for unbounded domains; the connection coefficients can simply be wrapped around the periodic domain as necessary [START_REF] Restrepo | Inner product computations using periodized Daubechies wavelets[END_REF].

Haar scaling functions

The simplest Daubechies wavelet family requires only two filter coefficients (p 0 D p 1 D 1) and is commonly known as the Haar wavelet family [START_REF] Haar | Zur theorie der orthogonalen funktionensysteme erste mitteilung[END_REF]. The Haar scaling functions are rectangular tophat-type functions; the father scaling function is

defined on t 2 R as .t / D ( 1 .0 Ä t < 1/ 0 .otherwise/: (33)
Since the compact support of the father scaling function is S 1 , the periodized function is equivalent. Example Haar scaling functions for J D 0; 1; : : : ; 5 are shown in Fig. 5. The values of a c and ˛in Eq. ( 5) were selected to simulate a reasonably rigid penalty function while maintaining good conditioning of the system of nonlinear equations. Table 1 lists the two forcing frequency cases considered: 150 Hz and 1275 Hz. The maximum tip displacements without contact for this system are approximately 1.4 mm at 150 Hz and 29 mm at 1275 Hz.

To act as a comparison solution, the unilateral contact finiteelement equations detailed above are solved using a variableorder numerical differentiation formula (NDF) time-stepping algorithm [START_REF] Shampine | The matlab ode suite[END_REF]. As the time-stepping solution does not solve for the periodic response directly, the solution is deemed to have converged to its periodic state when the rms of the relative error between the tip displacements for the i and i 1 periods is below 10 7 . This is possible due to the structural damping terms listed in Table 1 causing any transient response to decay relatively quickly.

Weighted residual formulation

Table 2 lists the combination of trial and weighting functions employed in the current investigation; where wavelet families are listed (B-spline, Haar, Daubechies), it is implied the scaling functions are used according to Eq. ( 24). The integrals involved in Eq. ( 17) are computed as discrete inner products using 256 points for all but the Galerkin formulations using Daubechies functions. For the Daubechies Galerkin cases, the inner products are derived using the connection coefficients as mentioned in Section 3.4.

Results

The six combinations of trial and weighting functions listed in Table 2 are employed to solve the nonlinear weighted residual formulation for forcing frequencies of 150 Hz and 1275 Hz.

Tip displacement

Samples of the approximate tip displacement responses at 150 Hz at 1275 Hz using 64 (i.e J D 6) basis functions are provided in Fig. 6 and Fig. 7, respectively.

For this number of basis functions HBM (Fourier:Fourier) approximates the tip displacements well at both 150 Hz and 1275 Hz compared to the time-stepping solution. However, close examination of the contact plateau for the 150 Hz case shows oscillations due to Gibbs phenomenon and rounding off of the sharp gradient changes. This will be further discussed in Section 5.3. Alternatively, the Haar:Fourier combination appears to approximate the tip displacement outside the contact zone less well due to the blocky nature of the Haar scaling functions. However, the contact plateau is well represented by the constant Haar scaling functions. The Haar basis appears to be less suited for the response at 1275 Hz as the contact plateau is smaller and not as flat; the tip penetration through the contact boundary is a result of the inexact modeling of hard contact using the exponential spring penalty function.

Tip velocity

Consideration is also given to the accuracy with which the trial functions considered can determine velocities. For all displacement approximations where the trial basis functions can be differentiated pointwise to approximate velocities using Eq. 11, the following equality is used:

P u i .t / D M X kD1 a i k P q k .t /: (34)
The only cases where Eq. ( 34) is not valid are those involving Haar trial functions; Haar scaling functions are piecewise constant functions thus they cannot be directly differentiated in the usual sense (i.e. pointwise) to determine velocities. Velocities for these cases can be approximated by assuming the derivative of a Haar scaling function is the combination of a positive delta function at t 1 and a negative delta function at t 2 , where t 1 corresponds to the positive jump discontinuity in the Haar scaling function and t 2 corresponds to the negative jump. Using a set of these functions as P q k .t / in Eq. ( 34) provides a reasonable approximation to the corresponding velocity field for cases involving Haar functions.

The approximate tip velocity response at 150 Hz is plotted in Fig. 8 for all cases using 64 basis functions. It is visible from these approximations that there is a sharp jump in velocity due to the contact condition. As expected, this results in ringing due to Gibbs phenomenon when the Fourier trial functions are used. Interestingly, the ringing is more pronounced in the Galerkin Bspline case. The DB6 trial functions appear to do the best job of approximating the tip response. This is attributed to the compact support of the scaling functions allowing accurate representation of rapid changes in gradient. For cases involving Haar trial functions the velocity function envelope is reasonably approximated by the delta function representation discussed above.

Tip contact force

The tip contact force is calculated using the penalty function provided in Eq. ( 5) in conjunction with the predicted tip displacement and presented in Fig. 9; the contact force has been normalized by the magnitude of the external force input. Notice the time duration of the contact forces coincide with the tip displacement plateaus of Fig. 6 as this is the only period during which the tip is in contact with the rigid foundation.

The effect of Gibbs phenomenon can be seen in the HBM and Galerkin B-spline cases. The cases where Haar trial functions are used approximate the contact force relatively well. It is hypothesized that when a rigid contact law is enforced the Haar scaling functions will perform even better relative to the other functions; this will be investigated in future work. Again, the DB6 trial functions appear to most accurately approximate the contact forces compared to the other cases.

Energy norm convergence

Convergence of the normalized energy norms for each basis combination at 150 Hz and 1275 Hz are given in Fig. 10(a) and 10(b), respectively. For these results the distributional form of the governing equations are used in all cases.

The system energy is calculated as

E D 1 2 P u > M P u C 1 2 u > Ku C E contact : (35)
The scalar energy norm is calculated as the rms value of the system energy E over a single period; this is then normalized by the energy norm from the time-stepping solution. The velocity terms P u are calculated as detailed in Section 5.2. As shown in Fig. 10, all combinations of trial and weighting functions converge to values similar to the time-stepping solution (i.e. energy norms are approximately unity using 128 basis functions). At 150 Hz, the HBM and Galerkin B-spline approaches show energy norms close to unity even using relatively few basis functions. The Petrov-Galerkin combinations utilizing Haar trial functions overpredict the system energy when few trial functions are used. This is evident in Fig. 6(e) where the blocky nature gives a coarse representation of the periodic response, thus affecting the resulting energy norm. An interesting comparison is between the DB6:DB6 and the DB6:Fourier curves. When Fourier functions are used as the weighting functions, the energy norm is better approximated. This implies projecting the residuals onto the Fourier basis, which has been shown to represent the solution relatively well, increases the accuracy of the prediction.

At 1275 Hz, the six combinations converge to a common energy norm more slowly than at the lower forcing frequency. In this case the HBM and Galerkin B-spline approaches again converge quickly to the time-stepping solution. The cases with Haar trial functions show larger difference with respect to the normalizing solution than that for the 150 Hz case. This is attributed to the displacement response at this frequency (see Fig. 7(c)); there is no distinct contact plateau thus the Haar trial functions do not approximate the solution as closely as the 150 Hz case. Finally, it is less obvious now that DB6:Fourier combination out performs the DB6:DB6 method at fewer basis functions.

Tip displacement relative error

It is also desirable to quantify the approximation error for the tip displacements predicted using the weighted residual combinations where the time-stepping solution is again used as a comparison solution. It provides a good relative measure to compare the accuracy of trial:weighting function combinations. The rms error over one period is plotted for the six basis function combinations in Fig. 11(a) for 150 Hz and Fig. 11(b) for 1275 Hz. As shown in Fig. 11(a), the HBM and Galerkin B-spline approaches have the lowest error when only 8 basis functions are used. This supports the convergence behaviour noted above for these combinations. As the number of basis functions is increased, the Galerkin Bspline case does not perform as well as the HBM method. This is attributed to a localized error which was found to develop at the end of the contact plateau (see Fig. 6(d)). The Galerkin DB6 case performs poorly using 16 trial functions, but rapidly improves to relative error values comparable to the HBM method. The ability of the nonsmooth DB6 scaling functions to accurately predict both the smooth portions of the response as well as the sharp gradient variation is attributed to both the compact support and vanishing moment properties of the the Daubechies functions. The relative error using Haar trial functions is highest out of the combinations investigated. This is expected due to the blocky nature of the Haar scaling functions; however the mid-points of the Haar functions match very well to the trend of the time-stepping solution. This suggests that if only a low resolution approximation of the response is required, the Haar scaling functions still perform well.

The relative error at 1275 Hz shows similar trends to those discussed at 150 Hz. One important difference is the Galerkin B-spline error as the number of basis functions increases: as the unknowns are increased from 32 to 128, the relative error stays reasonably constant. This suggests that this combination of basis functions has reached its limit of prediction accuracy.

Frequency response curve

Another feature of interest is the ability to construct the frequency response function for the unilateral contact condition. The penalty function simulating contact is effectively a stiffening condition, thus the resonant peaks for the system are expected to increase in frequency relative to the no-contact response. To capture this nonlinear response an arc-length continuation method [START_REF] Von Groll | The harmonic balance method with arc-length continuation in rotor/stator contact problems[END_REF] is utilized. The frequency response function over the system's fundamental frequency is shown in Fig. 12 using both the Fourier:Fourier and Haar:Fourier combinations with 64 basis functions.

As shown, the stiffening behaviour of the unilateral contact condition is captured by both trial:weighting function combinations; other combinations showed similar approximations. These results are promising as it implies the arc-length continuation method can be accurately applied to both Galerkin and Petrov-Galerkin approaches using the distributional formulation for a number of different bases.

Conclusions

The method of weighted residuals for capturing the periodic responses of unilateral contact problems has been investigated for a number of trial:weighting function combinations. The contact condition is simulated using an exponential penalty function approach. To extend the allowable trial function bases, a weak and a distributional formulations are presented which transfer trial function continuity requirements to the weighting functions. This allowed piecewise constant Haar scaling functions to be used as a trial basis in the current investigation. Results show that a number of trial:weighting function combinations produce accurate solutions which rapidly converge as the size of the discrete spaces is increased. As expected, Fourier functions perform well as a trial basis, though nonsmooth functions such as Haar and Daubechies scaling functions are also attractive since they provide comparable prediction accuracy and even out perform the Fourier functions in some measures. It is also shown that all the basis combinations considered can be used in an arc-length continuation framework to capture the nonlinear frequency response of the unilateral contact problem.

Research is continuing on this subject. The goal is to improve on the penalty approach by implementing an exact unilateral contact method (e.g Lagrange multiplier, LCP, etc.). This should result in perfectly flat contact plateaus within the periodic response. We also intend to apply this method to a standard rotor blade to validate the method for three-dimensional models involving more complicated contact interface geometries.

FIGURE 1 .

 1 FIGURE 1. ONE-DIMENSIONAL LINEAR BAR UNDERGOING UNILATERAL CONTACT CONDITIONS

FIGURE 4 .

 4 FIGURE 4. EXAMPLES OF THE FIRST SIX LEVELS OF PERI-ODIZED 6 COEFFICIENT DAUBECHIES SCALING FUNCTIONS (DB6)

17 FIGURE 5 .

 175 FIGURE 5. EXAMPLES OF THE FIRST SIX LEVELS OF PERI-ODIZED HAAR SCALING FUNCTIONS

FIGURE 6 .FIGURE 7 .

 67 FIGURE 6. TIP DISPLACEMENT VERSUS NORMALIZED TIME FOR CASE 1

FIGURE 8 .

 8 FIGURE 8. TIP VELOCITY VERSUS NORMALIZED TIME FOR CASE 1

FIGURE 9 .

 9 FIGURE 9. NORMALIZED CONTACT FORCE VERSUS NOR-MALIZED TIME FOR CASE 1

FIGURE 10 .

 10 FIGURE 10. SYSTEM ENERGY NORM CONVERGENCE USING SIX TRIAL:WEIGHTING FUNCTION COMBINATIONS

FIGURE 11 .

 11 FIGURE 11. RMS OF THE RELATIVE APPROXIMATION ER-ROR BETWEEN TRIAL:WEIGHTING FUNCTION COMBINATIONS AND THE TIME-STEPPING SOLUTION

FIGURE 12 .

 12 FIGURE 12. FREQUENCY RESPONSE FUNCTION OVER THE SECOND FUNDAMENTAL FREQUENCY OF THE UNILATERAL CONTACT SYSTEM

TABLE 1 .

 1 Test cases

		Case 1	Case 2
	forcing frequency	150 Hz	1275 Hz
	penalty constant, a c	1 10 8	1 10 8
	penalty exponent, ˛50	20
	damping factor, ˇ1 10 4 5 10 6
	initial gap, g 0	1.0 mm	14.5 mm

TABLE 2 .

 2 List of trial and weighting function combinations

	Trial function Weighting function
	Fourier	Fourier
	B-spline	B-spline
	Daubechies 6	Daubechies 6
	Haar	Fourier
	Daubechies 6	Fourier
	Haar	B-spline