A Privacy-preserving Disaggregation Algorithm for Non-intrusive Management of Flexible Energy
Paulin Jacquot, Olivier Beaude, Pascal Benchimol, Stéphane Gaubert, Nadia Oudjane
AA Privacy-preserving Disaggregation Algorithmfor Non-intrusive Management of Flexible Energy
Paulin Jacquot, Olivier Beaude, Pascal Benchimol, Stéphane Gaubert, Nadia Oudjane
Abstract — We consider a resource allocation problem involv-ing a large number of agents with individual constraints subjectto privacy, and a central operator whose objective is to opti-mizing a global, possibly non-convex, cost while satisfying theagents’ constraints. We focus on the practical case of the man-agement of energy consumption flexibilities by the operator of amicrogrid. This paper provides a privacy-preserving algorithmthat does compute the optimal allocation of resources, avoidingeach agent to reveal her private information (constraints andindividual solution profile) neither to the central operatornor to a third party. Our method relies on an aggregationprocedure: we maintain a global allocation of resources, andgradually disaggregate this allocation to enforce the satisfactionof private contraints, by a protocol involving the generation ofpolyhedral cuts and secure multiparty computations (SMC). Toobtain these cuts, we use an alternate projections method à laVon Neumann, which is implemented locally by each agent,preserving her privacy needs. Our theoretical and numericalresults show that the method scales well as the number of agentsgets large, and thus can be used to solve the allocation problemin high dimension, while addressing privacy issues.
I. INTRODUCTION
Motivation.
Consider an operator of an electricity mi-crogrid optimizing the joint production schedules of renew-able and thermal power plants in order to satisfy, at eachtime period, the consumption constraints of its consumers.To optimize the costs and the renewables integration, thisoperator relies on demand response techniques, that is, takingadvantage of the flexibilities of some of the consumerselectric appliances—those which can be controlled withoutimpacting the consumer’s confort, as electric vehicles orwater heaters [1]. However, for privacy reasons, consumersare not willing to provide neither their consumption con-straints nor their consumption profiles to a central operatoror any third party, as this information could be used to induceprivate information such as their presence at home.The global problem of the operator is to find an allocationof power (aggregate consumption) p = ( p t ) t at each timeperiod ( resource ) t ∈ T , such that p ∈ P (feasibilityconstraints of power allocation, induced by the power plantsconstraints). Besides, this aggregate allocation has to matchan individual comsumption profile x n = ( x n,t ) t ∈T for eachof the consumer (agent) n ∈ N considered. The problem can P. Jacquot, O. Beaude, P. Benchimol and N. Oudjane are with EDF R&D,OSIRIS, Palaiseau, France.P. Jacquot and S. Gaubert are with Inria Saclay and CMAP, Ecolepolytechnique, Palaiseau, France. [email protected] be written as follows: min x ∈ R N × T , p ∈P f ( p ) (1a) x n ∈ X n , ∀ n ∈ N (1b) (cid:88) n ∈N x n,t = p t , ∀ t ∈ T , (1c)The (aggregate) allocation p can be made public , that is,revealed to all agents. However, the individual constraint set X n and individual profiles x n constitute private informationof agent n , and should not be reavealed to the operator orany third party.The approach adopted in our paper is to deal with theproblem (1) as two kinds of interdependent subproblems.The firsts are optimal resource allocation problems, or mas-ter problems , min p ∈P ( s ) f ( p ) which consists in finding an aggregate allocation over T resources ( P ( s ) ⊂ P ⊂ R T ) .The second kind is problems of finding the disaggregation of a given aggregate allocation p , that is, to find an individualprofile x n for each agent (consumer) n satisfying her individ-ual constraint (1b), such that the aggregate of the profiles isthe optimal allocation (1c) determined in a master problem.Aside from the example above, ressource allocation prob-lems (optimizing common resources shared by multipleagents) find many applications in energy [1, 2], logistics[3], distributed computing [4], health care [5] and telecom-munications [6]. In these applications, several entities oragents (e.g. consumers, stores, tasks) share a common re-source (energy, products, CPU time, broadband) which hasa global cost for the system. For large systems composedof multiple agents, the dimension of the overall problemcan be prohibitive and one can rely on decomposition anddistributed approaches [7–9] to answer to this issue. Besides,agents’ individual constraints are often subject to privacyissues [10]. These considerations have paved the way to thedevelopment of privacy-preserving, or non-intrusive methodsand algorithms, e.g. [11, 12].In this work, we consider that each agent has a globaldemand constraint (e.g. energy demand or product quantity),which confers to the disaggregation problem the particularstructure of a transportation polytope [13]: the sum over theagents is fixed by the aggregate solution p , while the sumover the T resources are fixed by the agent global demandconstraint. Besides, individual constraints can also includeminimal and maximal levels on each resource, as for instanceelectricity consumers require, through their appliances, aminimal and maximal power at each time period. a r X i v : . [ c s . M A ] M a r ain Results. The major contribution of the paper isto provide a non-intrusive and distributed algorithm (Algo-rithm 4) that computes an aggregated resource allocation p ,optimal solution of the—possibly nonconvex—optimizationproblem (1), along with feasible individual profiles x foragents, without revealing the individual constraints of eachagent to a third party, either another agent or a centraloperator. The algorithm solves iteratively intances of masterproblem min p ∈P ( s ) f ( p ) by constructing successive approx-imations P ( s ) ⊂ P of the aggregate feasible set of (1) forwhich a disaggregation exists, by adding a new constraint on p to P ( s ) , before solving the next master problem.To identify whether or not disaggregation is feasible andto add a new constraint in the latter case, our algorithmrelies on the alternating projections method (APM) [14, 15]for finding a point in the intersection of convex sets. Here,we consider the two following sets: on the one hand, theaffine space defined by the aggregation to a given resourceprofile, and on the other hand, the set defined by all agentsindividual constraints (demands and bounds). As the latteris defined as a Cartesian product of each agent’s feasibilityset, APM can operate in a distributed fashion. The sequenceconstructed by the APM converges to a single point if theintersection of the convex sets is nonempty, and it convergesto a periodic orbit of length otherwise. Our key resultis the following: if the APM converges to a periodic orbit,meaning that the disaggregation is not feasible, we constructfrom this orbit a polyhedral cut , i.e. a linear inequalitysatisfied by all feasible solutions p of the global problem (1),but violated from the current resource allocation (Thm. 4).Adding this cut to the master problem, we can recomputea new resource allocation and repeat this procedure untildisaggregation is possible. Another major result stated inthis paper is the explicit upper bound on the convergencespeed of APM in our framework (Thm. 2), which is obtainedby spectral graph theory methods, exploiting also geometricproperties of transportation polytopes. This explicit speedshows a linear impact of the number of agents, which is astrong argument for the applicability of the method in largedistributed systems. Related Work.
A standard approach to solve resourceallocation problems in a distributed way is to use a La-grangian (dual) decomposition technique [8, 16, 17]. Thosetechniques are generally used to decompose a large problemsinto several subproblems of small dimension. They mayalso be implemented in a way which preserve privacy (seeRemark 2 in Sec. IV). However, Lagrangian decompositionmethods are based on strong duality property, requiringglobal convexity hypothesis which are not satisfied in manypractical problems (e.g. MILP, see Sec. V). On the contrary,our method can be used when the master allocation problemis not convex. In [2], the authors study a disaggregationproblem similar to the one considered in this paper. Theirresults concern zonotopic sets, which is different from thestructure we described in Sec. II. The APM has been thesubject of several works in itself [15, 18, 19]. The authors of[20] provide general results on the convergence rate of APM for semi-algebraic sets. They show that the convergence isgeometric for polyhedra. However, it is generally hard tocompute explicitly the geometric convergence rate of APM,as this requires to bound the singular values of certainmatrices arising from the polyhedral constraints. In [21],the authors provide an explicit convergence rate for APMon a class of polyhedra arising in submodular optimization.The sets they consider differ from the present transportationpolytopes.
Structure.
In Sec. II, we describe the master resourceallocation problem and formulate the associated disaggre-gation problem. In Sec. III, we focus on the APM andstate our main results. In Sec. IV, we apply these resultsto describe a non-intrusive version of APM (NI-APM) thatis used to describe our non-intrusive algorithm for computingan optimal resource allocation. Finally, in Sec. V, we providea concrete numerical example based on a MILP to model themanagement of a local electricity system (microgrid), andstudy numerically the influence of the number of agents onthe time needed for convergence of our algorithm.
Notation.
Vectors and matrices are denoted by bold fonts, v (cid:62) denotes the transpose of v , K denotes the vector (1 . . . (cid:62) of size K , U ([ a, b ]) stands for the uniform dis-tribution on [ a, b ] . We use (cid:107) x (cid:107) to denote the Frobeniusnorm (cid:107) x (cid:107) = (cid:80) n,t x n,t , and P C ( . ) to denote the Euclideanprojection on a convex set C .II. M ASTER PROBLEM AND DISAGGREGATIONSTRUCTURE
In this work, we suppose an operator wishes to determinean allocation of resources, represented by a T -dimensionalvector p , in order to minimize a global cost function f ,for instance, an electricity power economic dispatch (or theallocation of different types of merchandise in warehouses inlogistics applications) subject to a set of constraints describedby a feasibility set P . This problem can be nonconvex eitherbecause of nonconvex costs f or because of a nonconvexfeasible set P (see Sec. V). In the proposed method, theoperator will consider master problems of the form: min p ∈ R T f ( p ) (2a)s.t. p ∈ P ( s ) , (2b)where the set P ( s ) ⊂ P is an aggregate approximation ofdisaggregation constraints. Indeed, the resource allocation p has to be shared between N agents (e.g. consumers). Eachagent has a global demand (total energy needed) E n andsome lower and upper bounds on each of the resource t ∈ T .The admissible set of profiles of agent n is therefore: X n def = { x n ∈ R T | x (cid:62) n T = E n and ∀ t, x n,t ≤ x n,t ≤ x n,t } . (3)The disaggregation problem consists in finding individualprofiles x = ( x n,t ) n,t ∈ R NT of a given aggregatedallocation p such that x n is feasible for each agent n :F IND x ∈ Y p ∩ X (4)where Y p def = { y ∈ R NT | y (cid:62) N = p } and X def = (cid:89) n ∈N X n . ollowing (4), the disaggregated profile refers to x , whilethe aggregated profile refers to the allocation p .Problem (4) may not always be feasible. Some neces-sary conditions for a disaggreagation to exist, obtained bysumming the individual constraints on N , are the following aggregated constraints: p (cid:62) T = E (cid:62) N and x (cid:62) N ≤ p ≤ x (cid:62) N . (5)However, (II) are not sufficient conditions, as shown in Fig. 1where the problem (4) is represented as a flow or circulation problem from source nodes t ∈ T to sink nodes n ∈ N .Indeed, with this circulation representation of the disaggre-gation problem (4), an immediate consequence of Hoffmanntheorem [22, Thm. 3.18][23] is the following characterizationof the disaggregation feasibility, which involves an exponen-tial number of inequalities: Theorem 1.
The disaggregation problem (4) is feasible (i.e.
X ∩ Y p (cid:54) = ∅ ) iff for any T in ⊂ T , N in ⊂ N : (cid:88) t/ ∈T in p t ≤ (cid:88) t/ ∈T in ,n ∈N in x n,t − (cid:88) t ∈T in ,n/ ∈N in x n,t + (cid:88) n/ ∈N in E n . (6)The inequality (6) has a simple interpretation: the residualdemand (the left hand side composed of demand and exportsminus production) in T in ∪ N in cannot exceed the importcapacity (right hand side of the inequality). One can see that,in the example of Fig. 1, inequality (6) does not hold whenusing the cut composed of the dashed nodes p and E . p = 0 p = 3 E = 2 E = 0 . E = 0 . Fig. 1. Example of disaggregation structure ( T = 2 , N = 3) , with x = and x := . Although the aggregate constraints (II) are satisfied,the disaggregation (4) of p is not feasible in this example (see Thm. 1). There are two main reasons for which solving (1) is harderthan solving (2) and (4) separately:i) the dimension of (1) can be huge, as the number ofagents N can be really important, for instance in the exampleof individual consumers;ii) also, and this is the main motivation of this work, theinformation related to ( x n ) n , ( x n ) n and ( E n ) n might not beavailable to the centralized operator in charge of optimizingresources p , as this information may be confidential and keptby each agent n , not willing to reveal it to any third party.In the next sections, we provide a method that addressesthose two issues, by considering subproblems (2) and (4) in-dependently and iteratively, and exploiting the decomposablestructure of problem (4).III. A LTERNATE P ROJECTION M ETHOD (APM)
A. Convergence of APM on Transportation Polytopes
In this section, we consider a fixed aggregated profile p and present the Von Neumann Alternate Projections Method (APM) [14] which solves the problem Eq. (4) of finding apoint in the intersection X ∩ Y p . In the remaining, we willoften ommit p and just write Y to denote Y p . The key ideaof the method proposed in this paper is to use results ofAPM to generate a cut in the form of (6) and to add it as anew constraint in the master problem (2) to “improve” theaggregated profile p for the next iteration. As described inAlgorithm 1, APM can be used to decompose (4) and onlyinvolves local operations. Algorithm 1
Alternate Projections Method (APM)
Require:
Start with y (0) , k = 0 , ε cvg , a norm (cid:107) . (cid:107) on R NT repeat x ( k +1) ← P X ( y ( k ) ) y ( k +1) ← P Y ( x ( k +1) ) k ← k + 1 until (cid:13)(cid:13) y ( k ) − y ( k − (cid:13)(cid:13) < ε cvg The convergence of Algorithm 1 is proved by Thm. 2:
Theorem 2 ([15]) . Let X and Y be two convex sets with X bounded, and let ( x ( k ) ) k and ( y ( k ) ) k be the two infinitesequences generated by Algorithm 1 with ε cvg = 0 . Thenthere exists x ∞ ∈ X and y ∞ ∈ Y such that: x ( k ) −→ k →∞ x ∞ , y ( k ) −→ k →∞ y ∞ ; (7a) (cid:107) x ∞ − y ∞ (cid:107) = min x ∈X , y ∈Y (cid:107) x − y (cid:107) . (7b) In particular, if
X ∩ Y (cid:54) = ∅ , then ( x ( k ) ) k and ( y ( k ) ) k converge to a same point x ∞ ∈ X ∩ Y . If disaggregation is not feasible, Thm. 2 states that APMwill “converge” to an orbit ( x ∞ , y ∞ ) of period 2.The convergence rate of APM has been the subject ofseveral works [18, 20], and it strongly depends on thestructure of the sets on which the projections are done: forinstance, if the sets are polyhedral, [20, Prop. 4.2] shows thatthe convergence is geometric. However, there are very fewcases in which an explicit upper bound on the convergencerate has been proved. In our case, we are able to obtain sucha bound, as shown in the following theorem: Theorem 3.
For the sets X and Y defined in (3-4), thetwo subsequences of alternate projections converge at ageometric rate to x ∞ ∈ X , y ∞ ∈ Y , with: (cid:107) x ( k ) − x ∞ (cid:107) ≤ (cid:107) x (0) − x ∞ (cid:107) × ρ kNT where ρ NT def = 1 − (cid:0) N ( T + 1) ( T − (cid:1) < , Same inequalities hold for the convergence of y ( k ) to y ∞ .Proof. Appendix II provides a sketch of the proof.Thm. 3 shows that the APM is efficient in our case ofbounded transport polytopes. It shows that the number ofiterations for a given accuracy grows linearly in the numberof agents N .As stated in (4), the set X is a Cartesian product (cid:81) n X n ,so that the projection (13) can be computed by N projectionsn ( X n ) n , which can be executed in parallel. Now, insteadof solving the quadratic program by standard interior pointmethods and due to its particular structure, we can use thealgorithm of Brucker [24], which has a complexity in O ( T ) .On the other hand, P Y ( . ) is a projection on an affine space,and the solution can be obtained explicitly as: ∀ n, t, y n,t = x n,t + ν t and ν = N ( p − x (cid:62) N ) . (8) B. Generation of a cut from APM iterates
Our key result is the following: in the case where APMconverges to a periodic orbit ( x ∞ , y ∞ ) with x ∞ (cid:54) = y ∞ (seeThm. 2), we obtain from ( x ∞ , y ∞ ) an inequality (6) that isviolated by p : Theorem 4.
For the sets X and Y defined in (3-4) and if X ∩ Y = ∅ , the following sets given by the limit orbit ( x ∞ , y ∞ ) defined in Thm. 2: T = { t | p t > (cid:80) n ∈N x ∞ n,t } (9a) N = { n | E n − (cid:80) t/ ∈T x n,t − (cid:80) t ∈T x n,t < } (9b) define a Hoffman cut of form (6) violated by p , that is: (cid:88) n ∈N E n − (cid:88) t ∈T p t + (cid:88) t ∈T ,n/ ∈N x n,t − (cid:88) t/ ∈T ,n ∈N x n,t < . (10) This cut can be reformulated in terms of (cid:62) N x ∞ as: A T < (cid:88) t ∈T p t with A T def = (cid:88) t ∈T (cid:88) n ∈N x ∞ n,t . (11) Proof.
Appendix I gives the sketch of the proof of Thm. 4.The complete proofs will be given elsewhere.One can see that, intuitively, N is the subset associatedto T that minimizes the right hand side of (6). Note thatThm. 4 gives an alternative constructive proof of Hoffmancirculation’s theorem (Thm. 1) in the case of a bipartitegraph of the form of Fig. 1. Moreover, in the case where thedisaggregation problem (4) is not feasible, the negation ofequation (11) provides a new valid constraint as a conditionfor the existence of a disaggreagated profile of p . Thisconstraint can be used in the master problem (2) to update thevector of resources p for the next iteration. This constraintonly involves the aggregate information (cid:62) N x ∞ on the usersprofile. To make the process fully non-intrusive , we explainin Sec. IV-A how the operator can compute this constraintwithout making the agents reveal their profiles ( x ∞ n ) n ∈N .IV. N ON - INTRUSIVE P ROJECTIONS AND C OMPUTATIONOF D ISAGGREGATED O PTIMAL R ESOURCES
A. Non-Intrusive Alternate Projections Method (NI-APM)
Because of the particular structure of the problem, theprojections in APM can be computed separately by theoperator and the agents. The projection P Y is made by theoperator, which only requires to know p and the aggregateprofile x (cid:62) N according to from (8). The projection P X on X = (cid:81) n X n is executed in parallel by each agent: n computes P X n which only needs her private information E n and x n , x n . However, in the way APM is described in Algorithm 1, the operator and the agents still need toexchange the iterates x ( k ) , y ( k ) at each step. To avoid thetransmission of agents’ profiles to the operator, we use asecure multiparty computation (SMC) technique (see [25])which enables the operator to obtain the aggregate profile S ( k ) := (cid:62) N x ( k ) in a non-intrusive manner, as described inAlgorithm 2.The main idea of SMC is that, instead of sending herprofile x n , agent n splits x n,t for each t into N randomparts ( s n,t,m ) m , according to an uniform distribution andsumming to x n,t (Lines 2-3). Thus, each part s n,t,m takenindividually does not reveal any information on x n nor on X n , and can be sent to agent m . Once all exchanges ofparts are completed (Line 5), and n has herself received theparts from other agents, agent n computes a new aggregatequantity σ n (Line 7), which does not contain either anyinformation about any of the agents, and sends it to theoperator (Line 8). The operator can finally compute thequantity S = x (cid:62) N = σ (cid:62) N . Algorithm 2
SMC of Aggregate (SMCA) (cid:80) n ∈N x n Require:
Each agent has a profile ( x n ) n ∈N for each agent n ∈ N do Draw ∀ t, ( s n,t,m ) N − m =1 ∈U ([0 , A ] N − ) and set ∀ t, s n,t,N def = x n,t − (cid:80) N − m =1 s n,t,m Send ( s n,t,m ) t ∈T to agent m ∈ N done for each agent n ∈ N do Compute ∀ t, σ n,t = (cid:80) m ∈N s m,t,n Send ( σ n,t ) t ∈T to operator done Operator computes S = (cid:80) n ∈N σ n Remark 1. As σ n , and s n are random by construction, aneavesdropper aiming to learn the profile x n of n has nochoice but to intercept all the communications of n to allother agents (to learn ( s n,t,m ) m (cid:54) = n and ( s m,t,n ) m (cid:54) = n ) andto the operator (to learn σ n ). To increase the confidentialityof the procedure, one could use any encryption scheme (suchas RSA [26]) for all communications involved in Algorithm 2. We can use this non-intrusive computation of aggregate S in APM to obtain a non-intrusive algorithm NI-APM (Algorithm 3) in which agents do not reveal neither theirprofiles nor their constraints to the operator.One can see that x and y computed in Lines 3 and 8in Algorithm 3 correspond to the projections computed inthe original APM Algorithm 1. In Algorithm 3, the operatorobtains the aggregate profile S ( k ) (Line 5), computes andsends the corrections ν ( k ) to all agents (Line 6). Then, eachagent can compute locally the projection y ( k ) n = P Y ( x ( k ) n ) by applying the correction ν ( k ) (Line 8).Using (8), we get ν ( k ) → ν ∞ def = N ( p − (cid:62) N x ∞ ) . Thm. 4uses this limit value through T ∞ def = { t ∈ T | < ν ∞ t } .Yet, from APM, one can only access to ν ( k ) and thus to lgorithm 3 Non-intrusive APM (NI-APM)
Require:
Start with y (0) , k =0 , ε cvg , ε dis , norm (cid:107) . (cid:107) on R NT repeat for each agent n ∈ N do x ( k ) n ← P X n ( y ( k − n ) done Operator obtains S ( k ) ← SMCA( x ( k ) ) ( cf Algo.2) and sends ν ( k ) := N ( p − S ( k ) ) ∈ R T to agents N for each agent n ∈ N do Compute y ( k ) n ← x ( k ) n + ν ( k ) done k ← k + 1 until (cid:13)(cid:13) x ( k ) − x ( k − (cid:13)(cid:13) < ε cvg if (cid:13)(cid:13) x ( k ) − y ( k ) (cid:13)(cid:13) ≤ ε dis then (cid:46) found a ε dis -solution of the disaggregation problem Each agent adopts profile x ( k ) n return D ISAG ← T RUE15: else (cid:46) have to find a valid constraint violated by p Operator computes T ← { t ∈ T | Bε cvg < ν ( k ) t } Operator computes A T def = (cid:80) t ∈T (cid:62) N σ ( k ) t if A T − (cid:80) t ∈T p t < then return D ISAG ← F ALSE , A T else (cid:46) need to run APM with higher precision Return to Line 1 with ε cvg ← ε cvg / end end the approximation T , computed on Line 16), where B is apre-defined constant. However, we show that for ε cvg smallenough and a well-chosen value of B , we obtain T = T ∞ ,so that we get the termination result: Proposition 1.
For
B > (1 − ρ NT ) − , Algorithm 3 termi-nates in finite time. The termination of the loop Lines 1-11 is ensured byThm. 3. In the case where (cid:13)(cid:13) x ( k ) − y ( k ) (cid:13)(cid:13) ≤ ε dis , Algorithm 3terminates. Otherwise, if (cid:107) x ∞ − y ∞ (cid:107) > ε dis , then Algo-rithm 3 terminates (i.e. Line 18 is True and a new cut isfound) as soon as Bε cvg < min (cid:110) (cid:107) x ∞ − y ∞ (cid:107) √ N , ν (cid:111) , where ν def = min {| ν ∞ t | > } and with (cid:107) . (cid:107) = (cid:107) . (cid:107) . The completeproof is ommited here.In practice, we can start with a large ε cvg to obtain the firstconstraints while avoiding useless computation, and then half ε cvg if needed (Line 21) until the termination condition holds. Remark 2.
SMC techniques could also be used to implementnon-intrusive Lagrangian decomposition methods. However,these methods rely on a convexity hypothesis that we do notneed in the proposed method.B. Non-intrusive Disaggregation of Optimal Allocation
In this section, we describe a method to compute a solutionof the global problem (1), that is, an optimal resourceallocation p for which a disaggregation exists, along withan associated disaggregated profile x n for each agent n . This computation is done in a non-intrusive manner: theoperator in charge of p does not have access neither to thebounding constraints x and x of the agents nor to the agentsdisaggregated profile x , as detailed in Algorithm 4 below. Algorithm 4
Non-intrusive Optimal Disaggregation
Require: s = 0 , P (0) = P ; D ISAG = F
ALSE1: while
Not D
ISAG do Compute p ( s ) = arg min p ∈P ( s ) cs p D ISAG ← NI-APM( p ( s ) ) if D ISAG then Operator adopts p ( s ) else Obtain T ( s )0 , A ( s ) T from NI-APM( p ( s ) ) P ( s +1) ← P ( s ) ∩ { p | (cid:80) t ∈T ( s )0 p t ≤ A ( s ) T } end s ← s + 1 done Algorithm 4 iteratively calls NI-APM (Algorithm 3) and incase disaggregation is not possible (Line 6), a new constraintis added (Line 8), obtained from the quantity A T definedin (11), to the resource problem (2). This constraint isan inequality on p and thus does not reveal significantindividual information to the operator. The algorithm stopswhen disaggregation is possible (Line 4). The termination ofAlgorithm 4 is ensured by the following property and theform of the constraints added (10): Proposition 2.
Algorithm 4 stops after a finite number ofiterations, as at most T constraints (Line 8) can be addedto the master problem (Line 2). Although there exist some instances with an exponentialnumber of independent constraints, this does not jeopardizethe proposed method: in practice, the algorithm stops after avery small number of constraints added (see the example ofSec. V). Intuitively, we will only add constraints “support-ing” the optimal allocation p .Thus, Algorithm 4 is a method which enables the operatorto compute a resource allocation p and the N agents to adoptprofiles ( x n ) n , such that ( x , p ) solves the global problem(1), and the method ensures that both:1) the information relative to each agent constraints (upperbounds x n , lower bounds x n , demand E n );2) the final disaggregated profile x n (as well as the iterates ( x ( k ) ) k and ( y ( k ) ) k in NI-APM)are kept confidential by agent n and can not be induced bya third party (either the operator or any other agent m (cid:54) = n ).V. A PPLICATION TO M ANAGEMENT OF A M ICROGRID
We apply the proposed method to solve a nonconvexdistributed problem in the energy field. We consider a micro-grid [27] composed of N electricity consumers with flexibleappliances (such as electric vehicles or water heaters), aphotovoltaic (PV) power plant and a conventional generator. . Mixed Integer Problem Formulation
The operator responsible of the microgrid aims at satisfy-ing the demand constraints of consumers over a set of timeperiods T = { , . . . , T } , while minimizing the energy costfor the community. We have the following characteristics: • the PV plant generates a nondispatchable power profile ( p PV t ) t ∈T at marginal cost zero; • the conventional generator has a starting cost C ST , min-imal and maximal power production p g , p g , and piecewise-linear and continuous generation cost function p g (cid:55)→ f ( p g ) : f ( p g ) = α k + c k p g , if p g ∈ I k def = [ θ k − , θ k [ , k = 1 . . . K, where θ = 0 and θ K def = p g ; • each agent n ∈ N has some flexible appliances whichrequire a global energy demand E n on T , and has consump-tion constraints on the total household consumption, on eachtime period t ∈ T , that are formulated with x n , x n . Theseparameters are confidential because they could for instancecontain some information on agent n habits.The master problem (2) can be written as a MILP (12): min p , p g , ( p gk ) , ( b k ) , b ON , b ST (cid:88) t ∈T (cid:16) α b ON t + (cid:88) k c k p gkt + C ST b ST t (cid:17) (12a) p gt = (cid:80) Kk =1 p gk,t , ∀ t ∈ T (12b) b ,t θ ≤ p g ,t ≤ θ , ∀ t ∈ T (12c) b ,t ( θ − θ ) ≤ p g ,t ≤ b ,t ( θ − θ ) , ∀ t ∈ T (12d)... ... (12e) ≤ p gK,t ≤ b K − ,t ( θ K − θ K − ) , ∀ t ∈ T (12f) b ST t ≥ b ON t − b ON t − , ∀ t ∈ { , . . . , T } (12g) p g b ON t ≤ p gt ≤ p g b ON t , ∀ t ∈ T (12h) b ON t , b ST t , b ,t , . . . , b K − ,t ∈ { , } , ∀ t ∈ T (12i) p ≤ p PV + p g (12j) p (cid:62) T = E (cid:62) N (12k) x (cid:62) N ≤ p ≤ x (cid:62) N . (12l)In this formulation (12b-12f) are a mixed integer formu-lation of the generation cost function f : one can show thatthe boolean variable b k,t is equal to one iff p gk,t ≥ θ k foreach k ∈ { , . . . , K − } . Note that only α appears in (12a)because of the continuity assumption on f .Constraints (12g-12h) ensure the on/off and starting con-straints of the power plant, (12j) ensures that the power allo-cated to consumption is not above the total production, and(12k-12l) are the aggregated feasibility conditions alreadyreferred to in (II). Note that more complex and realisticMILP models exist for power plants (e.g. [28]), but with thesame structure than (12). The nonconvexity of (12) comesfrom the existence of starting costs and on minimal thepower constraint, which makes necessary to use boolean statevariables b ST , b ON . B. Parameters
We simulate the problem described above for differentvalues of N ∈ { , , , , } and one hundred instances with random parameters for each value of N . A scalingfactor κ N def = N/ is applied on parameters to ensurethat production capacity is large enough to meet consumersdemand. The parameters are chosen as follows: • T = 24 (hours of a day); • production costs: K = 3 , θ = [0 , , , κ N , c =[0 . , . , . , p g =50 κ N , p g =300 κ N , α =4 and C ST = 15 ; • photovoltaic: p PV t = (cid:104) − cos( ( t − π )+ U ([0 , (cid:105) κ N for t ∈ { , . . . , } , p PV t = 0 otherwise (see Fig. 3); • for consumption parameters, we used x n,t ∼ U ([0 , , x n,t ∼ U ([0 , x n,t and E n ∼ U ([ (cid:62) T x n , (cid:62) T x n ]) , so thatindividual feasibility ( X n (cid:54) = ∅ ) is ensured. N = 2 NUMBER OF SUBPROBLEMS SOLVED ( AVERAGE ON
INSTANCES )Fig. 2. Total number of computed projections for different values of N . We observe that the number of agents N has a sublinear impact on thetotal number of projections needed. C. A limited impact of the number N of agents We implement Algorithm 4 using Python 3.5. The MILP(12) is solved using Cplex Studio 12.6 and Pyomo inter-face. Simulations are run on a single core of a cluster at3GHz. For the convergence criteria (see Lines 11 and 12 ofAlgorithm 3), we use ε dis = 0 . with the operator normdefined by (cid:107)| x |(cid:107) = max n ∈N (cid:80) t | x n,t | (to avoid the √ N factor in the convergence criteria appearing with (cid:107) . (cid:107) ), andstarts with ε cvg = 0 . . The largest instances took around10 minutes to be solved in this configuration and withoutparallel implementation. As the CPU time needed depend onthe cluster load, it is not a reliable indicator of the influenceon N on the complexity of the problems. Moreover, oneadvantage of the proposed method is that the projections inAPM can be computed locally by each agent in parallel,which could not be implemented here for practical reasons.Tab. I gives two robust indicators of the influence of N on the problem complexity: the number of master problemssolved and the total number of projections computed, onaverage over the hundred instance for each value of N : • one observes that the number of master problems solved(MILP (12) ), which corresponds to the number of constraintsr “cuts” added to the master problem, remains almostconstant when N increases; • in all instances, this number is way below the upperbound of > , × possible constraints (see proof ofProp. 2), which suggests that only a polynomial number ofconstraints are added in practice; • the average total number of projections computed foreach instance (total number of iterations of the while loopof Algorithm 3, Line 1 over all calls of APM in the instance)increases in a sublinear way, as illustrated in Fig. 2, whichis even better that one could expect from the upper boundgiven in Thm. 3. Fig. 3. Example of optimal resource allocation (power production p = p PV + p g ) in the example of Sec. IV, with N = 20 agents. Exploiting consumption flexibilities, the consumption is higher during thePV production periods.
VI. C
ONCLUSION
We provided a non-intrusive algorithm that enables tocompute an optimal resource allocation, solution of a–possibly nonconvex–optimization problem, and affect to eachagent an individual profile satisfying a global demand andlower and upper bounds constraints. Our method uses localprojections and works in a distributed fashion. Hence, it en-sures that the problem is not affected by the high dimensionrelative to the large number of agents, and that it is privacy-preserving, as agents do not need to reveal any informationon their constraints or their individual profile to a third party.Several extensions and generalizations can be considered.First, we could generalize the abstract circulation problem onthe bipartite graph depicted in Fig. 1 to an arbitrary network,where the set of nodes is partitioned in K parts defining K sets on which we could make alternating projections. Second,our method relies on the particular structure obtained fromthe form of constraints. Although these kind of constraintsare widely used to model many practical situations, it wouldalso be useful to obtain similar results for arbitrary (orpolyhedral) agents constraints. Last, a deeper complexityanalysis, with a thinner upper bound on the maximal numberof constraints (cuts) added in the algorithm (see Prop. 2 andTab. I) would constitute interesting results. R EFERENCES[1] P. Jacquot, O. Beaude, S. Gaubert, and N. Oudjane, “Analysis andimplementation of an hourly billing mechanism for demand responsemanagement,”
IEEE Trans. Smart Grid , 2018.[2] F. L. Müller, J. Szabó, O. Sundström, and J. Lygeros, “Aggregationand disaggregation of energetic flexibility from distributed energyresources,”
IEEE Trans. Smart Grid , 2017.[3] K. K. Lai, K. Lam, and W. K. Chan, “Shipping container logistics andallocation,”
J. Oper. Res. Soc. , vol. 46, no. 6, pp. 687–697, 1995.[4] P.-Y. R. Ma et al. , “A task allocation model for distributed computingsystems,”
IEEE Trans. Computers , vol. 100, no. 1, pp. 41–47, 1982.[5] A. Rais and A. Viana, “Operations research in healthcare: a survey,”
Int. Trans. Oper. Res. , vol. 18, no. 1, pp. 1–31, 2011.[6] M. Zulhasnine, C. Huang, and A. Srinivasan, “Efficient resource allo-cation for device-to-device communication underlaying lte network,”in
WiMob, 2010 IEEE 6th Int. Conference . IEEE, 2010, pp. 368–375.[7] D. P. Bertsekas and J. N. Tsitsiklis,
Parallel and distributed computa-tion: numerical methods . Prentice hall Englewood Cliffs, NJ, 1989,vol. 23.[8] D. P. Palomar and M. Chiang, “A tutorial on decomposition methodsfor network utility maximization,”
IEEE J. Sel. Areas Commun. ,vol. 24, no. 8, pp. 1439–1451, 2006.[9] L. Xiao and S. Boyd, “Optimal scaling of a gradient method fordistributed resource allocation,”
J. Optim. Theory. Appl. , vol. 129,no. 3, pp. 469–488, 2006.[10] B. A. Huberman, E. Adar, and L. R. Fine, “Valuating privacy,”
IEEEsecurity & privacy , vol. 3, no. 5, pp. 22–25, 2005.[11] A. Zoha, A. Gluhak, M. A. Imran, and S. Rajasegarar, “Non-intrusiveload monitoring approaches for disaggregated energy sensing: Asurvey,”
Sensors , vol. 12, no. 12, pp. 16 838–16 866, 2012.[12] G. Jagannathan, K. Pillaipakkamnatt, and R. N. Wright, “A newprivacy-preserving distributed k-clustering algorithm,” in
Proc. of the2006 SIAM Int. Conf. on Data Mining . SIAM, 2006, pp. 494–498.[13] E. D. Bolker, “Transportation polytopes,”
Journal of CombinatorialTheory, Series B , vol. 13, no. 3, pp. 251–262, 1972.[14] J. Von Neumann,
Functional operators: Measures and integrals .Princeton University Press, 1950, vol. 1.[15] L. Gubin, B. Polyak, and E. Raik, “The method of projections forfinding the common point of convex sets,”
USSR Comput. Math. &Math. Phys. , vol. 7, no. 6, pp. 1 – 24, 1967.[16] L. Xiao, M. Johansson, and S. P. Boyd, “Simultaneous routing andresource allocation via dual decomposition,”
IEEE Trans. Comm. ,vol. 52, no. 7, pp. 1136–1144, 2004.[17] K. Seong, M. Mohseni, and J. M. Cioffi, “Optimal resource allocationfor ofdma downlink systems,” in
Information Theory, 2006 IEEE Int.Sym.
IEEE, 2006, pp. 1394–1398.[18] H. H. Bauschke and J. M. Borwein, “On the convergence of vonneumann’s alternating projection algorithm for two sets,”
Set-ValuedAnalysis , vol. 1, no. 2, pp. 185–212, 1993.[19] H. H. Bauschke, J. Chen, and X. Wang, “A bregman projectionmethod for approximating fixed points of quasi-bregman nonexpansivemappings,”
Applicable Analysis , vol. 94, no. 1, pp. 75–84, 2015.[20] J. M. Borwein, G. Li, and L. Yao, “Analysis of the convergence rate forthe cyclic projection algorithm applied to basic semialgebraic convexsets,”
SIAM J. Optim. , vol. 24, no. 1, pp. 498–527, 2014.[21] R. Nishihara, S. Jegelka, and M. I. Jordan, “On the convergence rateof decomposable submodular function minimization,” in
Advances inNeural Information Processing Systems , 2014, pp. 640–648.[22] W. J. Cook, W. Cunningham, W. Pulleyblank, and A. Schrijver,
Combinatorial optimization . Springer, 2009.[23] A. J. Hoffman, “Some recent applications of the theory of linearinequalities to extremal combinatorial analysis,” in
Proc. of Symposiaon Applied Mathematics , 1960, pp. 113–127.[24] P. Brucker, “An O(n) algorithm for quadratic knapsack problems,”
Oper. Res. Lett. , vol. 3, no. 3, pp. 163–166, 1984.[25] A. C. Yao, “How to generate and exchange secrets,” in , Oct 1986, pp. 162–167.[26] R. L. Rivest, A. Shamir, and L. Adleman, “A method for obtainingdigital signatures and public-key cryptosystems,”
Communications ofthe ACM , vol. 21, no. 2, pp. 120–126, 1978.[27] F. Katiraei, R. Iravani, N. Hatziargyriou, and A. Dimeas, “Microgridsmanagement,”
IEEE power and energy magazine , vol. 6, no. 3, 2008.[28] M. Carrión and J. M. Arroyo, “A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem,”
IEEE Trans. Pow. Sys. , vol. 21, no. 3, pp. 1371–1378, 2006.
PPENDIX IP ROOF OF P ROP . 3To show Thm. 4, we formulate the projections P X and P Y as the solutions of the constrained quadratic programs: min x ∈ R NT (cid:107) x − y (cid:107) (13a) x T = E ( λ ) (13b) x ≤ x ≤ x ( µ , µ ) (13c)and: min y ∈ R NT (cid:107) x − y (cid:107) (14a) y (cid:62) N = p ( ν ) , (14b)where λ , µ , µ , ν are the Lagrangian multipliers associatedto the constraints. Although there is no such explicit char-acterization of the solution of (13) as the one (8) given for(14), we can obtain the following properties: Proposition 3.
Suppose that
X ∩ Y = ∅ and consider thesets T and N given by (9) . Then we have the following:(i) ∀ t ∈ T , ∀ n / ∈ N , y ∞ n,t ≥ x n,t and x ∞ n,t = x n,t ;(ii) ∀ n ∈ N , λ n < ;(iii) ∀ t / ∈ T , ∀ n ∈ N , x ∞ n,t = x n,t ;(iv) the sets T , T c , N and N c are nonempty. The proof of Prop. 3 relies on the KKT optimality condi-tions associated to (13) and (14). Then we use Prop. 3 andthe convergence condition to prove Thm. 4: (cid:88) n ∈N E n + (cid:88) t ∈T ,n/ ∈N x n,t − (cid:88) t/ ∈T ,n ∈N x n,t − (cid:88) t ∈T p t = (cid:88) t ∈T (cid:32) (cid:88) n ∈N x ∞ n,t − (cid:88) n ∈N y ∞ n,t (cid:33) = (cid:88) t ∈T (cid:32) (cid:88) n ∈N − ν ∞ t (cid:33) = (cid:88) t ∈T (cid:32) − (cid:88) n ∈N (cid:12)(cid:12) x ∞ n,t − y ∞ n,t (cid:12)(cid:12)(cid:33) = − (cid:107) x ∞ − y ∞ (cid:107) < , where the last equality comes from (cid:80) t ∈T (cid:80) n ∈N ( x En,t − x P n,t ) = 0 . The compact form (11) also follow from Prop. 3: A T ( x ) = (cid:88) n ∈N E n + (cid:88) t ∈T ,n/ ∈N x n,t − (cid:88) t/ ∈T ,n ∈N x n,t = (cid:88) n ∈N ,t ∈T x ∞ n,t + (cid:88) t ∈T ,n/ ∈N x ∞ n,t = (cid:88) t ∈T (cid:88) n ∈N x ∞ n,t . A PPENDIX
IIP
ROOF OF T HM . 2For this analysis, we use the space R NT = R T ×· · ·× R T ,where the ( n − T +1 to nT coordinates correspond to agent n , for ≤ n ≤ N . We make use of the following results: Lemma 1 ([21]) . For APM on polyhedra X and Y , thesequences ( x ( k ) k and ( y ( k ) k converge at a geometric rate,where the rate is bounded by the maximal value of the squareof the cosine of the Friedrichs angle c F ( U, V ) between a face U of X and a face V of Y , where c F ( U, V ) is given by: c F ( U, V ) = sup { u T v | (cid:107) u (cid:107) ≤ , (cid:107) v (cid:107) ≤ u ∈ U ∩ ( U ∩ V ) ⊥ , v ∈ V ∩ ( U ∩ V ) ⊥ } . Lemma 2 ([21]) . Let A and B be matrices with or-thonormal rows and with equal numbers of columns and Λ sv ( AB (cid:62) ) the set of singular values of AB (cid:62) . Then if Λ sv ( AB (cid:62) ) = { } , then c F (Ker(A) , Ker(B)) = 0 . Other-wise, c F (Ker(A) , Ker(B)) = max λ< { λ ∈ Λ sv (AB (cid:62) ) } . In our case, the polyhedra Y is an affine subspace Y = { x ∈ R NT | A x = √ N − T } with A def = √ N − J ,N ⊗ I T ,where ⊗ denotes the Kronecker product. The matrix A hasorthonormal rows and the direction of Y is Ker ( A ) .Describing the faces X is more complex. We have apolyhedral description of X , and the faces of X are subsets ofthe collection of affine subspaces indexed by ( T n , T n ) n ⊂T N (with T ∩ T = ∅ ): A ( T n , T n ) n def = (cid:110) ( x ) nt |∀ n, x (cid:62) n T = E n and ∀ t ∈ T n , x n,t = x n,t , and ∀ t ∈ T n , x n,t = x n,t (cid:111) . The associated linear subspace is given by