# Privacy-Preserving Publication of Time-Series Data in Smart Grid

Authors: Franklin Leukam Lako; Paul Lajoie-Mazenc; Maryline Laurent
Journal: Security and Communication Networks
 (2021)

---

## Abstract

The collection of fine-grained consumptions of users in the smart grid enables energy suppliers and grid operators to propose new services (e.g., consumption forecasts and demand-response protocols) allowing to improve the efficiency and reliability of the grid. These services require the knowledge of aggregate consumption of users. However, an aggregate can be vulnerable to reidentification attacks which allow revealing the users’ individual consumption. Revealing an aggregate data is a key privacy concern. This paper focuses on publishing an aggregate of time-series data such as fine-grained consumptions, without indirectly disclosing individual consumptions. We propose novel algorithms which guarantee differential privacy, based on the discrete Fourier transform and the discrete wavelet transform. Experimental results using real data from the Irish Commission for Regulation of Utilities (CRU) demonstrate that our algorithms achieve better utility than previously proposed algorithms.


---

## Body

## 1. Introduction

A smart city is a designation given to a city that incorporates information and communication technologies to enhance the quality and performance of urban services such as energy, transportation, and utilities in order to reduce resource consumption, wastage, and overall costs. The overarching aim of a smart city is to enhance the quality of living for its citizens through smart technology [1–3].The smart grid is an important part of the smart city. Indeed, the smart grid allows greater penetration of highly variable renewable energy sources such as solar and wind power in the smart city.The smart grid modernizes the traditional electricity grid by establishing a communication infrastructure in parallel to the energy delivery network. This infrastructure is used by the grid operators and suppliers to remotely collect fine-grained consumptions from household smart meters and to provide new energy services such as consumption forecasts or demand-response. These services are suitable for improving the efficiency and reliability of the grid, saving energy and, more generally, for optimizing energy usage. In particular, forecasting enables the supplier to predict future consumptions based on past aggregate data in order to improve the grid and retail operations and enhance energy trading [4], while demand-response (DR) aims to shift the users’ consumption from peak to off-peak periods in order to avoid consumption peaks in the smart city.However, aggregates are vulnerable to reidentification attacks, such as set difference attacks [5] in which two aggregates that differ by a single consumer allow learning this individual consumption. Since the individual consumption data collected by smart meters reflect the use of all electric appliances by inhabitants in a household over time and enable to deduce the behaviors, activities, age, or preferences of the inhabitants [6–11], revealing an aggregate is a key privacy concern.Differential privacy (DP) [12] allows publishing an aggregate data while guaranteeing that an attacker does not learn any individual inputs from the aggregate. However, the noise added by DP often leads to a loss of utility. Moreover, publishing time-series data such as users’ consumption, which are correlated, by using DP, results in more noise added than publishing a single aggregate for the same privacy guarantee. Thus, disclosing time-series data leads to more loss of utility. Utility can be improved by increasing the size of the aggregate. Eibl and Engel [13] showed that for real-world smart metering, the aggregation group size must be of the order of thousands of smart meters in order to have reasonable utility. This paper shows how to obtain good utility with a group size smaller than 600. We obtain a mean relative error lower than 10% between the original data and the published one, which is considered practically suitable by energy experts for consumption forecasts.The Laplace mechanism [14] is a popular mechanism to enable DP, by adding independent and identically distributed (IID) Laplace noise to each component of the time-series. However, adding IID noise for correlated time-series is not appropriate. In fact, an adversary can use refinement methods, such as filtering, to sanitize the IID noise and improve the probability of disclosing individual data [15, 16].This paper focuses on disclosing an aggregate of users’ consumption data without learning individual data and proposes methods with improved utility. We summarize our contributions as follows:(i)
We revisit the Fourier perturbation algorithm (FPA) [17] in order to correct some mistakes leading to poor users’ privacy protection. We show that, in order to ensure the desired budget of privacy ϵ, a factor 2T must be added to the noise, where T is the size of the time-series. However, this reduces the utility of FPA.(ii)
We propose the “clamping Fourier perturbation algorithm (CFPA)” using the clamping mechanism proposed in [18], for reducing the sensitivity, and thus the noise introduced in FPA. This new algorithm is an improvement of the Fourier perturbation algorithm (FPA). Experimental results show a utility improvement by a factor more than 6.(iii)
We also propose the “clamping wavelet perturbation algorithm” (CWPA), a similar adaptation of wavelet perturbation algorithm (WPA) [19], with a utility improvement by a factor 2.(iv)
We compare FPA, CFPA, WPA, and CWPA by analyzing their relative errors on a real dataset, and we explain why CFPA obtains the best utility.The remainder of this paper is structured as follows. Section2 provides an overview of the literature, while Section 3 presents preliminaries. Section 4 correctly computes the sensitivity of DFT in order to make FPA ϵ-differentially private. Section 5 details our privacy-preserving publication techniques using clamping mechanism, DFT, and DWT. Section 6 reports our experimental results. Section 7 concludes the paper.Table1 lists the acronyms used in this paper.Table 1
List of acronyms.
AcronymMeaningCFPAClamping Fourier perturbation algorithmCWPAClamping wavelet perturbation algorithmDFTDiscrete Fourier transformDPDifferential privacyDWTDiscrete wavelet transformFPAFourier perturbation algorithmMREMean relative estimation errorWPAWavelet perturbation algorithm


## 2. Related Work

Demand-response protocols [20–23], and secure aggregation protocols [24–30] aim to protect the privacy of users while supporting energy services such as demand-response, smart metering, billing, or forecasting.In this paper, we investigate tools enabling forecasting and demand-response. In particular, we are interested in publishing an aggregate of individual consumptions, while preserving privacy.Differential privacy (DP), introduced by Dwork in 2006, guarantees that the publication of an aggregate does not indirectly reveal the individual data [12]. Moreover, DP guarantees that two aggregates that differ by a single consumer are almost indistinguishable. DP has evolved overtime [31] and was adopted by organizations such as the US Census Bureau [32], Google [33], Apple [34], and Microsoft [35]. The Laplace mechanism [14] is a popular mechanism that allows guaranteeing DP by adding a noise drawn from the Laplace distribution ℒ⋅ to the aggregate.The Laplace mechanism takes as input two parameters: the privacy budgetϵ and the sensitivity of the function to publish (in our case, the sum of users’ consumption). Smaller values of ϵ lead to a better protection, but add a bigger noise to the aggregate.Utility can be improved by increasing the size of the aggregate in order that the effect of noise is small enough that the result can be utilized. Eibl and Engel [13] showed that for real-world smart metering, the aggregation group size must be of the order of thousand smart meters in order to have reasonable utility. This paper shows how to obtain good utility with a group size smaller than a thousand.DP is typically applied to static data, i.e., to a single query. In this paper, we consider time-series consumption, which is equivalent to multiple queries on correlated data. Applying the Laplace mechanism independently to each data point of the time-series is not appropriate. Indeed, an adversary can use refinement methods, such as filtering, to sanitize the Laplace noise and improve the probability disclosing individual data [15, 16]. Thus, the data points of the time-series are correlated. The composition theorem [14] states that the privacy budget ϵ of T correlated queries adds up, i.e., setting the privacy budget for a single query to εq=0.5, the privacy budget of T=48 single queries (corresponding to a day profile with a time interval of 30 min) is ε=0.5×48=24. In order to guarantee a global privacy budget of ε, one solution is to set the privacy budget of each query to ε/T. Of course, this leads to more noise in the aggregate and a loss of utility.One method to guarantee DP for correlated time-series data publishing consists in transforming the original correlated time-series into another representation while maintaining its major characteristics before adding the Laplace noise. Rastogi and Nath [17] proposed the Fourier perturbation algorithm (FPA) that combines discrete Fourier transform (DFT) with DP to support time-series of count queries while not disclosing any individual data and ensuring good utility. We note that the sensitivity of count queries is 1, and the global sensitivity is T for a time-series of length T. Ács et al. [36] proposed an optimization of the FPA allowing to release histograms, where the global sensitivity is 1. They show through experimental evaluation that their scheme improves the utility of the initial FPA by a factor 10. Lyu et al. [19] applied FPA to time-series consumptions and proposed wavelet perturbation algorithm (WPA) by replacing DFT by discrete wavelet transform (DWT). The authors show through experimental results that WPA ensures better utility than FPA.We apply these approaches to time-series of consumption data and refine them by reducing the sensitivity of the queries in order to reduce the relative error of the final result.


## 3. Preliminaries

### 3.1. System and Threat Model

The entities involved in this paper are as follows:(i)
Trustworthy homes, which smart meter (SM) enables to collect their true individual time-series consumption.(ii)
A honest aggregator, which collects users’ individual consumption, and which publishes an aggregate time-series consumption of users to a forecaster in a privacy-preserving way for the forecaster not to be able to deduce any individual consumption of users.(iii)
A forecaster, which predicts future consumptions based on the aggregate consumption received in order to improve the grid and retail operations and enhance energy trading. The forecaster is considered honest-but-curious as it provides appropriate forecasts, but it may attempt to infer the users’ individual consumption from the aggregate in order to deduce the behaviors, activities, age, or preferences of the inhabitants.Figure1 depicts the system model. In a real scenario, the aggregator can be an energy distributor, and the forecaster can be a municipality that seeks to find out the total consumption of the inhabitants of the municipality.Figure 1
System model.Considering the case where the aggregator and the forecaster belong to two entities of the same energy provider, the publication of aggregate users’ consumption to forecasters in a privacy-preserving way reduces the risk of disclosing users’ individual consumption. Moreover, this avoids the need for forecasters to ask for explicit consent from customers in accordance with the GDPR [37] to process their personal data.LetN be the number of smart meters (SMs) in a district. Let Xj=x1j,x2j,…,xTj be the time-series of energy consumptions collected by SM j, where xtj is the consumption at time slot t (t=1,…,T) collected by SM j (j=1,…,N), with T being the time period considered. Each time-series consumption Xj is sent to an aggregator who computes the following aggregate:(1)S=S1,…,ST=∑j=1Nx1j,∑j=1Nx2j,…,∑j=1NxTj.To revealS to a forecaster without indirectly disclosing individual consumptions Xj (j=1,…,N), the aggregator can use differential privacy (DP).


### 3.2. Differential Privacy

Differential privacy is a framework introduced by Dwork allowing quantifying the privacy guarantees of a request on a database [38]. This request can be the publication of a database, or a more precise one, such as “what is the sum of energy consumptions of users in this database?”.A request on databases is said to be differentially private if this request makes two similar databases indistinguishable from looking only at the output of the request. Differential privacy relies on a parameter, notedε, called the privacy budget. The formal definition of a differentially private algorithm is given as follows.Definition 1 (ε−differentially private).
A requestA:D⟶S is ε−differentially private if and only if for all databases D1,D2∈D differing by at most one record, and for all subsets O⊂S,(2)PrAD1∈O≤expεPrAD2∈O.This definition can be applied not only to requests on databases but also to any function, by considering the domain of the function as a database format.Dwork also proposes the Laplace mechanism, which allows making any (vectors of) real-valued functionε-differentially private [38]. This mechanism relies on the notion of sensitivity of a function, which represents how a single record of the database can influence the output of the function.Definition 2 (sensitivity).
Letf:D⟶ℝd be a function; the sensitivity of f is(3)Δ1f=maxD1,D2∈Ds.t.dD1,D2≤1fD1−fD21.This sensitivity is also calledL1-sensitivity due to the L1−norm used in its definition and is denoted by Δ1f. Similarly, the L2-sensitivity used later and denoted by ε is computed using the L2−norm (the L1−norm and the L2 norm of a vector S=s1,…,sT are respectively equal to S1=∑j=1Tsj and S2=∑j=1Tsj2).The Laplace mechanism consists of adding a random value to the original result of the query, where the random value follows the Laplace distribution, where the parameter depends on the chosenϵ and on the sensitivity of the function, as follows.Theorem 1 (Laplacian mechanism).
For all functionsf:D⟶ℝd, the algorithm AD=fD+ℒ1Δ1f/ε,…,ℒdΔ1f/ε is ϵ−differentially private, where ℒ⋅ is the distribution of Laplace and Δ1f is the sensitivity of f.DP introduces noise in order to guarantee privacy. This noise can decrease the utility of the function. We quantify this loss using mean relative estimation error (MRE), defined as follows.Definition 3 (mean relative estimation error).
The mean relative estimation error (MRE) between two vectorsa and b of size T is 1/T⋅∑j=1Taj−bj/aj+12 (we add 1 to the denominator in order to avoid dividing by zero. This definition is also used in [27]).Consider the aggregateS=S1,…,ST defined in (1). Let M be the maximum consumption in the domain. One naive solution to publish S without revealing any individual consumption is to use the Laplace mechanism to add independent Laplace noise to each component of S and to release the results: S^=S1+ℒM⋅T/ε,…,ST+ℒM⋅T/ε, where the sensitivity of the sum of time-series consumption is M⋅T. However, this simple approach leads to excessive noise rendering the aggregate useless [13].Example 1.
Figures2(a) and 2(b), respectively, show the aggregated consumption of 250 homes from December 30th, 2009, to January 5th, 2010, taken from the CER dataset [39], and its noisy version using the naively applied Laplace mechanism, with ε=1 per day. Figure 2(a) shows two consumption peaks at 12 am and 6 pm which respectively correspond to lunch and dinner time. We also observe that in the night (from 12 pm to 6 am) the consumption decreases. Figure 2(b) shows that the noisy version is completely different from the original aggregate (Figure 2(a)). In this example, the MRE between the aggregate consumption and the noisy version is 141%, which is not usable. Moreover, the noisy version has inconsistent values such as negative consumptions.Figure 2
Aggregated time-series consumption of 250 homes from December 30th, 2009, to January 5th, 2010, taken from the CER dataset [39], and its noisy version using the naively applied Laplace mechanism, with ε=1 per day. (a) Aggregated consumption. (b) Noisy version using the naive solution.
(a)(b)Rastogi and Nath [17] introduce the Fourier perturbation algorithm (FPA) and show that is an effective tool for reducing the noise introduced by the Laplace mechanism for time-series. Section 3.3 presents the FPA. However, there are some mistakes in this version relying on the estimation of the FPA sensitivity. These mistakes are presented in Section 4, along with the corrected FPA.Table2 lists the symbols used in the rest of the paper.Table 2
List of symbols.
NotationDescriptionNNumber of smart meters (SMs) in the districtMMaximum consumption in the datasetTTime period during the collection of time-series consumptionXj=x1j,x2j,…,xTjTime-series of energy consumptions collected by SMj, where xtj is the consumption at time slot t (t=1,…,T) collected by SM j (j=1,…,N)S=S1,…,STSum of users’ time-series consumptions to be published, whereSt=∑j=1Nxtj for t=1,…,TεBudget of privacyS^Noisy version ofSkNumber of the first DFT or DWT coefficients conserved in the Fourier perturbation algorithm (FPA), wavelet perturbation algorithm (WPA), clamping Fourier perturbation algorithm (CFPA), and clamping wavelet perturbation algorithm (CWPA)Δ1L1-sensitivityΔ2L2-sensitivity


### 3.3. Fourier Perturbation Algorithm

The Fourier perturbation algorithm (FPA) presented in [17, 19, 36] takes as input a time-series S=S1,…,ST and an integer k≪T and returns the noisy time-series S^=S^1,…,S^T, as shown in Algorithm 1.Algorithm 1: Fourier perbutation algorithm [17].
Inputs:S=S1,…,ST, k, the maximum consumption M of the domain, and the privacy budget ϵ.(1)
Compute the discrete Fourier transform ofS: F=DFTS.(2)
Keep only the firstk coefficients of F, denoted by Fk.(3)
Generate the noisy version ofFk, denoted by F^kby adding a Laplace noise ℒMTk/ε to each coefficient in Fk.(4)
PadF^k to a T-dimensional vector, denoted by PADTF^k by appending T−k zeroes.(5)
Apply the inverse DFT toPADTF^k to obtain a noisy version of S denoted by S^.Rastogi and Nath [17] show that FPA is ε−differentially private. However, there are some mistakes in their proof of Theorem 4.1 of [17] which justified that FPA is ε-differentially private. These mistakes rely on the estimation of the FPA sensitivity and are presented in Section 4.


### 3.4. Wavelet Perturbation Algorithm

By replacing the DFT with the discrete Haar wavelet transform (DWT), Lyu et al. [19] proposed the wavelet perturbation algorithm (WPA) and showed that WPA guarantees better utility than DFT. Algorithm 2 describes WPA.Figure3 shows the same aggregated consumption presented in Example 1 and its noisy version using WPA (Algorithm 2) with Haar wavelet, ε=1 per day and k=5. In Figure 3, the MRE is however higher than 10% (18%). In the noisy aggregate, the first peak of the morning is masked and the peak of the evening is truncated, as well as the trough of the night.Algorithm 2: Wavelet perbutation algorithm [19].
Inputs:S=S1,…,ST, k, the maximum consumption M of the domain, and the privacy budget: ε(1)
Compute the DWT coefficients ofS : W=DWTS.(2)
Keep only the firstk coefficients of W, denoted by Wk.(3)
Generate the noisy version ofWk, denoted by W^k by adding a Laplace noise ℒMTk/ε to each coefficient in Wk.(4)
PadW^k to a T-dimensional vector, denoted by PADTW^k by appending T−k zeroes.(5)
Apply the inverse DWT toPADTW^k to obtain a noisy version of S denoted by S^.Figure 3
Aggregated consumption of 250 homes from 30th December 2009 to 5th January 2010 of dataset from CER [39] and its noisy version using WPA (Algorithm 2), with Haar wavelet, ε=1 per day and k=5.Theorem 2.
Wavelet perturbation algorithm (WPA) isε-differentially private.Proof.
DWT is orthonormal [40], i.e., W has the same L2 norm as S, that is, Δ2W=Δ2S. Furthermore, Δ2Wk≤Δ2S (because T−k DWT coefficients of W are set to 0). With the inequality of norm, Δ1Wk≤kΔ2Wk. Then, Δ1Wk≤kΔ2S≤M⋅kT. Thus, the noise introduced in Step 3 is justified and WPA guarantees ε−differential privacy.


## 3.1. System and Threat Model

The entities involved in this paper are as follows:(i)
Trustworthy homes, which smart meter (SM) enables to collect their true individual time-series consumption.(ii)
A honest aggregator, which collects users’ individual consumption, and which publishes an aggregate time-series consumption of users to a forecaster in a privacy-preserving way for the forecaster not to be able to deduce any individual consumption of users.(iii)
A forecaster, which predicts future consumptions based on the aggregate consumption received in order to improve the grid and retail operations and enhance energy trading. The forecaster is considered honest-but-curious as it provides appropriate forecasts, but it may attempt to infer the users’ individual consumption from the aggregate in order to deduce the behaviors, activities, age, or preferences of the inhabitants.Figure1 depicts the system model. In a real scenario, the aggregator can be an energy distributor, and the forecaster can be a municipality that seeks to find out the total consumption of the inhabitants of the municipality.Figure 1
System model.Considering the case where the aggregator and the forecaster belong to two entities of the same energy provider, the publication of aggregate users’ consumption to forecasters in a privacy-preserving way reduces the risk of disclosing users’ individual consumption. Moreover, this avoids the need for forecasters to ask for explicit consent from customers in accordance with the GDPR [37] to process their personal data.LetN be the number of smart meters (SMs) in a district. Let Xj=x1j,x2j,…,xTj be the time-series of energy consumptions collected by SM j, where xtj is the consumption at time slot t (t=1,…,T) collected by SM j (j=1,…,N), with T being the time period considered. Each time-series consumption Xj is sent to an aggregator who computes the following aggregate:(1)S=S1,…,ST=∑j=1Nx1j,∑j=1Nx2j,…,∑j=1NxTj.To revealS to a forecaster without indirectly disclosing individual consumptions Xj (j=1,…,N), the aggregator can use differential privacy (DP).


## 3.2. Differential Privacy

Differential privacy is a framework introduced by Dwork allowing quantifying the privacy guarantees of a request on a database [38]. This request can be the publication of a database, or a more precise one, such as “what is the sum of energy consumptions of users in this database?”.A request on databases is said to be differentially private if this request makes two similar databases indistinguishable from looking only at the output of the request. Differential privacy relies on a parameter, notedε, called the privacy budget. The formal definition of a differentially private algorithm is given as follows.Definition 1 (ε−differentially private).
A requestA:D⟶S is ε−differentially private if and only if for all databases D1,D2∈D differing by at most one record, and for all subsets O⊂S,(2)PrAD1∈O≤expεPrAD2∈O.This definition can be applied not only to requests on databases but also to any function, by considering the domain of the function as a database format.Dwork also proposes the Laplace mechanism, which allows making any (vectors of) real-valued functionε-differentially private [38]. This mechanism relies on the notion of sensitivity of a function, which represents how a single record of the database can influence the output of the function.Definition 2 (sensitivity).
Letf:D⟶ℝd be a function; the sensitivity of f is(3)Δ1f=maxD1,D2∈Ds.t.dD1,D2≤1fD1−fD21.This sensitivity is also calledL1-sensitivity due to the L1−norm used in its definition and is denoted by Δ1f. Similarly, the L2-sensitivity used later and denoted by ε is computed using the L2−norm (the L1−norm and the L2 norm of a vector S=s1,…,sT are respectively equal to S1=∑j=1Tsj and S2=∑j=1Tsj2).The Laplace mechanism consists of adding a random value to the original result of the query, where the random value follows the Laplace distribution, where the parameter depends on the chosenϵ and on the sensitivity of the function, as follows.Theorem 1 (Laplacian mechanism).
For all functionsf:D⟶ℝd, the algorithm AD=fD+ℒ1Δ1f/ε,…,ℒdΔ1f/ε is ϵ−differentially private, where ℒ⋅ is the distribution of Laplace and Δ1f is the sensitivity of f.DP introduces noise in order to guarantee privacy. This noise can decrease the utility of the function. We quantify this loss using mean relative estimation error (MRE), defined as follows.Definition 3 (mean relative estimation error).
The mean relative estimation error (MRE) between two vectorsa and b of size T is 1/T⋅∑j=1Taj−bj/aj+12 (we add 1 to the denominator in order to avoid dividing by zero. This definition is also used in [27]).Consider the aggregateS=S1,…,ST defined in (1). Let M be the maximum consumption in the domain. One naive solution to publish S without revealing any individual consumption is to use the Laplace mechanism to add independent Laplace noise to each component of S and to release the results: S^=S1+ℒM⋅T/ε,…,ST+ℒM⋅T/ε, where the sensitivity of the sum of time-series consumption is M⋅T. However, this simple approach leads to excessive noise rendering the aggregate useless [13].Example 1.
Figures2(a) and 2(b), respectively, show the aggregated consumption of 250 homes from December 30th, 2009, to January 5th, 2010, taken from the CER dataset [39], and its noisy version using the naively applied Laplace mechanism, with ε=1 per day. Figure 2(a) shows two consumption peaks at 12 am and 6 pm which respectively correspond to lunch and dinner time. We also observe that in the night (from 12 pm to 6 am) the consumption decreases. Figure 2(b) shows that the noisy version is completely different from the original aggregate (Figure 2(a)). In this example, the MRE between the aggregate consumption and the noisy version is 141%, which is not usable. Moreover, the noisy version has inconsistent values such as negative consumptions.Figure 2
Aggregated time-series consumption of 250 homes from December 30th, 2009, to January 5th, 2010, taken from the CER dataset [39], and its noisy version using the naively applied Laplace mechanism, with ε=1 per day. (a) Aggregated consumption. (b) Noisy version using the naive solution.
(a)(b)Rastogi and Nath [17] introduce the Fourier perturbation algorithm (FPA) and show that is an effective tool for reducing the noise introduced by the Laplace mechanism for time-series. Section 3.3 presents the FPA. However, there are some mistakes in this version relying on the estimation of the FPA sensitivity. These mistakes are presented in Section 4, along with the corrected FPA.Table2 lists the symbols used in the rest of the paper.Table 2
List of symbols.
NotationDescriptionNNumber of smart meters (SMs) in the districtMMaximum consumption in the datasetTTime period during the collection of time-series consumptionXj=x1j,x2j,…,xTjTime-series of energy consumptions collected by SMj, where xtj is the consumption at time slot t (t=1,…,T) collected by SM j (j=1,…,N)S=S1,…,STSum of users’ time-series consumptions to be published, whereSt=∑j=1Nxtj for t=1,…,TεBudget of privacyS^Noisy version ofSkNumber of the first DFT or DWT coefficients conserved in the Fourier perturbation algorithm (FPA), wavelet perturbation algorithm (WPA), clamping Fourier perturbation algorithm (CFPA), and clamping wavelet perturbation algorithm (CWPA)Δ1L1-sensitivityΔ2L2-sensitivity


## 3.3. Fourier Perturbation Algorithm

The Fourier perturbation algorithm (FPA) presented in [17, 19, 36] takes as input a time-series S=S1,…,ST and an integer k≪T and returns the noisy time-series S^=S^1,…,S^T, as shown in Algorithm 1.Algorithm 1: Fourier perbutation algorithm [17].
Inputs:S=S1,…,ST, k, the maximum consumption M of the domain, and the privacy budget ϵ.(1)
Compute the discrete Fourier transform ofS: F=DFTS.(2)
Keep only the firstk coefficients of F, denoted by Fk.(3)
Generate the noisy version ofFk, denoted by F^kby adding a Laplace noise ℒMTk/ε to each coefficient in Fk.(4)
PadF^k to a T-dimensional vector, denoted by PADTF^k by appending T−k zeroes.(5)
Apply the inverse DFT toPADTF^k to obtain a noisy version of S denoted by S^.Rastogi and Nath [17] show that FPA is ε−differentially private. However, there are some mistakes in their proof of Theorem 4.1 of [17] which justified that FPA is ε-differentially private. These mistakes rely on the estimation of the FPA sensitivity and are presented in Section 4.


## 3.4. Wavelet Perturbation Algorithm

By replacing the DFT with the discrete Haar wavelet transform (DWT), Lyu et al. [19] proposed the wavelet perturbation algorithm (WPA) and showed that WPA guarantees better utility than DFT. Algorithm 2 describes WPA.Figure3 shows the same aggregated consumption presented in Example 1 and its noisy version using WPA (Algorithm 2) with Haar wavelet, ε=1 per day and k=5. In Figure 3, the MRE is however higher than 10% (18%). In the noisy aggregate, the first peak of the morning is masked and the peak of the evening is truncated, as well as the trough of the night.Algorithm 2: Wavelet perbutation algorithm [19].
Inputs:S=S1,…,ST, k, the maximum consumption M of the domain, and the privacy budget: ε(1)
Compute the DWT coefficients ofS : W=DWTS.(2)
Keep only the firstk coefficients of W, denoted by Wk.(3)
Generate the noisy version ofWk, denoted by W^k by adding a Laplace noise ℒMTk/ε to each coefficient in Wk.(4)
PadW^k to a T-dimensional vector, denoted by PADTW^k by appending T−k zeroes.(5)
Apply the inverse DWT toPADTW^k to obtain a noisy version of S denoted by S^.Figure 3
Aggregated consumption of 250 homes from 30th December 2009 to 5th January 2010 of dataset from CER [39] and its noisy version using WPA (Algorithm 2), with Haar wavelet, ε=1 per day and k=5.Theorem 2.
Wavelet perturbation algorithm (WPA) isε-differentially private.Proof.
DWT is orthonormal [40], i.e., W has the same L2 norm as S, that is, Δ2W=Δ2S. Furthermore, Δ2Wk≤Δ2S (because T−k DWT coefficients of W are set to 0). With the inequality of norm, Δ1Wk≤kΔ2Wk. Then, Δ1Wk≤kΔ2S≤M⋅kT. Thus, the noise introduced in Step 3 is justified and WPA guarantees ε−differential privacy.


## 4. Correctly Estimating the Sensitivity of FPA

In [17], authors show that FPA, as described in Section 3, guarantees ε−differential privacy. The authors estimated the sensitivity of DFT to be MTk, while it should be MT2k, with T being the size of the time-series and M being the maximum consumption in the domain. Thus, for a given privacy budget ε, the utility of FPA is worse than presented in [17].This section correctly computes the sensitivity of DFT, which allows to make render FPAε−differential private. Before that we recall the definition of DFT.


### 4.1. Discrete Fourier Transform (DFT)

LetS=S1,…,ST be a time-series. DFT takes S as input and returns a time-series of T complex numbers F=F1,…,FT such that(4)Fk=1T∑j=1TSje−2πij−1k−1/Tfork=1,…,T,where i2=−1. The inverse of the DFT is computed as follows:(5)Sk=1T∑j=1TFje−2πij−1k−1/T,fork=1,…,T.This version of the DFT is normalized, that is,DFTS2=S2.DFT can be defined in other ways, for instance, the1/T, present in both the DFT and the inverse definitions above, can be replaced by a factor 1 in the DFT and 1/T in the inverse DFT. In that case, the DFT is not normalized.In [17, 19, 36], the authors use the latter version of DFT, which is not normalized. However, the sensitivity computation relies on the equality DFTS2=S2, while it should be DFTS2=T⋅S2. Thus, the correct total privacy budget is T⋅ε instead of ε. This is the first mistake in this approach and can be resolved by using the normalized DFT.Another error lies in the fact that the Laplacian mechanism is only applied to the real part of the Fourier coefficients, which are complex numbers. This mistake can be resolved by applying the Laplace mechanism to both real and imaginary parts of the Fourier coefficients.The following section computes the sensitivity of the DFT, and thus of FPA, and takes into account those two errors.


### 4.2. Sensitivity of the DFT

LetDFTk be the function which takes a time-series S=S1,…,ST as input and returns the first k DFT coefficients of S. This function can be seen as a DFTk:ℝT⟶ℝ2k, the function which returns the real and imaginary parts of the first k Fourier coefficients. This function is a real-valued function, we can thus use the Laplace mechanism on it. First, we need to compute the L1-sensitivity of DFTk.Lemma 1.
LetDFTk be defined as follows:(6)DFTk:ℝT⟶ℝ2kS⟼DFTkS=a1,b1,…,ak,bk.We denotecj=aj+ibj the j-th coefficient of DFTS, with i2=−1 and j=1,…,k. aj and bj respectively represent the real and imaginary parts of cj.TheL1-sensitivity of DFTk, DFTkS1, is M⋅2Tk when the DFT is normalized (respectively, MT⋅2k when the DFT is not normalized), with M as the maximum value in the dataset.Proof.
LetDFTk be defined as in Lemma 1.(7)DFTkS1=a1,b1,…,ak,bk1=∑j=1kaj,bj1≤2∑j=1kaj,bj2Minkowski inequality≤2∑j=1kcj≤2c1,…,ck1≤2kc1,…,ck2Minkowski inequality≤2ks1,…,sT2asTFourier coefficients have the sameL2norm asS.Then,Δ1DFTk≤2kΔ2S=M2Tk.This result is true when the DFT is normalized (2) as in our case. In [17, 19, 36], the L2 norm of Fourier coefficients equals to T times the L2 norm of S (Parvesal’s theorem). This result is valid when the normalized DFT (2) is used as in our case. When the DFT is not normalized, as is the case in [17, 19, 36], the sensitivity of the first k DFT coefficients should be Δ1DFTk=T×M2Tk=MT2k instead of (MTk). Thus, using the normalized DFT, the function then becomes(8)DFTk˜:ℝT⟶ℝ2kS⟼DFTS˜=a1,b1,…,ak,bk+y1,1,y1,2,…,yk,1,yk,2,which is ε−DP, with yj,ℓ=ℒM2Tk/ε, for all j=1,…,k and ℓ=1,2.For simplicity, in the following, we writecj+ℒM2Tk/ε instead of aj,bj+ℒM2Tk/ε,ℒM2Tk/ε, meaning that two independent Laplace noises ℒM2Tk/ε are added to the real and imaginary parts of cj.Algorithm3 shows the Fourier perturbation algorithm (FPA) revisited.Algorithm 3: Fourier perturbation algorithm (FPA) revisited.
Inputs: S=S1,…,ST, k, the maximum consumption M of the domain, and the privacy budget ε.(1)
Compute the normalized DFT coefficients ofS : F=DFTS.(2)
Keep only the firstk coefficients of F, denoted by Fk.(3)
Generate the noisy version ofFk, denoted by F^k by adding a Laplace noise ℒM2Tk/ε to each coefficient in Fk.(4)
PadF^k to a T-dimensional vector, denoted by PADTF^k by appending T−k zeroes.(5)
Apply the inverse DFT toPADTF^k to obtain a noisy version of S denoted by S^.


### 4.3. Differences between the Initial, yet Incorrect, FPA, and the Corrected FPA

For a budget of privacyε, the differences between the initial incorrect FPA and the corrected one can be highlighted as follows:(1)
The DFT used in the initial incorrect FPA [17] is not normalized, while it is normalized in the corrected FPA. Thus, a factor 2T is missing in the Laplace noise in Algorithm 1.(2)
In the initial incorrect FPA [17], Laplace noises are only added to the real part of the DFT coefficients, while they should be added to the real and imaginary parts of the DFT coefficients as in the corrected FPA (Algorithm 3). Thus, k imaginary coefficients are not noised in Algorithm 1.Figure4 shows the same aggregated consumption presented in Example 1 and its noisy version using the corrected FPA (Algorithm 3) with ε=1 per day and k=5. Figure 4 shows that the corrected FPA obtains a large MRE (84%), making it useless. The noisy aggregate has negative consumptions and does not contain the peaks present in the original aggregate.Figure 4
Aggregated consumption of 250 homes from 30th December 2009 to 5th January 2010 of dataset from CER [39] and its noisy version using the corrected FPA (Algorithm 3), with ε=1 per day and k=5.For the sake of simplicity, in the following sections, we use FPA to talk about the corrected version.


## 4.1. Discrete Fourier Transform (DFT)

LetS=S1,…,ST be a time-series. DFT takes S as input and returns a time-series of T complex numbers F=F1,…,FT such that(4)Fk=1T∑j=1TSje−2πij−1k−1/Tfork=1,…,T,where i2=−1. The inverse of the DFT is computed as follows:(5)Sk=1T∑j=1TFje−2πij−1k−1/T,fork=1,…,T.This version of the DFT is normalized, that is,DFTS2=S2.DFT can be defined in other ways, for instance, the1/T, present in both the DFT and the inverse definitions above, can be replaced by a factor 1 in the DFT and 1/T in the inverse DFT. In that case, the DFT is not normalized.In [17, 19, 36], the authors use the latter version of DFT, which is not normalized. However, the sensitivity computation relies on the equality DFTS2=S2, while it should be DFTS2=T⋅S2. Thus, the correct total privacy budget is T⋅ε instead of ε. This is the first mistake in this approach and can be resolved by using the normalized DFT.Another error lies in the fact that the Laplacian mechanism is only applied to the real part of the Fourier coefficients, which are complex numbers. This mistake can be resolved by applying the Laplace mechanism to both real and imaginary parts of the Fourier coefficients.The following section computes the sensitivity of the DFT, and thus of FPA, and takes into account those two errors.


## 4.2. Sensitivity of the DFT

LetDFTk be the function which takes a time-series S=S1,…,ST as input and returns the first k DFT coefficients of S. This function can be seen as a DFTk:ℝT⟶ℝ2k, the function which returns the real and imaginary parts of the first k Fourier coefficients. This function is a real-valued function, we can thus use the Laplace mechanism on it. First, we need to compute the L1-sensitivity of DFTk.Lemma 1.
LetDFTk be defined as follows:(6)DFTk:ℝT⟶ℝ2kS⟼DFTkS=a1,b1,…,ak,bk.We denotecj=aj+ibj the j-th coefficient of DFTS, with i2=−1 and j=1,…,k. aj and bj respectively represent the real and imaginary parts of cj.TheL1-sensitivity of DFTk, DFTkS1, is M⋅2Tk when the DFT is normalized (respectively, MT⋅2k when the DFT is not normalized), with M as the maximum value in the dataset.Proof.
LetDFTk be defined as in Lemma 1.(7)DFTkS1=a1,b1,…,ak,bk1=∑j=1kaj,bj1≤2∑j=1kaj,bj2Minkowski inequality≤2∑j=1kcj≤2c1,…,ck1≤2kc1,…,ck2Minkowski inequality≤2ks1,…,sT2asTFourier coefficients have the sameL2norm asS.Then,Δ1DFTk≤2kΔ2S=M2Tk.This result is true when the DFT is normalized (2) as in our case. In [17, 19, 36], the L2 norm of Fourier coefficients equals to T times the L2 norm of S (Parvesal’s theorem). This result is valid when the normalized DFT (2) is used as in our case. When the DFT is not normalized, as is the case in [17, 19, 36], the sensitivity of the first k DFT coefficients should be Δ1DFTk=T×M2Tk=MT2k instead of (MTk). Thus, using the normalized DFT, the function then becomes(8)DFTk˜:ℝT⟶ℝ2kS⟼DFTS˜=a1,b1,…,ak,bk+y1,1,y1,2,…,yk,1,yk,2,which is ε−DP, with yj,ℓ=ℒM2Tk/ε, for all j=1,…,k and ℓ=1,2.For simplicity, in the following, we writecj+ℒM2Tk/ε instead of aj,bj+ℒM2Tk/ε,ℒM2Tk/ε, meaning that two independent Laplace noises ℒM2Tk/ε are added to the real and imaginary parts of cj.Algorithm3 shows the Fourier perturbation algorithm (FPA) revisited.Algorithm 3: Fourier perturbation algorithm (FPA) revisited.
Inputs: S=S1,…,ST, k, the maximum consumption M of the domain, and the privacy budget ε.(1)
Compute the normalized DFT coefficients ofS : F=DFTS.(2)
Keep only the firstk coefficients of F, denoted by Fk.(3)
Generate the noisy version ofFk, denoted by F^k by adding a Laplace noise ℒM2Tk/ε to each coefficient in Fk.(4)
PadF^k to a T-dimensional vector, denoted by PADTF^k by appending T−k zeroes.(5)
Apply the inverse DFT toPADTF^k to obtain a noisy version of S denoted by S^.


## 4.3. Differences between the Initial, yet Incorrect, FPA, and the Corrected FPA

For a budget of privacyε, the differences between the initial incorrect FPA and the corrected one can be highlighted as follows:(1)
The DFT used in the initial incorrect FPA [17] is not normalized, while it is normalized in the corrected FPA. Thus, a factor 2T is missing in the Laplace noise in Algorithm 1.(2)
In the initial incorrect FPA [17], Laplace noises are only added to the real part of the DFT coefficients, while they should be added to the real and imaginary parts of the DFT coefficients as in the corrected FPA (Algorithm 3). Thus, k imaginary coefficients are not noised in Algorithm 1.Figure4 shows the same aggregated consumption presented in Example 1 and its noisy version using the corrected FPA (Algorithm 3) with ε=1 per day and k=5. Figure 4 shows that the corrected FPA obtains a large MRE (84%), making it useless. The noisy aggregate has negative consumptions and does not contain the peaks present in the original aggregate.Figure 4
Aggregated consumption of 250 homes from 30th December 2009 to 5th January 2010 of dataset from CER [39] and its noisy version using the corrected FPA (Algorithm 3), with ε=1 per day and k=5.For the sake of simplicity, in the following sections, we use FPA to talk about the corrected version.


## 5. Clamping Transform Perturbation Algorithm

The intuition behind our approach, “Clamping transform perturbation algorithm,” lies in the perturbation error, caused by the Laplace mechanism, which depends on the sensitivity of the sum of consumptions. As such, by reducing the sensitivity, we expect to reduce the perturbation error.To estimate the sensitivity of consumptions, we split our database ofN users into two almost equal parts: D1 corresponding to the consumptions of the first half of users (a training dataset) and D2 containing the second half of users’ consumptions (a validation dataset). Using D1, we compute the distribution of users’ consumptions in the frequency domain. We denote by M=M1,…,Mk the maximum magnitude (by ignoring outliers) of the k first coefficients.For example, using the Irish consumption database [39], the distribution of the individual consumption of the first half customers (from 1 to 1818) in the frequency domain is given in Figure 5. In Figure 5, the maximum magnitudes (rounded) of the 5 first coefficients are M=M1,M2,M3,M4,M5=9,4,3,2,2.Figure 5
Distribution of users’ consumptions in the frequency domain using the discrete Fourier transform (DFT).The databaseD2 is used for testing our methodology. Let X=X1,…,Xn with Xj=x1i,…,xTi for all j=1,…,n be the users’ individual consumptions. To publish the sum of consumptions S=S1,…,ST=∑j=1nx1j,…,∑j=1nxTj, our methodology, which can be applied to either the Fourier transform or to wavelet transforms, is described as follows:(1)
For all individual consumptionsXj (j=1,…,n), compute the corresponding magnitude in the domain of the transform and keep the first k coefficients denoted by Cj=C1j,…,Ckj.(2)
If the modulus of coefficientCℓj is greater than Mℓ (1≤ℓ≤k), replace Cℓj with Cℓj⋅Mℓ/Cℓ so that all coefficients have a modulus smaller than Mℓ and their phase, if the coefficient is complex, is unchanged.(3)
Compute the sum of coefficientsC=∑j=1nC1j,…,∑j=1nCkj.(4)
Add a noise following the distribution of Laplaceℒ⋅, depending on the sensitivity of the transform, to each coefficient Cℓ1≤ℓ≤k of C. The result is denoted by C^. We note that the Laplace noise is added to the real and imaginary parts of each coefficient when the DFT is used.(5)
Pad the vectorC^ by n−k zeroes and compute the inverse transform to obtain the noisy version of the consumption S^.Section5.1 presents an adaptation of this methodology using the discrete Fourier transform.


### 5.1. Clamping Fourier Perturbation Algorithm

This section describes the clamping Fourier perturbation algorithm (CFPA) detailed in Algorithm4. This algorithm allows an aggregator to compute and publish an aggregate guaranteeing ε−differential privacy.Algorithm 4: Clamping Fourier perturbation algorithm (CFPA).
Inputs:(i)
Consumptions:X=X1,…,Xn with Xi=x1i,…,xTi for all i=1,…,n(ii)
k(iii)
The maximum magnitudes ofk first DFT coefficients: M=M1,…,Mk∈ℝ+k(iv)
Privacy budget:ε(1)
Clamping: for each individual time-series consumption Xj,(i)
compute thek first DFT coefficients of Xj: Fj=DFTXj1,…,DFTXjk=F1j,…,Fkj(ii)
ifFℓj>Mℓ, then replace Fℓj with Fℓj⋅Mℓ/Fℓj for all ℓ=1,…,k(2)
Laplacian mechanism: compute the sum of noisy consumptions of each DWT coefficient: F^ℓ=∑j=1nFℓj+ℒMℓ2/ε/k for all ℓ=1,…,k. We denote F^k=F^1,…,F^k. We note that the noise is added to the real and imaginary parts of the sum of coefficients.(3)
PadF^k with T−k zeros; the result is denoted by PADTF^k(4)
Compute the inverse DFT ofPADTF^k to get the noisy sum of consumptions denoted by S^=S^1,…,S^T of the initial sum S=∑j=1nx1j,…,∑j=1nxTj.CFPA takes as inputs the individual time-series consumptions ofn consumers, the maximum magnitudes of DFT coefficients of individual consumptions M (computed over database D1), the number k of DFT coefficients to be considered, and the privacy budget ε, and it returns the noisy time-series sum of consumptions of n consumers.Step 1, called clamping, computes the firstk DFT coefficients of each individual time-series consumption. If the magnitude of a coefficient Fℓj is greater than the maximum magnitude Mℓ, then this coefficient is clamped and replaced by Fℓj⋅Mℓ/Fℓj, in which magnitude is Fℓj⋅Mℓ/Fℓj=Mℓ. Thus, for all individual consumptions Xj, the maximum magnitude of the k first DFT coefficients Fj=F1j,…,Fkj is M=M1,…,Mk, i.e., the final values of the coefficients have the same phase as the initial values, but their magnitudes are bounded by M1,…,Mk.After computing the firstk DFT coefficients Fj=F1j,…,Fkj of each individual time-series consumption of consumers (j=1,…,n), Step 2 consists in computing the sum F1,…,Fk=∑j=1nF1j,…,∑j=1nFkj of these coefficients using the Laplacian mechanism. The result is denoted by F^k=F^1,…,F^k.Finally, the noisy sum of consumptions is equal to the inverse of the noisy DFT coefficients padded withT−k zeros.Theorem 3.
Algorithm CFPA isε-differentially private.Proof.
To prove that Algorithm4 is ε-differentially private, we need to prove that the sensitivity of the sum of DFT coefficients of users’ individual consumptions F1 (resp. F2,…,Fk) is 2⋅M1 (resp. 2⋅M2,…,2⋅Mk). This is done in Lemma 2.
Then, as a Laplacian noiseℒMℓ2/ε/k is added to each component Fℓ (ℓ=1,…,k), the resulting F^1 (resp. F^2,…,F^k) is ε/k-differentially private. Finally, the composition theorem [14] guarantees that any computation on the k components of F^1,…,F^k is ε-differentially private; thus, the inverse DFT of those coefficients is ε-DP.Lemma 2.
LetFj=DFTXj1,…,DFTXjk=F1j,…,Fkj be the first k DFT coefficients of the individual consumption of consumer j (j=1,…,n), obtained after the clamping mechanism. The sensitivity of the sum of each DFT coefficient Fℓj (ℓ=1,…,k) of n consumers’ individual consumptions is Mℓ⋅2.Proof.
LetFj=DFTXj1,…,DFTXjk=F1j,…,Fkj. After the clamping, the magnitude of each DFT coefficient Fℓj is smaller than Mℓ for ℓ=1,…,k, and the sensitivity of the function fℓ:Dn⟶ℂ≡ℝ2 defined by fℓ:X1,…,Xn⟼∑j=1nFℓj≡∑j=1naℓj,∑j=1nbℓj, with Fℓj=aℓj+ibℓj being equal to(9)Δ1fℓ=maxaℓj,bℓj∑j=1naℓj,∑j=1nbℓj−∑j=2naℓj,∑j=2nbℓj1=maxaℓ1,bℓ11≤max2aℓ1,bℓ12Minkowski inequality=max2aℓ12+bℓ12=max2Fℓ1=Mℓ2.Thus, Lemma2 proves Theorem 3, and algorithm CFPA guarantees ε-differential privacy.For example, Figure6 shows the same aggregated consumption presented in Example 1 and its noisy version using CFPA (Algorithm 4) with ε=1 per day and k=5. Figure 6 shows that CFPA obtains a good utility with an MRE equal to 9.7%. This good utility of CFPA can be explained by the fact that Laplace noise added in CFPA depends on the amplitude of each coefficient, while in FPA, the same noise ℒM⋅2Tk/ε is added to every DFT coefficients, where M is the maximum consumption in the dataset.Figure 6
Aggregated consumption of 250 homes from 30th December 2009 to 5th January 2010 of dataset from CER [39] and its noisy version using CFPA (Algorithm 4), with ε=1 per day and k=5.


### 5.2. Clamping Wavelet Perturbation Algorithm

The clamping wavelet perturbation algorithm (CWPA), as presented in Algorithm5, is obtained by replacing DFT with DWT in Algorithm 4. The computation of DWT is based on multiresolution analysis which determines the number of approximation coefficients (scaling functions) and detail coefficients (wavelet functions) [40]. DWT takes as input a time-series of length a power of 2. If the input’s length is not a power of 2, we can pad it with zeroes [41].Algorithm 5: Clamping wavelet perturbation algorithm (CWPA).
Inputs:(i)
Consumptions:X=X1,…,Xn with Xj=x1j,…,xTj for all j=1,…,n(ii)
k(iii)
The maximum magnitudes ofk first DWT coefficients: M=M1,…,Mk∈ℝ+k(iv)
Privacy budget:ε(1)
Clamping: for each individual time-series consumption Xj,(i)
compute the firstk DWT coefficients of Xj: Wj=DWTXj1,…,DWTXjk=W1j,…,Wkj(ii)
ifWℓj>Mℓ, then replace Wℓj with Wℓj⋅Mℓ/Wji for ℓ=1,…,k(2)
Laplacian Mechanism: compute the sum of noisy consumptions of each DWT coefficient:W^ℓ=∑j=1nWℓj+ℒMℓ/ε/k for all ℓ=1,…,k. We denote W^k=W^1,…,W^k(3)
PadW^k with T−k zeroes; the result is denoted by PADTW^k(4)
Compute the inverse DWT ofPADTW^k to get the noisy sum of consumptions denoted by S^=S^1,…,S^T of the initial sum S=∑j=1nx1j,…,∑j=1nxTjAlgorithm5 takes as inputs the maximum magnitudes of the first k DWT coefficients which are obtained in the training process on D1, by computing the distribution of DWT coefficients of individual consumptions.We note that there are multiple DWTs, such as Haar, Daubechies, Symlets, and Coiflets. In this paper, we use Haar and Daubechies wavelets as shown inSection 6, because they give a low reconstruction error, as will be discussed in Section 6.1.Theorem 4.
The clamping wavelet perturbation algorithm (CWPA), Algorithm5, is ε-differentially private.Proof.
The proof is similar to the one for Theorem3. We need to prove that the sensitivity of the sum of DWT coefficients of users’ individual consumptions W1 (resp. W2,…,Wk) is M1 (resp. M2,…,Mk). This is done in Lemma 3.
Then, as a Laplacian noiseℒMℓ/ε/k is added to each component Wℓ (ℓ=1,…,k), the resulting W^1 (resp. W^2,…,W^k) is ϵ/k-differentially private. Finally, the composition theorem [14] guarantees that any computation on the k components W^1,…,W^k is ϵ-differentially private; thus, the inverse DWT of those coefficients is ϵ-DP.Lemma 3.
LetWj=DWTXj1,…,DWTXjk=W1j,…,Wkj be the first k DWT coefficients of the individual consumption of consumer j (j=1,…,n), obtained after the clamping mechanism. The sensitivity of the sum of each DWT coefficient Wℓj (ℓ=1,…,k) of n consumers’ individual consumptions is Mℓ.Proof.
LetWj=DWTXj1,…,DWTXjk=W1j,…,Wkj. After the clamping, the magnitude of each DWT coefficient Wℓj is smaller than Mℓ for ℓ=1,…,k, and the sensitivity of the function wℓ:Dn⟶ℝ defined by wℓ:X1,…,Xn⟼∑j=1nWℓj is equal to(10)Δ1wℓ=maxWℓj∑j=1nWℓj−∑j=2nWℓj=maxWℓ1=Mℓ.For example, Figure7 shows the same aggregated consumption presented in Example 1 and its noisy version using CWPA (Algorithm 5) with Haar wavelet, ε=1 per day and k=5. However, Figure 3 shows that the MRE of CWPA is still higher than 10%. We explain this result in Section 6.Figure 7
Aggregated consumption of 250 homes from 30th December 2009 to 5th January 2010 of dataset from CER [39] and its noisy version using CWPA (Algorithm 5), with Haar wavelet, ε=1 per day and k=5.


## 5.1. Clamping Fourier Perturbation Algorithm

This section describes the clamping Fourier perturbation algorithm (CFPA) detailed in Algorithm4. This algorithm allows an aggregator to compute and publish an aggregate guaranteeing ε−differential privacy.Algorithm 4: Clamping Fourier perturbation algorithm (CFPA).
Inputs:(i)
Consumptions:X=X1,…,Xn with Xi=x1i,…,xTi for all i=1,…,n(ii)
k(iii)
The maximum magnitudes ofk first DFT coefficients: M=M1,…,Mk∈ℝ+k(iv)
Privacy budget:ε(1)
Clamping: for each individual time-series consumption Xj,(i)
compute thek first DFT coefficients of Xj: Fj=DFTXj1,…,DFTXjk=F1j,…,Fkj(ii)
ifFℓj>Mℓ, then replace Fℓj with Fℓj⋅Mℓ/Fℓj for all ℓ=1,…,k(2)
Laplacian mechanism: compute the sum of noisy consumptions of each DWT coefficient: F^ℓ=∑j=1nFℓj+ℒMℓ2/ε/k for all ℓ=1,…,k. We denote F^k=F^1,…,F^k. We note that the noise is added to the real and imaginary parts of the sum of coefficients.(3)
PadF^k with T−k zeros; the result is denoted by PADTF^k(4)
Compute the inverse DFT ofPADTF^k to get the noisy sum of consumptions denoted by S^=S^1,…,S^T of the initial sum S=∑j=1nx1j,…,∑j=1nxTj.CFPA takes as inputs the individual time-series consumptions ofn consumers, the maximum magnitudes of DFT coefficients of individual consumptions M (computed over database D1), the number k of DFT coefficients to be considered, and the privacy budget ε, and it returns the noisy time-series sum of consumptions of n consumers.Step 1, called clamping, computes the firstk DFT coefficients of each individual time-series consumption. If the magnitude of a coefficient Fℓj is greater than the maximum magnitude Mℓ, then this coefficient is clamped and replaced by Fℓj⋅Mℓ/Fℓj, in which magnitude is Fℓj⋅Mℓ/Fℓj=Mℓ. Thus, for all individual consumptions Xj, the maximum magnitude of the k first DFT coefficients Fj=F1j,…,Fkj is M=M1,…,Mk, i.e., the final values of the coefficients have the same phase as the initial values, but their magnitudes are bounded by M1,…,Mk.After computing the firstk DFT coefficients Fj=F1j,…,Fkj of each individual time-series consumption of consumers (j=1,…,n), Step 2 consists in computing the sum F1,…,Fk=∑j=1nF1j,…,∑j=1nFkj of these coefficients using the Laplacian mechanism. The result is denoted by F^k=F^1,…,F^k.Finally, the noisy sum of consumptions is equal to the inverse of the noisy DFT coefficients padded withT−k zeros.Theorem 3.
Algorithm CFPA isε-differentially private.Proof.
To prove that Algorithm4 is ε-differentially private, we need to prove that the sensitivity of the sum of DFT coefficients of users’ individual consumptions F1 (resp. F2,…,Fk) is 2⋅M1 (resp. 2⋅M2,…,2⋅Mk). This is done in Lemma 2.
Then, as a Laplacian noiseℒMℓ2/ε/k is added to each component Fℓ (ℓ=1,…,k), the resulting F^1 (resp. F^2,…,F^k) is ε/k-differentially private. Finally, the composition theorem [14] guarantees that any computation on the k components of F^1,…,F^k is ε-differentially private; thus, the inverse DFT of those coefficients is ε-DP.Lemma 2.
LetFj=DFTXj1,…,DFTXjk=F1j,…,Fkj be the first k DFT coefficients of the individual consumption of consumer j (j=1,…,n), obtained after the clamping mechanism. The sensitivity of the sum of each DFT coefficient Fℓj (ℓ=1,…,k) of n consumers’ individual consumptions is Mℓ⋅2.Proof.
LetFj=DFTXj1,…,DFTXjk=F1j,…,Fkj. After the clamping, the magnitude of each DFT coefficient Fℓj is smaller than Mℓ for ℓ=1,…,k, and the sensitivity of the function fℓ:Dn⟶ℂ≡ℝ2 defined by fℓ:X1,…,Xn⟼∑j=1nFℓj≡∑j=1naℓj,∑j=1nbℓj, with Fℓj=aℓj+ibℓj being equal to(9)Δ1fℓ=maxaℓj,bℓj∑j=1naℓj,∑j=1nbℓj−∑j=2naℓj,∑j=2nbℓj1=maxaℓ1,bℓ11≤max2aℓ1,bℓ12Minkowski inequality=max2aℓ12+bℓ12=max2Fℓ1=Mℓ2.Thus, Lemma2 proves Theorem 3, and algorithm CFPA guarantees ε-differential privacy.For example, Figure6 shows the same aggregated consumption presented in Example 1 and its noisy version using CFPA (Algorithm 4) with ε=1 per day and k=5. Figure 6 shows that CFPA obtains a good utility with an MRE equal to 9.7%. This good utility of CFPA can be explained by the fact that Laplace noise added in CFPA depends on the amplitude of each coefficient, while in FPA, the same noise ℒM⋅2Tk/ε is added to every DFT coefficients, where M is the maximum consumption in the dataset.Figure 6
Aggregated consumption of 250 homes from 30th December 2009 to 5th January 2010 of dataset from CER [39] and its noisy version using CFPA (Algorithm 4), with ε=1 per day and k=5.


## 5.2. Clamping Wavelet Perturbation Algorithm

The clamping wavelet perturbation algorithm (CWPA), as presented in Algorithm5, is obtained by replacing DFT with DWT in Algorithm 4. The computation of DWT is based on multiresolution analysis which determines the number of approximation coefficients (scaling functions) and detail coefficients (wavelet functions) [40]. DWT takes as input a time-series of length a power of 2. If the input’s length is not a power of 2, we can pad it with zeroes [41].Algorithm 5: Clamping wavelet perturbation algorithm (CWPA).
Inputs:(i)
Consumptions:X=X1,…,Xn with Xj=x1j,…,xTj for all j=1,…,n(ii)
k(iii)
The maximum magnitudes ofk first DWT coefficients: M=M1,…,Mk∈ℝ+k(iv)
Privacy budget:ε(1)
Clamping: for each individual time-series consumption Xj,(i)
compute the firstk DWT coefficients of Xj: Wj=DWTXj1,…,DWTXjk=W1j,…,Wkj(ii)
ifWℓj>Mℓ, then replace Wℓj with Wℓj⋅Mℓ/Wji for ℓ=1,…,k(2)
Laplacian Mechanism: compute the sum of noisy consumptions of each DWT coefficient:W^ℓ=∑j=1nWℓj+ℒMℓ/ε/k for all ℓ=1,…,k. We denote W^k=W^1,…,W^k(3)
PadW^k with T−k zeroes; the result is denoted by PADTW^k(4)
Compute the inverse DWT ofPADTW^k to get the noisy sum of consumptions denoted by S^=S^1,…,S^T of the initial sum S=∑j=1nx1j,…,∑j=1nxTjAlgorithm5 takes as inputs the maximum magnitudes of the first k DWT coefficients which are obtained in the training process on D1, by computing the distribution of DWT coefficients of individual consumptions.We note that there are multiple DWTs, such as Haar, Daubechies, Symlets, and Coiflets. In this paper, we use Haar and Daubechies wavelets as shown inSection 6, because they give a low reconstruction error, as will be discussed in Section 6.1.Theorem 4.
The clamping wavelet perturbation algorithm (CWPA), Algorithm5, is ε-differentially private.Proof.
The proof is similar to the one for Theorem3. We need to prove that the sensitivity of the sum of DWT coefficients of users’ individual consumptions W1 (resp. W2,…,Wk) is M1 (resp. M2,…,Mk). This is done in Lemma 3.
Then, as a Laplacian noiseℒMℓ/ε/k is added to each component Wℓ (ℓ=1,…,k), the resulting W^1 (resp. W^2,…,W^k) is ϵ/k-differentially private. Finally, the composition theorem [14] guarantees that any computation on the k components W^1,…,W^k is ϵ-differentially private; thus, the inverse DWT of those coefficients is ϵ-DP.Lemma 3.
LetWj=DWTXj1,…,DWTXjk=W1j,…,Wkj be the first k DWT coefficients of the individual consumption of consumer j (j=1,…,n), obtained after the clamping mechanism. The sensitivity of the sum of each DWT coefficient Wℓj (ℓ=1,…,k) of n consumers’ individual consumptions is Mℓ.Proof.
LetWj=DWTXj1,…,DWTXjk=W1j,…,Wkj. After the clamping, the magnitude of each DWT coefficient Wℓj is smaller than Mℓ for ℓ=1,…,k, and the sensitivity of the function wℓ:Dn⟶ℝ defined by wℓ:X1,…,Xn⟼∑j=1nWℓj is equal to(10)Δ1wℓ=maxWℓj∑j=1nWℓj−∑j=2nWℓj=maxWℓ1=Mℓ.For example, Figure7 shows the same aggregated consumption presented in Example 1 and its noisy version using CWPA (Algorithm 5) with Haar wavelet, ε=1 per day and k=5. However, Figure 3 shows that the MRE of CWPA is still higher than 10%. We explain this result in Section 6.Figure 7
Aggregated consumption of 250 homes from 30th December 2009 to 5th January 2010 of dataset from CER [39] and its noisy version using CWPA (Algorithm 5), with Haar wavelet, ε=1 per day and k=5.


## 6. Experimental Results

This section compares FPA, CFPA, WPA, and CWPA and explains through experimentations why CFPA achieves a better utility than other publication techniques. After presenting the raw results, we explain them by decomposing the mean relative error into a perturbation error, caused by the clamping mechanism and the Laplace mechanism, and a reconstruction error, due to ignoringT−k coefficients of the transform. The analysis of the error is thus conducted in the next two Subsections 6.1 and 6.2. Section 6.1 analyzes the reconstruction error, while Section 6.2 analyzes the perturbation one.Conditions: the experiments rely on data originating from the Irish Commission for Energy Regulation (CER) [39]. This dataset contains real time-series consumptions. The achieved results are valid for this very specific case, for Irish consumptions with an Irish weather being never too hot or too cold. The results show that the approach is good, but will probably have to be adapted for other datasets, i.e., by computing the maximum magnitudes of the k first coefficients of the considered transform over a subpart of the dataset. Consumption data from the CER were collected every 30 minutes from 2009 to 2010 with the participation of more than 5,000 Irish homes and businesses. This experiment only considers homes. We divided the database in two parts: D1, corresponding to the first half of consumers (1 to 1818), and D2, corresponding to the second half (1819 to 3639). D1 is used to calibrate the algorithms by computing the maximum magnitudes M=M1,…,Mk of the first k coefficients in the frequency domain, and D2 is used to test the publication techniques FPA, CFPA, WPA, and CWPA.Notations: we noteN as the number of homes or smart meters considered in the district to compute the time-series of the sum of consumptions. For each day (48 time slots), we compute the sum of consumptions of 50 different districts of N random homes, and we execute FPA, WPA, CFPA, and CWPA with privacy budget ε∈1;3 for each day and k∈5;8;12. The discrete wavelet transforms used here are Haar transform (which represents the same wavelet as Daubechies with order 1, noted db1), Daubechies with order 2, and Daubechies with order 3, respectively, noted db2 and db3.Raw results and analysis: Figures8 and 9 show the distribution of the mean relative estimation error (MRE) according to the number of homes in the district (N from 50 to 450) and k from 5 to 12 for the budget of privacy ε=1 and ε=3, respectively. The boxes extend from the lower to upper quartile values of the MRE, with a line at the median and a triangle representing the mean. The whiskers extend from the box to show the range of the MRE. In order to make consumption forecasts, an MRE lower than 10% is required in practice by experts in the energy sector. In this section, an MRE of less than 10% is therefore considered useful.Figure 8
Mean relative estimation error (MRE) of FPA vs CFPA vs WPA vs CWPA, using DFT and DWT with Haar, Daubechies 2 (db2), and Daubechies 3 (db3) wavelets, according tok, and the number of smart meters (N) in the district, with ε=1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 9
Mean relative estimation error (MRE) of FPA vs CWPA vs WPA vs CWPA with Haar, Daubechies 2 (db2), and Daubechies 3 (db3) wavelets, according tok, and the number of smart meters (N) in the district, with ε=3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)In Figures8 and 9, the first column corresponds to the comparison between the FPA and the CFPA. The other columns correspond to the comparison between the WPA and the CWPA using Haar wavelet with 2 approximation coefficients, Daubechies 2 (db2) with 5 approximation coefficients, and Daubechies 3 (db3) with 10 approximation coefficients, respectively.Figures8 and 9 show that CFPA has a better utility than FPA. For example, for ε=1 (Figure 8), when k=5 and the number of homes N=350, the MRE of CFPA is 12%, while the MRE of FPA is 75%. In that configuration, the MRE of CFPA is 6.25 times lower than that of FPA. Similarly, the CWPA obtains a better utility than the WPA. For example, for k=5 and the number of homes N=350, the MRE of CWPA using Haar wavelet is 15%, while the MRE of the WPA is 30%. In that configuration, the MRE of CWPA is 2 times lower than that of WPA.Generally, the larger the size of the districtN, the smaller the MRE is. Similarly, the larger the budget of privacy ε, the smaller the MRE is; Figure 9 (ε=3) shows a better utility than Figure 8 (ε=1). However, Figures 8 and 9 show that WPA and CWPA using db3 are not useful for k=5 because, as shown in Section 6.1, the reconstruction error is high (between 70% and 80%).Figures8 and 9 show that for larger k, the MRE of WPA and CWPA using db3 is smaller. Moreover, in Figure 9, for k=12 and when the number of homes is higher than 250 and ε=3, CWPA using db3 has the median of MRE smaller than 11%. The CWPA using Haar wavelet obtains the second best utility, with the median of MRE smaller than 10% when N is greater than 250, while the CFPA gets the best utility, with the median of MRE decreasing to 5% when N=550 and k=8.However, the utility of FPA and WPA decreases whenk increases. This is caused by the perturbation error; indeed, the greater the k, the greater the Laplacian noise added to each coefficient is. This noise is attenuated by the clamping as shown by the CFPA. Indeed, when k goes from 5 to 8, the reconstruction error decreases and the clamping also decreases the perturbation error leading to the total error reduction. However, when k goes from 8 to 12, although the reconstruction error decreases, clamping does not reduce the perturbation error sufficiently. This explains why the MRE of CFPA is a little bigger when k=12 compared to k=8.In Figures8 and 9, we notice that the median of MRE of WPA and CWPA converge to a threshold and never goes below it. For example, for ε=3 and k=5, the median of MRE of WPA and CWPA using db2 converges to 23%. This is caused by the reconstruction error.


### 6.1. Reconstruction Error

The reconstruction error is due to considering only thek first transform coefficients, thus removing the precision brought by coefficients (k+1,k+2,…). To measure this error, a first solution consists in computing the cumulative distribution function (CDF) of the coefficients as a first assessment of the impact of the transform coefficients and, then, to get confirmation through some experimental reconstruction error measurements. Intuitively, if the CDF of some coefficients k is close to 1, it means that the coefficients after k (k+1,k+2,…) have less impact on the reconstruction, and thus, when set to zero, lead to a smaller reconstruction error.The CDF is computed for a district of 50 homes of several transformations: discrete Fourier transform (DFT), discrete wavelet transform (DWT) using Haar, and Daubechies 2 and Daubechies 3 wavelets. The closer to 1 the cumulative distribution function atk is, the smaller the reconstruction error is. Figure 10 compares the cumulative distribution function of DFT and DWT with different wavelet transforms. This figure shows that DFT has a higher cumulative distribution than DWT for the considered range value of k (k≤10).In order to analyze this error more precisely, we define formally the reconstruction error below, and we then compute it experimentally.Figure 10
Comparison of cumulative distribution function of DFT and DWT with Haar, and Daubechies 2 and 3 wavelets for a district of 50 homes.Definition 4 (reconstruction error).
LetS=S1,S2,…,ST be a sum of time-series consumptions and C=C1,C2,…,CT be the coefficients in the frequency domain of this time-series. We denote PADTCk=C1,…,Ck,0,…,0 as the first k coefficients padded with zeros and S˜=S˜1,S˜2,…,S˜T as the inverse of PADTCk (in the time domain). The reconstruction error (RE) of PADTCk is equal to the mean relative estimation error between S and S˜ given by (we add 1 to the denominator in order to avoid the division by zero)(11)REPADTCk=1T⋅∑j=1TSj−S˜jSj+1.Figure11 shows the reconstruction error for DFT and DWT with different wavelet transforms for a district of 50 and 450 homes. This figure shows that the DFT obtains the smallest relative error (lower than 10% when k is greater than 5) followed by Haar and Daubechies. We note that the reconstruction error of Daubechies 2 is higher than 23% when k=5, which leads to a total error higher than 23% and justifies the relative error obtained in Figures 8 and 9.Figure 11
Comparison of reconstruction error of DFT and DWT with Haar amd Daubechies 2 and Daubechies 3 wavelets for a district of 50 and 450 homes. (a)N = 50. (b) N = 450.
(a)(b)Moreover, whenk=5, the reconstruction error of Daubechies 3 is higher than 70%, which justifies why its total error is higher than 70% when k=5, according to Figures 8 and 9.According to the database from the Irish Commission for Energy Regulation (CER) [39], the discrete Fourier transform gets the smaller reconstruction error, followed respectively by Haar (which is the same as Daubechies 1) and Daubechies 2 and Daubechies 3 wavelets.


### 6.2. Perturbation Error

The perturbation error is caused by the Laplace mechanism, applied on the firstk transform coefficients. The higher the transform coefficients, the lower the impact of this perturbation in terms of relative error, and thus the lower the perturbation error.We note that the amplitude of the Laplace noise introduced by the Laplace mechanism is different for CFPA and CWPA; it is2 times greater for CFPA than for CWPA. Indeed, for all ℓ=1,…,k, the parameter for the Laplace noise is ℒMℓ2/ε/k for CFPA and ℒMℓ/ε/k for CWPA. Moreover, in the CFPA, 2k coefficients (the real and imaginary parts of the k DFT coefficients) are noisy while only k coefficients are noisy in the CWPA.For a district of 50 homes, we compute the distribution of the magnitude of DFT and DWT with Haar amd Daubechies 2 and 3 wavelets, and we compare their coefficient distribution median in Figure12.Figure 12
Comparison of the median of coefficients distribution of DFT and DWT with Haar and Daubechies 2 and Daubechies 3 wavelets for a district of 50 homes.Figure12 shows that the coefficient values vary according to the values of k and the considered transforms. For instance, when k is in the interval 7,10, Daubechies 3 obtains the highest magnitudes of coefficients, followed by Daubechies 2 and DFT.In clamping perturbation algorithms (CFPA, CWPA), the clamping mechanism allows to add a noise proportional to the modulus of the coefficients of the considered transform (DFT, DWT). This reduces the impact of noise compared to perturbation algorithms (FPA, WPA); however, at the price of a perturbation error induced by the clamping of the coefficients. Formally, the perturbation error of clamping perturbation algorithms (CFPA, CWPA) is defined as follows:Definition 5.
Perturbation error for clamping perturbation algorithms (CFPA, CWPA).
LetX1,…,XN be the individual time-series of energy consumptions of N homes, with Xi=x1i,…,xTi for i=1,…,N. The sum of time-series consumptions is noted as S=S1,…,ST=∑i=1Nx1i,…,∑i=1NxTi. For all i=1,…N, we note C¯i=c¯1i,…,c¯ki,ck+1i,…,cTi as the result of the considered transform of the time-series consumption Xi whose first k coefficients c¯1i,…,c¯ki have been clamped. We note M=M1,…,Mk as the maximum magnitude of the first k coefficients of the considered transform. Let C¯=∑i=1Nc¯1i+ℒδ1/ε/k,…,∑i=1Nc¯ki+ℒδk/ε/k,∑i=1Nck+1i,…,∑i=1NcTi be the sum of coefficients of the considered transform by perturbing only the first k coefficients, with δj=Mj2 (respectively δj=Mj) for CFPA (respectively for CWPA), for j=1,…,k.
LetS¯=S1¯,…,ST¯ be the inverse transform of C¯. The perturbation error of C¯ equals to the mean relative estimation error (MRE) between S and S¯, given by(we add 1 to the denominator in order to avoid the division by zero). For CWPA, the Laplace noise ℒMj2/ε/k must be replaced by ℒMj/ε/k for j=1,…,k and the DFT by the DWT,(12)PEC¯=1T∑ℓ=1TSℓ−S¯ℓSℓ+1.The perturbation error depends on the following parameters,k, Mj, ε, and N for j=1,…,k. k, Mj (j=1,…,k) and ε are parameters of the Laplace distribution, so they have a direct impact on the amplitude of the added noise. Let ε and Mj be fixed; the bigger the k, the smaller the Laplace distribution parameter δk/ε/k is, and thus, the bigger the noise added on the k first coefficients is. This makes the perturbation error more significant. The choice of Mj is important to define the clamping threshold and it directly impacts the perturbation of the Laplace mechanism. The greater the Mj, the bigger the Laplace noise is, and thus, the more the perturbation error is. The smaller the Mj (close to zero), the less the Laplace noise is, but the more the coefficients are clamped, and thus, the more the perturbation error is. The number of homes N indirectly plays a role in the perturbation error; the larger the N, the more diluted the added noise is. This leads to decrease the perturbation error.Figure13 (respectively, Figure 14) shows the distribution of the perturbation error of the clamping perturbation algorithms (CFPA and CWPA) according to k, N, with ε=1 (respectively, ε=3).Figure 13
Distribution of perturbation error according to clamping perturbation algorithms CFPA and CWPA, with Haar, Daubechies 2, and Daubechies 3, and according tok and the number of homes N, for a fixed privacy budget, ε=1. Note that the scales in (a)–(d) are different. (a) N = 50. (b) N = 250. (c) N = 450. (d). N = 550.
(a)(b)(c)(d)Figure 14
Distribution of perturbation error according to clamping perturbation algorithms CFPA and CWPA, with Haar, Daubechies 2, and Daubechies 3, and according tok and the number of homes N, for a fixed privacy budget, ε=3. Note that the scales in 14(a)–14(d) are different. (a) N = 50 (b) N = 250. (c) N = 450. (d). N = 550.
(a)(b)(c)(d)Figures13 and 14 show that the perturbation error of CFPA is higher than that of CWPA. This result is explained by(1)
The number of coefficients to be noised in CFPA is twice as many as the number of coefficients to be noised in CWPA. Indeed, in CFPA, the DFT coefficients are complex numbers, so both real and imaginary parts must be noised.(2)
The absolute value of the noise added in the CFPA is2 times greater than that in the CWPA.In addition, the greater theN, the more the added noise is diluted in the aggregate, causing the perturbation error to decrease. E.g., for ε=1 (Figure 13), when k=5, the median of the perturbation error of CFPA (respectively, CWPA with Haar) goes from 70% to 7% (respectively from 32% to 4%) when N goes from 50 to 550. Likewise, for ε=3 (Figure 14), when k=5, the median of perturbation error of CFPA (respectively, CWPA with Haar) goes from 25% to 2.7% (respectively from 11% to 2.1%) when N goes from 50 to 550. We notice that, the greater the N, the smaller the difference of the perturbation error between CFPA and CWPA is. This result is also true when ε increases. This can be explained by the decrease of the noise introduced on the coefficients of the transforms.Figures13 and 14 show that the perturbation error increases when k increases. The larger the k, the smaller the budget ε/k allocated to each coefficient is. This leads to a noise increase on each coefficient and thus on the perturbation error.


### 6.3. Summary of the Experimental Results

The combination of the reconstruction error (Figure11) and the perturbation error (Figures 13 and 14) enables to determine which transform is appropriate according to the number of homes N and the budget of privacy ε, for getting a total error as small as possible.Lemma 4.
The mean relative error (MRE) of CFPA (respectively, CWPA) is lower than or equal to the sum of the reconstruction error and the perturbation error of CFPA (respectively, CWPA).Proof.
The proof of the above lemma deferred to the appendix.Section6.1 shows that the reconstruction error of DFT is lower than that of the considered DWT. For example, when N=450 and k=5, the median of the reconstruction error is 6% for DFT, while it is 13% for the Haar and Daubechies 2 transforms.However, Section6.2 shows that algorithms based on DFT (e.g., CFPA) have a higher perturbation error than those based on DWT (e.g., CWPA).According to Lemma4, the total error (MRE) is less than or equal to the sum of reconstruction error and perturbation error. Thus, if the reconstruction error or perturbation error is greater than 10%, there is a high probability that the final error will not be less than this threshold.As the reconstruction error of the DWT is greater than9%, there is a high probability that the final error of CWPA will not be less than this threshold, even if the Laplace noise decreases, i.e., when the number of homes N or the privacy budget ε increases. However, as the reconstruction error of the DFT is small (the median is between 2% and 3% when k=7,8,9), then the total error of the CFPA may be lower than that of the CWPA when the impact of Laplace noise decreases. For example, the median of the perturbation error of CFPA is between 3% and 5% when k=7,8,9, N=550, and ε=3. This analysis explains why, for ε=1, the CWPA obtains a better utility than the CFPA when the number of homes N is less than 250. For example, when N=50 and k=5, the median of the perturbation error (respectively, the reconstruction error) of CFPA is 70% (respectively, 10%) against 32% (respectively, 18%) for CWPA using Haar. Thus, the median of MRE of CFPA is between 70% and 80% against 32% and 50% for CWPA.WhenN is higher than 250, CFPA gets a better utility than CWPA. For example, when N=450 and k=5, the median of the perturbation error (respectively, the reconstruction error) of CFPA is 8.5% (respectively, 6.5%) against 4.5% (respectively, 16%) for CWPA using Haar. Thus, the median of MRE of CFPA is between 8.5% and 15% against 16% and 20.5% for CWPA.In this use case, by comparing the different techniques for publishing time-series consumption, it appears that clamping perturbation algorithms (CFPA, CWPA) get a better utility than unbounded algorithms (FPA, WPA), which shows that the clamping mechanism reduces the total error. Furthermore, when the number of homes is greater than 250, CFPA obtains the best utility, with a mean relative error of less than10% when ε=3. When the budget of privacy ϵ=1, the mean relative error of CFPA is less than 10% for N=450 homes.The CWPA gets the best utility when the number of homesN is smaller than 150 and the budget of privacy ε is 1. This is justified by its low perturbation error.Table3 summarizes the publishing algorithm with the smallest MRE according to the budget of privacy ε and the number of homes in the district (N). Based on the dataset from the Irish Commission for Energy Regulation (CER) [39], Table 3 shows that the clamping Fourier perturbation algorithm (CFPA) achieves a lower MRE than the clamping wavelet perturbation algorithm (CWPA) for N>150. Hence, CFPA gets a better utility than CWPA for N>150.Table 3
Publishing algorithm with the smallest MRE according to the budget of privacyε and the number of homes in the district (N).
Number of homes (N)Budget of privacy,εBest algorithmCoefficients,kMedian of MRE (%)501CWPA-Haar, db25353CWPA-Haar211501CWPA-Haar5193CFPA112501CFPA516383501CFPA5123874501CFPA5103865501CFPA59385


## 6.1. Reconstruction Error

The reconstruction error is due to considering only thek first transform coefficients, thus removing the precision brought by coefficients (k+1,k+2,…). To measure this error, a first solution consists in computing the cumulative distribution function (CDF) of the coefficients as a first assessment of the impact of the transform coefficients and, then, to get confirmation through some experimental reconstruction error measurements. Intuitively, if the CDF of some coefficients k is close to 1, it means that the coefficients after k (k+1,k+2,…) have less impact on the reconstruction, and thus, when set to zero, lead to a smaller reconstruction error.The CDF is computed for a district of 50 homes of several transformations: discrete Fourier transform (DFT), discrete wavelet transform (DWT) using Haar, and Daubechies 2 and Daubechies 3 wavelets. The closer to 1 the cumulative distribution function atk is, the smaller the reconstruction error is. Figure 10 compares the cumulative distribution function of DFT and DWT with different wavelet transforms. This figure shows that DFT has a higher cumulative distribution than DWT for the considered range value of k (k≤10).In order to analyze this error more precisely, we define formally the reconstruction error below, and we then compute it experimentally.Figure 10
Comparison of cumulative distribution function of DFT and DWT with Haar, and Daubechies 2 and 3 wavelets for a district of 50 homes.Definition 4 (reconstruction error).
LetS=S1,S2,…,ST be a sum of time-series consumptions and C=C1,C2,…,CT be the coefficients in the frequency domain of this time-series. We denote PADTCk=C1,…,Ck,0,…,0 as the first k coefficients padded with zeros and S˜=S˜1,S˜2,…,S˜T as the inverse of PADTCk (in the time domain). The reconstruction error (RE) of PADTCk is equal to the mean relative estimation error between S and S˜ given by (we add 1 to the denominator in order to avoid the division by zero)(11)REPADTCk=1T⋅∑j=1TSj−S˜jSj+1.Figure11 shows the reconstruction error for DFT and DWT with different wavelet transforms for a district of 50 and 450 homes. This figure shows that the DFT obtains the smallest relative error (lower than 10% when k is greater than 5) followed by Haar and Daubechies. We note that the reconstruction error of Daubechies 2 is higher than 23% when k=5, which leads to a total error higher than 23% and justifies the relative error obtained in Figures 8 and 9.Figure 11
Comparison of reconstruction error of DFT and DWT with Haar amd Daubechies 2 and Daubechies 3 wavelets for a district of 50 and 450 homes. (a)N = 50. (b) N = 450.
(a)(b)Moreover, whenk=5, the reconstruction error of Daubechies 3 is higher than 70%, which justifies why its total error is higher than 70% when k=5, according to Figures 8 and 9.According to the database from the Irish Commission for Energy Regulation (CER) [39], the discrete Fourier transform gets the smaller reconstruction error, followed respectively by Haar (which is the same as Daubechies 1) and Daubechies 2 and Daubechies 3 wavelets.


## 6.2. Perturbation Error

The perturbation error is caused by the Laplace mechanism, applied on the firstk transform coefficients. The higher the transform coefficients, the lower the impact of this perturbation in terms of relative error, and thus the lower the perturbation error.We note that the amplitude of the Laplace noise introduced by the Laplace mechanism is different for CFPA and CWPA; it is2 times greater for CFPA than for CWPA. Indeed, for all ℓ=1,…,k, the parameter for the Laplace noise is ℒMℓ2/ε/k for CFPA and ℒMℓ/ε/k for CWPA. Moreover, in the CFPA, 2k coefficients (the real and imaginary parts of the k DFT coefficients) are noisy while only k coefficients are noisy in the CWPA.For a district of 50 homes, we compute the distribution of the magnitude of DFT and DWT with Haar amd Daubechies 2 and 3 wavelets, and we compare their coefficient distribution median in Figure12.Figure 12
Comparison of the median of coefficients distribution of DFT and DWT with Haar and Daubechies 2 and Daubechies 3 wavelets for a district of 50 homes.Figure12 shows that the coefficient values vary according to the values of k and the considered transforms. For instance, when k is in the interval 7,10, Daubechies 3 obtains the highest magnitudes of coefficients, followed by Daubechies 2 and DFT.In clamping perturbation algorithms (CFPA, CWPA), the clamping mechanism allows to add a noise proportional to the modulus of the coefficients of the considered transform (DFT, DWT). This reduces the impact of noise compared to perturbation algorithms (FPA, WPA); however, at the price of a perturbation error induced by the clamping of the coefficients. Formally, the perturbation error of clamping perturbation algorithms (CFPA, CWPA) is defined as follows:Definition 5.
Perturbation error for clamping perturbation algorithms (CFPA, CWPA).
LetX1,…,XN be the individual time-series of energy consumptions of N homes, with Xi=x1i,…,xTi for i=1,…,N. The sum of time-series consumptions is noted as S=S1,…,ST=∑i=1Nx1i,…,∑i=1NxTi. For all i=1,…N, we note C¯i=c¯1i,…,c¯ki,ck+1i,…,cTi as the result of the considered transform of the time-series consumption Xi whose first k coefficients c¯1i,…,c¯ki have been clamped. We note M=M1,…,Mk as the maximum magnitude of the first k coefficients of the considered transform. Let C¯=∑i=1Nc¯1i+ℒδ1/ε/k,…,∑i=1Nc¯ki+ℒδk/ε/k,∑i=1Nck+1i,…,∑i=1NcTi be the sum of coefficients of the considered transform by perturbing only the first k coefficients, with δj=Mj2 (respectively δj=Mj) for CFPA (respectively for CWPA), for j=1,…,k.
LetS¯=S1¯,…,ST¯ be the inverse transform of C¯. The perturbation error of C¯ equals to the mean relative estimation error (MRE) between S and S¯, given by(we add 1 to the denominator in order to avoid the division by zero). For CWPA, the Laplace noise ℒMj2/ε/k must be replaced by ℒMj/ε/k for j=1,…,k and the DFT by the DWT,(12)PEC¯=1T∑ℓ=1TSℓ−S¯ℓSℓ+1.The perturbation error depends on the following parameters,k, Mj, ε, and N for j=1,…,k. k, Mj (j=1,…,k) and ε are parameters of the Laplace distribution, so they have a direct impact on the amplitude of the added noise. Let ε and Mj be fixed; the bigger the k, the smaller the Laplace distribution parameter δk/ε/k is, and thus, the bigger the noise added on the k first coefficients is. This makes the perturbation error more significant. The choice of Mj is important to define the clamping threshold and it directly impacts the perturbation of the Laplace mechanism. The greater the Mj, the bigger the Laplace noise is, and thus, the more the perturbation error is. The smaller the Mj (close to zero), the less the Laplace noise is, but the more the coefficients are clamped, and thus, the more the perturbation error is. The number of homes N indirectly plays a role in the perturbation error; the larger the N, the more diluted the added noise is. This leads to decrease the perturbation error.Figure13 (respectively, Figure 14) shows the distribution of the perturbation error of the clamping perturbation algorithms (CFPA and CWPA) according to k, N, with ε=1 (respectively, ε=3).Figure 13
Distribution of perturbation error according to clamping perturbation algorithms CFPA and CWPA, with Haar, Daubechies 2, and Daubechies 3, and according tok and the number of homes N, for a fixed privacy budget, ε=1. Note that the scales in (a)–(d) are different. (a) N = 50. (b) N = 250. (c) N = 450. (d). N = 550.
(a)(b)(c)(d)Figure 14
Distribution of perturbation error according to clamping perturbation algorithms CFPA and CWPA, with Haar, Daubechies 2, and Daubechies 3, and according tok and the number of homes N, for a fixed privacy budget, ε=3. Note that the scales in 14(a)–14(d) are different. (a) N = 50 (b) N = 250. (c) N = 450. (d). N = 550.
(a)(b)(c)(d)Figures13 and 14 show that the perturbation error of CFPA is higher than that of CWPA. This result is explained by(1)
The number of coefficients to be noised in CFPA is twice as many as the number of coefficients to be noised in CWPA. Indeed, in CFPA, the DFT coefficients are complex numbers, so both real and imaginary parts must be noised.(2)
The absolute value of the noise added in the CFPA is2 times greater than that in the CWPA.In addition, the greater theN, the more the added noise is diluted in the aggregate, causing the perturbation error to decrease. E.g., for ε=1 (Figure 13), when k=5, the median of the perturbation error of CFPA (respectively, CWPA with Haar) goes from 70% to 7% (respectively from 32% to 4%) when N goes from 50 to 550. Likewise, for ε=3 (Figure 14), when k=5, the median of perturbation error of CFPA (respectively, CWPA with Haar) goes from 25% to 2.7% (respectively from 11% to 2.1%) when N goes from 50 to 550. We notice that, the greater the N, the smaller the difference of the perturbation error between CFPA and CWPA is. This result is also true when ε increases. This can be explained by the decrease of the noise introduced on the coefficients of the transforms.Figures13 and 14 show that the perturbation error increases when k increases. The larger the k, the smaller the budget ε/k allocated to each coefficient is. This leads to a noise increase on each coefficient and thus on the perturbation error.


## 6.3. Summary of the Experimental Results

The combination of the reconstruction error (Figure11) and the perturbation error (Figures 13 and 14) enables to determine which transform is appropriate according to the number of homes N and the budget of privacy ε, for getting a total error as small as possible.Lemma 4.
The mean relative error (MRE) of CFPA (respectively, CWPA) is lower than or equal to the sum of the reconstruction error and the perturbation error of CFPA (respectively, CWPA).Proof.
The proof of the above lemma deferred to the appendix.Section6.1 shows that the reconstruction error of DFT is lower than that of the considered DWT. For example, when N=450 and k=5, the median of the reconstruction error is 6% for DFT, while it is 13% for the Haar and Daubechies 2 transforms.However, Section6.2 shows that algorithms based on DFT (e.g., CFPA) have a higher perturbation error than those based on DWT (e.g., CWPA).According to Lemma4, the total error (MRE) is less than or equal to the sum of reconstruction error and perturbation error. Thus, if the reconstruction error or perturbation error is greater than 10%, there is a high probability that the final error will not be less than this threshold.As the reconstruction error of the DWT is greater than9%, there is a high probability that the final error of CWPA will not be less than this threshold, even if the Laplace noise decreases, i.e., when the number of homes N or the privacy budget ε increases. However, as the reconstruction error of the DFT is small (the median is between 2% and 3% when k=7,8,9), then the total error of the CFPA may be lower than that of the CWPA when the impact of Laplace noise decreases. For example, the median of the perturbation error of CFPA is between 3% and 5% when k=7,8,9, N=550, and ε=3. This analysis explains why, for ε=1, the CWPA obtains a better utility than the CFPA when the number of homes N is less than 250. For example, when N=50 and k=5, the median of the perturbation error (respectively, the reconstruction error) of CFPA is 70% (respectively, 10%) against 32% (respectively, 18%) for CWPA using Haar. Thus, the median of MRE of CFPA is between 70% and 80% against 32% and 50% for CWPA.WhenN is higher than 250, CFPA gets a better utility than CWPA. For example, when N=450 and k=5, the median of the perturbation error (respectively, the reconstruction error) of CFPA is 8.5% (respectively, 6.5%) against 4.5% (respectively, 16%) for CWPA using Haar. Thus, the median of MRE of CFPA is between 8.5% and 15% against 16% and 20.5% for CWPA.In this use case, by comparing the different techniques for publishing time-series consumption, it appears that clamping perturbation algorithms (CFPA, CWPA) get a better utility than unbounded algorithms (FPA, WPA), which shows that the clamping mechanism reduces the total error. Furthermore, when the number of homes is greater than 250, CFPA obtains the best utility, with a mean relative error of less than10% when ε=3. When the budget of privacy ϵ=1, the mean relative error of CFPA is less than 10% for N=450 homes.The CWPA gets the best utility when the number of homesN is smaller than 150 and the budget of privacy ε is 1. This is justified by its low perturbation error.Table3 summarizes the publishing algorithm with the smallest MRE according to the budget of privacy ε and the number of homes in the district (N). Based on the dataset from the Irish Commission for Energy Regulation (CER) [39], Table 3 shows that the clamping Fourier perturbation algorithm (CFPA) achieves a lower MRE than the clamping wavelet perturbation algorithm (CWPA) for N>150. Hence, CFPA gets a better utility than CWPA for N>150.Table 3
Publishing algorithm with the smallest MRE according to the budget of privacyε and the number of homes in the district (N).
Number of homes (N)Budget of privacy,εBest algorithmCoefficients,kMedian of MRE (%)501CWPA-Haar, db25353CWPA-Haar211501CWPA-Haar5193CFPA112501CFPA516383501CFPA5123874501CFPA5103865501CFPA59385


## 7. Conclusion

The large deployment of smart meters provides users and suppliers with the capacity to optimize the energy consumption through forecasting and demand-response services. This paper proposes an original and efficient approach to mitigate privacy leakages of users’ consumptions. This approach uses differential privacy and time-series transformations for supporting high privacy guarantees and utility. The clamping Fourier perturbation algorithm (CFPA) we propose achieves an error 6 times lower than the Fourier perturbation algorithm (FPA). Similarly, the clamping wavelet perturbation algorithm (CWPA) achieves an error 2 times lower than the wavelet perturbation algorithm (WPA). Thanks to our algorithm, the publication of aggregate time-series consumptions is now possible while guaranteeing that the aggregate does not reveal any individual consumptions and while achieving better utility than existing algorithms. These privacy-preserving aggregate time-series consumptions can then be used as a building block, enabling services such as forecasting and demand-response, which are suitable for improving the efficiency and reliability of the electric grid.In the future, we plan to investigate how to decentralize our clamping transform perturbation algorithm in order to resist to malicious aggregators. We plan to examine how to combine secure multiparty computation (SMC) with differential privacy (DP). SMC enables parties to compute a joint function without learning any individual inputs. SMC combined with DP could allow homes to compute and publish their aggregated consumptions without relying on an aggregator. However, SMC incurs a communication cost, which might have an impact on the running time performance.



---

*Source: 6643566-2021-03-26.xml*