A Privacy-preserving Method to Optimize Distributed Resource Allocation

TL;DR

This paper introduces a privacy-preserving algorithm for distributed resource allocation that ensures agents' private information remains confidential while achieving optimal global resource distribution, applicable to large-scale, high-dimensional problems.

Contribution

The paper proposes a novel aggregation-based method combining polyhedral cuts and secure multiparty computations to optimize resource allocation without revealing private data.

Findings

Algorithm converges efficiently in high-dimensional settings.

Effective in preserving privacy during the allocation process.

Numerical results validate practical applicability and convergence.

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 adress 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.

Tables

Table 1. Table 1: number of subproblems solved (average on 100 instances)

$N=$	$2^4$	$2^5$	$2^6$	$2^7$	$2^8$
# master	193.6	194.1	225.5	210.9	194.0
# projs.	9507	15367	24319	26538	26646

Equations

$$\underset{\mathbf{x}\in {\mathbb{R}}^{N\times T},\mathbf{p}\in \mathcal{P}}{\min }f(\mathbf{p})$$

$${\mathbf{x}}_n\in {\mathcal{X}}_n,\forall n\in [N]$$

$$\sum _{n\in [N]}{x}_{n,t}=p_t,\forall t\in [T],$$

$$\text{textsc}Find\mathbf{x}\in {\mathcal{Y}}_{\mathbf{p}}\cap \mathcal{X}$$

$${\mathcal{Y}}_{\mathbf{p}}:=\{\mathbf{y}\in {\mathbb{R}}^{NT}|{\mathbf{y}}^{\top }\text{mathds}{1}_N=\mathbf{p}\}\text{hfill}\text{ and }\text{hfill}\mathcal{X}:=\prod _{n\in [N]}{\mathcal{X}}_n.$$

$$\underset{\mathbf{p}\in \mathcal{Q}}{\min }f(\mathbf{p}).$$

$$\mathcal{P}_{\mathrm{D}}:=\big{\{}\bm{p}\in\mathcal{P}\ |\ \exists\bm{x}\in\mathcal{X}\ ;\ \bm{p}=\textstyle\sum_{n}\bm{x}_{n}\big{\}}\ .$$

$${\mathcal{X}}_n=\{{\mathbf{x}}_n\in {\mathbb{R}}^T:{\sum }_{t\in [T]}{x}_{n,t}=E_n\text{ and }{\underline{x}}_{n,t}⩽x_{n,t}⩽{\overline{x}}_{n,t}\}\ne \varnothing .$$

$${\mathbf{p}}^{\top }\text{mathds}{1}_T={\mathbf{E}}^{\top }\text{mathds}{1}_N$$

$${\underline{\mathbf{x}}}^{\top }\text{mathds}{1}_N⩽\mathbf{p}⩽{\overline{\mathbf{x}}}^{\top }\text{mathds}{1}_N.$$

$$\forall \mathcal{A}\subset \mathcal{V},\sum _{e\in {\delta }^{+}(\mathcal{A})}{u}_e⩾\sum _{v\in \mathcal{A}}{d}_v+\sum _{e\in {\delta }^{+}({\mathcal{A}}^c)}{\ell }_e,$$

$$\forall \mathcal{T}\subset [T],\forall \mathcal{N}\subset [N],\sum _{t\in \mathcal{T}}{p}_t-\sum _{n\in \mathcal{N}}{E}_n+\sum _{t\notin \mathcal{T},n\in \mathcal{N}}{\underline{x}}_{n,t}⩽\sum _{t\in \mathcal{T},n\notin \mathcal{N}}{\overline{x}}_{n,t}.$$

$$\sum _{t\in \mathcal{T}}{p}_t⩽\underset{\mathcal{N}\subset [N]}{\min }\left\{\sum _{n\in \mathcal{N}}{E}_n-\sum _{t\notin \mathcal{T},n\in \mathcal{N}}{\underline{x}}_{n,t}+\sum _{t\in \mathcal{T},n\notin \mathcal{N}}{\overline{x}}_{n,t}\right\},$$

$$\mathcal{X}:={\mathcal{X}}_1\times \cdots \times {\mathcal{X}}_N\text{ where }{\mathcal{X}}_n:=\left\{{\mathbf{x}}_n\in {\mathbb{R}}_{+}^T|{\sum }_{t\in [T]}{x}_{n,t}=E_n\text{ and }\forall t,{\underline{x}}_{n,t}⩽x_{n,t}⩽{\overline{x}}_{n,t}\right\}.$$

$$\displaystyle\mathcal{Y}_{\bm{p}}:=\big{\{}\bm{x}\in\mathbb{R}^{NT}\ |\ \forall t\in[T],\textstyle\sum_{n\in[N]}x_{n,t}=p_{t}\big{\}}\ .$$

$$\mathcal{P}={\mathcal{P}}^{(0)}\supset {\mathcal{P}}^{(1)}\supset \cdots \supset {\mathcal{P}}^{(S)}\supset {\mathcal{P}}_{\mathrm{D}}.$$

$$\underset{\mathbf{p}\in {\mathbb{R}}^T}{\min }f(\mathbf{p})$$

$$\mathbf{p}\in {\mathcal{P}}^{(s)}.$$

$${\mathbf{x}}^{(k)}\underset{k\to \infty }{⟶}{\mathbf{x}}^{\infty },{\mathbf{y}}^{(k)}\underset{k\to \infty }{⟶}{\mathbf{y}}^{\infty };$$

$$\Vert {\mathbf{x}}^{\infty }-{\mathbf{y}}^{\infty }{\Vert }_2=\underset{\mathbf{x}\in \mathcal{X},\mathbf{y}\in \mathcal{Y}}{\min }\Vert \mathbf{x}-\mathbf{y}{\Vert }_2.$$

$$\underset{\mathbf{x},\mathbf{y}}{\min }\frac{1}{2}\Vert \mathbf{x}-\mathbf{y}{\Vert }_2^2$$

$$\forall n\in [N],\sum _{t\in [T]}{x}_{n,t}=E_n$$

$$\forall n\in [N],\forall t\in [T],{\underline{x}}_{n,t}⩽x_{n,t}⩽{\overline{x}}_{n,t}$$

$$\forall t\in [T],\sum _{n\in [N]}{y}_{n,t}=p_t$$

$$\begin{array}{rl}\mathcal{L}(\mathbf{x},\mathbf{y},\mathbf{\lambda },\mathbf{\mu },\mathbf{\nu })& =\frac{1}{2}\Vert \mathbf{x}-\mathbf{y}{\Vert }_2^2-{\mathbf{\lambda }}^{\top }{\left({\sum }_t{x}_{n,t}-E_n\right)}_n-{\underline{\mathbf{\mu }}}^{\top }(\mathbf{x}-\underline{\mathbf{x}})\\ & -{\overline{\mathbf{\mu }}}^{\top }(\overline{\mathbf{x}}-\mathbf{x})-{\mathbf{\nu }}^{\top }\left({\sum }_n{\mathbf{y}}_n-\mathbf{p}\right).\end{array}$$

$$\forall n\in [N],\forall t\in [T],y_{n,t}=x_{n,t}+{\nu }_t.$$

$$\mathcal{T}$$

$$\mathcal{N}$$

$$\sum_{t\in\mathcal{T}}p_{t}\leqslant\min_{\mathcal{N}\subset[N]}\left\{\sum_{n\in\mathcal{N}}\Big{(}E_{n}-\sum_{t\notin\mathcal{T}}\underline{x}_{n,t}-\sum_{t\in\mathcal{T}}\overline{x}_{n,t}\Big{)}\right\}+\sum_{t\in\mathcal{T},n\in[N]}\overline{x}_{n,t}.$$

$$\sum _{n\in \mathcal{N}}{E}_n-\sum _{t\in \mathcal{T}}{p}_t+\sum _{t\in \mathcal{T},n\notin \mathcal{N}}{\overline{x}}_{n,t}-\sum _{t\notin \mathcal{T},n\in \mathcal{N}}{\underline{x}}_{n,t}<0.$$

$$A_{\mathcal{T}}({\mathbf{x}}^{\infty })<\sum _{t\in \mathcal{T}}{p}_t\text{ with }{A}_{\mathcal{T}}({\mathbf{x}}^{\infty }):=\sum _{t\in \mathcal{T}}\sum _{n\in [N]}{x}_{n,t}^{\infty }.$$

$$\sum _{n\in \mathcal{N}}{E}_n+\sum _{t\in \mathcal{T},n\notin \mathcal{N}}{\overline{x}}_{n,t}-\sum _{t\notin \mathcal{T},n\in \mathcal{N}}{\underline{x}}_{n,t}-\sum _{t\in \mathcal{T}}{p}_t$$

$$=\sum _{n\in \mathcal{N}}\sum _t{x}_{n,t}+\sum _{t\in \mathcal{T},n\notin \mathcal{N}}{\overline{x}}_{n,t}-\sum _{t\notin \mathcal{T},n\in \mathcal{N}}{\underline{x}}_{n,t}-\sum _{t\in \mathcal{T}}{p}_t(\text{from }\text{eqref}eq:minSquare:consE)$$

$$\displaystyle=\sum_{n\in\mathcal{N}}\Big{(}\sum_{t\notin\mathcal{T}}\underline{x}_{n,t}+\sum_{t\in\mathcal{T}}x_{n,t}\Big{)}+\sum_{t\in\mathcal{T},n\notin\mathcal{N}}x_{n,t}-\sum_{t\notin\mathcal{T},n\in\mathcal{N}}\underline{x}_{n,t}-\sum_{t\in\mathcal{T}}p_{t}\;\text{ (from Prop.\leavevmode\nobreak\ \ref{prop:factsOncut} \ref{claimFactcut1} and \ref{claimFactcut3})}$$

$$\displaystyle=\sum_{t\in\mathcal{T}}\Big{(}\sum_{n\in\mathcal{N}}x_{n,t}+\sum_{n\notin\mathcal{N}}x_{n,t}\Big{)}-\sum_{t\in\mathcal{T}}p_{t}=\sum_{t\in\mathcal{T}}\Big{(}\sum_{n\in[N]}x_{n,t}-\sum_{n\in[N]}y_{n,t}\Big{)}$$

$$\displaystyle=\sum_{t\in\mathcal{T}}\Big{(}\sum_{n\in[N]}-\nu_{t}\Big{)}=\sum_{t\in\mathcal{T}}\Big{(}-\sum_{n\in[N]}\left|x_{n,t}-y_{n,t}\right|\Big{)}$$

$$\sum _{t\in [T]}\sum _{n\in [N]}(x_{n,t}-y_{n,t})=\sum _{n\in [N]}{E}_n-\sum _{t\in [T]}{p}_t=0,\text{vspace}-0.2cm$$

Full text

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) $\bm{p}=(p_{t})_{t\in[T]}$, where $[T]:=\{1,\dots,T\}$ is the set of time periods, such that $\bm{p}\in\mathcal{P}$ (feasibility constraints of power allocation, induced by the power plants constraints), while minimizing the cost $f(\bm{p})$ (representing the production and operations costs of the centrally controlled power plants). Besides, this aggregate allocation has to match an individual consumption profile $\bm{x}_{n}=(x_{n,t})_{t\in[T]}$ for each of the consumer (agent) $n\in[N]$ considered. The problem can be written as follows:

[TABLE]

The (aggregate) allocation $\bm{p}$ can be made public, that is, revealed to all agents. However, the individual constraint set $\mathcal{X}_{n}$ and individual profiles $\bm{x}_{n}$ 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) as the combination of two interdependent subproblems:

1. i)

given an aggregate allocation $\bm{p}$, the disaggregation problem consists in finding, for each agent $n$, an individual profile $\bm{x}_{n}$ satisfying her individual constraint (1b), so that constraint (1c) is satisfied. This is equivalent to:

[TABLE]

When (2) has a solution, we say that a disaggregation exists for $\bm{p}$; 2. ii)

For a given subset $\mathcal{Q}\subset\mathcal{P}$, we define the master problem,

[TABLE]

When $\mathcal{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).

Aside from the example above, resource allocation problems (optimizing common resources shared by multiple agents) with the same structure as (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\in[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 $\bm{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.

1.2 Main Results

The main contribution of the paper is to provide a non-intrusive and distributed algorithm (Algorithm 4) that computes an aggregate resource allocation $\bm{p}$, optimal solution of the—possibly nonconvex—optimization problem (1), along with feasible individual profiles $\bm{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 $\min_{\bm{p}\in\mathcal{P}^{(s)}}f(\bm{p})$, obtained by constructing a decreasing sequence of successive approximations $\mathcal{P}^{(s)}$ of the set of aggregate consumptions $\bm{p}$ for which a disaggregation exists (see (1)), starting from $\mathcal{P}^{(0)}=\mathcal{P}$. At each step, we reduce the set $\mathcal{P}^{(s)}$ by incorporating a new constraint on $\bm{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 (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 $\bm{x}\in\mathbb{R}^{NT}$ aggregating to a given resource allocation $\bm{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 $\bm{p}$ of the global problem (1), but violated from the current resource allocation (Theorem 3.3). Adding this cut to the master problem (3) by updating $\mathcal{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) 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.7 for more details).

The problem of the aggregation of constraints has been studied in the field of energy, in 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 Section 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 Section 3.2, we show how to guarantee the privacy of the proposed procedure by using secure multiparty computation techniques. In Section 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: Section 5.1 gives an illustrative toy example in dimension $T=4$, while in Section 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 “$\bm{x}$”, are used to denote vectors, while normal fonts, like “$x$”, refer to a scalar quantities. We denote by $\bm{v}^{\top}$ the transpose of a vector $\bm{v}$. Recall that we denote by $[N]$ the set $\{1,\dots,N\}$. Calligraphic letters such as $\mathcal{T},\mathcal{N},\mathcal{X}$ are used to denote sets, and if $\mathcal{T}\subset[T]$, $\mathcal{T}^{c}:=\{t\in[T]\setminus\mathcal{T}\}$ denotes the complementary set of $\mathcal{T}$. The notation $\mathcal{U}([a,b])$ stands for the uniform distribution on $[a,b]$. The notation $P_{\mathcal{C}}(.)$ refers to the Euclidean projection onto a convex set $\mathcal{C}$. For $d\in\mathbb{N},\mathds{1}_{d}$ denotes the vector of ones $(1\dots 1)^{\top}\in\mathbb{R}^{d}$.

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 $\bm{p}\mapsto f(\bm{p})$, where $\bm{p}\in\mathbb{R}^{T}$ is the aggregate allocation of $T$ dimensional resources (for example power production over $T$ time periods). This resource allocation $\bm{p}$ is to be shared between $N$ individual agents, each agent obtaining a part $\bm{x}_{n}\in\mathcal{X}_{n}$, where $\mathcal{X}_{n}$ denotes the individual feasibility set of agent $n$.

The global problem the operator wants to solve is described in (1). We assume that in Problem (1), the constraints set $\mathcal{X}_{n}$ and individual profile $\bm{x}_{n}$ are confidential to agent $n$ and should not be disclosed to the central operator or to another agent.

Let us define the set $\mathcal{P}_{\mathrm{D}}$ of feasible aggregate allocations that are disaggregable as:

[TABLE]

Feasibility of problem (1) is equivalent to having $\mathcal{P}_{\mathrm{D}}$ 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\in[N]$, there exists fixed parameters $E_{n}>0$, $\bm{\overline{x}}_{n}\in\mathbb{R}^{T}$, $\bm{\underline{x}}_{n}\in\mathbb{R}^{T}$ such that :

[TABLE]

In particular, $\mathcal{X}_{n}$ is convex and compact. Given an allocation $\bm{p}$, the structure obtained on the matrix $(\bm{x}_{n,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 $\mathcal{X}_{n}$ is a polyhedron.

Given a particular allocation $\bm{p}\in\mathcal{P}$, the operator will be interested to know if this allocation is disaggregable, that is, if there exists individual profiles $(\bm{x}_{n})_{n\in[N]}\in\prod_{n}\mathcal{X}_{n}$ summing to $\bm{p}$, or equivalently if the disaggregation problem (2) has a solution.

Following (2), the disaggregate profile refers to $\bm{x}$, while the aggregate profile refers to the allocation $\bm{p}$. Problem (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:

[TABLE]

Because they are necessary, we assume that those aggregate conditions (6) hold for vectors of $\mathcal{P}$, that is:

Assumption 2**.**

*All vectors $\bm{p}\in\mathcal{P}$ satisfy (6), that is, $\mathcal{P}\subset\{\bm{p}\in\mathbb{R}^{T}\ |\ \eqref{eq:agg-cond}\text{ hold }\}$. *

However, conditions (6) 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 $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with vertex set $\mathcal{V}$, edge set $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$, demands $d:\mathcal{V}\rightarrow\mathbb{R}$ (where $d_{v}<0$ means that $v$ is a production node), edge lower capacities $\ell:\mathcal{E}\rightarrow\mathbb{R}_{+}$ and upper capacities $u:\mathcal{E}\rightarrow\mathbb{R}_{+}$. A flow on $\mathcal{G}$ is a function $\bm{x}:\mathcal{E}\rightarrow\mathbb{R}_{+}$ such that $\bm{x}$ satisfies the capacity constraints, that is $\forall e\in\mathcal{E},\ell_{e}\leqslant x_{e}\leqslant u_{e}$, and Kirchoff’s law, that is, $\forall v\in\mathcal{V},\sum_{e\in\delta_{v}^{+}}x_{e}=d_{v}+\sum_{e\in\delta_{v}^{-}}x_{e}$, where $\delta^{+}_{v}$ (resp. $\delta^{-}_{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 $\mathcal{G}$ whose set of vertices is the disjoint union $\mathcal{V}=[T]\sqcup[N]$ and whose set of edges is $\mathcal{E}=\{(t,n)\}_{t\in[T],n\in[N]}$. Define demands on nodes $[T]$ by $d_{t}=-p_{t}$ and demands on nodes $[N]$ by $d_{n}=E_{n}$. Assign to each edge $(t,n)$ an upper capacity $u_{n,t}=\overline{x}_{n,t}$ and lower capacity $\ell_{n,t}=\underline{x}_{n,t}$. Then, finding a solution $\bm{x}$ to (2) is equivalent to finding a feasible flow in $\mathcal{G}$.

Hoffman [22] gave a necessary and sufficient condition for the flow problem to be feasible in a graph with balanced demands, that is $d(\mathcal{V}):=\sum_{{v\in\mathcal{V}}}d_{v}=0$ (total production matches total positive demand, in our case (6a)). 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 $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with demand $d\in\mathbb{R}^{\mathcal{V}}$ such that $d(\mathcal{V})=0$ and capacities $\ell\in(\mathbb{R}\cup\{-\infty\})^{\mathcal{E}}$ and $u\in(\mathbb{R}\cup{\infty})^{\mathcal{E}}$ with $\ell\leqslant u$, there exists a feasible flow $\bm{x}\in\mathcal{E}\rightarrow\mathbb{R}_{+}$ on $\mathcal{G}$ if and only if:

[TABLE]

where $\delta_{+}(\mathcal{A}):=\{(u,v)\in\mathcal{E}|u\in{\mathcal{A}^{c}},v\in\mathcal{A}\}$ is the set of edges coming to set $\mathcal{A}$ and ${\mathcal{A}^{c}}:=\mathcal{V}\setminus\mathcal{A}$.

Proposition 2.4 translates Theorem 2.3 in our framework:

Proposition 2.4**.**

Given an aggregate resource allocation $\bm{p}=(p_{t})_{t\in[T]}\in\mathcal{P}$, the disaggregation problem is feasible, meaning that $\bm{p}\in\mathcal{P}_{\mathrm{D}}$, iff:

[TABLE]

Proof 2.5**.**

We apply (7) with $\mathcal{A}:=\mathcal{T}^{c}\cup\mathcal{N}^{c}$ and use the equality $d(\mathcal{V})=0=\sum_{v\in\mathcal{A}}d_{v}+\sum_{v\in\mathcal{A}^{c}}d_{v}$.

From Theorem 2.3 or Proposition 2.4 above, one can see that the aggregate constraints (6) are in general not sufficient to ensure that the disaggregation problem has a solution.

For a given set $\mathcal{T}$, there is a choice of $\mathcal{N}$ which leads to the strongest inequality (8), namely:

[TABLE]

In this way, we get $2^{T}-2$ inequalities corresponding to the proper subsets $\mathcal{T}\subset[T]$. Moreover, in general, these $2^{T}-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) while preserving the privacy of the agent’s constraints $\mathcal{X}_{n}$ and individual profile $\bm{x}_{n}\in\mathbb{R}^{T}$. 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 (2).

Let us consider the polyhedron enforcing the agents constraints:

[TABLE]

Besides, given an allocation $\bm{p}\in\mathcal{P}$, we consider the set of profiles aggregating to $\bm{p}$ :

[TABLE]

Note that $\mathcal{Y}_{\bm{p}}$ is an affine subspace of $\mathbb{R}^{NT}$ (to be distinguished from $\mathcal{P}$ which is a subset of $\mathbb{R}^{T}$), and that $\mathcal{Y}_{\bm{p}}\cap\mathcal{X}$ is empty iff $\bm{p}\notin\mathcal{P}_{\mathrm{D}}$, according to the definition of $\mathcal{P}_{\mathrm{D}}$ in (4). The idea of the proposed algorithm is to build a finite sequence of decreasing subsets $(\mathcal{P}^{(s)})_{0\leqslant s\leqslant S}$ such that:

[TABLE]

At each iteration, a new aggregate resource allocation $\bm{p}^{(s)}$ is obtained by solving an instance of the master problem introduced in (3) with $\mathcal{Q}=\mathcal{P}^{(s)}$:

[TABLE]

In the sequel, we will refer to (10) as the master problem. Our procedure relies on the following immediate observation:

Proposition 3.1**.**

If $\bm{p}^{(s)}$ is a solution of (10) and $\bm{x}\in\mathcal{Y}_{\bm{p}^{(s)}}\cap\mathcal{X}$, then $(\bm{p}^{(s)},\bm{x})$ is an optimal solution of the initial problem (1). ■

Having in hands a solution $\bm{p}^{(s)}$, we can check whether $\mathcal{Y}_{\bm{p}^{(s)}}\cap\mathcal{X}\neq\emptyset$, and then find a vector $\bm{x}$ in this set, using APM on $\mathcal{X}$ and $\mathcal{Y}_{\bm{p}^{(s)}}$, as described in Algorithm 1 below (where $\mathcal{Y}=\mathcal{Y}_{\bm{p}}$).

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 $\mathcal{X}$ and $\mathcal{Y}$ be two closed convex sets with $\mathcal{X}$ bounded, and let $({\bm{x}}^{(k)})_{k}$ and $({\bm{y}}^{(k)})_{k}$ be the two infinite sequences generated by the APM on $\mathcal{X}$ and $\mathcal{Y}$ (Algorithm 1) with $\varepsilon_{\rm{cvg}}=0$. Then there exists $\bm{x}^{\infty}\in\mathcal{X}$ and $\bm{y}^{\infty}\in\mathcal{Y}$ such that:

[TABLE]

In particular, if $\mathcal{X}\cap\mathcal{Y}\neq\emptyset$, then $({\bm{x}}^{(k)})_{k}$ and $({\bm{y}}^{(k)})_{k}$ converge to a same point $\bm{x}^{\infty}\in\mathcal{X}\cap\mathcal{Y}$.

The convergence theorem is illustrated in Figure 2 in the case where $\mathcal{X}\cap\mathcal{Y}=\emptyset$, that is, when the disaggregation problem (2) is not feasible.

The idea of the algorithm proposed in this paper is, when $\mathcal{Y}_{\bm{p}^{(s)}}\cap\mathcal{X}=\emptyset$, to use the resulting vectors ${{\bm{x}}^{\infty}}$ and ${{\bm{y}}^{\infty}}$ to construct a new subset $\mathcal{P}^{(s+1)}$ by adding a constraint of type (8) to $\mathcal{P}^{(s)}$: indeed, from Proposition 2.4, we know that, if $\mathcal{Y}_{\bm{p}^{(s)}}\cap\mathcal{X}=\emptyset$, there exists at least one violated inequality (8).

The difficulty is to guess such a violated inequality among the $2^{T}$ possible inequalities. It turns out that, using the output of APM, we can build such an inequality.

Suppose that we obtain $\bm{x}^{\infty}\neq\bm{y}^{\infty}$ as defined in Theorem 3.2: we get a periodic cycle of APM, that is, we have $\bm{x}^{\infty}=P_{\mathcal{X}}(\bm{y}^{\infty})$ and $\bm{y}^{\infty}=P_{\mathcal{Y}}(\bm{x}^{\infty})$, and the couple $(\bm{x}^{\infty},\bm{y}^{\infty})$ is the solution of the following optimization problem, with given parameters $(E_{n})_{n}$ and $(p_{t})_{t}$.

[TABLE]

where $\lambda_{n}\in\mathbb{R}$, $\underline{\mu}_{n,t},\overline{\mu}_{n,t}\in\mathbb{R}_{+}$ and $\nu_{t}\in\mathbb{R}$ are the Lagrangian multipliers associated to the constraints (12b),(12c),(12d), with the associated Lagrangian function:

[TABLE]

The stationarity condition of the Lagrangian with respect to the variable $y_{n,t}$ yields:

[TABLE]

Let us consider the sets $\mathcal{T}\subset[T]$ and $\mathcal{N}\subset[N]$ defined from the output $(\bm{x}^{\infty},\bm{y}^{\infty})$ of APM on $\mathcal{X}$ and $\mathcal{Y}_{\bm{p}}$ as follows:

[TABLE]

In Theorem 3.3 below, we show that applying the inequality (8) with the sets $\mathcal{T}$ and $\mathcal{N}$ defined in (14) defines a valid inequality for the disaggregation problem violated by the current allocation $\bm{p}$.

The intuition behind the definition of $\mathcal{T}$ and $\mathcal{N}$ in (14) is the following: $\mathcal{T}$ is the subset of resources for which there is an over supply (which overcomes the upper bound for at least one agent). Once $\mathcal{T}$ is defined, $\mathcal{N}$ is the associated subset of $[N]$ minimizing the right hand side of (9). Indeed, (9) can be rewritten as:

[TABLE]

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 $({\bm{x}}^{(k)},{\bm{y}}^{(k)})_{k\in\mathbb{N}}$ generated by APM on $\mathcal{X}$ and $\mathcal{Y}_{\bm{p}}$ (see Algorithm 1). Then one of the following holds:

1. (i)

if $\mathcal{X}\cap\mathcal{Y}_{\bm{p}}\neq\emptyset$, then ${\bm{x}}^{(k)},{\bm{y}}^{(k)}\underset{k\rightarrow\infty}{\longrightarrow}\bm{x}^{\infty}\in\mathcal{X}\cap\mathcal{Y}_{\bm{p}}$; 2. (ii)

otherwise, if $\mathcal{X}\cap\mathcal{Y}_{\bm{p}}=\emptyset$, then ${\bm{x}}^{(k)}\underset{k\rightarrow\infty}{\longrightarrow}\bm{x}^{\infty}\in\mathcal{X}$ and ${\bm{y}}^{(k)}\underset{k\rightarrow\infty}{\longrightarrow}{{\bm{y}}^{\infty}}\in\mathcal{Y}_{\bm{p}}$. Then, considering the sets $\mathcal{T}$ and $\mathcal{N}$ defined in (14) gives an inequality of Hoffman (8) violated by $\bm{p}$, that is:

[TABLE]

Moreover, this Hoffman inequality can be written as a function of ${{\bm{x}}^{\infty}}$, as follows:

[TABLE]

Before giving the proof of Theorem 3.3, we need to show some technical properties on the sets $\mathcal{T}$, $\mathcal{N}$. For simplicity of notation, we use $\bm{x}$ and $\bm{y}$ to denote $\bm{x}^{\infty}$ and $\bm{y}^{\infty}$ in Proposition 3.4 and the proof of Theorem 3.3.

Proposition 3.4**.**

With $\bm{x}\neq\bm{y}$ solutions of problem (12) (outputs of APM on $\mathcal{X}$ and $\mathcal{Y}_{\bm{p}}$ with $\varepsilon_{\rm{cvg}}=0$),

1. (i)

$\forall t\in\mathcal{T},\forall n\notin\mathcal{N},\ {x}_{n,t}=\overline{x}_{n,t}$ and ${y}_{n,t}>\overline{x}_{n,t}$; 2. (ii)

$\mathcal{T}=\{t\ |\ \nu_{t}>0\}=\{t\ |\ p_{t}>\sum_{n}x_{n,t}\}$, where $\nu_{t}$ is the optimal Lagrangian multiplier associated to (12d); 3. (iii)

$\forall n\in\mathcal{N},\lambda_{n}<0$ ; 4. (iv)

$\forall t\notin\mathcal{T},\forall n\in\mathcal{N},\ {x}_{n,t}=\underline{x}_{n,t}$  ; 5. (v)

the sets $\mathcal{T},\mathcal{T}^{c}$, $\mathcal{N}$ and $\mathcal{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 3.5** (Proof of Theorem 3.3).**

We have:

[TABLE]

using the stationarity conditions (13) and for all $t\in\mathcal{T},\nu_{t}>0$ by Prop.3.4 (ii). Moreover, using:

[TABLE]

obtained from 2, we see that:

[TABLE]

which shows (15). We now show that inequality (16) is obtained by a rewriting of (15), indeed:

[TABLE]

Remark 3.6**.**

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 $\mathcal{X}$ and $\mathcal{Y}$ converge to two distinct points $\bm{x}^{\infty}$ and $\bm{y}^{\infty}$. Then, at each round $k$, we can define from (13) and considering any $n\in[N]$, the multiplier $\bm{\nu}^{(k)}=\bm{y}_{n}^{(k)}-\bm{x}_{n}^{(k)}$ tends to $\bm{y}_{n}^{\infty}-\bm{x}^{\infty}_{n}:=\bm{\nu}^{\infty}$. The set $\mathcal{T}$ of Theorem 3.3 is:

[TABLE]

which raises an issue for practical computation, as $\bm{\nu}^{\infty}$ is only obtained ultimately by APM, possibly in infinite time. To have access to $\mathcal{T}^{\infty}$ in finite time, that is, from one of the iterates $(\bm{\nu}^{(k)})_{k}$, we consider the set:

[TABLE]

where $\varepsilon_{\rm{cvg}}$ is the tolerance for convergence of APM as defined in Algorithm 1, $B>0$ is a constant, and $K$ (depending on $\varepsilon_{\rm{cvg}}$) is the first integer such that $\|\bm{x}^{(K)}-\bm{x}^{(K-1)}\|<\varepsilon_{\rm{cvg}}$.

We next show that we can choose $B$ to ensure that $\mathcal{T}^{(K)}=\mathcal{T}^{\infty}$ for $\varepsilon_{\rm{cvg}}$ small enough. We rely on the geometric convergence rate of APM on polyhedra [12, 33]:

Proposition 3.7** ([33]).**

If $\mathcal{X}$ and $\mathcal{Y}$ are polyhedra, there exists $\rho\in(0,1)$ such that the sequences $(\bm{x}^{(k)})_{k}$ and $(\bm{y}^{(k)})_{k}$ generated by APM verify for all $k\geqslant 1$:

[TABLE]

Proposition 3.7 applies to any polyhedra $\mathcal{X}$ and $\mathcal{Y}$. In Section 3.4 we shall give an explicit upper bound on the constant $\rho$ in the specific transportation case given by (1c) and (5).

From the previous proposition, we can quantify the distance to the limits in terms of $\rho$:

Lemma 3.8**.**

Consider an integer $K$ such that the sequence $(\bm{x}^{(k)})_{k\geqslant 0}$ generated by APM satisfies $\|\bm{x}^{(K)}-\bm{x}^{(K-1)}\|\leqslant\varepsilon_{\rm{cvg}}$, then we have for any $K^{\prime}\geqslant K-1$:

[TABLE]

Proof 3.9**.**

From Proposition 3.7, we have for any $k\geqslant K$:

[TABLE]

so that, by taking the limit $k\rightarrow\infty$, one obtains $\|{{\bm{x}}^{\infty}}-{\bm{x}}^{(K^{\prime})}\|\leqslant\tfrac{\varepsilon_{\rm{cvg}}}{1-\rho}$.

With this previous lemma, we can state the condition on $B$ ensuring the desired property:

Proposition 3.10**.**

Define $\underline{\nu}:=\min\{|\nu^{\infty}_{t}|>0\}$ (least nonzero element of $\bm{\nu}^{\infty}$). If the constants $B$ and $\varepsilon_{\rm{cvg}}>0$ are chosen such that $B>\frac{1}{1-\rho}\text{ and }\ \varepsilon_{\rm{cvg}}\times 2B<\underline{\nu},$ and Algorithm 1 stops at iteration $K$, then we have:

[TABLE]

Proof 3.11**.**

Let $t\in\mathcal{T}^{\infty}$, that is $\nu_{t}^{\infty}>0$ which is equivalent to $\nu_{t}^{\infty}\geqslant\underline{\nu}$ by definition of $\underline{\nu}$. We have:

[TABLE]

and this last quantity is greater than $B\varepsilon_{\rm{cvg}}$ as soon as $B\geqslant\tfrac{1}{1-\rho}$, thus $t\in\mathcal{T}^{(K)}$.

Conversely, if $t\in\mathcal{T}^{(K)}$, then:

[TABLE]

so that $t\in\mathcal{T}^{\infty}$. Furthermore, the “approximated” cut $\sum_{t\in\mathcal{T}}\left(\sum_{n\in[N]}x^{(K)}_{n,t}-p_{t}\right)\geqslant 0$ is violated by the current value of $\bm{p}$ (or $\bm{p}^{(s)}$ at iteration $s$) in the algorithm as:

[TABLE]

using (17) and (18). This last quantity is negative as soon as $\|{\bm{x}}^{(K)}-{{\bm{x}}^{\infty}}\|_{1}<\frac{1}{2}\|{{\bm{x}}^{\infty}}-{{\bm{y}}^{\infty}}\|_{1}$, which holds in particular if $B\varepsilon_{\rm{cvg}}<\frac{1}{2}\|{{\bm{x}}^{\infty}}-{{\bm{y}}^{\infty}}\|_{1}$.

This second proposition shows a surprising result: even if we do not have access to the limit $\bm{x}^{\infty}$, we can compute in finite time the exact left hand side term $A_{\mathcal{T}}(\bm{x}^{\infty})$ of the cut (16):

Proposition 3.12**.**

Under the hypotheses of Proposition 3.10, we have:

[TABLE]

Proof 3.13**.**

We start by showing some technical properties similar to Proposition 3.4:

Lemma 3.14**.**

The iterate $\bm{x}^{(K)}$ satisfies the following properties:

1. (i)

$\forall t\in\mathcal{T},\forall n\notin\mathcal{N},{x}_{n,t}^{(K)}=x^{\infty}_{n,t}=\overline{x}_{n,t}$ ; 2. (ii)

$\forall t\notin\mathcal{T},\forall n\in\mathcal{N},\ {x}_{n,t}^{(K)}=x^{\infty}_{n,t}=\underline{x}_{n,t}$  .

The proof of Lemma 3.14 is similar to Proposition 3.4 and is given in Appendix B. Then, having in mind that $\mathcal{T}^{(K)}=\mathcal{T}^{\infty}$ from Proposition 3.10, and $\mathcal{N}$ is obtained from $\mathcal{T}^{\infty}$ by (14), we obtain:

[TABLE]

which equals to $A_{\mathcal{T}}(\bm{x}^{\infty})$ as we have $\sum_{t\in[T]}x_{n,t}^{(K)}=E_{n}$ for each $n\in[N]$.

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.

3.2 Privacy-Preserving Projections through SMC

APM, as described in Algorithm 1, enables a distributed implementation in our context, by the structure of the algorithm itself: the operator computes the projection on $\mathcal{Y}_{\bm{p}}$ while each agent $n$ can compute, possibly in parallel, the projection on $\mathcal{X}_{n}$ 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 ${\bm{x}}$ to the operator.

For this, we use the fact that $\mathcal{Y}_{\bm{p}}$ is an affine subspace and thus the projection on $\mathcal{Y}_{\bm{p}}$ can be obtained explicitly component-wise. Indeed, as $\mathcal{P}_{\mathcal{Y}_{\bm{p}}}(\bm{x})$ is the solution $\bm{y}$ of the quadratic program $\min_{\bm{y}\in\mathbb{R}^{NT}|\sum_{n}\bm{y}_{n}=\bm{p}}\tfrac{1}{2}\|\bm{x}-\bm{y}\|_{2}^{2}$, we obtain optimality conditions similar to (13), and summing these equalities on $[N]$, we get $N\nu_{t}=\bm{p}_{t}-\sum_{n\in[N]}x_{n,t}$ for each $t\in[T]$, and thus:

[TABLE]

Thus, having access to the aggregate profile $\SS:=\sum_{n\in[N]}\bm{x}_{n}$, each agent can compute locally the component of the projection on $\mathcal{Y}_{\bm{p}}$ 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 $\SS$ can be computed in a non-intrusive manner and by several communications between agents and the operator, as described in Algorithm 2. The main idea of the SMC summation protocol [4, 38] below is that, instead of sending her profile ${\bm{x}}_{n}$, agent $n$ splits $x_{n,t}$ for each $t$ into $N$ random parts $(s_{n,t,m})_{m}$, according to an uniform distribution and summing to ${\bm{x}}_{n,t}$ (Lines 3-4). Thus, each part $s_{n,t,m}$ taken individually does not reveal any information on ${\bm{x}}_{n}$ nor on $\mathcal{X}_{n}$, 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 $\bm{\sigma}_{n}$ (Line 8), which does not contain either any information about any of the agents, and sends it to the operator (Line 9). The operator can finally compute the quantity $\SS={\bm{x}}^{\top}\mathds{1}_{N}=\bm{\sigma}^{\top}\mathds{1}_{N}$.

Remark 3.15**.**

A drawback of the proposed SMC summation (Algorithm 2) is that each agent need to exchange with all other nodes, therefore the communication overhead for each iteration and each node is $\mathcal{O}(N\ell)$ where $\ell$ 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 Section 3.3). However, as detailed in [4, Protocol I] the SMC summation can be adapted by splitting secret numbers among only $k\in\{1,\dots,N\}$ members. In that case, the communication overhead would be reduced to $\mathcal{O}(k\ell)$, 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 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.10) and in a privacy-preserving way using SMC.

To choose $B$ and $\varepsilon_{\rm{cvg}}$ satisfying the conditions of Proposition 3.10 a priori, one has to know the value of $\underline{\nu}$. 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 $\varepsilon_{\rm{cvg}}$. Indeed, one can start with $\varepsilon_{\rm{cvg}}$ arbitrary large so that APM will converge quickly, and then check if the cut obtained is violated by the current value of $\bm{p}$ (19): if it is not the case, we can continue the iterations of APM with convergence precision improved to $\varepsilon_{\rm{cvg}}/2$ (22). Proposition 3.10 ensures that this loop terminates in finite time.

The parameter $\varepsilon_{\rm{dis}}>0$ (Line 13 of Algorithm 3) has to be chosen a priori by the operator, depending on the precision required. In general in APM, $\bm{x}^{\infty}=\bm{y}^{\infty}$ will only be achieved in infinite time, so choosing $\varepsilon_{\rm{dis}}$ strictly positive is required.

We end this section by summarizing in Algorithm 4 the global iterative procedure to compute an optimal and disaggregable resource allocation $\bm{p}$, solution of the initial problem (1), using iteratively NI-APM (Algorithm 3) and adding constraints as stated in Theorem 3.3.

Algorithm 4 iteratively calls NI-APM (Algorithm 3) and in case disaggregation is not possible (Line 12), a new constraint is added (Line 14), obtained from the quantity $A_{\mathcal{T}}$ defined in (16), to the feasible set of resource allocations $\mathcal{P}^{(s)}$ in problem (10). This constraint is an inequality on $\bm{p}$ and thus does not reveal significant individual information to the operator. The algorithm stops when disaggregation is possible (Line 10). The termination of Algorithm 4 is ensured by the following property and the form of the constraints added (15):

Proposition 3.16**.**

Algorithm 4 stops after a finite number of iterations, as at most $2^{T}-2$ constraints (Line 14) can be added to the master problem (Line 3).

The following Proposition 3.17 shows the correctness of our Algorithm 4.

Proposition 3.17**.**

Let $B$ and $\varepsilon_{\rm{cvg}}$ satisfy the conditions of Proposition 3.10 and 3.12. Then:

1. (i)

if the problem (1) has no solution, Algorithm 4 exits at Line 5 after at most $2^{T}-2$ iterations; 2. (ii)

otherwise, Algorithm 4 computes, after at most $s\leqslant 2^{T}-2$ iterations, an aggregate solution $\bm{p}^{(s)}\in\mathcal{P}$, associated to individual profiles $(\bm{x}^{*})_{n}=$ NI-APM($\bm{p}^{(s)}$) such that:

[TABLE]

where $f^{*}$ is the optimal value of problem (1).

Proof 3.18**.**

The proof is immediate from Theorem 3.3, Proposition 3.10 and Proposition 3.12.

Remark 3.19**.**

The upper bound on the number of constraints added has no dependence on $N$ because, as stated in (9), 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 $\bm{p}$. Thus, Algorithm 4 is a method which enables the operator to compute a resource allocation $\bm{p}$ and the $N$ agents to adopt profiles $(\bm{x}_{n})_{n}$, such that $(\bm{x},\bm{p})$ solves the global problem (1), and the method ensures that both agent constraints (upper bounds $\bm{\overline{x}}_{n}$, lower bounds $\bm{\underline{x}}_{n}$, demand $E_{n}$); and disaggregate (individual) profile $\bm{x}_{n}$ (as well as the iterates $({\bm{x}}^{(k)})_{k}$ and $({\bm{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\neq n$).

Remark 3.20**.**

A natural approach to address problem (1) in a distributed way, assuming that both the cost function $\bm{p}\mapsto f(\bm{p})$ and the feasibility set $\mathcal{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 techniques, where Lagrangian multipliers associated to the (relaxed) aggregation constraint $\sum_{n}\bm{x}_{n}=\bm{p}$ would be updated using the SMC technique as described in Algorithm 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 4 proposed in this paper is that convergence is ensured (see Proposition 3.16) even if the cost function $\bm{p}\mapsto f(\bm{p})$ and the feasibility set $\mathcal{P}$ are not convex, which is the case in many practical situations (see Section 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 4, the operator obtains the following information for each iteration $s\in\{1,\dots,s^{\text{cv}}\}$ (where $s^{\text{cv}}<\infty$ is the number of iterations needed for convergence):

• •

the set $\mathcal{T}^{(s)}$ and the scalar $A_{\mathcal{T}}^{(s)}$ defining the cut;

• •

whenever ni-apm is called, the serie $(S^{(s,k)})_{k=1}^{K_{s}}$ of the aggregated projections.

Let us denote by \mathbb{I}^{\text{op}}:=\big{\{}(S^{(0,k)})_{k=1}^{K_{0}},\mathcal{T}^{(0)},A_{\mathcal{T}}^{(0)},\dots,(S^{(s^{\text{cv}},k)})_{k=1}^{K_{s^{\text{cv}}}},\mathcal{T}^{(s^{\text{cv}})},A_{\mathcal{T}}^{(s^{\text{cv}})}\big{\}} the information 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.21**.**

Consider two different instances of the resource allocation problem, differing only by a permutation of the agents, meaning that for all $n\in[N]$, $\mathcal{X}_{n}$ is replaced by $\mathcal{X}_{\sigma(n)}$, for some permutation $\sigma$ of the elements of $[N]$. Then, the information $\mathbb{I}^{\text{op}}$ gotten by the operator (as well as the information $(\bm{\nu}^{(k)})_{k}$ sent by the operator to agents in Algorithm 3), and the sequence of cuts added in the master problem (Algorithm 4), are precisely the same in both instances.

Proof 3.22**.**

The only information transmitted by the agents to the central operator consists of the aggregate profiles $\SS^{(k)}=$ SMCA$(\bm{x}^{(k)})$ and the sums $A_{\mathcal{T}}$ (line 5 and line 17 of Algorithm 3), for different values of the subset $\mathcal{T}$, corresponding to different Hoffman cuts. Each aggregate profile $\SS^{(k)}$ is invariant by permutation of the agents. Each sum $A_{\mathcal{T}}$ is of the form $A_{\mathcal{T}}=\sum_{n\in[N],t\in\mathcal{T}}x_{nt}^{(K)}$, 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.21, 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 $A_{\mathcal{T}}$, for different values of $\mathcal{T}\subset[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.23**.**

With $N>1$, a malicious (or semi-honest) operator using Algorithm 4 cannot infer the constraints $\mathcal{X}_{n}$ or profile $\bm{x}_{n}$ of a particular agent $n\in[N]$ with probability 1.

Remark 3.24**.**

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 $(E_{1},\overline{x}_{1,1},\overline{x}_{1,2}):=(1,2,0)$, $(E_{2},\overline{x}_{2,1},\overline{x}_{2,2}):=(3,1,2)$ and $\bm{\underline{x}}:=0$ (all unknown to the operator). The operator knows a priori the aggregate quantities $E_{1}+E_{2}=4,\ \bm{\overline{x}}_{1}+\bm{\overline{x}}_{2}=(3,2),\ \bm{\underline{x}}_{1}+\bm{\underline{x}}_{2}=(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 $(E^{*}_{1},\overline{x}^{*}_{1,1},\overline{x}^{*}_{1,2},\underline{x}^{*}_{1,1},\underline{x}^{*}_{1,2}):=(3,2,2,0,0)$.

Now suppose that the first aggregate profile proposed by the operator is $\bm{p}^{(0)}:=(3,1)$. From Algorithm 4, he obtains the cut defined by $\mathcal{T}^{(0)}=\{1\}$ and $A_{\{1\}}^{(0)}=2$. Thus, assuming the operator knows that $A_{\{1\}}^{(0)}$ is of the form, $A_{\{1\}}^{(0)}=\min_{\mathcal{N}\subset[N]}\{\sum_{n\in\mathcal{N}}E_{n}+\sum_{t\in\mathcal{T},n\notin\mathcal{N}}\overline{x}_{n,t}-\sum_{t\notin\mathcal{T},n\in\mathcal{N}}\underline{x}_{n,t}\}$ for a subset $\mathcal{N}$, he infers the following conditions on the parameters:

[TABLE]

With these conditions, the parameters values $(E^{*}_{1},\overline{x}^{*}_{1,1},\overline{x}^{*}_{1,2},\underline{x}^{*}_{1,1},\underline{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 $(E^{*}_{2},\overline{x}^{*}_{2,1},\overline{x}^{*}_{2,2},\underline{x}^{*}_{2,1},\underline{x}^{*}_{2,2}):=(1,1,0,0,0)$, which is incompatible with (26).

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.24 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 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 $\underline{x}_{n,t}$ and $\bar{x}_{n,t}$ are encoded by fixed (say double) precision numbers. Then, the total number of bits needed to encode the instance is $O(NT)$. However, it follows from Proposition 3.17 that the number of bits received by a malicious central operator is $O(2^{T})$, hence, in the regime $N\gg 2^{T}$ (large number of agents), the information leak is asymptotically negligible. Moreover, we give in Section 5.2 a real example with $N=2^{8}$ 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, $NT+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 $\bm{x}_{n}$ of agent $n$ has to intercept all the communications between $n$ and all $N-1$ other agents (to learn $(s_{n,t,m})_{m\neq n}$ and $(s_{m,t,n})_{m\neq n}$) and to the operator (to learn $\sigma_{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 2, the messages $s_{n,t,j}$ sent by agent $n$ to agent $j$, encoding shares of the consumption $x_{nt}$ 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 $\sigma_{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 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 2 exist to compute the sum $\sum_{n\in[N]}y_{n}$, while keeping $y_{n}$ 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 4: indeed, the only additional information that the collusion obtains is composed of the sequences $\{\bm{\nu}^{(k)}\}$ obtained from each execution of Algorithm 3. Thus, the symmetry argument given in Proposition 3.21 applies to give a privacy guarantee similar to Corollary 3.23:

Corollary 3.25**.**

Using Algorithm 4 with $N>1$, a collusion of less than $(N-1)$ malicious (or semi-honest) agents cannot infer the constraints $\mathcal{X}_{n}$ or profile $\bm{x}_{n}$ 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 4 will not converge to the global optimum. However, we still have the privacy guarantee given in Corollary 3.23: 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 1 on the sets $\mathcal{X}$ and $\mathcal{Y}_{\bm{p}}$ defined in Section 2.

A general result in [12] gives an upper bound of the sequences generated by APM on $\mathcal{X}$ and $\mathcal{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 $\mathcal{X}$ and $\mathcal{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 $\mathcal{X}$ and $\mathcal{Y}_{\bm{p}}$ defined in (5) and (2b):

Theorem 3.26**.**

For the two sets $\mathcal{X}$ and $\mathcal{Y}_{\bm{p}}$, the sequence of alternate projections converges to ${\bm{x}}^{*}\in\mathcal{X}$, ${\bm{y}}^{*}\in\mathcal{X}^{P}$ satisfying $\|{\bm{x}}^{*}-{\bm{y}}^{*}\|=\inf_{\bm{x}\in\mathcal{X},\bm{y}\in\mathcal{Y}_{\bm{p}}}\|\bm{x}-\bm{y}\|$, at the geometric rate:

[TABLE]

and the analogous inequalities hold for $({\bm{y}}^{(k)})_{k}$.

For the remaining of this section, we will just use $\mathcal{Y}$ to denote $\mathcal{Y}_{\bm{p}}$, as $\bm{p}$ remains fixed during APM. For the result stated in Theorem 3.26 above, we use several lemmas from [33].

Proof 3.27**.**

First, we use the fact stated in [33] that APM on subspaces $U$ and $V$ converge with geometric rate $c_{F}(U,V)^{2}$, where the rate is given by the square of the cosine of the Friedrichs angle between $U$ and $V$, given by:

[TABLE]

An intuitive generalization of this result for polyhedra $\mathcal{X}$ and $\mathcal{Y}$, considering all affine subspaces supporting the faces of $\mathcal{X}$ and $\mathcal{Y}$ is given in [33]:

Lemma 3.28** ([33]).**

For APM on polyhedra $\mathcal{X}$ and $\mathcal{Y}$ in $\mathbb{R}^{D}$, the convergence is geometric with rate bounded by the square of the maximal cosine of Friedrichs angle between subspaces supporting faces of $\mathcal{X}$ and $\mathcal{Y}$:

[TABLE]

where, for any $\bm{x}\in\mathbb{R}^{D}$, $\mathcal{X}_{\bm{x}}:=\arg\max_{\bm{v}\in\mathcal{X}}\bm{x}^{\top}\bm{v}$ is the face of $\mathcal{X}$ generated by direction $\bm{x}$ and $\operatorname{aff}_{0}(C)=\operatorname{aff}(C)-\bm{c}$ for some $\bm{c}\in C$ denotes the subspace supporting the affine hull of $C$, for $C=\mathcal{X}_{\bm{x}}$ or $C=\mathcal{Y}_{\bm{y}}$.

In the remaining of the proof, we bound the quantity (27) for our polyhedra $\mathcal{X}$ and $\mathcal{Y}$.

For this, we use the space $\mathbb{R}^{NT}=\mathbb{R}^{T}\times\dots\times\mathbb{R}^{T}$, where the $(n-1)T+1$ to $nT$ entries correspond to the profile of agent $n$, for $1\leqslant n\leqslant 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.29** ([33]).**

If $A$ and $B$ are matrices with orthonormal rows with same number of columns, then:

• – if all the singular values of $AB^{\top}$ are equal to one, then c_{F}\big{(}\operatorname{Ker}A,\ \operatorname{Ker}B\big{)}=0;
• – otherwise, c_{F}\big{(}\operatorname{Ker}A,\ \operatorname{Ker}B\big{)} is equal to the largest singular value of $AB^{\top}$ among those that are smaller than one.

We are left with finding a matrix representation of the faces of polyhedra $\mathcal{X}$ and $\mathcal{Y}$ and, then, bounding the corresponding singular values.

In our case, the polyhedra $\mathcal{Y}$ is an affine subspace $\mathcal{Y}=\{\bm{x}\in\mathbb{R}^{NT}\ |A\bm{x}=\sqrt{N}^{-1}\bm{p}\}$ where:

[TABLE]

where $\otimes$ denotes the Kronecker product. The matrix $A$ has orthonormal rows and the linear subspace associated to $\mathcal{Y}$ is equal to $\operatorname{Ker}(A)$.

Obtaining a matrix representation of the faces of $\mathcal{X}$ is more complex. The faces of $\mathcal{X}$ are obtained by considering, for each $n\in[N]$, subsets of the time periods that are at lower or upper bound (respectively $\underline{\mathcal{T}}_{n}$ and $\overline{\mathcal{T}}_{n}$, with $\underline{\mathcal{T}}_{n}\cap\overline{\mathcal{T}}_{n}=\emptyset$). Considering a collection of such subsets, a face of $\mathcal{X}$ can be written as:

[TABLE]

For some particular collection of subsets $(\overline{\mathcal{T}}_{n},\underline{\mathcal{T}}_{n})_{n}$, the set $\mathcal{A}_{(\overline{\mathcal{T}}_{n},\underline{\mathcal{T}}_{n})_{n}}$ might be empty. The linear subspace associated to $\mathcal{A}_{(\overline{\mathcal{T}}_{n},\underline{\mathcal{T}}_{n})_{n}}$ is given by $\{\bm{x}\in\mathbb{R}^{NT}|\ B\bm{x}=0\}=\operatorname{Ker}(B)$, where the $N$ first rows of $B$, corresponding to the constraints $\sum_{t}x_{n,t}=E_{n}$, are given before orthonormalization by:

[TABLE]

and the matrix $B$ has $b:=\sum_{n}|\mathcal{T}_{n}|$ more rows, where $\mathcal{T}_{n}:=\underline{\mathcal{T}_{n}}\cup\overline{\mathcal{T}_{n}}$, corresponding to the saturated bounds. Each of this row is given by the unit vector $\bm{e}_{n,t}^{\top}\in\mathbb{R}^{NT}$ for $n\in[N]$, $t\in\mathcal{T}_{n}$, which gives already an orthonormalized family of (unit) vectors. Therefore, a simple orthonormalized matrix $B\in\mathcal{M}_{N+b,NT}(\mathbb{R})$ giving the direction of $\mathcal{A}_{(\overline{\mathcal{T}}_{n},\underline{\mathcal{T}}_{n})_{n}}$ is given by:

[TABLE]

where $\bm{1}_{\mathcal{T}_{n}^{c}}\in\mathbb{R}^{T}$ is the vector where the indices in $\mathcal{T}_{n}^{c}$ are equal to 1 and 0 otherwise, and $B_{\mathcal{T}_{n}}:=\sum_{\begin{subarray}{c}1\leqslant k\leqslant|\mathcal{T}_{n}|\\ \mathcal{T}_{n}=\{t_{1},\dots,t_{|\mathcal{T}_{n}|}\}\end{subarray}}E_{kt_{k}}$ is the matrix $|\mathcal{T}_{n}|\times T$ with indices of $\mathcal{T}_{n}$. We obtain the double product:

[TABLE]

We observe that:

• – if $t_{0}\in\bigcap_{i=1}^{N}\mathcal{T}_{n}$, then $\bm{e}_{t_{0}}$ is an eigenvector associated to eigenvalue $\lambda_{t_{0}}=1$ ;
• – the vector $\bm{1}_{\bar{\mathcal{T}}}:=(\mathds{1}_{t\notin\cap_{n}\mathcal{T}_{n}})_{t\in[T]}\in\mathbb{R}^{T}$, where $\bar{\mathcal{T}}:=\cup_{n}\mathcal{T}_{n}^{c}$, is an eigenvector associated to eigenvalue $\lambda=1$. Indeed, if we denote by $\mathcal{N}_{\theta}=\{n\in[N]|\theta\in\mathcal{T}_{n}\}$, then $[\bm{1}_{\bar{\mathcal{T}}}]_{\theta}=1\Leftrightarrow\mathcal{N}_{\theta}^{c}\neq\emptyset$, and for each $\theta\in\bar{\mathcal{T}}$:

[TABLE]

To bound the other eigenvalues of the matrix $(AB^{\top})(BA^{\top})$, we rely on spectral graph theory arguments. Consider the weighted graph $\mathcal{G}=([T],\mathcal{E})$ whose vertices are the time periods $[T]$ and each edge $(k,\ell)\in[T]\times[T]$ with $k\neq\ell$ has weight $S_{k,\ell}=\tfrac{1}{N}\sum_{n}\tfrac{\mathds{1}_{\{k,\ell\}\subset\mathcal{T}_{n}^{c}}}{T-|\mathcal{T}_{n}|}$ (if this quantity is zero, then there is no edge between $k$ and $\ell$).

The matrix $P:=I_{T}-(AB^{\top})(BA^{\top})$ verifies for each $k\in[T]$:

[TABLE]

which shows that $P$ is the Laplacian matrix of graph $\mathcal{G}$. As Sp$(AB^{\top}BA^{\top})=1-$Sp$(P)$, we want to have a lower bound on the least eigenvalue of $P$ greater than 0, that we denote by $\lambda_{1}$.

By rearranging the indices of $[T]$ in two blocs $\bar{\mathcal{T}}$ and $\bar{\mathcal{T}}^{c}$, we observe that $P$ can be written as a block diagonal matrix $P=\text{diag}(P_{\bar{\mathcal{T}}},0_{\bar{\mathcal{T}}^{c}})$. As we are only interested in the positive eigenvalues of $\mathcal{P}$, we can therefore study the linear application associated to $P$ restricted to the subspace $\operatorname{Vect}(e_{t})_{t\in\bar{\mathcal{T}}}$.

As $\bm{1}_{\bar{\mathcal{T}}}$ is an eigenvector of $P$ associated to $\lambda_{0}=0$, from the minmax theorem, we have:

[TABLE]

Let us consider an eigenvector $u$ realizing (30). Let $u_{t^{*}}:=\max_{t}u_{t}$ and $u_{s^{*}}:=\min_{t}u_{t}$ and let $d_{s^{*},t^{*}}$ be the distance between $s^{*}$ and $t^{*}$ in $\mathcal{G}$, and let $(s^{*}\text{-}t^{*})$ denote a shortest path from $s^{*}$ to $t^{*}$ in $\mathcal{G}$. As $P$ is a Laplacian matrix, we have:

[TABLE]

where the last inequality is obtained from Cauchy-Schwartz inequality.

Let us write the path $(s^{*}\text{-}t^{*})=(t_{0},t_{1},\dots,t_{d})$. As $(s^{*}\text{-}t^{*})$ is a shortest path, for each $k\in\{0,d-1\}$, the edge $(t_{k},t_{k+1})$ exists so there exists $n\in[N]$ such that $\{t_{k},t_{k+1}\}\subset\mathcal{T}_{n}^{c}$. Moreover, for each $n$, we have $\mathcal{T}_{n}^{c}\cap\{t_{0},\dots,t_{k-1},t_{k+2},\dots,t_{d}\}=\emptyset$, otherwise we could “shortcut” the path $(s^{*}\text{-}t^{*})$, thus we have $|\mathcal{T}_{n}|\geqslant d-1$. We obtain:

[TABLE]

On the other hand, we have $(u_{t^{*}}-u_{s^{*}})\geqslant u_{t^{*}}+\tfrac{u_{t^{*}}}{T-1}=\tfrac{T}{T-1}u_{t^{*}}\geqslant\tfrac{T}{(T-1)\sqrt{T}}\|u\|_{2}$.

Using these bounds and (31), we obtain:

[TABLE]

Therefore, $\lambda_{1}\geqslant\frac{4}{N(T+1)^{2}(T-1)}:=\kappa_{N,T}$ and the greatest singular value lower than one of $(AB^{\top})(BA^{\top})$ is $1-\kappa_{N,T}$. We conclude by applying successively Lemma 3.29 and Lemma 3.28, to obtain the convergence rate stated in Theorem 3.26.

4 Generalization to Polyhedral Agents Constraints

In this section, we extend our results to a more general framework where for each $n\in[N]$, $\mathcal{X}_{n}$ is an arbitrary polyhedron, instead of having the particular structure given in (5). Let us now consider that $(\mathcal{X}_{n})_{n}$ are polyhedra with, for each $n$:

[TABLE]

with $A_{n}\in\mathcal{M}_{k_{n},T}(\mathbb{R})$ with $k_{n}\in\mathbb{N}$. The disaggregation problem (2), with $\bm{p}\in\mathcal{P}$ fixed, writes:

[TABLE]

where $A_{0}=J_{1,N}\otimes I_{T}$, $B=I_{T}$, (such that (33b) corresponds to the aggregation constraint $\sum_{n}\bm{x}_{n}=\bm{p}$) and $\bm{\lambda}_{0}\in\mathbb{R}^{T}$, $(\bm{\lambda}_{n})_{n\in[N]}\in\mathbb{R}_{+}^{\sum_{n}k_{n}}$ are the Lagrangian multipliers associated to (33b) and (33c).

With the polyhedral constraints (32), the graph representation of the disaggregation problem, as illustrated in Figure 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 $\bm{p}$. However, using duality theory, Proposition 4.1 below also gives a characterization of disaggregation:

Proposition 4.1**.**

A profile $\bm{p}\in\mathcal{P}$ is disaggregable iff:

[TABLE]

where

[TABLE]

with $A:=\operatorname{diag}(A_{n})_{n\in[N]}$.

Proof 4.2**.**

From strong duality, we have:

[TABLE]

If the polytope $\mathcal{Y}_{\bm{p}}\cap\mathcal{X}$ given by the constraints of (33) is empty, then there is an infeasibility certificate $\bm{\lambda}^{\top}=(\bm{\lambda}_{0}^{\top}\ \bm{\lambda}_{1}^{\top}\dots\bm{\lambda}_{N}^{\top})\in\mathbb{R}^{T}\times\prod_{n}\mathbb{R}_{+}^{k_{n}}$ such that:

[TABLE]

On the other hand, if $\mathcal{Y}_{\bm{p}}\cap\mathcal{X}$ is nonempty, then a solution to the dual problem (38) is bounded, which implies that \forall\bm{\lambda}:=(\bm{\lambda}_{0},(\bm{\lambda}_{n})_{n})\in\Lambda,\ \ \bm{\lambda}_{0}^{\top}B\bm{p}+\sum_{n}\bm{\lambda}_{n}^{\top}\bm{b}_{n}\geqslant 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 $\Lambda$ can be represented by a finite number of generators (edges), that is, there exists $\Lambda^{*}:=\{\bm{\lambda}^{*(1)},\dots,\bm{\lambda}^{*(d)}\}$ such that:

[TABLE]

Thus, we obtain the following corollary to Proposition 4.1:

Corollary 4.3**.**

There exists a finite set $\Lambda^{*}\subset\Lambda$ such that, for any $\bm{p}\in\mathcal{P}$, $\bm{p}$ is disaggregable iff:

[TABLE]

Remark 4.4**.**

In the transportation case (5), we can write each agent constraints in the form $A_{n}\bm{x}_{n}\leqslant\bm{b}_{n}$ (writing the equality $\sum_{t}x_{n,t}=E_{n}$ is written as two inequalities), and Hoffman conditions (8) can be written in the form (41). Moreover, Theorem 3.3 ensures that one possibility for $\Lambda^{*}$ of Corollary 4.3 is to consider the collection of $2^{T}$ multipliers corresponding to the subsets $\mathcal{T}\subset[T]$ and $\mathcal{N}$ minimizing (9). We skip the details here for brevity.

As in the first part of the paper, we want to use APM to decompose problem (1) and, in the case where disaggregation is not possible, use the result of APM to generate an inequality (34) violated by the current profile $\bm{p}$.

In the case of impossible disaggregation, APM converges to the orbit $(\bm{y}^{\infty},\bm{x}^{\infty})$, and $\bm{\mu}:=\bm{y}^{\infty}-\bm{x}^{\infty}$ defines a separating hyperplane $\bar{\bm{x}}+\bm{\mu}^{\perp}$ , where $\bar{\bm{x}}=\frac{\bm{y}^{\infty}+\bm{x}^{\infty}}{2}$, that satisfies, with $a:=\bar{\bm{x}}.\bm{\mu}$ (note that $\bar{\bm{x}}$ can be replaced by any $\bm{y}\in[\bm{y}^{\infty},\bm{x}^{\infty}]$):

[TABLE]

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):

[TABLE]

Strong duality on these problems implies that there exist $\bm{\lambda}_{0}$ and $\bm{\lambda}$ such that:

[TABLE]

Thus, we obtain $(\bm{\lambda}_{0},\bm{\lambda})\in\Lambda$ satisfying (39), that is, $\bm{\lambda}_{0}^{\top}B\bm{p}+\bm{b}^{\top}\bm{\lambda}<0$, and we can use this to add a new valid additional inequality on $\bm{p}$ of form (34), that will change the current profile $\bm{p}$:

[TABLE]

In Algorithm 5 below, we summarize the proposed decomposition of problem (1). This is a generalization of the decomposition principle used for Algorithm 4.

Remark 4.5**.**

We use the fact that $\bm{\mu}=\bm{y}^{\infty}-\bm{x}^{\infty}$ 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 $\bm{\lambda}_{0}^{\top}B\bm{p}+\bm{b}^{\top}\bm{\lambda}<0$. In practice, one can proceed as in the transportation case and Algorithm 3 use a large $\varepsilon_{\rm{cvg}}$ as stopping criteria in APM, then compute $(\bm{\lambda}_{0},\bm{\lambda})\in\Lambda$ and check if $\bm{\lambda}_{0}^{\top}B\bm{p}+\bm{b}^{\top}\bm{\lambda}<0$. If this is not the case, restart with $\varepsilon_{\rm{cvg}}=\varepsilon_{\rm{cvg}}/2$.

Remark 4.6**.**

When $\mathcal{Y}_{\bm{p}}=\{\bm{x}\in\mathbb{R}^{NT}|A_{0}\bm{x}=B_{0}\bm{p}\}=\{\bm{x}|\sum_{n}\bm{x}_{n}=\bm{p}\}$, we can obtain a non-intrusive version of APM on $\mathcal{Y}_{\bm{p}}$ and $\mathcal{X}$, similar to Algorithm 3. In this case, (52) ensures that we have $\bm{\mu}_{n,t}=-[{\bm{\lambda}_{0}}]_{t}$ for each $n\in[N]$, and $\bm{\lambda}_{0}$ is fixed by $\bm{\mu}$. 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 $\bm{p}$. Thus, Lines 16 to 19 of Algorithm 3 have to be replaced by Algorithm 6.

Remark 4.7** (Link with Benders’ Decomposition).**

In this generalized case, we obtain an algorithm related to Benders’ decomposition [8] (recall that in our specific case (33), the cost function does not involve the variable $\bm{x}$ but only variable $\bm{p}$).

The difference between the proposed Algorithm 5 and Benders’ decomposition lies in the way of generating the new cut. Benders’ decomposition would directly solve the dual problem (38): $\max_{\bm{\lambda}}\{-\bm{\lambda}_{0}^{\top}B_{0}\bm{p}-(\bm{\lambda}_{n})_{n}^{\top}\bm{b}\ |\ \bm{\lambda}_{0}A_{0}+(\bm{\lambda}_{n})_{n}A=0\}$ and obtain a cut if it is unbounded. However, this problem involves the constraints of all users (through $A$ and $\bm{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 example with $T=4$ and $N=3$. Assuming that we have to satisfy the aggregate constraint $\sum_{t}p_{t}=\sum_{n}E_{n}$, we can use the projections on this affine space of solutions of master problems $(\bm{p}^{(s)})_{s}$ to visualize them in dimension 3.

One can wonder if, in the transportation case, applying Algorithm 4 or Algorithm 5 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 4 converges in 3 iterations and Algorithm 5 needs 4 iterations.

We consider the problem (1) with agents constraints (5) with parameters $\bm{\underline{x}}:=0$ and:

[TABLE]

Considering the aggregate equality constraint $\sum_{1\leqslant t\leqslant 4}p_{t}=\sum_{1\leqslant n\leqslant 3}E_{n}=3.3$, we use the canonical projection of $4$ dimensional vectors into the 3 dimensional space $(p_{1},p_{2},p_{3})$ to visualize the cuts and solutions. In this example, there exist $2^{T}-2=14$ nontrivial Hoffman inequalities characterizing disaggregation from Theorem 2.3. The projection of the obtained polytope $\mathcal{P}_{\mathrm{D}}$, as defined in (4), is represented in Figure 4a. One can remark that this polytope has only 6 facets. Our Algorithm 4 applied on this instance with $\varepsilon_{\rm{dis}}=10^{-3}$ and $\varepsilon_{\rm{cvg}}=10^{-5}$ converges in 3 iterations, with successive solutions of the master problem (10) and cuts added:

[TABLE]

Figure 4b 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 5 with the same precision parameters$(\varepsilon_{\rm{dis}},\varepsilon_{\rm{cvg}})$, there are 3 cuts generated and 4 resolutions of master problem needed for convergence, given by (we refer the reader to Remark 4.5 for the numerical precision obtained in the values):

[TABLE]

The 4 successive solutions and the 3 added cuts are represented in the three dimensional space on Figure 4c: we observe that the last cut needed to obtain the convergence of Algorithm 5, corresponds to the first one added with Algorithm 4.

Due to the strict convexity of the cost function $\bm{p}\mapsto f(\bm{p})$, the final solution obtained is the same, unique aggregated optimal solution of (1).

The 4 successive solutions and the 3 added cuts are represented in the three dimensional space on Figure 4c: we observe that the last cut needed to obtain the convergence of Algorithm 5, corresponds to the first one added with Algorithm 4.

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 $\mathcal{T}=\{1,\dots,T\}$, while minimizing the energy cost for the community. We have the following characteristics:

• – the PV plant generates a nondispatchable power profile $(p^{\textsc{pv}}_{t})_{t\in[T]}$ at marginal cost zero;
• – the conventional generator has a starting cost $C^{\textsc{st}}$, minimal and maximal power production $\underline{p}^{g},\overline{p}^{g}$, and piecewise-linear and continuous generation cost function $p^{g}\mapsto f(p^{g})$:

[TABLE]

• – each agent $n\in[N]$ has some flexible appliances which require a global energy demand $E_{n}$ on $[T]$, and has consumption constraints on the total household consumption, on each time period $t\in[T]$, that are formulated with $\bm{\underline{x}}_{n},\bm{\overline{x}}_{n}$. These parameters are confidential because they could for instance contain some information on agent $n$ habits.

The master problem (10) can be written as the following MILP (55):

[TABLE]

In this formulation (55b-55c), where $b_{0,t}:=1$ and $b_{K,t}:=0$, are a mixed integer formulation of the generation cost function $f$. One can show that the Boolean variable $b_{k,t}$ is equal to one iff $p^{g}_{t}\geqslant\theta_{k}$ for each $k\in\{1,\dots,K-1\}$. Note that only $\alpha_{1}$ appears in (55a) because of the continuity of $f$.

Constraints (55d-55e) ensure the on/off and starting constraints of the power plant, (55g) ensures that the power allocated to consumption is not above the total production, and (55h-55i) are the aggregated feasibility conditions already referred to in (6). The nonconvexity of (55) comes from the existence of starting costs and constraints of minimal power, which makes necessary to use Boolean state variables $b^{\textsc{st}},b^{\textsc{on}}$.

We simulate the problem described above for different values of $N\in\{2^{4}$,$2^{5}$, $2^{6}$,$2^{7}$,$2^{8}\}$ and one hundred instances with random parameters for each value of $N$. A scaling factor $\kappa_{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);
• – production costs: $K=3$ , $\theta=[0,70,100,300]\kappa_{N},\bm{c}=[0.2,0.4,0.5]$, $\underline{p}^{g}=50\kappa_{N},\overline{p}^{g}=300\kappa_{N}$, $\alpha_{1}=4$ and $C^{\textsc{st}}=15$;
• – photovoltaic: $p^{\textsc{pv}}_{t}=\left[50(1-\mathrm{cos}(\frac{(t-6)2\pi}{16})+\mathcal{U}([0,10])\right]\kappa_{N}$ for $t\in\{6,\dots,20\}$, $p^{\textsc{pv}}_{t}=0$ otherwise;
• – consumption parameters are drawn randomly with: $\underline{x}_{n,t}\sim\mathcal{U}([0,10])$, $\overline{x}_{n,t}\sim\mathcal{U}([0,5])+\underline{x}_{n,t}$ and $E_{n}\sim\mathcal{U}([\mathds{1}_{T}^{\top}\bm{\underline{x}}_{n},\mathds{1}_{T}^{\top}\bm{\overline{x}}_{n}])$, so that individual feasibility ($\mathcal{X}_{n}\neq\emptyset$) is ensured.

We implement Algorithm 4 using Python 3.5. The MILP (55) 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 12 and 13 of Algorithm 3), we use $\varepsilon_{\rm{dis}}=0.01$ with the operator norm defined by $\left\||\bm{x}|\right\|=\max_{n\in[N]}\sum_{t}|x_{n,t}|$ (to avoid the $\sqrt{N}$ factor in the convergence criteria appearing with $\|.\|_{2}$), and starts with $\varepsilon_{\rm{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 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 (55) solved (number of “cuts” added), remains almost constant when $N$ increases. In all instances, this number is way below the upper bound of $2^{24}>1,6\times 10^{7}$ possible constraints (see Proposition 3.16), 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, Line 2 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.26.

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 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 A.1** (Proof of Item (i)).**

Let us write the stationarity conditions associated to problem (12):

[TABLE]

By summing the latter equalities on $t\in[T]$ and $n\in[N]$, we obtain the equalities:

[TABLE]

Moreover, by summing over $t\in\mathcal{T}^{\circ}_{n}$ the first series of equalities in (56), and by using the complementary slackness conditions, which entail that $\underline{\mu}_{n,t}=\overline{\mu}_{n,t}=0$ for $t\in\mathcal{T}^{\circ}_{n}$, we get

[TABLE]

where we define for each $n\in[N]$:

[TABLE]

From (57) and 2 giving the equality $\sum_{n}E_{n}=\sum_{t}p_{t}$, we obtain $\sum_{t}\nu_{t}=0$ and:

[TABLE]

Let us show that $\forall t\in\mathcal{T},\nu_{t}>0$. For $t\in\mathcal{T}$, there exists $n$ s.t. $y_{n,t}>\overline{x}_{n,t}$. From (56) we get:

[TABLE]

Suppose by contradiction that there exists $n\notin\mathcal{N}$ and $\hat{t}\in\mathcal{T}$ such that $x_{n,\hat{t}}<\overline{x}_{n,\hat{t}}$. We obtain from complementary slackness conditions and the first equality in (56):

[TABLE]

From this, we obtain $\overline{\mathcal{T}}_{n}\subset\mathcal{T}$: indeed, for $t\in\overline{\mathcal{T}}_{n}$, we have $y_{n,t}+\lambda_{n}\geqslant\overline{x}_{n,t}$ which gives:

[TABLE]

From the condition (59) and as $\nu_{t}>0$ for each $t\in\overline{\mathcal{T}}_{n}$ because $\overline{\mathcal{T}}_{n}\subset\mathcal{T}$ and (60), we get:

[TABLE]

which is strictly negative: this implies that there exists $t^{\prime}\in\underline{\mathcal{T}}_{n}$ such that $y_{n,t^{\prime}}<\underline{x}_{n,t^{\prime}}$. Necessarily, $t^{\prime}\notin\mathcal{T}$ because $\nu_{t^{\prime}}=y_{n,t^{\prime}}-x_{n,t^{\prime}}<\underline{x}_{n,t^{\prime}}-\underline{x}_{n,t^{\prime}}=0$. Then, as we have $\sum_{m\in[N]}y_{m,t^{\prime}}=p_{t^{\prime}}\geqslant\sum_{m}\underline{x}_{m,t^{\prime}}$, there exists $m\in[N]$ such that $y_{m,t^{\prime}}>\underline{x}_{m,t^{\prime}}$. If $\lambda_{m}\leqslant 0$, and as $x_{m,t^{\prime}}=y_{m,t^{\prime}}-\nu_{t^{\prime}}>\underline{x}_{m,t^{\prime}}$, we get:

[TABLE]

which is impossible, thus $\lambda_{m}>0$. Now, we observe that $\mathcal{T}^{\circ}_{n}\subset\mathcal{T}$. Indeed, otherwise, if $t"\in\mathcal{T}^{\circ}_{n}\cap\mathcal{T}^{c}$, we have $\nu_{t"}=-\lambda_{n}>0$ and $x_{m,t"}=y_{m,t"}-\nu_{t"}<y_{m,t"}<\overline{x}_{m,t"}$, thus we get:

[TABLE]

which is impossible, thus $\mathcal{T}^{\circ}_{n}\subset\mathcal{T}$.

Finally, since $\overline{\mathcal{T}}_{n}^{c}\neq\emptyset$, consider $t_{0}\in\arg\min_{t\notin\overline{\mathcal{T}}_{n}}\{\overline{x}_{n,t}-y_{n,t}\}$. By (58), we obtain:

[TABLE]

and thus:

[TABLE]

which contradicts $n\notin\mathcal{N}$. Thus we have shown that $\forall n\notin\mathcal{N},\forall t\in\mathcal{T},\ x_{n,t}=\overline{x}_{n,t}$. From (56) and as $\forall t\in\mathcal{T},\nu_{t}>0$ from (60), we obtain $\forall n\notin\mathcal{N},\forall t\in\mathcal{T},\ y_{n,t}=x_{n,t}+\nu_{t}>\overline{x}_{n,t}$. This terminates the proof for Item (i).

Proof A.2** (Proof of Item (ii)).**

Let us consider $\mathcal{T}^{\prime}:=\{t|\nu_{t}>0\}$. We already showed that $\mathcal{T}\subset\mathcal{T}^{\prime}$ in (60). To prove the other inclusion and obtain Item (ii), we observe that, considering $n\notin\mathcal{N}$ (nonempty as shown independently in Item (v)) and considering $\mathcal{T}^{\prime}$ instead of $\mathcal{T}$, all the facts established in the proof of Item (i) above also hold: by contradiction of $\forall t\in\mathcal{T}^{\prime}$, $x_{n,t}=\overline{x}_{n,t}$, we first obtain $\lambda_{n}<0$ as in (61). Then, as for $t"\in\mathcal{T}^{\circ}_{n}$, we have $\nu_{t"}=-\lambda_{n}>0$, we get $\mathcal{T}^{\circ}_{n}\cap\mathcal{T}^{\prime c}=\emptyset$. Finally the same sequence of inequalities as (63-64) show a contradiction. Consequently, for each $t\in\mathcal{T}^{\prime}$, $x_{n,t}=\overline{x}_{n,t}$ and $y_{n,t}=x_{n,t}+\nu_{t}>\overline{x}_{n,t}$, thus $t\in\mathcal{T}$ and $\mathcal{T}^{\prime}\subset\mathcal{T}$.

Proof A.3** (Proof of Item (iii)).**

Suppose on the contrary that there exists $n\in\mathcal{N}$ such that $\lambda_{n}\geqslant 0$. For $t\in\mathcal{T}^{\circ}_{n}$, we have $\nu_{t}=-\lambda_{n}\leqslant 0$, thus, $\mathcal{T}^{\circ}_{n}\subset\mathcal{T}^{c}$. Then, if $t\in\mathcal{T}$ and if $x_{n,t}<\overline{x}_{n,t}$, we would have:

[TABLE]

which is impossible, thus $x_{n,t}=\overline{x}_{n,t}$, and $\mathcal{T}\subset\overline{\mathcal{T}}_{n}$. As we show independently in Item (v) that $\mathcal{T}\neq\emptyset$, we know $\overline{\mathcal{T}}_{n}\neq\emptyset$. Let us consider $t_{0}\in\arg\min_{t\notin\underline{\mathcal{T}}_{n}}\{y_{n,t}-\underline{x}_{n,t}\}$. By (58), we obtain:

[TABLE]

and thus:

[TABLE]

which contradicts $n\in\mathcal{N}$ and terminates the proof of Item (iii).

Proof A.4** (Proof of Item (iv)).**

From (ii), we know that $\mathcal{T}^{c}=\{t|\nu_{t}\leqslant 0\}$, thus, if $t\notin\mathcal{T}$ and $n\in\mathcal{N}$, if $x_{n,t}>\underline{x}_{n,t}$ then we would have $x_{n,t}\leqslant y_{n,t}+\lambda_{n}=x_{n,t}+\nu_{t}+\lambda_{n}<x_{n,t}$, which is a contradiction.

Proof A.5** (Proof of Item (v)).**

From $\sum_{t}\nu_{t}=0$, we see that if $\mathcal{T}=\emptyset$, then this means that $\nu_{t}=0$ for all $t\in[T]$, and thus ${\bm{y}}={\bm{x}}$ which is a contradiction. Thus there exists $t_{0}$ such that $\nu_{t_{0}}>0$ and because $\sum_{t\in\mathcal{T}}\nu_{t}=0$, there exists $t^{\prime}_{0}$ such that $\nu_{t^{\prime}_{0}}<0$.

If $\mathcal{N}=\emptyset$, then using (i), we would have for all $n$, $y_{n,t_{0}}>\overline{x}_{n,t_{0}}$ and thus $p_{t_{0}}>\sum_{n\in[N]}\overline{x}_{n,t_{0}}$, which contradicts the aggregate upper bound constraint $\forall t,\ p_{t}\leqslant\sum_{n\in[N]}\overline{x}_{n,t}$.

If $\mathcal{N}^{c}=\emptyset$, then using (iv), we would have for all $n$, $y_{n,t_{0}^{\prime}}<\underline{x}_{n,t_{0}^{\prime}}$ and thus $p_{t_{0}^{\prime}}<\sum_{n\in[N]}\underline{x}_{n,t_{0}^{\prime}}$, which contradicts the aggregate lower bound constraint $\forall t,\ p_{t}\geqslant\sum_{n\in[N]}\underline{x}_{n,t}$.

Appendix B Proof of Lemma 3.14

Proof B.1** (Proof of Item (i)).**

From $\bm{x}^{(K)}=P_{\mathcal{X}}(\bm{y}^{(K-1)})$ and $\bm{y}^{(K)}=P_{\mathcal{Y}}(\bm{x}^{(K)})$, we obtain, similarly to (56):

[TABLE]

where the Lagrangian multipliers $\lambda_{n}^{(K)},\underline{\mu}_{n,t}^{(K)},\overline{\mu}_{n,t}^{(K)}$ (resp. $\nu_{t}^{(K)}$) are associated to the quadratic problem characterizing the projections $P_{\mathcal{X}}(\bm{y}^{(K-1)})$ (resp. $\bm{y}^{(K)}=P_{\mathcal{Y}}(\bm{x}^{(K)})$). We obtain equalities similar to (57, 58). We proceed as for Proposition 3.4(i) and suppose that there exists $n\notin\mathcal{N}$ and $\hat{t}\in\mathcal{T}$ such that $x_{n,\hat{t}}<\overline{x}_{n,\hat{t}}$. Then, as $\|\bm{y}^{(K)}-\bm{y}^{\infty}\|\leqslant\tfrac{\varepsilon_{\rm{cvg}}}{1-\rho}$ and $\sum_{t\in[T]}\bm{y}_{t}^{(K)}=\sum_{t\in[T]}\bm{y}^{\infty}_{t}$, we have for each $n\in[N],t\in[T]$, $|y_{n,t}^{(K)}-y_{n,t}^{\infty}|\leqslant\tfrac{\varepsilon_{\rm{cvg}}}{2(1-\rho)}$, and thus we get:

[TABLE]

as $\nu_{\hat{t}}^{\infty}\geqslant\underline{\nu}>2B\varepsilon_{\rm{cvg}}$. Let us now consider $t^{\prime}\in{\mathcal{T}^{\circ}_{n}}^{(K)}\cup\overline{\mathcal{T}}_{n}^{(K)}$, then:

[TABLE]

which shows that $t^{\prime}\in\mathcal{T}^{(K)}=\mathcal{T}^{\infty}$ and thus ${\mathcal{T}^{\circ}_{n}}^{(K)}\cup\overline{\mathcal{T}}_{n}^{(K)}\subset\mathcal{T}$. Then, the same sequence of inequalities as (62, 63, 64) applied to $\bm{y}^{(K-1)}$ gives a contradiction to $n\notin\mathcal{N}$.

Proof B.2** (Proof of Item (ii)).**

The proof of Item (ii) is symmetric to the one of Item (i): if we suppose that there exists $n\in\mathcal{N}$ and $\hat{t}\notin\mathcal{T}$ such that $x_{n,\hat{t}}^{(K)}>\underline{x}_{n,\hat{t}}$, we obtain, symmetrically to (69), that $\lambda_{n}^{(K)}\geqslant-\tfrac{\varepsilon_{\rm{cvg}}}{2(1-\rho)}$. Then, considering $t^{\prime}\in\underline{\mathcal{T}}_{n}^{(K)}\cup{\mathcal{T}^{\circ}_{n}}^{(K)}$, we show, symmetrically to (70), that $\nu_{t^{\prime}}^{(K)}<B\varepsilon_{\rm{cvg}}$ i.e.  $t^{\prime}\notin\mathcal{T}$ and thus $\underline{\mathcal{T}}_{n}^{(K)}\cup{\mathcal{T}^{\circ}_{n}}^{(K)}\subset\mathcal{T}^{c}$. We conclude by obtaining a contradiction to $n\in\mathcal{N}$ by the same sequence of inequalities as (65, 66, 67).