A privacy-preserving method to optimize distributed resource allocation ∗

arXiv:1908.03080v3 [math.OC] 22 Jun 2020

Olivier Beaude † , Pascal Benchimol‡ , Stéphane Gaubert § , Paulin Jacquot¶, and Nadia
Oudjane†

Abstract. We consider a resource allocation problem involving a large number of agents with individual constraints subject to privacy, and a central operator whose objective is to optimize a global, possibly
nonconvex, cost while satisfying the agents’ constraints, for instance an energy operator in charge
of the management of energy consumption flexibilities of many individual consumers. We provide a
privacy-preserving algorithm that does compute the optimal allocation of resources, avoiding each
agent to reveal her private information (constraints and individual solution profile) neither to the
central operator nor to a third party. Our method relies on an aggregation procedure: we compute
iteratively a global allocation of resources, and gradually ensure existence of a disaggregation, that
is individual profiles satisfying agents’ private constraints, by a protocol involving the generation of
polyhedral cuts and secure multiparty computations (SMC). To obtain these cuts, we use an alternate projection method, which is implemented locally by each agent, preserving her privacy needs.
We address especially the case in which the local and global constraints define a transportation polytope. Then, we provide theoretical convergence estimates together with numerical results, showing
that the algorithm can be effectively used to solve the allocation problem in high dimension, while
addressing privacy issues.

1. Introduction.
1.1. Motivation. Consider an operator of an electricity microgrid optimizing the joint
production schedules of renewable and thermal power plants in order to satisfy, at each time
period, the consumption constraints of its consumers. To optimize power generation or market costs and the integration of renewable energies, this operator relies on demand response
techniques, taking advantage of the flexibilities of some of the consumers electric appliances—
those which can be controlled without impacting the consumer’s comfort, as electric vehicles
or water heaters [25]. However, for confidentiality reasons, consumers may not be willing to
provide their consumption constraints or their consumption profiles to a central operator or
to any third party, as this information could be used to infer private information such as their
presence at home.
The global problem of the operator is to find an allocation of power (aggregate flexible
consumption) p = (pt )t∈[T ] , where [T ] := {1, . . . , T } is the set of time periods, such that
p ∈ P (feasibility constraints of power allocation, induced by the power plants constraints),
while minimizing the cost f (p) (representing the production and operations costs of the centrally controlled power plants). Besides, this aggregate allocation has to match an individual
consumption profile xn = (xn,t )t∈[T ] for each of the consumer (agent) n ∈ [N ] considered. The
∗
Initial version submitted on August 2, 2019. Revision of June 24, 2020. Some of the present results have been
announced in the conference proceedings [24].
†
EDF R&D, OSIRIS, Palaiseau, France ([olivier.beaude,nadia.oudjane]@edf.fr)
‡
CMAP,
École
polytechnique,
CNRS
and
Triscale
Innov,
Palaiseau,
France
(pascal.benchimol@polytechnique.edu).
§
Inria and CMAP, École polytechnique, CNRS, Palaiseau, France (stephane.gaubert@inria.fr)
¶
EDF R&D, OSIRIS, Inria and CMAP, École polytechnique,
CNRS, Palaiseau,
France,
paulin.jacquot@polytechnique.edu

1

problem can be written as follows:
(1.1a)
(1.1b)
(1.1c)

min

x∈RN ×T, p∈P

f (p)

xn ∈ Xn , ∀n ∈ [N ]
X
xn,t = pt , ∀t ∈ [T ] ,
n∈[N ]

The (aggregate) allocation p can be made public, that is, revealed to all agents. However, the
individual constraint set Xn and individual profiles xn constitute private information of agent
n, and should not be revealed to the operator or any third party.
It will be helpful to think of Problem (1.1) as the combination of two interdependent
subproblems:
i) given an aggregate allocation p, the disaggregation problem consists in finding, for each agent
n, an individual profile xn satisfying her individual constraint (1.1b), so that constraint (1.1c)
is satisfied. This is equivalent to:
(1.2a)
(1.2b)

Find x ∈ Yp ∩ X
where Yp := {y ∈ RN T |y > 1N = p} and X :=

Y

Xn .

n∈[N ]

When (1.2) has a solution, we say that a disaggregation exists for p;
ii) For a given subset Q ⊂ P, we define the master problem,
(1.3)

min f (p) .
p∈Q

When Q is precisely the set of aggregate allocations for which a disaggregation exists, the
optimal solutions of the master problem correspond to the optimal solutions of (1.1).
Aside from the example above, resource allocation problems (optimizing common resources
shared by multiple agents) with the same structure as (1.1), find many applications in energy
[31, 25], logistics [28], distributed computing [30], health care [35] and telecommunications
[45]. In these applications, several entities or agents (e.g. consumers, stores, tasks) share
a common resource (energy, products, CPU time, broadband) which has a global cost for
the system. For large systems composed of multiple agents, the dimension of the overall
problem can be prohibitive. Hence, a solution is to rely on decomposition and distributed
approaches [10, 34, 40]. Besides, agents’ individual constraints are often subject to privacy
issues [23]. These considerations have paved the way to the development of privacy-preserving,
or non-intrusive methods and algorithms, e.g. [44, 26].
In this work, except in Section 4, we consider that each agent n ∈ [N ] has a global demand
constraint (e.g. energy demand or product quantity), which confers to the disaggregation
problem the particular structure of a transportation polytope [11]: the sum over the agents
is fixed by the aggregate solution p, while the sum over the T resources is fixed by the
agent global demand constraint. Besides, individual constraints can also include minimal and
maximal levels on each resource. For instance, electricity consumers require, through their
appliances, a minimal and maximal power at each time period.
2

1.2. Main Results. The main contribution of the paper is to provide a non-intrusive and
distributed algorithm (Algorithm 3.4) that computes an aggregate resource allocation p, optimal solution of the—possibly nonconvex—optimization problem (1.1), along with feasible
individual profiles x for agents, without revealing the individual constraints of each agent to
a third party, either another agent or a central operator. The algorithm solves iteratively
instances of master problems minp∈P (s) f (p), obtained by constructing a decreasing sequence
of successive approximations P (s) of the set of aggregate consumptions p for which a disaggregation exists (see (1.1)), starting from P (0) = P. At each step, we reduce the set P (s) by
incorporating a new constraint on p (a cutting plane), before solving the next master problem. We shall see that this cutting plane can be computed and added to the master problem
without revealing any individual information on the agents.
More precisely, to identify whether or not disaggregation (1.2) is feasible and to add a
new constraint in the latter case, our algorithm relies on the alternate projections method
(APM) [39, 20] for finding a point in the intersection of convex sets. Here, we consider the
two following sets: on the one hand, the affine space of profiles x ∈ RN T aggregating to a
given resource allocation p, and on the other hand, the set defined by all agents individual
constraints (demands and bounds). As the latter is defined as a Cartesian product of each
agent’s feasibility set, APM can operate in a distributed fashion. The sequence constructed
by APM converges to a single point if the intersection of the convex sets is nonempty, and
it converges to a periodic orbit of length 2 otherwise. If APM converges to a periodic orbit,
meaning that the disaggregation is not feasible, we construct from this orbit a polyhedral
cut, i.e. a linear inequality satisfied by all feasible solutions p of the global problem (1.1),
but violated from the current resource allocation (Theorem 3.3). Adding this cut to the
master problem (1.3) by updating Q to a specific subset, we can recompute a new resource
allocation and repeat this procedure until disaggregation is possible. At this stage, the use of
a cryptographic protocol, secure multiparty computation, allows us to preserve the privacy of
agents. Another main result stated in this paper is the explicit upper bound on the convergence
speed of APM in our framework (Theorem 3.2), which is obtained by spectral graph theory
methods, exploiting also geometric properties of transportation polytopes. This explicit speed
shows a linear impact of the number of agents, which is a strong argument for the applicability
of the method in large distributed systems.
1.3. Related Work. A standard approach (e.g. [34, 41, 37]) to solve resource allocation
problems in a distributed way is to rely on a Lagrangian based decomposition technique: for
instance dual subgradient methods [9, Ch. 6] or ADMM [18]. Such techniques are generally
used to decompose a large problems into several subproblems of small dimension. However,
those methods often require global convexity hypothesis, which are not satisfied in many
practical problems (e.g. MILP). We refer the reader to [9, Chapter 6] for more background.
On the contrary, our method can be used when the allocation problem (1.1) is not convex.
As explained in Section 4, the method proposed here can be related to Bender’s decomposition [8]. It differs from Bender’s approach in the way cuts are generated: instead of solving
linear programs, we use APM and our theoretical results, which provides a decentralized,
privacy-preserving and scalable procedure. In contrast, at each stage, Benders’ algorithm requires to solve a linear program requiring the knowledge of the private constraints of each
individual agent (see Remark 4.6 for more details).
The problem of the aggregation of constraints has been studied in the field of energy, in
3

the framework of smart grids [31, 3]. In [31], the authors study the management of energy
flexibilities and propose to approximate individual constraints by zonotopic sets to obtain an
aggregate feasible set. A centralized aggregated problem is solved via a subgradient method,
and a disaggregation procedure of a solution computes individual profiles. In [3], the authors
propose to solve the economic power dispatch of a microgrid, subject to several agents private
constraints, by using a Dantzig-Wolfe decomposition method.
APM has been the subject of several works in itself [20, 5, 7]. The authors of [12] provide
general results on the convergence rate of APM for semi-algebraic sets. They show that the
convergence is geometric for polyhedra. However, it is generally hard to compute explicitly the
geometric convergence rate of APM, as this requires to bound the singular values of certain
matrices arising from the polyhedral constraints. A remarkable example where an explicit
convergence rate for APM has been obtained is [33]. The authors consider there a different
class of polyhedra arising in submodular optimization. A common point with our results is
the use of spectral graph theory arguments to estimate singular values.
1.4. Structure. In Section 2, we describe the class of resource allocation problems that we
address. We formulate the idea of decomposition via disaggregation subproblems. In Section 3,
we focus on APM, the subroutine used to solve the disaggregation subproblems. In Subsection 3.1, after recalling basic properties of APM, we establish the key result on which relies
the proposed decomposition: how to generate a new cut to add in the master problem, from
the output of APM. In Subsection 3.2, we show how to guarantee the privacy of the proposed
procedure by using secure multiparty computation techniques. In Subsection 3.4 , we establish an explicit upper bound on the rate of convergence of APM in our case. In Section 4, we
generalize part of our results and propose a modified algorithm which applies to polyhedral
agents constraints. Finally, in Section 5, we present numerical examples: Subsection 5.1 gives
an illustrative toy example in dimension T = 4, while in Subsection 5.2, we consider a larger
scale, nonconvex example, coming from the microgrid application presented at the beginning
of the introduction.
Notation. Bold fonts, like “x”, are used to denote vectors, while normal fonts, like “x”,
refer to a scalar quantities. We denote by v > the transpose of a vector v. Recall that we
denote by [N ] the set {1, . . . , N }. Calligraphic letters such as T , N , X are used to denote
sets, and if T ⊂ [T ], T c := {t ∈ [T ] \ T } denotes the complementary set of T . The notation
U([a, b]) stands for the uniform distribution on [a, b]. The notation PC (.) refers to the Euclidean
projection onto a convex set C. For d ∈ N, 1d denotes the vector of ones (1 . . . 1)> ∈ Rd .
2. Resource Allocation and Transportation Structure.
2.1. A Decomposition Method Based on Disaggregation. As stated in the introduction,
we consider a centralized entity (e.g. an energy operator) interested in minimizing a possibly
nonconvex cost function p 7→ f (p), where p ∈ RT is the aggregate allocation of T dimensional
resources (for example power production over T time periods). This resource allocation p is
to be shared between N individual agents, each agent obtaining a part xn ∈ Xn , where Xn
denotes the individual feasibility set of agent n.
The global problem the operator wants to solve is described in (1.1). We assume that in
Problem (1.1), the constraints set Xn and individual profile xn are confidential to agent n and
should not be disclosed to the central operator or to another agent.
4

Let us define the set PD of feasible aggregate allocations that are disaggregable as:

P
(2.1)
PD := p ∈ P | ∃x ∈ X ; p = n xn .
Feasibility of problem (1.1) is equivalent to having PD not empty.
Constraints for each agent are composed of a global demand over the resources and lower
and upper bounds over each resource, as given below:
Assumption 1. For each n ∈ [N ], there exists fixed parameters En > 0, xn ∈ RT , xn ∈ RT
such that :
P
6 ∅.
(2.2)
Xn = {xn ∈ RT : t∈[T ] xn,t = En and xn,t ⩽ xn,t ⩽ xn,t } =
In particular, Xn is convex and compact. Given an allocation p, the structure obtained on
the matrix (xn,t )n,t , where sums of coefficients along columns and along rows are fixed, is
often referred to as a transportation problem. The latter has many applications (see e.g. [2,
32]). We focus on transportation problems in Sections 2 and 3, while in Section 4, we give a
generalization of some of our results in the general case where Xn is a polyhedron.
Given a particular allocation p ∈ P, the operator will be interested to Q
know if this allocation is disaggregable, that is, if there exists individual profiles (xn )n∈[N ] ∈ n Xn summing to
p, or equivalently if the disaggregation problem (1.2) has a solution.
Following (1.2), the disaggregate profile refers to x, while the aggregate profile refers to
the allocation p. Problem (1.2) may not always be feasible. Some necessary conditions for
a disaggregation to exist, obtained by summing the individual constraints on [N ], are the
following aggregate constraints:
(2.3a)

p> 1T = E > 1N

(2.3b)

and x> 1N ⩽ p ⩽ x> 1N .

Because they are necessary, we assume that those aggregate conditions (2.3) hold for vectors
of P, that is:
Assumption 2. All vectors p ∈ P satisfy (2.3), that is, P ⊂ {p ∈ RT | (2.3) hold }.
However, conditions (2.3) are not sufficient in general, as explained in the the following section.
2.2. Equivalent Flow Problem and Hoffman Conditions. Owing to its special structure,
the problem under study can be rewritten as as a flow problem in a graph, as stated in Proposition 2.2. We refer the reader to the book [15, Chapter 3] for terminology and background.
Definition 2.1. Consider a directed graph G = (V, E) with vertex set V, edge set E ⊂ V × V,
demands d : V → R (where dv < 0 means that v is a production node), edge lower capacities
` : E → R+ and upper capacities u : E → R+ . A flow on G is a function x : E → R+ such that
x satisfies
that is ∀e ∈ E, `e ⩽ xe ⩽ ue , and Kirchoff ’s law, that is,
P the capacity constraints,
P
∀v ∈ V, e∈δv+ xe = dv + e∈δv− xe , where δv+ (resp. δv− ) is the set of edges ending at (resp.
departing from) vertex v.
The following proposition is immediate:
Proposition 2.2. Consider the bipartite graph G whose set of vertices is the disjoint union
V = [T ]t[N ] and whose set of edges is E = {(t, n)}t∈[T ],n∈[N ] . Define demands on nodes [T ] by
dt = −pt and demands on nodes [N ] by dn = En . Assign to each edge (t, n) an upper capacity
5

p1 = 0

p2 = 3

1

E1 = 2

2

E2 = 0.5

3
[N ]

E3 = 0.5

t=1

t=2
[T ]

Figure 2.1: Example of a flow representation of the disaggregation problem (T = 2 and N = 3,
x = 0, x = 1). Here, the aggregate constraints (2.3) are verified, but condition (2.4) written
with A = {t1 , n1 } (dashed nodes) does not hold.
un,t = xn,t and lower capacity `n,t = xn,t . Then, finding a solution x to (1.2) is equivalent to
finding a feasible flow in G.
Hoffman [22] gave a necessary and sufficient condition
for the flow problem to be feasible
P
in a graph with balanced demands, that is d(V) := v∈V dv = 0 (total production matches
total positive demand, in our case (2.3a)). This generalizes a result of Gale (1957). The stated
condition is intuitive: there cannot be a subset of nodes whose demand exceeds its “import
capacity”.
Theorem 2.3 ([22]). Given a digraph G = (V, E) with demand d ∈ RV such that d(V) = 0
and capacities ` ∈ (R ∪ {−∞})E and u ∈ (R ∪ ∞)E with ` ⩽ u, there exists a feasible flow
x ∈ E → R+ on G if and only if:
X
X
X
(2.4)
∀A ⊂ V,
ue ⩾
dv +
`e ,
e∈δ + (A)

v∈A

e∈δ + (Ac )

where δ+ (A) := {(u, v) ∈ E|u ∈ Ac , v ∈ A} is the set of edges coming to set A and Ac := V \A.
Proposition 2.4 translates Theorem 2.3 in our framework:
Proposition 2.4. Given an aggregate resource allocation p = (pt )t∈[T ] ∈ P, the disaggregation problem is feasible, meaning that p ∈ PD , iff:
X
X
X
X
(2.5)
∀T ⊂ [T ], ∀N ⊂ [N ],
pt −
En +
xn,t ⩽
xn,t .
t∈T

n∈N

t∈T
/ ,n∈N

t∈T ,n∈N
/

c
c
P Proof. We apply (2.4) with A := T ∪ N and use the equality d(V) = 0 =
v∈Ac dv .

P

v∈A dv +

From Theorem 2.3 or Proposition 2.4 above, one can see that the aggregate constraints (2.3)
are in general not sufficient to ensure that the disaggregation problem has a solution.
For a given set T , there is a choice of N which leads to the strongest inequality (2.5),
namely:


X

X
X
X
xn,t ,
(2.6)
pt ⩽ min
En −
xn,t +

N ⊂[N ] 
t∈T

n∈N

t∈T
/ ,n∈N

6

t∈T ,n∈N
/

In this way, we get 2T − 2 inequalities corresponding to the proper subsets T ⊂ [T ]. Moreover,
in general, these 2T − 2 inequalities are not redundant. Although this is not stated in [22],
this is a classical result whose proof is elementary.
3. Disaggregation Based on APM.
3.1. Generation of Hoffman’s Cuts by APM. In this section, we propose an algorithm
that solves (1.1) while preserving the privacy of the agent’s constraints Xn and individual
profile xn ∈ RT . To do this, the proposed algorithm is implemented in a decentralized manner
and relies on the method of alternate projections (APM) to solve the disaggregation problem
(1.2).
Let us consider the polyhedron enforcing the agents constraints:
n
o
P
X := X1 × · · · × XN where Xn := xn ∈ RT+ |
x
=
E
and
∀t,
x
⩽
x
⩽
x
.
n,t
n
n,t
n,t
n,t
t∈[T ]
Besides, given an allocation p ∈ P, we consider the set of profiles aggregating to p :

P
Yp := x ∈ RN T | ∀t ∈ [T ], n∈[N ] xn,t = pt .
Note that Yp is an affine subspace of RN T (to be distinguished from P which is a subset of
RT ), and that Yp ∩ X is empty iff p ∈
/ PD , according to the definition of PD in (2.1). The idea
of the proposed algorithm is to build a finite sequence of decreasing subsets (P (s) )0⩽s⩽S such
that:
P = P (0) ⊃ P (1) ⊃ · · · ⊃ P (S) ⊃ PD .
At each iteration, a new aggregate resource allocation p(s) is obtained by solving an instance
of the master problem introduced in (1.3) with Q = P (s) :
min f (p)

(3.1a)

p∈RT

s.t. p ∈ P (s) .

(3.1b)

In the sequel, we will refer to (3.1) as the master problem. Our procedure relies on the following
immediate observation:
Proposition 3.1. If p(s) is a solution of (3.1) and x ∈ Yp(s) ∩X , then (p(s) , x) is an optimal
solution of the initial problem (1.1).
Having in hands a solution p(s) , we can check whether Yp(s) ∩ X =
6 ∅, and then find a
vector x in this set, using APM on X and Yp(s) , as described in Algorithm 3.1 below (where
Y = Yp ).
The idea of using cyclic projections to compute a point in the intersection of two sets
comes from von Neumann [39], who applied it to the case of affine subspaces. We next recall
the basic convergence result concerning the APM.
Theorem 3.2 ([20]). Let X and Y be two closed convex sets with X bounded, and let (x(k) )k
and (y (k) )k be the two infinite sequences generated by the APM on X and Y (Algorithm 3.1)
with εcvg = 0. Then there exists x∞ ∈ X and y ∞ ∈ Y such that:
(3.2a)

x(k) −→ x∞ ,

(3.2b)

kx

k→∞
∞
∞

− y k2 =

y (k) −→ y ∞ ;
k→∞

min

x∈X ,y∈Y

7

kx − yk2 .

Algorithm 3.1 Alternate Projections Method (APM)
Require: Start with y (0) , k = 0 , εcvg , a norm k · k on RN T
1: repeat
2:
x(k+1) ← PX (y (k) )
3:
y (k+1) ← PY (x(k+1) )
4:
k ←k+1
5: until kx(k) − x(k−1) k < εcvg
In particular, if X ∩ Y =
6 ∅, then (x(k) )k and (y (k) )k converge to a same point x∞ ∈ X ∩ Y.
The convergence theorem is illustrated in Figure 3.1 in the case where X ∩ Y = ∅, that is,
when the disaggregation problem (1.2) is not feasible. The idea of the algorithm proposed in
y (1)y

y∞

x∞

(0)

Y

•
x(1) x(0)

X

Figure 3.1: Alternate projections method (APM) on two sets X and Y. When X ∩ Y = ∅,
APM cycles over two points x∞ and y ∞ .
this paper is, when Yp(s) ∩ X = ∅, to use the resulting vectors x∞ and y ∞ to construct a new
subset P (s+1) by adding a constraint of type (2.5) to P (s) : indeed, from Proposition 2.4, we
know that, if Yp(s) ∩ X = ∅, there exists at least one violated inequality (2.5).
The difficulty is to guess such a violated inequality among the 2T possible inequalities. It
turns out that, using the output of APM, we can build such an inequality.
Suppose that we obtain x∞ 6= y ∞ as defined in Theorem 3.2: we get a periodic cycle of
APM, that is, we have x∞ = PX (y ∞ ) and y ∞ = PY (x∞ ), and the couple (x∞ , y ∞ ) is the
solution of the following optimization problem, with given parameters (En )n and (pt )t .
(3.3a)
(3.3b)

1
min kx − yk22
x,y 2
X
∀n ∈ [N ],
xn,t = En

(λn )

t∈[T ]

(3.3c)
(3.3d)

∀n ∈ [N ], ∀t ∈ [T ], xn,t ⩽ xn,t ⩽ xn,t
X
∀t ∈ [T ],
yn,t = pt

(µn,t , µn,t )
(νt ) ,

n∈[N ]

where λn ∈ R, µn,t , µn,t ∈ R+ and νt ∈ R are the Lagrangian multipliers associated to the
constraints (3.3b),(3.3c),(3.3d), with the associated Lagrangian function:
P
1
L(x, y, λ, µ, ν) = kx − yk22 − λ> ( t xn,t − En )n − µ> (x − x)
2
P
− µ> (x − x) − ν > ( n yn − p) .
8

The stationarity condition of the Lagrangian with respect to the variable yn,t yields:
∀n ∈ [N ], ∀t ∈ [T ], yn,t = xn,t + νt .

(3.4)

Let us consider the sets T ⊂ [T ] and N ⊂ [N ] defined from the output (x∞ , y ∞ ) of APM
on X and Yp as follows:

∞
(3.5a)
T := t ∈ [T ] | ∃n ∈ [N ], yn,t
> xn,t ,

P
P
N := n ∈ [N ] | En − t∈T
(3.5b)
/ xn,t −
t∈T xn,t < 0 .
In Theorem 3.3 below, we show that applying the inequality (2.5) with the sets T and N
defined in (3.5) defines a valid inequality for the disaggregation problem violated by the current
allocation p.
The intuition behind the definition of T and N in (3.5) is the following: T is the subset
of resources for which there is an over supply (which overcomes the upper bound for at least
one agent). Once T is defined, N is the associated subset of [N ] minimizing the right hand
side of (2.6). Indeed, (2.6) can be rewritten as:
(
)

X
X
X
X
X
En −
xn,t −
pt ⩽ min
xn,t
+
xn,t .
t∈T

N ⊂[N ]

n∈N

t∈T

t∈T
/

t∈T ,n∈[N ]

The Theorem 3.3 below is the key result on which relies the algorithm proposed in this paper.
Theorem 3.3. Consider the sequence of iterates (x(k) , y (k) )k∈N generated by APM on X
and Yp (see Algorithm 3.1). Then one of the following holds:
(i) if X ∩ Yp 6= ∅, then x(k) , y (k) −→ x∞ ∈ X ∩ Yp ;
k→∞

(ii) otherwise, if X ∩ Yp = ∅, then x(k) −→ x∞ ∈ X and y (k) −→ y ∞ ∈ Yp . Then,
k→∞

k→∞

considering the sets T and N defined in (3.5) gives an inequality of Hoffman (2.5) violated by
p, that is:
X
X
X
X
(3.6)
En −
xn,t −
xn,t < 0 .
pt +
n∈N

t∈T

t∈T ,n∈N
/

t∈T
/ ,n∈N

Moreover, this Hoffman inequality can be written as a function of x∞ , as follows:
X
X X
(3.7)
AT (x∞ ) <
pt with AT (x∞ ) :=
x∞
n,t .
t∈T

t∈T n∈[N ]

Before giving the proof of Theorem 3.3, we need to show some technical properties on the sets
T , N . For simplicity of notation, we use x and y to denote x∞ and y ∞ in Proposition 3.4
and the proof of Theorem 3.3.
Proposition 3.4. With x 6= y solutions of problem (3.3) (outputs of APM on X and Yp
with εcvg = 0),
yn,t > xn,t ;
(i) ∀t ∈ T , ∀n ∈
/ N , xn,t = xn,t and P
(ii) T = {t | νt > 0} = {t | pt > n xn,t }, where νt is the optimal Lagrangian multiplier
associated to (3.3d);
(iii) ∀n ∈ N , λn < 0 ;
9

(iv) ∀t ∈
/ T , ∀n ∈ N , xn,t = xn,t ;
(v) the sets T , T c , N and N c are nonempty.
The proof of Proposition 3.4 is technical and given in Appendix A. With Proposition 3.4, we
are now ready to prove Theorem 3.3.
Proof of Theorem 3.3. We have:
X
X
X
X
xn,t −
xn,t −
En +
pt
n∈N

=

X

xn,t +

t∈T ,n∈N
/

=

X X

X

t∈T
/

t∈T

=

X X

X

=

X X

n∈N

t∈T

t∈T

xn,t +
xn,t +

n∈N

−νt =

X

xn,t +


X

−

X

xn,t −

xn,t −

pt

(from (3.3b))

X

xn,t −

X

pt =

t∈T

xn,t −

n∈[N ]

|xn,t − yn,t |

pt (from Prop. 3.4 (i) and (iv))

t∈T

t∈T
/ ,n∈N

X X

t∈T

X

X
t∈T

t∈T ,n∈N
/

t∈T

n∈[N ]

xn,t −

t∈T
/ ,n∈N



n∈N
/



X

xn,t −

t

n∈N

t∈T

t∈T
/ ,n∈N

t∈T ,n∈N
/

XX

X

yn,t



n∈[N ]



n∈[N ]

using the stationarity conditions (3.4) and for all t ∈ T , νt > 0 by Prop.3.4 (ii). Moreover,
using:
X X
X
X
(3.8)
(xn,t − yn,t ) =
En −
pt = 0 ,
t∈[T ] n∈[N ]

n∈[N ]

t∈[T ]

obtained from Assumption 2, we see that:

X
X
(3.9)
−
|xn,t − yn,t | = −(kx − yk1 )/2 < 0 ,
t∈T

n∈[N ]

which shows (3.6). We now show that inequality (3.7) is obtained by a rewriting of (3.6),
indeed:
X
X
X
xn,t −
xn,t
En +
n∈N

=
=

t∈T
/ ,n∈N

t∈T ,n∈N
/

X X

xn,t +

X

xn,t −

X

n∈N t∈[T ]

t∈T ,n∈N
/

t∈T
/ ,n∈N

X

X

X X

t∈T ,n∈N

xn,t +

t∈T ,n∈N
/

xn,t =

xn,t

(from Prop.3.4 (i) and (iv))

xn,t = AT (x).

t∈T n∈[N ]

Remark 3.5. An alternative to APM is to use Dykstra’s projections algorithm [17]. It
is shown in [6] that the outputs of Dykstra’s algorithm also satisfy the conditions given in
Theorem 3.2, thus Theorem 3.3 will also hold for this algorithm. However, in this paper we
focus on APM as it is simpler and, to our knowledge, there is no guarantee that Dykstra’s
algorithm would be faster in this framework.
Suppose, as before, that the two sequences generated by APM on X and Y converge to two
distinct points x∞ and y ∞ . Then, at each round k, we can define from (3.4) and considering
10

(k)

(k)

any n ∈ [N ], the multiplier ν (k) = yn − xn
Theorem 3.3 is:

∞
tends to yn∞ − x∞
n := ν . The set T of

T ∞ := {t ∈ [T ] | 0 < νt∞ },

(3.10)

which raises an issue for practical computation, as ν ∞ is only obtained ultimately by APM,
possibly in infinite time. To have access to T ∞ in finite time, that is, from one of the iterates
(ν (k) )k , we consider the set:
(K)

T (K) := {t ∈ [T ] | Bεcvg < νt

(3.11)

},

where εcvg is the tolerance for convergence of APM as defined in Algorithm 3.1, B > 0 is a
constant, and K (depending on εcvg ) is the first integer such that kx(K) − x(K−1) k < εcvg .
We next show that we can choose B to ensure that T (K) = T ∞ for εcvg small enough. We
rely on the geometric convergence rate of APM on polyhedra [12, 33]:
Proposition 3.6 ([33]). If X and Y are polyhedra, there exists ρ ∈ (0, 1) such that the
sequences (x(k) )k and (y (k) )k generated by APM verify for all k ⩾ 1:
kx(k+1) − x(k) k2 ⩽ ρkx(k) − x(k−1) k2

ky (k+1) − y (k) k2 ⩽ ρky (k) − y (k−1) k2 .

and

Proposition 3.6 applies to any polyhedra X and Y. In Subsection 3.4 we shall give an explicit
upper bound on the constant ρ in the specific transportation case given by (1.1c) and (2.2).
From the previous proposition, we can quantify the distance to the limits in terms of ρ:
Lemma 3.7. Consider an integer K such that the sequence (x(k) )k⩾0 generated by APM
satisfies kx(K) − x(K−1) k ⩽ εcvg , then we have for any K 0 ⩾ K − 1:
0

kx∞ − x(K ) k ⩽

εcvg
.
1−ρ

Proof. From Proposition 3.6, we have for any k ⩾ K:
(k)

kx

−x

(K 0 )

k⩽

0
k−K
X

(K+s+1)

kx

(K+s)

−x

k⩽

s=0

0
k−K
X

0

0

ρs kx(K +1) − x(K ) k ⩽

s=0
0

1
εcvg ,
1−ρ

ε

cvg
so that, by taking the limit k → ∞, one obtains kx∞ − x(K ) k ⩽ 1−ρ
.

With this previous lemma, we can state the condition on B ensuring the desired property:
Proposition 3.8. Define ν := min{|νt∞ | > 0} (least nonzero element of ν ∞ ). If the con1
stants B and εcvg > 0 are chosen such that B > 1−ρ
and εcvg × 2B < ν, and Algorithm 3.1
stops at iteration K, then we have:
T (K) = T ∞ .
Proof. Let t ∈ T ∞ , that is νt∞ > 0 which is equivalent to νt∞ ⩾ ν by definition of ν. We
have:
X (K)
X
1
1
1 X ∞ X (K)
(K)
νt = (pt −
xn,t ) = (pt −
x∞
(
x −
xn,t )
n,t ) +
N
N
N n n,t
n
n
n
ε

ε

cvg
cvg
1
> νt∞ − 1−ρ
⩾ ν − 1−ρ
> εcvg (2B − 1−ρ
),

11

1
and this last quantity is greater than Bεcvg as soon as B ⩾ 1−ρ
, thus t ∈ T (K) .
Conversely, if t ∈ T (K) , then:

νt∞ =

X
X (K)
1
1
1 X ∞ X (K)
(pt −
x∞
(pt −
xn,t ) − (
x −
xn,t )
n,t ) =
N
N
N n n,t
n
n
n
(K)

⩾ νt

(K)

B
> νt
− 1−ρ

− Bεcvg ⩾ (B − B)εcvg ⩾ 0 ,

so that t ∈ T ∞ . Furthermore, the “approximated” cut

P

P

(K)

n∈[N ] xn,t − pt



t∈T
violated by the current value of p (or p(s) at iteration s) in the algorithm as:


X


(K)

X

xn,t − pt  ⩽


t∈T


X

(K)



+
xn,t − x∞
n,t


t∈T

n∈[N ]


X

X


X


t∈T

n∈[N ]

⩾ 0 is


x∞
n,t − pt

n∈[N ]

1
⩽ kx(K) − x∞ k1 − kx∞ − y ∞ k1
2
using (3.8) and (3.9). This last quantity is negative as soon as kx(K) − x∞ k1 < 21 kx∞ − y ∞ k1 ,
which holds in particular if Bεcvg < 21 kx∞ − y ∞ k1 .
This second proposition shows a surprising result: even if we do not have access to the
limit x∞ , we can compute in finite time the exact left hand side term AT (x∞ ) of the cut
(3.7):
Proposition 3.9. Under the hypotheses of Proposition 3.8, we have:
X X

AT (x(K) ) =

(K)

xn,t = AT (x∞ ) .

t∈T n∈[N ]

Proof. We start by showing some technical properties similar to Proposition 3.4:
Lemma 3.10. The iterate x(K) satisfies the following properties:
(K)
(i) ∀t ∈ T , ∀n ∈
/ N , xn,t = x∞
n,t = xn,t ;
(K)

(ii) ∀t ∈
/ T , ∀n ∈ N , xn,t = x∞
n,t = xn,t .
The proof of Lemma 3.10 is similar to Proposition 3.4 and is given in Appendix B. Then,
having in mind that T (K) = T ∞ from Proposition 3.8, and N is obtained from T ∞ by (3.5),
we obtain:
X (K) 
X
X X
X
(K)
AT (x(K) ) =
xn,t −
xn,t +
xn,t
xn,t +
n∈N

⩽

t∈T
/

X X

(K)
xn,t −

n∈N t∈[T ]

which equals to AT (x∞ ) as we have

t∈T

t∈T
/ ,n∈N

X

xn,t +

t∈T
/ ,n∈N

X

t∈T ,n∈N
/

xn,t

(from Lemma 3.10)

t∈T ,n∈N
/

(K)
t∈[T ] xn,t = En for each n ∈ [N ].

P

Before presenting our algorithm using this last result, we focus on the technique of secure multiparty computation (SMC) which will be used here to ensure the privacy of agent’s
constraints and profiles while running APM.
12

3.2. Privacy-Preserving Projections through SMC. APM, as described in Algorithm 3.1,
enables a distributed implementation in our context, by the structure of the algorithm itself:
the operator computes the projection on Yp while each agent n can compute, possibly in
parallel, the projection on Xn of the new profile transmitted by the operator. This enables
each agent (as well as the operator) to keep her individual constraint and not reveal it to the
operator or other agents. However, if this agent had to transmit back her newly computed
individual profile to the operator for the next iteration, privacy would not be preserved. We
next show that, using a secure multi-party computation (SMC) summation protocol [42] we
can avoid this communication of individual profiles and implement APM without revealing
the sequence of agents profiles x to the operator.
For this, we use the fact that Yp is an affine subspace and thus the projection on Yp can
be obtained explicitly component-wise. Indeed, as PYp (x) is the solution y of the quadratic
program miny∈RN T | Pn yn =p 12 kx − yk22 , we obtain optimality
conditions similar to (3.4), and
P
summing these equalities on [N ], we get N νt = pt − n∈[N ] xn,t for each t ∈ [T ], and thus:
P
(3.12)
∀n ∈ [N ], [PYp (x)]n = xn + N1 (p − m∈[N ] xm ) .
P
Thus, having access to the aggregate profile S := n∈[N ] xn , each agent can compute locally
the component of the projection on Yp of her profile, instead of transmitting the profile to the
operator for computing the projection in a centralized way.
Using SMC principle, introduced by [42] and [19], the sum S can be computed in a
non-intrusive manner and by several communications between agents and the operator, as
described in Algorithm 3.2. The main idea of the SMC summation protocol [4, 38] below is
that, instead of sending her profile xn , agent n splits xn,t for each t into N random parts
(sn,t,m )m , according to an uniform distribution and summing to xn,t (Lines 2-3). Thus, each
part sn,t,m taken individually does not reveal any information on xn nor on Xn , and can be
sent to agent m. Once all exchanges of parts are completed (Line 5), and n has herself received
the parts from other agents, agent n computes a new aggregate quantity σn (Line 7), which
does not contain either any information about any of the agents, and sends it to the operator
(Line 8). The operator can finally compute the quantity S = x> 1N = σ > 1N .
Algorithm 3.2 SMC of Aggregate (SMCA)

P

n∈[N ] xn

P
Require: profile xn for each n ∈ [N ], parameter R > maxt max{ n xn,t } given to each agent
1: for each agent n ∈ [N ] do
−1
2:
Draw ∀t, (sn,t,m )N
R]N −1 )
m=1 ∈U([0,
PN −1
3:
and set ∀t, sn,t,N :=xn,t − m=1 sn,t,m mod R
4:
Send (sn,t,m )t∈[T ] to agent m ∈ [N ]
5: done
6: for each agent n ∈ [N ] do
P
7:
Compute ∀t, σn,t = m∈[N ] sm,t,n mod R
8:
Send (σn,t )t∈[T ] to operator
9: done
P
10: Operator computes S =
n∈[N ] σn mod R (and broadcasts it to agents)
Remark 3.11. A drawback of the proposed SMC summation (Algorithm 3.2) is that each
agent need to exchange with all other nodes, therefore the communication overhead for each
13

iteration and each node is O(N `) where ` is an upper bound on the length of the message
(e.g number of bits encoding constant R). Splitting the messages with the other N − 1 nodes
is required to obtain privacy guarantees against collusion of strictly less than N − 1 nodes
(see Subsection 3.3). However, as detailed in [4, Protocol I] the SMC summation can be
adapted by splitting secret numbers among only k ∈ {1, . . . , N } members. In that case, the
communication overhead would be reduced to O(k`), but privacy would only be protected
against collusion of less than k − 1 agents. This choice has to be made as a tradeoff between
computation overhead and privacy guarantees.
We sum up in Algorithm 3.3 below the procedure of generating a new constraint as stated in
Theorem 3.3 from the output of APM in finite time (see Proposition 3.8) and in a privacypreserving way using SMC.
Algorithm 3.3 Non-intrusive APM (NI-APM)
Require: Start with y (0) , k = 1, εcvg , εdis , norm k.k on RN T
1: repeat
2:
for each agent n ∈ [N ] do
(k)
(k−1)
3:
xn ← PXn (yn
)
4:
done
5:
Operator obtains S (k) ←SMCA(x(k) ) (cf Algo. 3.2)
6:
and sends ν (k) := N1 (p − S (k) ) ∈ RT to agents [N ]
7:
for each agent n ∈ [N ] do
(k)
(k)
8:
Compute yn ← xn + ν (k) . from (3.4) and (3.12), y (k) = PYp (x(k) )
9:
done
10:
k ←k+1
11: until kx(k) − x(k−1) k < εcvg
12: if kx(k) − y (k) k ⩽ εdis then
. found a εdis -solution of the disaggregation problem
(k)
13:
Each agent adopts profile xn
14:
return Disag ← True
15: else . have to find a valid constraint violated by p
(k)
16:
Operator computes T ← {t ∈ [T ] | Bεcvg < νt }
(k)
17:
OperatorP
computes AT ← SMCA( (xt )t∈T )
18:
if AT − t∈T pt < 0 then
19:
return Disag ← False, T , AT
20:
else . need to run APM with higher precision
21:
Return to Line 1 with εcvg ← εcvg /2
22:
end
23: end
To choose B and εcvg satisfying the conditions of Proposition 3.8 a priori, one has to
know the value of ν. Although a conservative lower bound could be obtained by Diophantine
arguments if we consider rationals as inputs of the algorithm, in practice it is easier and more
efficient to proceed in an iterative manner for the value of εcvg . Indeed, one can start with εcvg
arbitrary large so that APM will converge quickly, and then check if the cut obtained is violated
by the current value of p (Subsection 3.2): if it is not the case, we can continue the iterations of
APM with convergence precision improved to εcvg /2 (Subsection 3.2). Proposition 3.8 ensures
14

that this loop terminates in finite time.
The parameter εdis > 0 (Line 12 of Algorithm 3.3) has to be chosen a priori by the operator,
depending on the precision required. In general in APM, x∞ = y ∞ will only be achieved in
infinite time, so choosing εdis strictly positive is required.
We end this section by summarizing in Algorithm 3.4 the global iterative procedure to compute an optimal and disaggregable resource allocation p, solution of the initial problem (1.1),
using iteratively NI-APM (Algorithm 3.3) and adding constraints as stated in Theorem 3.3.
Algorithm 3.4 Non-intrusive Optimal Disaggregation
Require: s = 0 , P (0) = P ; Disag= False
1: while Not Disag do
2:
Solve minp∈P (s) f (p)
3:
if problem infeasible then
4:
Exit
5:
else
6:
Compute p(s) = arg minp∈P (s) f (p)
7:
end
8:
Disag← NI-APM(p(s) ) (Algo. 3.3)
9:
if Disag then
10:
Operator adopts p(s)
11:
else
(s)
12:
Obtain T (s) , AT from NI-APM(p(s) )
P
(s)
13:
P (s+1) ← P (s) ∩ {p| t∈T (s) pt ⩽ AT }
14:
end
15:
s←s+1
16: done
Algorithm 3.4 iteratively calls NI-APM (Algorithm 3.3) and in case disaggregation is not
possible (Line 11), a new constraint is added (Line 13), obtained from the quantity AT defined
in (3.7), to the feasible set of resource allocations P (s) in problem (3.1). This constraint is an
inequality on p and thus does not reveal significant individual information to the operator. The
algorithm stops when disaggregation is possible (Line 9). The termination of Algorithm 3.4 is
ensured by the following property and the form of the constraints added (3.6):
Proposition 3.12. Algorithm 3.4 stops after a finite number of iterations, as at most 2T − 2
constraints (Line 13) can be added to the master problem (Line 2).
The following Proposition 3.13 shows the correctness of our Algorithm 3.4.
Proposition 3.13. Let B and εcvg satisfy the conditions of Proposition 3.8 and 3.9. Then:
(i) if the problem (1.1) has no solution, Algorithm 3.4 exits at Line 4 after at most 2T − 2
iterations;
(ii) otherwise, Algorithm 3.4 computes, after at most s ⩽ 2T − 2 iterations, an aggregate
solution p(s) ∈ P, associated to individual profiles (x∗ )n = NI-APM(p(s) ) such that:
p(s) ∈ P,

, ∀n ∈ [N ], x∗n ∈ Xn ,

∗
(s)
n∈[N ] xn − p

P

where f ∗ is the optimal value of problem (1.1).
15

⩽ εdis

and f (p(s) ) ⩽ f ∗ ,

Proof. The proof is immediate from Theorem 3.3, Proposition 3.8 and Proposition 3.9.
Remark 3.14. The upper bound on the number of constraints added has no dependence
on N because, as stated in (2.6), once a subset of [T ] is chosen, the constraint we add in the
algorithm is found by taking the minimum over the subsets of [N ].
Although there exist some instances with an exponential number of independent constraints,
this does not jeopardize the proposed method: in practice, the algorithm stops after a very
small number of constraints added. Intuitively, we will only add constraints “supporting” the
optimal allocation p. Thus, Algorithm 3.4 is a method which enables the operator to compute
a resource allocation p and the N agents to adopt profiles (xn )n , such that (x, p) solves the
global problem (1.1), and the method ensures that both agent constraints (upper bounds xn ,
lower bounds xn , demand En ); and disaggregate (individual) profile xn (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 by a
third party (either the operator or any other agent m 6= n).
Remark 3.15. A natural approach to address problem (1.1) in a distributed way, assuming
that both the cost function p 7→ f (p) and the feasibility set P are convex, is to rely on
Lagrangian based decomposition techniques. Examples of such methods are Dual subgradient
methods [9, Chapter 6], auxiliary problem principle method [14], ADMM [18],[43] or bundle
methods [29].
It is conceivable to develop a privacy-preserving implementation of those
P techniques, where
Lagrangian multipliers associated to the (relaxed) aggregation constraint n xn = p would be
updated using the SMC technique as described in Algorithm 3.2. However, those techniques
usually ask for strong convexity hypothesis: for instance, in ADMM, in order to keep the
decomposition structure in agent by agent, a possibility is to use multi-blocs ADMM with
N + 1 blocs (N agents and the operator), which is known to converge in the condition that
strong convexity of the cost function in at least N of the N + 1 variables holds [16]. The study
of privacy-preserving implementations of Lagrangian decomposition methods is left for further
work.
The advantage of Algorithm 3.4 proposed in this paper is that convergence is ensured (see
Proposition 3.12) even if the cost function p 7→ f (p) and the feasibility set P are not convex,
which is the case in many practical situations (see Subsection 5.2).
3.3. Privacy guarantees. We next analyze the protection of privacy of the present optimization algorithm, considering different situations. For this analysis, we consider the paradigm of honest, semi-honest and malicious agents. The term semi-honest (also called honest
but curious [1],[36]) refers to a user who obeys the protocol, but wishes to use the data obtained through the algorithm to infer private information. In contrast, a malicious user may
even disobey the protocol.
Malicious operator and honest agents. As a first step, we consider an ideal model, in which
we only analyze the loss of privacy caused by the algorithm, through communications between
honest agents and a malicious (or semi-honest) central operator.
We neglect, in this situation, the information leaks induced by inter-agent communications
through the secure multiparty protocol.
Through Algorithm 3.4, the operator obtains the following information for each iteration
s ∈ {1, . . . , scv } (where scv < ∞ is the number of iterations needed for convergence):
(s)
• the set T (s) and the scalar AT defining the cut;
16

s
• whenever ni-apm is called, the serie (S (s,k) )K
k=1 of the aggregated projections.

cv
(0) , A(0) , . . . , (S (scv ,k) )Kscv , T (scv ) , A(s ) the informa0
Let us denote by Iop := (S (0,k) )K
T
T
k=1
k=1 , T
tion obtained by the operator. The next proposition implies that, in this ideal model, the
operator cannot discern the individual profiles or individual constraints.

Proposition 3.16. Consider two different instances of the resource allocation problem, differing only by a permutation of the agents, meaning that for all n ∈ [N ], Xn is replaced by
Xσ(n) , for some permutation σ of the elements of [N ]. Then, the information Iop gotten by the
operator (as well as the information (ν (k) )k sent by the operator to agents in Algorithm 3.3),
and the sequence of cuts added in the master problem (Algorithm 3.4), are precisely the same
in both instances.
Proof. The only information transmitted by the agents to the central operator consists
of the aggregate profiles S (k) = SMCA(x(k) ) and the sums AT (line 5 and line 17 of Algorithm 3.3), for different values of the subset T , corresponding to different Hoffman cuts. Each
aggregate profile S (k) is invariant by permutation of the agents. Each sum AT is of the form
P
(K)
AT = n∈[N ],t∈T xnt , where K is the iteration at which the violated inequality is found.
Hence it is also invariant by a permutation of the agents.
To illustrate Proposition 3.16, suppose Alice and Bob live in a small town, in which there is a
sports Hall, with two kind of lessons every day, boxing, from 12h00 to 13h00, and dancing, from
18h00 to 19h00. Suppose Alice goes to the boxing lesson, and that Bob goes to the dancing
lesson. Suppose further that Alice and Bob do not consume any electric energy during their
lessons. Then, by analyzing the successive aggregates AT , for different values of T ⊂ [T ], the
central operator may deduce that the global consumption decreases at the time of the lessons,
implying that one agent boxes whereas another one dances. However, whether it is Alice or
Bob who does boxing or dancing remains inaccessible to the operator, owing to the symmetry
argument. This protection persists even with a malicious operator. This privacy property is
stated formally as follows:
Corollary 3.17. With N > 1, a malicious (or semi-honest) operator using Algorithm 3.4
cannot infer the constraints Xn or profile xn of a particular agent n ∈ [N ] with probability 1.
Remark 3.18. Even if the identity of users is fully protected, in some specific cases the
central operator could infer limited information about constraints or profile that an agent
within the population may bear. Indeed, let us consider an illustrative example with T = 2,
N = 2, and constraints parameters (E1 , x1,1 , x1,2 ) := (1, 2, 0), (E2 , x2,1 , x2,2 ) := (3, 1, 2) and
x := 0 (all unknown to the operator). The operator knows a priori the aggregate quantities
E1 + E2 = 4, x1 + x2 = (3, 2), x1 + x2 = (0, 0).
Before the algorithm, from the operator viewpoint, any parameter values satisfying the
aggregate conditions above are possible. For instance, it is possible that agent 1 has the
parameters (E1∗ , x∗1,1 , x∗1,2 , x∗1,1 , x∗1,2 ) := (3, 2, 2, 0, 0).
Now suppose that the first aggregate profile proposed by the operator is p(0) := (3, 1).
(0)
From Algorithm 3.4, he obtains the cut defined by T (0) = {1} and A{1} = 2. Thus, assuming
P
P
(0)
(0)
the operator knows that A{1} is of the form, A{1} = minN ⊂[N ] { n∈N En + t∈T ,n∈N
/ xn,t −
17

P

t∈T
/ ,n∈N xn,t } for a subset N , he infers the following conditions on the parameters:

(3.13)

E2 + x1,1 − x2,2 = 2
E1 + x2,1 − x1,2 = 2
OR
and E2 + x1,1 − x2,2 ⩾ 2
and E1 + x2,1 − x1,2 ⩾ 2 .

With these conditions, the parameters values (E1∗ , x∗1,1 , x∗1,2 , x∗1,1 , x∗1,2 ) are no longer possible for
agent 1 (neither for agent 2 as the information the operator obtains is necessarily symmetric).
Indeed, from the initial aggregate conditions, this implies that (E2∗ , x∗2,1 , x∗2,2 , x∗2,1 , x∗2,2 ) :=
(1, 1, 0, 0, 0), which is incompatible with (3.13).
It seems difficult to have an algorithm where the operator will not learn any information on
the distribution of parameters. The kind of information leak illustrated in Remark 3.18 is also
observed in other privacy-preserving algorithms. One such example is the privacy-preserving
consensus method to compute an average value of numbers secretly owned by agents of [36].
Even if the secret number of a particular agent cannot be learned by another agent (see [36,
Th.3]), at the end of the procedure, each agent has obtained new information on the joint
distribution of initial numbers of other agents, as each agent obtains the average value Avg[0]
of those numbers.
Fortunately, in most cases, an operator using Algorithm 3.4 will only learn very limited
information on the distribution of constraints and profiles of users: for this, it is instructive
to estimate the information leak caused by the algorithm by comparing the information sent
to the central operator with the one needed to encode the instance. Suppose, for simplicity,
that all the numbers manipulated by the algorithm (in particular the bounds xn,t and x̄n,t
are encoded by fixed (say double) precision numbers. Then, the total number of bits needed
to encode the instance is O(N T ). However, it follows from Proposition 3.13 that the number
of bits received by a malicious central operator is O(2T ), hence, in the regime N  2T
(large number of agents), the information leak is asymptotically negligible. Moreover, we
give in Subsection 5.2 a real example with N = 28 agents, T = 24, in which the algorithm
terminates after 194 iterations. Here, the information leak consists of 194 doubles, to be
compared with the encoding size of the instance, N T + N = 12544 doubles.
Semi-honest agents. Let us now consider agents that are only semi-honest, meaning that
they still obey the protocol, but wish to exploit the data received in the algorithm to infer
information about the other users.
A semi-honest agent (or an external eavesdropper) aiming to learn the profile xn of agent
n has to intercept all the communications between n and all N − 1 other agents (to learn
(sn,t,m )m6=n and (sm,t,n )m6=n ) and to the operator (to learn σn ) during the SMC protocol
in order to succeed (otherwise, this semi-honest agent has just access to random numbers).
Besides, if we think of collusion of semi-honest agents aiming to learn the profile of a specific
agent, we can observe that the collusion must involve N − 1 agents (all except one) to succeed
in learning anything.
Indeed, in the SMC summation protocol (Algorithm 3.2, the messages sn,t,j sent by agent
n to agent j, encoding shares of the consumption xnt of agent n at different times, are all
uniformly distributed random variables, because the sum of independent uniform random
variables modulo R is still uniform. Similarly, the messages σn,t sent by agent n to the central
operator are uniformly distributed random variables. Furthermore, there are no correlations
between the random variables received by an agent at different iterations. This entails that,
even if some agents keep the same consumption profile over time, a semi-honest agent cannot
18

learn the profiles of the other agents by comparing the information received during different
iterations of the SMC protocol.
We refer the reader to [4, Protocol I], where a finer variant of SMC (secure split protocol),
splitting secret numbers only among k < n randomly chosen players, is analyzed.
Alternative
SMC methods to the proposed protocol Algorithm 3.2 exist to compute the
P
sum n∈[N ] yn , while keeping yn secret to user n, for instance [13]. Other techniques could
also be considered such as consensus-based aggregation algorithms [21].
From the above paragraph, one deduces that a collusion of less than (N − 1) semi-honest
agents cannot obtain more information than the operator from an execution of Algorithm 3.4:
indeed, the only additional information that the collusion obtains is composed of the sequences
{ν (k) } obtained from each execution of Algorithm 3.3. Thus, the symmetry argument given
in Proposition 3.16 applies to give a privacy guarantee similar to Corollary 3.17:
Corollary 3.19. Using Algorithm 3.4 with N > 1, a collusion of less than (N − 1) malicious
(or semi-honest) agents cannot infer the constraints Xn or profile xn of a particular agent n
(which is not in the collusion) with probability 1.
Malicious agents and Robustness. The proposed procedure is not robust against malicious
agents, i.e., agents lying about their consumption to jeopardize the system. Indeed, if agents
lie about their consumptions or feasibility constraints, Algorithm 3.4 will not converge to
the global optimum. However, we still have the privacy guarantee given in Corollary 3.17:
malicious agents lying to the operator in order to learn information about other agents can
only learn global information, and will not be able to learn individual information (profile or
constraints) about a specific agent.
A more detailed privacy analysis of our method and of its possible refinements can be an
avenue for further work.
In the next section, we focus on the convergence rate of APM in the particular case
of transportation constraints, and give an explicit bound on the geometric rate stated in
Theorem 3.2.
3.4. Complexity Analysis of APM in the Transportation Case. In this section we analyze
the speed of convergence of the alternate projections method (APM) described in Algorithm 3.1
on the sets X and Yp defined in Section 2.
A general result in [12] gives an upper bound of the sequences generated by APM on X and
Y if these two sets are semi-algebraic. In particular, it establishes the geometric convergence
for polyhedral sets. However, as stated in [33], given two particular polyhedral sets X and Y,
it is not straightforward to deduce an explicit rate of convergence from their result.
The authors in [33] established in a particular case a geometric convergence with an explicit
upper bound on the convergence rate. They consider APM on two sets P and Q, where P is
a linear subspace and Q is a product of base polytopes of submodular functions.
In this section, we also establish an explicit upper bound on the convergence rate of APM
in the transportation case, that is with X and Yp defined in (2.2) and (1.2b):
Theorem 3.20. For the two sets X and Yp , the sequence of alternate projections converges
to x∗ ∈ X , y ∗ ∈ X P satisfying kx∗ − y ∗ k = inf x∈X ,y∈Yp kx − yk, at the geometric rate:

k
kx(k) − x∗ k ⩽ 2kx(0) − x∗ k × 1 − N (T +1)42 (T −1) ,
19

and the analogous inequalities hold for (y (k) )k .
For the remaining of this section, we will just use Y to denote Yp , as p remains fixed during
APM. For the result stated in Theorem 3.20 above, we use several lemmas from [33].
Proof. First, we use the fact stated in [33] that APM on subspaces U and V converge with
geometric rate cF (U, V )2 , where the rate is given by the square of the cosine of the Friedrichs
angle between U and V , given by:
cF (U, V ) = sup{uT v | u ∈ U ∩ (U ∩ V )⊥ , v ∈ V ∩ (U ∩ V )⊥ , kuk ⩽ 1, kvk ⩽ 1}.
An intuitive generalization of this result for polyhedra X and Y, considering all affine subspaces
supporting the faces of X and Y is given in [33]:
Lemma 3.21 ([33]). For APM on polyhedra X and Y in RD , the convergence is geometric
with rate bounded by the square of the maximal cosine of Friedrichs angle between subspaces
supporting faces of X and Y:

max cF aff 0 (Xx ), aff 0 (Yy ) ,

(3.14)

x,y

where, for any x ∈ RD , Xx := arg maxv∈X x> v is the face of X generated by direction x and
aff 0 (C) = aff(C) − c for some c ∈ C denotes the subspace supporting the affine hull of C, for
C = Xx or C = Yy .
In the remaining of the proof, we bound the quantity (3.14) for our polyhedra X and Y.
For this, we use the space RN T = RT × · · · × RT , where the (n − 1)T + 1 to nT entries
correspond to the profile of agent n, for 1 ⩽ n ⩽ N . As in [33], we use a result connecting
angles between subspaces and the eigenvalues of matrices giving the directions of these spaces:
Lemma 3.22 ([33]). If A and B are matrices with orthonormal rows with same number of
columns, then:

> are equal to one, then c
– if all the singular values of AB
F Ker A, Ker B = 0;

– otherwise, cF Ker A, Ker B is equal to the largest singular value of AB > among those that
are smaller than one.
We are left with finding a matrix representation of the faces of polyhedra X and Y and, then,
bounding the corresponding singular values.
√ −1
In our case, the polyhedra Y is an affine subspace Y = {x ∈ RN T |Ax = N p} where:
√
A :=

N

−1

J1,N ⊗ IT ,

where ⊗ denotes the Kronecker product. The matrix A has orthonormal rows and the linear
subspace associated to Y is equal to Ker(A).
Obtaining a matrix representation of the faces of X is more complex. The faces of X are
obtained by considering, for each n ∈ [N ], subsets of the time periods that are at lower or
upper bound (respectively T n and T n , with T n ∩ T n = ∅). Considering a collection of such
subsets, a face of X can be written as:
n
o
P
A(T n ,T )n:= (x)n,t |∀n,
x
=
E
and
∀t
∈
T
,
x
=
x
,
and
∀t
∈
T
,
x
=
x
.
n,t
n
n,t
n
n,t
n,t
n
n,t
t
n

20

For some particular collection of subsets (T n , T n )n , the set A(T n ,T )n might be empty. The
n
linear subspace associated to A(T n ,T )n is given by {x ∈ RN T | Bx = 0} = Ker(B), where the
n
P
N first rows of B, corresponding to the constraints t xn,t = En , are given before orthonormalization by:
√
T

−1

IN ⊗ J1,T ,

P

and the matrix B has b := n |Tn | more rows, where Tn := Tn ∪ Tn , corresponding to the
N T for n ∈ [N ],
saturated bounds. Each of this row is given by the unit vector e>
n,t ∈ R
t ∈ Tn , which gives already an orthonormalized family of (unit) vectors. Therefore, a simple
orthonormalized matrix B ∈ MN +b,N T (R) giving the direction of A(T n ,T )n is given by:
n

p

p
−1
−1 >
B := diag
T − |T1 | 1>
,
.
.
.
,
T
−
|T
|
1
c
c
N
T
T


1

N

>
,
diag (BT1 , . . . , BTN )

|

T
c
where 1P
Tnc ∈ R is the vector where the indices in Tn are equal to 1 and 0 otherwise, and
Ektk is the matrix |Tn | × T with indices of Tn . We obtain the double
BTn :=
1⩽k⩽|Tn |
Tn ={t1 ,...,t|Tn | }

product:


P 1k∈T
1 X >
/ n ∧`∈T
/ n
+
BTn BTn
n
T − |Tn | 1⩽k,`⩽T N
n


1 X P
1 P 1{k,`}⊂Tnc
+
( n 1t∈Tn ) Et,t .
=
n
T − |Tn | 1⩽k,`⩽T N
N

1
(AB )(BA ) =
N
>

>

1⩽t⩽T

We observe
T that:
– if t0 ∈ N
i=1 Tn , then et0 is an eigenvector associated to eigenvalue λt0 = 1 ;
T
c
– the vector 1T̄ := (1t∈∩
/ n Tn )t∈[T ] ∈ R , where T̄ := ∪n Tn , is an eigenvector associated to
eigenvalue λ = 1. Indeed, if we denote by Nθ = {n ∈ [N ]|θ ∈ Tn }, then [1T̄ ]θ = 1 ⇔ Nθc 6= ∅,
and for each θ ∈ T̄ :


X X 1t∈T
X
1
/ n

[(AB > )(BA> )]θ 1T̄ =
[1 ]t +
1θ∈Tn [1T̄ ]θ 
N
T
−
|Tn | T̄
c
n
t
i∈Nθ


X
|N c | + |Nθ |[1T̄ ]θ
1  X T − |Tn |
=
1+
1 × [1T̄ ]θ  = θ
= [1T̄ ]θ .
N
T − |Tn |
N
c
i∈Nθ

i∈Nθ

To bound the other eigenvalues of the matrix (AB > )(BA> ), we rely on spectral graph
theory arguments. Consider the weighted graph G = ([T ], E) whose vertices are the time
P 1
c
n
periods [T ] and each edge (k, `) ∈ [T ] × [T ] with k 6= ` has weight Sk,` = N1 n {k,`}⊂T
T −|Tn | (if
this quantity is zero, then there is no edge between k and `).
The matrix P := IT − (AB > )(BA> ) verifies for each k ∈ [T ]:
(3.15)

X

−Pk,` =

`6=k

(3.16)

= N1

X
`6=k

1
N

X 1{k,`}⊂T c
n

n

T − |Tn |

= N1

X 1k∈T c (T − |Tn | − 1)
n

n

X
X 1k∈T c
n
(1 − 1k∈Tn ) − N1
= Pkk ,
T
−
|T
n|
n
n
21

T − |Tn |

which shows that P is the Laplacian matrix of graph G. As Sp(AB > BA> ) = 1−Sp(P ), we
want to have a lower bound on the least eigenvalue of P greater than 0, that we denote by λ1 .
By rearranging the indices of [T ] in two blocs T̄ and T̄ c , we observe that P can be written
as a block diagonal matrix P = diag(PT̄ , 0T̄ c ). As we are only interested in the positive
eigenvalues of P, we can therefore study the linear application associated to P restricted to
the subspace Vect(et )t∈T̄ .
As 1T̄ is an eigenvector of P associated to λ0 = 0, from the minmax theorem, we have:
λ1 =

(3.17)

u> P u
.
u⊥1T̄ ,u6=0 u> u
min

Let us consider an eigenvector u realizing (3.17). Let ut∗ := maxt ut and us∗ := mint ut and
let ds∗ ,t∗ be the distance between s∗ and t∗ in G, and let (s∗ -t∗ ) denote a shortest path from
s∗ to t∗ in G. As P is a Laplacian matrix, we have:
(3.18)
u> P u =

1
1 X
−Pk,` (uk − u` )2 ⩾
2
2

X

−Pk,` (uk − u` )2 ⩾

{k,`}∈(s∗ -t∗ )

k,`∈T̄

min

(−Pk,` )

k,`∈(s∗ -t∗ )

(ut∗ − us∗ )2
,
ds∗ ,t∗

where the last inequality is obtained from Cauchy-Schwartz inequality.
Let us write the path (s∗ -t∗ ) = (t0 , t1 , . . . , td ). As (s∗ -t∗ ) is a shortest path, for each k ∈
{0, d−1}, the edge (tk , tk+1 ) exists so there exists n ∈ [N ] such that {tk , tk+1 } ⊂ Tnc . Moreover,
for each n, we have Tnc ∩{t0 , . . . , tk−1 , tk+2 , . . . , td } = ∅, otherwise we could “shortcut” the path
(s∗ -t∗ ), thus we have |Tn | ⩾ d − 1. We obtain:
−Ptk ,tk+1 =

1 X 1{tk ,tk+1 }⊂Tnc
1
⩾
.
N n
T − |Tn |
N (T − d + 1)

T √
t∗
kuk2 .
= T T−1 ut∗ ⩾ (T −1)
On the other hand, we have (ut∗ − us∗ ) ⩾ ut∗ + Tu−1
T

Using these bounds and (3.18), we obtain:
u> P u ⩾

(ut∗ − us∗ )2
4
4T
⩾
kuk22 ⩾
kuk22 .
2
2
N (T − ds∗ ,t∗ + 1)ds∗ ,t∗
N (T + 1) (T − 1)
N (T + 1)2 (T − 1)

Therefore, λ1 ⩾ N (T +1)42 (T −1) := κN,T and the greatest singular value lower than one of
(AB > )(BA> ) is 1−κN,T . We conclude by applying successively Lemma 3.22 and Lemma 3.21,
to obtain the convergence rate stated in Theorem 3.20.
4. Generalization to Polyhedral Agents Constraints. In this section, we extend our
results to a more general framework where for each n ∈ [N ], Xn is an arbitrary polyhedron,
instead of having the particular structure given in (2.2). Let us now consider that (Xn )n are
polyhedra with, for each n:
(4.1)

Xn = {xn ∈ RT |An xn ⩽ bn } ,
22

Figure 3.2: Evolution of the convergence rate, given as λ1 (P ) (lowest nonzero eigenvalue of
P ), with N = 6 and T ∈ {4, 6, 8, 12, 20, 60}. The worst convergence rate is evaluated by taking
100 × T random draws of the sets Tn ⊂ [T ] for each n, and evaluating the eigenvalue of the
matrix. The slope is around -0.93, which indicates that in practice the convergence rate is
O(T −1 ), faster than the upper bound in O(T −3 ) established in Theorem 3.20.
with An ∈ Mkn ,T (R) with kn ∈ N. The disaggregation problem (1.2), with p ∈ P fixed, writes:

min 0

(4.2a)

x∈RN T

s.t. A0 x = Bp (λ0 )

(4.2b)

An xn ⩽ bn , ∀n ∈ [N ] (λn ) .

(4.2c)

where A0 = J1,N ⊗ IT , B = IT , (such thatP (4.2b) corresponds to the aggregation constraint
P
T
n kn
are the Lagrangian multipliers associated to
n xn = p) and λ0 ∈ R , (λn )n∈[N ] ∈ R+
(4.2b) and (4.2c).
With the polyhedral constraints (4.1), the graph representation of the disaggregation problem, as illustrated in Figure 2.1 is no longer valid. Consequently, one can not directly apply
Hoffman’s theorem (Theorem 2.3) to obtain a characterization of disaggregation feasibility by
inequalities on p. However, using duality theory, Proposition 4.1 below also gives a characterization of disaggregation:
Proposition 4.1. A profile p ∈ P is disaggregable iff:
(4.3)


P
>
∀ λ0 , λ1 , . . . λN ∈ Λ, λ>
0 Bp +
n∈[N ] λn bn ⩾ 0 ,

where
>
Λ := {λ0 ∈ Rk0 , ∀n ∈ [N ], λn ∈ Rk+n | A>
0 λ0 + A (λn )n = 0},

with A := diag(An )n∈[N ] .
23

Proof. From strong duality, we have:
(4.4)

min

max

x∈RN T λ0 ∈Rk0 ,λn ∈Rkn

X

λ>
0 (A0 x − Bp) +

n

+

max

(4.5)

=

λ>
n (An xn − bn )

−λ>
0 Bp −

X

λ0 ∈Rk0 ,λn ∈Rk+n
n
>
s.t. λ>
0 A0 + (λn )n A = 0

λ>
n bn

.

If the polytope Yp ∩ X given by the constraints
Q knof (4.2) is empty, then there is an infeasibility
>
>
>
>
T
certificate λ = (λ0 λ1 . . . λN ) ∈ R × n R+ such that:
>
>
>
λ>
0 A0 + (λn )n A = 0 and λ0 Bp + (λn )n b < 0 .

(4.6)

On the other hand, if Yp ∩X is nonempty, then a solution
P to>the dual problem (4.5) is bounded,
>
which implies that ∀λ := (λ0 , (λn )n ) ∈ Λ, λ0 Bp + n λn bn ⩾ 0 .
As opposed to Hoffman circulation’s theorem where disaggregation is characterized by a finite
number of inequalities, Proposition 4.1 involves a priori an infinite number of inequalities.
However, we know that the polyhedral cone Λ can be represented by a finite number of
generators (edges), that is, there exists Λ∗ := {λ∗(1) , . . . , λ∗(d) } such that:
nP
o
∗(i) | (α ) ∈ Rd
(4.7)
Λ=
i i
+ .
1⩽i⩽d αi λ
Thus, we obtain the following corollary to Proposition 4.1:
Corollary 4.2. There exists a finite set Λ∗ ⊂ Λ such that, for any p ∈ P, p is disaggregable
iff:
∀(λ0 , (λn )n ) ∈ Λ∗ , λ0 Bp +

(4.8)

P

n λn b n ⩾ 0 .

Remark 4.3. In the transportation case
P (2.2), we can write each agent constraints in the
form An xn ⩽ bn (writing the equality
t xn,t = En is written as two inequalities), and
Hoffman conditions (2.5) can be written in the form (4.8). Moreover, Theorem 3.3 ensures
that one possibility for Λ∗ of Corollary 4.2 is to consider the collection of 2T multipliers
corresponding to the subsets T ⊂ [T ] and N minimizing (2.6). We skip the details here for
brevity.
As in the first part of the paper, we want to use APM to decompose problem (1.1) and, in
the case where disaggregation is not possible, use the result of APM to generate an inequality
(4.3) violated by the current profile p.
In the case of impossible disaggregation, APM converges to the orbit (y ∞ , x∞ ), and µ :=
∞
∞
∞
y − x∞ defines a separating hyperplane x̄ + µ⊥ , where x̄ = y +x
, that satisfies, with
2
a := x̄.µ (note that x̄ can be replaced by any y ∈ [y ∞ , x∞ ]):
∀x ∈ Yp ; µ> x > a

∀x ∈ X ; a > µ> x,

which give lower bounds on the linear problems (the second one is decomposed because A is
a block-diagonal matrix, but it can also be written in one problem):
min µ> x

x∈RN T

(4.9)

=

A0 x = Bp (λ0 )
(4.10)

and

∀n ∈ [N ],

max µn xn

x∈RN T

An xn ⩽ bn (λn )
24

=

max −λ>
0 Bp

λ0 ∈Rk0

µ = −A>
0 λ0
minλn ∈Rkn b>
n λn
+

µn = λ>
n An

.

Strong duality on these problems implies that there exist λ0 and λ such that:
(4.11)

>
µ = −A>
0 λ0 and a < −λ0 Bp

µ = λ> A and a > b> λ .

>
Thus, we obtain (λ0 , λ) ∈ Λ satisfying (4.6), that is, λ>
0 Bp + b λ < 0, and we can use this to
add a new valid additional inequality on p of form (4.3), that will change the current profile
p:
>
λ>
0 Bp + λ b ⩾ 0

(4.12)

In Algorithm 4.1 below, we summarize the proposed decomposition of problem (1.1). This
is a generalization of the decomposition principle used for Algorithm 3.4.
Algorithm 4.1 Non-intrusive optimal disaggregation with polyhedral constraints
Require: Start with Λ(0) = {}, k = 0, Disag= false
1: while not Disag do
>
(k)
2:
get solution p(k) of problem minp∈P {f (p) | λ>
0 Bp + λ b ⩾ 0, ∀λ ∈ Λ }
3:
get µ(k) = y ∞ − x∞ ← APM(Yp(k) , X )
4:
if µ(k) 6= 0 then
(k)
(k) | µ(k) = −A> λ }
5:
obtain λ0 ← maxλ0 ∈Rk0 {−λ>
0 Bp
0 0
(k)

(k)

>
obtain for each n, λn ← maxλn ⩾0 {b>
n λn | µn = λn An }
(k)
7:
add Λ(k+1) = Λ(k) ∪ {(λ0 , λ(k) )}
8:
else
9:
Return Disag = true, p(k) as optimal solution
10:
end
11:
k ←k+1
12: done

6:

Remark 4.4. We use the fact that µ = y ∞ − x∞ although, as before, we only have an
approximation of this quantity. The approximation has to be precise enough to ensure that
>
the solution obtained verifies λ>
0 Bp + b λ < 0. In practice, one can proceed as in the
transportation case and Algorithm 3.3 use a large εcvg as stopping criteria in APM, then
>
compute (λ0 , λ) ∈ Λ and check if λ>
0 Bp + b λ < 0. If this is not the case, restart with
εcvg = εcvg /2.
P
Remark 4.5. When Yp = {x ∈ RN T |A0 x = B0 p} = {x| n xn = p}, we can obtain a
non-intrusive version of APM on Yp and X , similar to Algorithm 3.3. In this case, (4.11)
ensures that we have µn,t = −[λ0 ]t for each n ∈ [N ], and λ0 is fixed by µ. The only difference
with the transportation case for a non-intrusive APM in the general polyhedral case, is in the
way of computing the valid constraint violated by p. Thus, Lines 16 to 19 of Algorithm 3.3
have to be replaced by Algorithm 4.2.
Remark 4.6 (Link with Benders’ Decomposition). In this generalized case, we obtain an
algorithm related to Benders’ decomposition [8] (recall that in our specific case (4.2), the cost
function does not involve the variable x but only variable p).
The difference between the proposed Algorithm 4.1 and Benders’ decomposition lies in the
way of generating the new cut. Benders’ decomposition would directly solve the dual problem
25

Algorithm 4.2 Modification of Lines 16-19 of Algorithm 3.3 for NI-APM with polyhedral
constraints
16: for each agent n ∈ [N ] do
17:
compute Mn optimal value of (4.10).
18: done
19: Operator computes M ← SMCA((Mn )n )
20: if −ν.p + M < 0 then
21:
return Disag ← False, −ν, M
>
(4.5): maxλ {−λ>
0 B0 p − (λn )n b | λ0 A0 + (λn )n A = 0} and obtain a cut if it is unbounded.
However, this problem involves the constraints of all users (through A and b), and is it not
straightforward to obtain a method to solve this problem in a decentralized and efficient way.

5. Numerical Examples.
5.1. An Illustrative Example with T=4. In this section we illustrate the iterations of the
method proposed in this paper on an
Pwith T = 4 and N = 3. Assuming that we have
Pexample
to satisfy the aggregate constraint t pt = n En , we can use the projections on this affine
space of solutions of master problems (p(s) )s to visualize them in dimension 3.
One can wonder if, in the transportation case, applying Algorithm 3.4 or Algorithm 4.1
will always lead to the same cuts and solutions: the answer is no, as shown by the instance
considered in this section, for which Algorithm 3.4 converges in 3 iterations and Algorithm 4.1
needs 4 iterations.
We consider the problem (1.1) with agents constraints (2.2) with parameters x := 0 and:

(5.1)

[0.8, 0.2, 0.7, 0.1]
E1 = 1.8
x := [0.5, 0.1, 0.3, 0.6] , E2 = 0.4 ,
[0.1, 0.1, 0.7, 0.2]
E3 = 1.1

∀p ∈ R4 , f (p) :=

X

0.8 × pt + 0.1 × p2t .

1⩽t⩽4

P
P
Considering the aggregate equality constraint
1⩽t⩽4 pt =
1⩽n⩽3 En = 3.3, we use the
canonical projection of 4 dimensional vectors into the 3 dimensional space (p1 , p2 , p3 ) to visualize the cuts and solutions. In this example, there exist 2T − 2 = 14 nontrivial Hoffman
inequalities characterizing disaggregation from Theorem 2.3. The projection of the obtained
polytope PD , as defined in (2.1), is represented in Figure 5.1a. One can remark that this
polytope has only 6 facets. Our Algorithm 3.4 applied on this instance with εdis = 10−3 and
εcvg = 10−5 converges in 3 iterations, with successive solutions of the master problem (3.1)
and cuts added:
cut

p(1) = [1., 0.4, 1., 0.9]

−→

p(2) = [0.75, 0.4, 1.4, 0.75]

cut

−→

p1 + p2 + p4 ⩽ 1.9
p2 + p3 + p4 ⩽ 2.4

p(3) = [0.9, 0.4, 1.4, 0.6] .
Figure 5.1b represents in the projection space the three successive solutions and the two
generated cuts (in red for each iteration).
On the other hand, applying Algorithm 4.1 with the same precision parameters(εdis , εcvg ),
there are 3 cuts generated and 4 resolutions of master problem needed for convergence, given
26

by (we refer the reader to Remark 4.4 for the numerical precision obtained in the values):
cut

p̂(1) = [1., 0.4, 1., 0.9]

−→

p̂(2) = [0.8097, 0.4, 1.3984, 0.6919]

cut

−→
cut

p̂(3) = [0.9062, 0.4, 1.3823, 0.6115] −→
p̂

(4)

−0.25p1 − 0.25p2 + 1.0p3 − 0.5p4 ⩾ 0.75

1.0p1 − 0.509p2 + 0.018p3 − 0.509p4 ⩾ 0.4161
−0.333p1 − 0.333p2 + 1.0p3 − 0.333p4 ⩾ 0.7666

= [0.9, 0.4, 1.4, 0.6] .

The 4 successive solutions and the 3 added cuts are represented in the three dimensional space
on Figure 5.1c: we observe that the last cut needed to obtain the convergence of Algorithm 4.1,
corresponds to the first one added with Algorithm 3.4.
Due to the strict convexity of the cost function p 7→ f (p), the final solution obtained is
the same, unique aggregated optimal solution of (1.1). The 4 successive solutions and the 3

p̂(2) p

(2)

p(2)

p(2)

p(3)= p̂(4)
p̂(3)

p(3)

p3

p3
p(1)

p2

p(1)
p(1) = p̂(1)

p2
p1

p1

(a) Projected Polyhedron PD

(b) Solutions and cuts (projected) obtained with Algorithm 3.4

p̂(2)

p̂(2)

p̂(4) p̂(2)

p̂(3)

p3

p̂(3)

p3

p3
p2
p̂(1)

p1

p2
p̂(1)

p1

p̂(1)

p2
p1

(c) Solutions and cuts (projected) obtained with Algorithm 4.1

Figure 5.1: Illustration of the iterations of the proposed decomposition method. The cut
p3 ⩾ 1.4, which is added at first for Algorithm 3.4 , is only added at the third iteration of
Algorithm 4.1
added cuts are represented in the three dimensional space on Figure 5.1c: we observe that the
last cut needed to obtain the convergence of Algorithm 4.1, corresponds to the first one added
with Algorithm 3.4.
27

5.2. A Nonconvex Example: Management of a Microgrid. In this section, we illustrate
the proposed method on a larger scale practical example from energy. We consider an electricity microgrid [27] composed of N electricity consumers with flexible appliances (such as electric
vehicles or water heaters), a photovoltaic (PV) power plant and a conventional generator. The
operator of the microgrid aims at satisfying the demand constraints of consumers over a set of
time periods T = {1, . . . , T }, while minimizing the energy cost for the community. We have
the following characteristics:
– the PV plant generates a nondispatchable power profile (ppv
t )t∈[T ] at marginal cost zero;
– the conventional generator has a starting cost C st , minimal and maximal power production
pg , pg , and piecewise-linear and continuous generation cost function pg 7→ f (pg ):
f (pg ) = αk + ck pg , if pg ∈ Ik := [θk−1 , θk [, k = 1 . . . K, where θ0 := 0 and θK := pg ;
– each agent n ∈ [N ] has some flexible appliances which require a global energy demand En
on [T ], and has consumption constraints on the total household consumption, on each time
period t ∈ [T ], that are formulated with xn , xn . These parameters are confidential because
they could for instance contain some information on agent n habits.
The master problem (3.1) can be written as the following MILP (5.2):

X
X
g
on
st st
(5.2a)
min
α
b
+
c
p
+
C
b
1 t
k kt
t
g
p,pg ,(pk ),(bk ),bon ,bst

(5.2b)
(5.2c)
(5.2d)
(5.2e)
(5.2f)
(5.2g)

t∈[T ]

k

g
k=1 pk,t , ∀t ∈ [T ]
bk,t (θk − θk−1 ) ⩽ pgk,t ⩽ bk−1,t (θk − θk−1 ), ∀1 ⩽ k ⩽ K, ∀t ∈ [T ]
on
on
bst
t ⩾ bt − bt−1 , ∀t ∈ {2, . . . , T }
g
g on
pg bon
t ⩽ pt ⩽ p bt , ∀t ∈ T
st
bon
t , bt , b1,t , . . . , bK−1,t ∈ {0, 1}, ∀t ∈ [T ]
pv
g

pgt =

PK

p⩽p

+p

>

(5.2h)

p 1T = E > 1N

(5.2i)

x> 1N ⩽ p ⩽ x> 1N .

In this formulation (5.2b-5.2c), where b0,t := 1 and bK,t := 0, are a mixed integer formulation
of the generation cost function f . One can show that the Boolean variable bk,t is equal to one
iff pgt ⩾ θk for each k ∈ {1, . . . , K − 1}. Note that only α1 appears in (5.2a) because of the
continuity of f .
Constraints (5.2d-5.2e) ensure the on/off and starting constraints of the power plant, (5.2g)
ensures that the power allocated to consumption is not above the total production, and (5.2h5.2i) are the aggregated feasibility conditions already referred to in (2.3). The nonconvexity
of (5.2) comes from the existence of starting costs and constraints of minimal power, which
makes necessary to use Boolean state variables bst , bon .
We simulate the problem described above for different values of N ∈ {24 ,25 , 26 ,27 ,28 }
and one hundred instances with random parameters for each value of N . A scaling factor
κN = N/20 is applied on parameters to ensure that production capacity is large enough to
meet consumers demand. The parameters are chosen as follows:
– T = 24 (hours of a day);
28

– production costs: K = 3 , θ = [0, 70, 100, 300]κN , c = [0.2, 0.4, 0.5], pg = 50κN , pg = 300κN ,
α1 =4 and C st = 15; h
i
(t−6)2π
pv
– photovoltaic: ppv
t = 50(1−cos( 16 )+U([0, 10]) κN for t ∈ {6, . . . , 20}, pt = 0 otherwise;
– consumption parameters are drawn randomly with: xn,t ∼ U([0, 10]), xn,t ∼ U([0, 5]) + xn,t
and En ∼ U([1T> xn , 1T> xn ]), so that individual feasibility (Xn 6= ∅) is ensured.
N=

24

25

26

27

28

# master
# projs.

193.6
9507

194.1
15367

225.5
24319

210.9
26538

194.0
26646

Table 5.1: number of subproblems solved (average on 100 instances)
We implement Algorithm 3.4 using Python 3.5. The MILP (5.2) is solved using Cplex Studio
12.6 and Pyomo interface. Simulations are run on a single core of a cluster at 3GHz. For the
convergence criteria (see Lines 11 and
√we use εdis = 0.01 with the operator
P 12 of Algorithm 3.3),
norm defined by k|x|k = maxn∈[N ] t |xn,t | (to avoid the N factor in the convergence criteria
appearing with k.k2 ), and starts with εcvg = 0.1. The largest instances took around 10 minutes
to be solved in this configuration and without parallel implementation. As the CPU time
needed depends on the cluster load, it is not a reliable indicator of the influence on N on
the complexity of the problems. Moreover, one advantage of the proposed method is that the
projections in APM can be computed locally by each agent in parallel, which could not be
implemented here for practical reasons.
Table 5.1 gives the number of master problems solved and the total number of projections
computed, on average over the hundred instance for each value of N .
One observes that the number of master problems (5.2) solved (number of “cuts” added),
remains almost constant when N increases. In all instances, this number is way below the
upper bound of 224 > 1, 6 × 107 possible constraints (see Proposition 3.12), which suggests
that only a limited number of constraints are added in practice. The average total number
of projections computed for each instance (total number of iterations of the while loop of
Algorithm 3.3, Line 1 over all calls of APM in the instance) increases in a sublinear way which
is even better that one could expect from the upper bound given in Theorem 3.20.
6. Conclusion. We provided a non-intrusive algorithm that enables to compute an optimal
resource allocation, solution of a–possibly nonconvex–optimization problem, and affect to each
agent an individual profile satisfying a global demand and lower and upper bounds constraints.
Our method uses local projections and works in a distributed fashion. Hence, the resolution
of the problem is still efficient even in the case of a very large number of agents. The method
is also privacy-preserving, as agents do not need to reveal any information on their constraints
or their individual profile to a third party.
Several extensions and generalizations can be considered for this work. Section 4 generalizes the procedure to arbitrary polyhedral constraints for agents. However, the number
of constraints (cuts) added to the master problem is not proved to be finite as done in the
transportation case. Proving that only a finite number of constraints can be added (maybe up
to a refinement procedure of the current constraint obtained) will enable to have a termination
result for the algorithm in the general polyhedral case. In the transportation case, we showed
the geometric convergence of APM with a rate linear in the number of agents. Moreover, the
29

number of cuts added in the procedure is finite but the upper bound that we have remains
exponential. In practice however, the number of constraints to consider remains small, as seen
in Section 5. A thinner upper bound on the number of cuts added in the algorithm in this
case would constitute an interesting result.
Acknowledgments. We thank the referees for their comments, leading to improvements of
this paper. We thank Daniel Augot for very helpful discussions concerning secure multiparty
computations and cryptographic protocols.
Appendix A. Proof of Proposition 3.4.
Proof of Item (i). Let us write the stationarity conditions associated to problem (3.3):
(A.1)

∀n ∈ [N ], ∀t ∈ [T ],

0 = (xn,t − yn,t ) − λn − µn,t + µn,t

and

yn,t = xn,t + νt .

By summing the latter equalities on t ∈ [T ] and n ∈ [N ], we obtain the equalities:
X
X
X
(A.2)
νt =
yn,t − En , ∀n ∈ [N ]
pt =
xn,t + N νt , ∀t ∈ [T ] .
t

n

t

Moreover, by summing over t ∈ Tn◦ the first series of equalities in (A.1), and by using the
complementary slackness conditions, which entail that µn,t = µn,t = 0 for t ∈ Tn◦ , we get
(A.3)

X

|Tn◦ |λn = En −

xn,t −

t∈T n

X
t∈Tn◦

yn,t −

X

xn,t , ∀n ∈ [N ] ,

t∈T n

where we define for each n ∈ [N ]:
Tn◦ := {t | xn,t < xn,t < xn,t },

T n = {t | xn,t = xn,t } and T n = {t | xn,t = xn,t } .
P
P
P
From (A.2) and Assumption 2 giving the equality n En = t pt , we obtain t νt = 0 and:
P
(A.4)
∀n ∈ [N ],
t∈[T ] yn,t = En .
Let us show that ∀t ∈ T , νt > 0. For t ∈ T , there exists n s.t. yn,t > xn,t . From (A.1) we get:
νt = yn,t − xn,t ⩾ yn,t − xn,t > 0 .

(A.5)

Suppose by contradiction that there exists n ∈
/ N and t̂ ∈ T such that xn,t̂ < xn,t̂ . We obtain
from complementary slackness conditions and the first equality in (A.1):
(A.6)

xn,t̂ ⩾ yn,t̂ + λn = xn,t̂ + νt̂ + λn which implies that λn < 0 .

From this, we obtain T n ⊂ T : indeed, for t ∈ T n , we have yn,t + λn ⩾ xn,t which gives:
yn,t ⩾ xn,t − λn > xn,t which implies t ∈ T .
From the condition (A.4) and as νt > 0 for each t ∈ T n because T n ⊂ T and (A.5), we get:
X
X
X
X
0=
(yn,t − xn,t ) =
(yn,t − xn,t ) +
(−λn ) +
νt i.e.,
t∈[T ]

t∈Tn◦

t∈T n

30

t∈T n

X

(yn,t − xn,t ) =

X

λn −

t∈Tn◦

t∈T n

X

νt ,

t∈T n

which is strictly negative: this implies that there exists t0 ∈ T n such that yn,t0 < xn,t0 .
Necessarily, t0 ∈
/ T because
νt0 = yn,t0 − xn,t0 < xn,t0 − xn,t0 = 0. Then, as we have
P
P
0 = pt0 ⩾
y
x
m,t
m∈[N ]
m m,t0 , there exists m ∈ [N ] such that ym,t0 > xm,t0 . If λm ⩽ 0,
and as xm,t0 = ym,t0 − νt0 > xm,t0 , we get:
xm,t0 = min(xm,t0 , ym,t0 + λm ) ⩽ ym,t0 + λm ⩽ ym,t0 = xm,t0 + νt0 < xm,t0 ,
which is impossible, thus λm > 0. Now, we observe that Tn◦ ⊂ T . Indeed, otherwise, if
t” ∈ Tn◦ ∩ T c , we have νt” = −λn > 0 and xm,t” = ym,t” − νt” < ym,t” < xm,t” , thus we get:
xm,t” = max(xm,t” , ym,t” + λm ) ⩾ ym,t” + λm > ym,t” = xm,t” + νt” + λm > xm,t” ,
which is impossible, thus Tn◦ ⊂ T .
c
Finally, since T n 6= ∅, consider t0 ∈ arg mint∈T
/ n {xn,t − yn,t }. By (A.3), we obtain:
(A.7)

yn,t0 + λn < xn,t0 ⇐⇒ En −

X

xn,t −

X

yn,t −

t∈Tn◦

t∈T n

X

xn,t < |Tn◦ |(xn,t0 − yn,t0 )

t∈T n

and thus:
X
X
X
X
X
X
En −
xn,t = En −
xn,t −
xn,t
xn,t −
xn,t +
xn,t −
t∈T c

t∈T

t∈T n

t∈T n ∩T

t∈T n ∪Tn◦

t∈T n ∩T

(as T n ∪ Tn◦ ⊂ T )
(A.8)

<

X

(xn,t0 − yn,t0 ) − (xn,t − yn,t ) +

t∈Tn◦

(A.9)

⩽0

X

(xn,t − xn,t )

(from (A.7))

T n ∩T

(from the definition of t0 and xn,t ⩽ xn,t ),

which contradicts n ∈
/ N . Thus we have shown that ∀n ∈
/ N , ∀t ∈ T , xn,t = xn,t . From (A.1)
and as ∀t ∈ T , νt > 0 from (A.5), we obtain ∀n ∈
/ N , ∀t ∈ T , yn,t = xn,t + νt > xn,t . This
terminates the proof for Item (i).
Proof of Item (ii). Let us consider T 0 := {t|νt > 0}. We already showed that T ⊂ T 0 in
(A.5). To prove the other inclusion and obtain Item (ii), we observe that, considering n ∈
/N
(nonempty as shown independently in Item (v)) and considering T 0 instead of T , all the facts
established in the proof of Item (i) above also hold: by contradiction of ∀t ∈ T 0 , xn,t = xn,t ,
we first obtain λn < 0 as in (A.6). Then, as for t” ∈ Tn◦ , we have νt” = −λn > 0, we get
Tn◦ ∩ T 0c = ∅. Finally the same sequence of inequalities as (A.8-A.9) show a contradiction.
Consequently, for each t ∈ T 0 , xn,t = xn,t and yn,t = xn,t + νt > xn,t , thus t ∈ T and T 0 ⊂ T .
Proof of Item (iii). Suppose on the contrary that there exists n ∈ N such that λn ⩾ 0.
For t ∈ Tn◦ , we have νt = −λn ⩽ 0, thus, Tn◦ ⊂ T c . Then, if t ∈ T and if xn,t < xn,t , we would
have:
xn,t = max(xn,t , yn,t + λn ) ⩾ xn,t + 0 + νt > xn,t ,
31

which is impossible, thus xn,t = xn,t , and T ⊂ T n . As we show independently in Item (v) that
T 6= ∅, we know T n 6= ∅. Let us consider t0 ∈ arg mint∈T
/ n {yn,t − xn,t }. By (A.3), we obtain:
(A.10)

yn,t0 + λn > xn,t0 ⇐⇒ En −

xn,t −

En −

X

X

yn,t −

xn,t > |Tn◦ |(xn,t0 − yn,t0 )

t∈T n

X
X
X
X
X
X
xn,t −
xn,t = En −
xn,t −
xn,t −
xn,t +
xn,t
xn,t −

t∈T c

>

T

X

t∈T n

t∈Tn◦

t∈T c ∩T n

t∈T c ∩T

t∈T c ∩T n

t∈T n

(as T ⊂ T n )
X

(yn,t − xn,t ) − (yn,t0 − xn,t0 ) +
xn,t − xn,t

t∈Tn◦

(A.12)

X
t∈Tn◦

t∈T n

and thus:

(A.11)

X

(from (A.10))

n

(from the definition of t0 and xn,t ⩽ xn,t ),

⩾0

which contradicts n ∈ N and terminates the proof of Item (iii).
Proof of Item (iv). From (ii), we know that T c = {t|νt ⩽ 0}, thus, if t ∈
/ T and n ∈ N ,
if xn,t > xn,t then we would have xn,t ⩽ yn,t + λn = xn,t + νt + λn < xn,t , which is a
contradiction.
P
Proof of Item (v). From t νt = 0, we see that if T = ∅, then this means that νt = 0 for
all t ∈ [T ], and
P thus y = x which is a contradiction. Thus there exists t0 such that νt0 > 0
and because t∈T νt = 0, there exists t00 such that νt00 < 0.
P
If N = ∅, then using (i), we would have for all n, yn,t0 > xn,t0 and
P thus pt0 > n∈[N ] xn,t0 ,
which contradicts the aggregate upper bound constraint ∀t, pt ⩽ n∈[N ] xn,t .
c
0
0
P If N = ∅, then using (iv), we would have for all n, yn,t0 < xn,t00 and
P thus pt0 <
n∈[N ] xn,t00 , which contradicts the aggregate lower bound constraint ∀t, pt ⩾
n∈[N ] xn,t .
Appendix B. Proof of Lemma 3.10.
Proof of Item (i). From x(K) = PX (y (K−1) ) and y (K) = PY (x(K) ), we obtain, similarly
to (A.1):
(K)

(K−1)

(B.1) ∀n ∈ [N ], ∀t ∈ [T ], 0 = (xn,t − yn,t
(K)

(K)

(K)

(K)

(K)

(K)
) − λ(K)
n − µn,t + µn,t and yn,t = xn,t + νt
(K)

(K)

where the Lagrangian multipliers λn , µ(K)
, µn,t (resp. νt
n,t

.

) are associated to the quadratic

(y (K−1) ) (resp.

problem characterizing the projections PX
y (K) = PY (x(K) )). We obtain
equalities similar to (A.2, A.3). We proceed as for Proposition 3.4(i) and suppose that there
P
(K)
εcvg
exists n ∈
/ N and t̂ ∈ T such that xn,t̂ < xn,t̂ . Then, as ky (K) − y ∞ k ⩽ 1−ρ
and t∈[T ] yt =
P
(K)
εcvg
∞
∞
t∈[T ] yt , we have for each n ∈ [N ], t ∈ [T ], |yn,t − yn,t | ⩽ 2(1−ρ) , and thus we get:
(K)

(B.2)

(K−1)

xn,t̂ ⩾ xn,t̂ ⩾ yn,t̂
ε

ε

ε

cvg
cvg
∞
+ λ(K)
⩾ yn,
− 2(1−ρ)
+ λ(K)
= x∞
+ νt̂∞ − 2(1−ρ)
+ λ(K)
n
n
n
t̂
n,t̂

cvg
=⇒ λ(K)
< 2(1−ρ)
v − νt̂∞ <
n

Bεcvg
− 2Bεcvg = − 32 Bεcvg ,
2
(K)

as νt̂∞ ⩾ ν > 2Bεcvg . Let us now consider t0 ∈ Tn◦ (K) ∪ T n , then:
(B.3)
(K)
(K)
(K)
(K)
(K)
ε
εcvg
3
νt0 = yn,t0 − xn,t0 ⩾ yn,t0 − yn,t0 − λ(K)
> − cvg
n
2 + 2 Bεcvg > Bεcvg + 2 (B − 1) ⩾ Bεcvg ,
32

(K)

which shows that t0 ∈ T (K) = T ∞ and thus Tn◦ (K) ∪ T n ⊂ T . Then, the same sequence of
inequalities as (A.7, A.8, A.9) applied to y (K−1) gives a contradiction to n ∈
/ N.
Proof of Item (ii). The proof of Item (ii) is symmetric to the one of Item (i): if we suppose
(K)
that there exists n ∈ N and t̂ ∈
/ T such that xn,t̂ > xn,t̂ , we obtain, symmetrically to (B.2),
(K)

that λn

(K)

ε

cvg
⩾ − 2(1−ρ)
. Then, considering t0 ∈ T n

(K)

∪ Tn◦ (K) , we show, symmetrically to (B.3),

(K)

that νt0 < Bεcvg i.e. t0 ∈
/ T and thus T n ∪ Tn◦ (K) ⊂ T c . We conclude by obtaining a
contradiction to n ∈ N by the same sequence of inequalities as (A.10, A.11, A.12).
References.
[1] E. A. Abbe, A. E. Khandani, and A. W. Lo, Privacy-preserving methods for sharing financial
risk exposures, American Economic Review, 102 (2012), pp. 65–70.
[2] Y. P. Aneja and K. P. Nair, Bicriteria transportation problem, Management Science, 25
(1979), pp. 73–78.
[3] M. F. Anjos, A. Lodi, and M. Tanneau, A decentralized framework for the optimal coordination of distributed energy resources, IEEE Trans. Pow. Sys., 34 (2018), pp. 349–359.
[4] M. Atallah, M. Bykova, J. Li, K. Frikken, and M. Topkara, Private collaborative forecasting and benchmarking, in Proceedings of the 2004 ACM workshop on Privacy in the electronic
society, 2004, pp. 103–114.
[5] H. H. Bauschke and J. M. Borwein, On the convergence of Von Neumann’s alternating
projection algorithm for two sets, Set-Valued Analysis, 1 (1993), pp. 185–212.
, Dykstra’s alternating projection algorithm for two sets, Journal of Approximation Theory,
[6]
79 (1994), pp. 418–443.
[7] H. H. Bauschke, J. Chen, and X. Wang, A Bregman projection method for approximating
fixed points of quasi-Bregman nonexpansive mappings, Applicable Analysis, 94 (2015), pp. 75–84.
[8] J. F. Benders, Partitioning procedures for solving mixed variables programming problems, Numerische mathematik, 4 (1962), pp. 238–252.
[9] D. P. Bertsekas, Nonlinear programming, Athena Scientific, 1999.
[10] D. P. Bertsekas and J. N. Tsitsiklis, Parallel and distributed computation: numerical methods, vol. 23, Prentice hall Englewood Cliffs, NJ, 1989.
[11] E. D. Bolker, Transportation polytopes, Journal of Combinatorial Theory, Series B, 13 (1972),
pp. 251–262.
[12] J. M. Borwein, G. Li, and L. Yao, Analysis of the convergence rate for the cyclic projection
algorithm applied to basic semialgebraic convex sets, SIAM J. Optim., 24 (2014), pp. 498–527.
[13] C. Clifton, M. Kantarcioglu, J. Vaidya, X. Lin, and M. Y. Zhu, Tools for privacy
preserving distributed data mining, ACM Sigkdd Explorations Newsletter, 4 (2002), pp. 28–34.
[14] G. Cohen and D. L. Zhu, Decomposition and coordination methods in large scale optimization
problems: The nondifferentiable case and the use of augmented lagrangians, Adv. in Large Scale
Systems 1, (1984), pp. 203–266.
[15] W. J. Cook, W. Cunningham, W. Pulleyblank, and A. Schrijver, Combinatorial optimization, Springer, 2009.
[16] W. Deng, M.-J. Lai, Z. Peng, and W. Yin, Parallel multi-block admm with o (1/k) convergence, Journal of Scientific Computing, 71 (2017), pp. 712–736.
[17] R. L. Dykstra, An algorithm for restricted least squares regression, Journal of the American
Statistical Association, 78 (1983), pp. 837–842.
[18] R. Glowinski and A. Marroco, Sur l’approximation, par éléments finis d’ordre un, et la
résolution, par pénalisation-dualité d’une classe de problèmes de dirichlet non linéaires, ESAIM,
9 (1975), pp. 41–76.
[19] O. Goldreich, S. Micali, and A. Wigderson, How to play any mental game, in Proceedings
of the Nineteenth Annual ACM Symposium on Theory of Computing, STOC ’87, New York, NY,
USA, 1987, Association for Computing Machinery, p. 218–229.
33

[20] L. Gubin, B. T. Polyak, and E. Raik, The method of projections for finding the common
point of convex sets, USSR Comput. Math. & Math. Phys., 7 (1967), pp. 1 – 24.
[21] J. He, L. Cai, P. Cheng, J. Pan, and L. Shi, Consensus-based data-privacy preserving data
aggregation, IEEE Trans. Autom. Control, (2019), pp. 1–1.
[22] A. J. Hoffman, Some recent applications of the theory of linear inequalities to extremal combinatorial analysis, in Proc. of Symposia on Applied Mathematics, 1960, pp. 113–127.
[23] B. A. Huberman, E. Adar, and L. R. Fine, Valuating privacy, IEEE security & privacy, 3
(2005), pp. 22–25.
[24] P. Jacquot, O. Beaude, P. Benchimol, S. Gaubert, and N. Oudjane, A privacypreserving disaggregation algorithm for non-intrusive management of flexible energy, in IEEE
58th Conference on Decision and Control (CDC), IEEE, 2019.
[25] P. Jacquot, O. Beaude, S. Gaubert, and N. Oudjane, Analysis and implementation of an
hourly billing mechanism for demand response management, IEEE Trans. Smart Grid, 10 (2019),
pp. 4265–4278.
[26] G. Jagannathan, K. Pillaipakkamnatt, and R. N. Wright, A new privacy-preserving
distributed k-clustering algorithm, in Proc. of the 2006 SIAM Int. Conf. on Data Mining, SIAM,
2006, pp. 494–498.
[27] F. Katiraei, R. Iravani, N. Hatziargyriou, and A. Dimeas, Microgrids management, IEEE
power and energy magazine, 6 (2008).
[28] K. K. Lai, K. Lam, and W. K. Chan, Shipping container logistics and allocation, J. Oper.
Res. Soc., 46 (1995), pp. 687–697.
[29] C. Lemaréchal, A. Nemirovskii, and Y. Nesterov, New variants of bundle methods, Math.
Program., 69 (1995), pp. 111–147.
[30] P.-Y. R. Ma et al., A task allocation model for distributed computing systems, IEEE Trans.
Computers, 100 (1982), pp. 41–47.
[31] F. L. Müller, J. Szabó, O. Sundström, and J. Lygeros, Aggregation and disaggregation
of energetic flexibility from distributed energy resources, IEEE Trans. Smart Grid, (2017).
[32] J. Munkres, Algorithms for the assignment and transportation problems, SIAM J. App. Math.,
5 (1957), pp. 32–38.
[33] R. Nishihara, S. Jegelka, and M. I. Jordan, On the convergence rate of decomposable
submodular function minimization, in NIPS, 2014, pp. 640–648.
[34] D. P. Palomar and M. Chiang, A tutorial on decomposition methods for network utility
maximization, IEEE J. Sel. Areas Commun., 24 (2006), pp. 1439–1451.
[35] A. Rais and A. Viana, Operations research in healthcare: a survey, Int. Trans. Oper. Res., 18
(2011), pp. 1–31.
[36] M. Ruan, H. Gao, and Y. Wang, Secure and privacy-preserving consensus, IEEE Transactions
on Automatic Control, 64 (2019), pp. 4035–4049.
[37] K. Seong, M. Mohseni, and J. M. Cioffi, Optimal resource allocation for ofdma downlink
systems, in Information Theory, 2006 IEEE Int. Sym., IEEE, 2006, pp. 1394–1398.
[38] R.-h. Shi, Y. Mu, H. Zhong, J. Cui, and S. Zhang, Secure multiparty quantum computation
for summation and multiplication, Scientific reports, 6 (2016), pp. 1–9.
[39] J. Von Neumann, Functional operators: Measures and integrals, vol. 1, Princeton University
Press, 1950.
[40] L. Xiao and S. Boyd, Optimal scaling of a gradient method for distributed resource allocation,
J. Optim. Theory. Appl., 129 (2006), pp. 469–488.
[41] L. Xiao, M. Johansson, and S. P. Boyd, Simultaneous routing and resource allocation via
dual decomposition, IEEE Trans. Comm., 52 (2004), pp. 1136–1144.
[42] A. C. Yao, How to generate and exchange secrets, in 27th Annual Symp. Found. of Comp. Sci.
(SFCS), 1986, pp. 162–167.
[43] H. Yu and M. J. Neely, A simple parallel algorithm with an o(1/t) convergence rate for general
convex programs, SIAM J. Optim., 27 (2017), pp. 759–783.
[44] A. Zoha, A. Gluhak, M. A. Imran, and S. Rajasegarar, Non-intrusive load monitoring
34

approaches for disaggregated energy sensing: A survey, Sensors, 12 (2012), pp. 16838–16866.
[45] M. Zulhasnine, C. Huang, and A. Srinivasan, Efficient resource allocation for device-todevice communication underlaying lte network, in WiMob, 2010 IEEE 6th Int. Conference, IEEE,
2010, pp. 368–375.

35

