Assessing the Scalability and Privacy of Energy Communities by using a Large-Scale Distributed and Parallel Real-Time Optimization

In the context of the energy transition, energy communities are gaining increasing attention all over the world, in recent years. By participating in an energy community, prosumers may take a leading role in the energy transition and improve the self-consumption of renewable energy produced inside the community. Prosumers can carry out energy exchanges inside the energy community and provide ancillary services to the system operators, thus contributing to improve the efficiency and stability of the grid. A novel scalable, privacy-preserving, and real-time distributed parallel optimization is proposed to manage a large-scale energy community, considering energy exchanges inside the community according to the model of virtual self-consumption and the provision of ancillary services. The proposed method preserves the privacy of prosumers and allows the assessment of the impact of energy exchanges on the ancillary services provided by an energy community. Simulation results confirmed that the proposed method is superior in terms of privacy if compared with the equivalent centralized optimization and that it has a convergence rate higher than that of the splitting conic solver (SCS).


A. Background
The recent evolution of the European regulatory framework promotes the aggregation of prosumers in Energy Communities (ECs) to maximize the self-consumption of renewable energy produced inside the community and to provide services to system operators [1], [2]. By participating in an EC, consumers with the production capability (prosumers) can improve the exploitation of the renewable energy produced inside the community and carry out energy exchanges in the community also according to the model of virtual self-consumption [3]. Prosumers inside an EC may also increase their economic benefits by providing ancillary services to the Transmission System Operator (TSO) or Distribution System Operator (DSO).

B. Solutions for the management of an energy community
Centralized optimization methods for the management of a community of prosumers and the provision of energy flexibility in the ancillary services market are proposed in [4][5][6][7][8][9][10][11]. A large number of residential prosumers are optimized using a Mixed-Integer Linear Programming (MILP) optimization model by the authors of [4]. The optimization aims at providing up and down-regulation without considering energy exchanges inside the community. A method for an energy community, consisting of 48 households and 42 commercial buildings, is proposed in [5] for the delivery of the frequency restoration service in the Ancillary Service Market (ASM). The authors of [6] propose a two-stage stochastic method based on a Mixed-Integer Linear Programming (MILP) optimization model for the optimal management of an energy community consisting of 50 households. The optimal operation of the aggregators is achieved by the authors in [7] and [8]. A new framework to assist the interactions between the system operator and aggregators in a balancing market is designed considering the transactive energy concept [7]. The framework prevents voltage level violations and congestions in the distribution network. An optimization model based on convexification of the mixed-integer formulation is proposed by the authors of [8] to manage the operation of aggregators with the grid operators, considering the uncertainties due to prediction errors. The authors of [9] propose a robust virtual battery model to manage the flexibility of aggregated prosumers.
Distribution locational marginal prices are used to optimally schedule the resources of prosumers and avoid network congestions. A three-level optimization problem is proposed by the authors of [10] to coordinate the operation of prosumers considering Volt/VAR optimization in a mediumvoltage distribution network. A coordinated energy scheduling for a microgrid with prosumers is proposed by the authors in [11] and a genetic algorithm for the scheduling of the batteries of prosumers to reduce energy losses and the energy exchange with the main grid.
Centralized methods, however, involve a central coordinator and may suffer in some cases as in [4], from the issue of privacy. Considering the benefits, objectives, and challenging issues related to the management of the energy community and the scalability and privacy requirements, distributed optimization methods for energy management have recently gained a growing interest [12][13]. Column generation and Dantzig-Wolfe decomposition is considered in [12] where a Mixed-Integer Linear Programming (MILP) optimization model is proposed to optimally coordinate Distributed Energy Resources (DERs). Distributed optimization methods exhibit, indeed, some improvements if compared to centralized approaches since they allow solving large-scale problems by distributing them into simpler subproblems to be solved by many processors operating in parallel, thus reducing the computational efforts [13]. The prerequisite to keeping the private information of prosumers preserved is also considered in [13] where a scalable and privacy-preserving distributed parallel optimization is proposed to manage prosumers with residential Photovoltaic (PV) and battery systems and allow them to participate in the ASM. Although the proposed method, based on Linear Programming (LP), exhibits better performances in terms of scalability if compared to previous methods based on integer and Quadratic programming (QP), electrical energy exchanges in an energy community are not considered. A distributed online control algorithm for the coordination and decentralized optimization of flexible energy resources is proposed in [14]. The method is applied to optimize the energy management for eight households in a neighborhood, however, the issue of scalability is not assessed. The authors of [15] consider the issue of scalability of a large number of storage devices and propose the alternating direction method of multipliers (ADMM) to solve a MILP problem to provide flexibility to the Distribution System Operator (DSO) considering only 100 users. Unlike the basic version of ADMM, described in [16], which is designed to be solved serially, parallelizable versions need carefully chosen hyperparameters such as step sizes which may be hard to tune for large-scale instances.
To be able to scale up to a very large number of prosumers, LP problems should be solved extremely fast. To this end, an ADMM-based interior method [17], which can be parallelized, is proposed in [18] and [19]. Another parallelizable method, which is even more scalable than the one described in [18] and [19], is presented in [20] and tested for problems with trillions of variables. The brilliant idea of the authors of [20] is the use of a quadratic regularization term which allows applying an accelerated gradient. As a result, O(1/k^2) convergence rate is reached while as far as we know parallel ADMM has a convergence rate equal to O(1/k) [21,22,23].
It is worth noting that, especially in the data science community, parallel computing, using increasingly ubiquitous multi-threaded CPUs and Graphics Processing Units (GPUs), is used [23]. This approach is valid when the privacy of prosumers is not the main concern while, in the real-time optimization considered in this paper, two important goals should be achieved. The first one is related to the scalability of the approach, requiring that the real-time optimization is carried out efficiently, with an acceptable computational time regardless of the number of prosumers. In this regard, in the proposed method, similar to other firstorder methods (FOMs) [20], the bottleneck is matrix-vector multiplications, which can be easily distributed between clusters [24], involving an acceptable number of communications between prosumers and the virtual layer as it will be described in Section III. The second goal is related to the privacy of the approach. Differently from what has been considered in [20], calculations should be carried out privately since the spatial-temporal context information of prosumers is vulnerable to third-party exploitations [25].
Recent research on privacy issues in smart grids mainly concentrates on the load profiles obtained by smart meters and indicates some technological methods to impose privacy [26]. Detailed surveys on privacy-preserving methods can be found in [27,28,29]. Common approaches used to preserve privacy include anonymizing or pseudonymization [30] and aggregating smart metering data [31]. Indeed, data aggregation for securing the consumers' data combines data from the meter reading at the gateway so that the adversarial agent or cyber attacker cannot identify an individual user's information. Also, this approach can be combined with the notion of differential privacy [26] as it is suggested at the end of Section III-B. Further methods for preserving privacy consider the relation between privacy and load data resolution and propose to split the load data into groups with different resolutions and distinct authorization levels [26]. Restricted disclosure of data in a decomposition framework is proposed by the authors in [32], while in [33] the privacy is protected by limiting the revelation of the individual consumption using encryption of electricity usage. The addition of numerical artifacts, including perturbing noises or random sequences to the original load data, is proposed in by the authors in [34] to obscure different contributions. Secure signal processing in the encrypted domain can be also used for privacy preservation as discussed by the authors in [35].

C. The novel contribution of the proposed method
Considering that in the previous research in the field of large-scale energy communities the issues of scalability and privacy of prosumers have been only partially explored, with methods and simulation studies considering only a limited This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ number of prosumers, in this paper a novel extreme-scale distributed parallel optimization is proposed.
The method evaluates the flexibility contribution of an energy community in the ASM. Differently from [13], in the proposed method the influence of energy exchanges inside the community on the provision of ancillary services is assessed considering the model of virtual self-consumption in the EC. Due to the inclusion of energy exchanges in the model, and the requirement of privacy preservation, a new and more challenging optimization problem is considered and solved with a new optimization approach that is distributed but not decentralized since it involves a virtual layer as a means of communication between all prosumers as it will be explained in the following sections.
To make the proposed real-time optimization flexible enough to be applicable even in the case of fast frequency reserve service and considering that higher communication delays can happen during the optimization in practice, we believe that a time window of a few seconds is desirable for the real-time optimization [36].
The obtained numerical results show that the real-time optimization can be carried out with an acceptable amount of communication time which makes it feasible for real-world applications. The designed warm-start strategy, in which the optimal solution of the previous optimization is taken as the initial starting point for the current optimization, allows reducing the number of communications. The warm-start strategy is a common theme in many real-time optimization scenarios such as Model Predictive Control (MPC) and Reinforcement Learning (RL) [37,38].
The main innovative contributions of this paper are as follows: -a novel scalable, privacy-preserving, and real-time distributed parallel optimization is proposed to manage a large scale energy community considering energy exchanges inside the community and the provision of ancillary services; -a new method is proposed to preserve the privacy of prosumers of the energy community by carrying out parallel matrix-vector multiplications in a privacypreserving mode, prosumers can execute some local calculations coded up and executed locally on their microcontrollers; -as a result of the warm-start strategy in which we take advantage of the optimal solution of the previous optimization, the proposed distributed-parallel optimization can be executed in real-time with a reduced number of communications between prosumers and the virtual layer; -the proposed method allows the assessment of the impact of energy exchanges on the ancillary services provided by an energy community. The paper is structured as follows: Section II described the mathematical model of the optimization. The novel scalable and privacy-preserving extreme-scale and real-time distributed parallel optimization is described in Section III. Some numerical results are presented in Section IV and discussed in Section V. Conclusions are given in Section VI. The Appendix describes the local constraints of the proposed method and presents an approach to evaluate its robustness with some results.

A. Considered energy community and the role of prosumers
An energy community of residential prosumers and consumers has been considered in the proposed approach as presented in [6,13,40]. The term prosumer is a combination of the word "producer" and "consumer" and represents consumers, equipped with PV generation and Energy Storage Systems (ESSs), that can also generate energy. Their role in the energy community is fundamental since they can exchange some excess energy produced locally by their PV systems with other prosumers at a price lower than that of the energy absorbed from the grid.
The proposed approach is based on the virtual model of self-consumption, where the self-consumption is virtual so that the energy balance inside the community is achieved by taking into account the energy exchanges evaluated at each Point Of Delivery (POD) [40]. The number of prosumers in the energy community is limited by the market zone as described in [40]. Even if according to this model energy exchanges inside the community can be extended towards a market zone, with prosumers having PODs in different feeders, in this paper the virtual self-consumption is limited within the same MV-LV transformer substation. This is achieved in the mathematical formulation by considering the constraint (4) valid for each substation. The distribution network hosting capacity for connecting the prosumers is considered. The distribution network operator only allows the connection to the network of new prosumers that do not violate the security constraints of the network. Each prosumer can exchange a maximum power with the grid limited by its contractual agreement. The maximum powers that can be imported/exported from/to the grid are denoted by , or , as detailed in the Appendix (in the case studies it has been assumed that they coincide and can assume the values of 3 kW, 4.5 kW, or 6 kW). For each prosumer, the sum of the powers exchanged with the grid and with other prosumers is limited by as , , , as shown by constraints (27) and (28) described in the Appendix. The interactions of each prosumer with other ones are thus limited and the grid hosting capacity constraints are considered.
An energy community manager (ECM), by scheduling and operating the batteries of prosumers, tries to maximize the profit of the energy community by minimizing the energy costs and maximizing the revenue deriving from the provision of flexibility to the TSO. Depending on their need and on the choices of the ECM, prosumers can absorb the energy from the local network and use the electrical energy produced by their PV plant by self-consumption or sharing with other prosumers or storing it in their local batteries. When the excess of the produced electrical energy is shared with other prosumers, the local public distribution network is used without necessitating any physical modification in the structure of the distribution network since the configuration has been also defined virtually.
By carrying out energy exchanging inside the community, prosumers not only aim at reducing their energy costs by improving the self-consumption of the energy produced by the local energy sources or stored in their batteries but also to obtain revenue by providing ancillary services to the TSO by trading electrical energy with the grid in the ancillary services market (ACM) [41]. Prosumers, by exchanging some excess energy produced locally by their PV systems with other prosumers can contribute to reducing the dependency on the grid by efficiently allocating available energy within the community, thus improving the grid efficiency.

B. Mathematical models of the day-ahead and realtime optimizations
It is assumed that the ECM, based on the day-ahead forecasts of the photovoltaic (PV) production and electrical load demand, can carry out a day ahead optimization to decide if it is convenient to offer services of flexibility to the TSO assuming its baseline as a reference [4,13]. The TSO, if the flexibility offer is accepted, can ask the ECM to provide the services during the real-time phase. An imbalance cost can be applied if the agreed flexibility is not respected, therefore a real-time optimization based on shortterm forecasts is carried out by the energy community.
The energy community objective function shown in (1) aims at maximizing the revenue coming from the flexibility provision and minimizing the costs of the community. A day is divided into 96-time intervals of fifteen minutes (∆t), accordingly, the continuous variables are the mean powers in the quarter of an hour t in the range [1,96]. For the day ahead optimization, to take into account the uncertain behaviors of PV facilities and electrical loads, scenarios determined by day-ahead forecasts are generated as in [4] and [42]. In the following equations, a single scenario of electrical load and photovoltaic power is assumed only to simplify the notation. More details related to the way the day-ahead stochastic optimization is carried out are given in [42]. The objective function to minimize is defined in Eq. 1.
where is the revenue coming from the flexibility provision, is the cost/revenue related to the energy exchange with the grid.
is the power injected/absorbed into/from the grid by the prosumer . , ( )/ , ( ) is the cost/revenue of the energy absorbed/injected into the network.
is the offer decided by the ECM for down and upregulation. For each prosumer , , ( ) is the baseline, , ℎ ( )/ , ℎ ( ) is the charging/discharging power of the battery [13]. and ∀ ∈ . , , and are fixed parameters for the real-time optimization since they are determined during the day-ahead optimization. They can vary in the range [1,96] to define the time ranges [ , ] for up-regulation and [ , ] for down-regulation during which the provision of energy flexibility in the ASM should be provided. In realtime optimization, there are some fixed parameters determined by using the variable neighborhood search method as described in [13]. , 2 ( )/ , 2 ( ) is the power imported/exported for each prosumer. The global equality constraint is related to the energy exchanges among all prosumers as shown in Eq. (4), and the global inequality constraints as shown in Eq. (5) and (6), respectively in the case of up and down-regulation are introduced to respect the minimum power flexibility / that the community should provide. ∑ , The real-time optimization minimizes the objective function presented in (7) for each of the intervals in the time ranges [ , ] and/or [ , ] where , , and are determined during the day-ahead phase. It is worth noting that, since optimization is based on real-time forecasts of electrical loads and PV production, also considering measurements in the time-step, can be assumed to be quite accurate. It is then assumed that the forecast errors in each quarter are normally distributed variables with a variance such that we have a maximum percentage error of 30% for the electrical load forecast and 20% for PV production forecast. Even if in the real-time optimization considered here the computational time required for the optimization should comply with the requirement to achieve an optimal solution every quarter, the proposed approach is also applicable in the case of fast frequency reserve service where a computational time of at most few seconds is required.

∀ ∈ [ , ] and ∀ ∈ [ , ]
The problem is subject to previously described local constraints and the following global equality constraint for energy exchange: (8) where , 2 ( )/ , 2 is the mean power imported/exported by each prosumer from/to the EC. In the real-time optimization, instead of (5) and (6), the global equality constraints (9) and (10) should be respected according to the power flexibilities, ( ) and/or ( ), accepted for every quarter by the TSO during the day-ahead optimization.
Both the day-ahead and real-time problems have local and global equality and inequality constraints. The local constraints for each prosumer include the equality constraints related to power balance considering energy exchanges inside the EC, the inequality constraints related to the state of charge (SC) of the battery, and the equation modeling the SC as described in the Appendix. Local inequality constraints are also considered to limit the power from/to the grid and the discharging/charging power of the battery as described in the Appendix.
The proposed method has been implemented by using the layered architecture shown in Fig. 1 and presented in [13] [43]. In the architecture, a physical layer is essential to provide each prosumer with the infrastructure and devices necessary to connect to the grid and to communicate with other prosumers, while the virtual layer allows all prosumers to have a privacy-preserving exchange of only some information to carry out parallel computing which involves the local computational resources of all users during the realtime operation.

A. Description of the proposed method and comparison with ADMM
In this section, a more abstract version of the optimization problem introduced in the previous section is developed to take advantage of recent theoretical progress which makes the solution of large-scale linear programming problems feasible in a distributed private way.
This section has two parts. At first, the basic version of a new parallelizable LP solver is described which has been successful to tackle LPs with trillions of variables. It is worth noting that, in the original context for which the algorithm has been proposed, the issue of privacy has not been considered. The novel contribution of the proposed method compared to ECLIPSE [20] is related to a new mathematical formulation that allows not only to maintain very high scalability of the original method but also to keep the privacy of prosumers preserved.
To better work with the resulting LP problem presented in the previous section, it is conventional among the optimization community to write every equality constraint = as two inequality ones ≤ and − ≤ − . So it is possible to consider (8), (9), and (10) as inequality constraints. Therefore, the LP is written as follows: . . , are vectors of unknown variables. All of the global constraints are re-formulated as 1 1 +. . . + ≤ , the local constraints presented in [4] are re-formulated as ̂≤̂, = 1, . . . , and the bound constraints are represented by ≤ ≤ . Now solving the resulting LP is challenging in very largescale problems, therefore the idea is to approximately solve the problem in parallel. The splitting methods, such as ADMM, have been widely used by researchers in so many areas, including smart girds [15] in which quadratic programming (QP) needs to be solved in each iteration by each processor (user). Solving each QP boils down to solving a linear system which reduces the scalability of the method. An example of this line of approaches is the one presented in [17] where an ADMM-based interior method is presented which can be parallelized, but the issue is that each step involves solving a linear system that is computationally demanding. Also, the numerical test cases are far from what these days is considered a "large-scale LP". To address the issue of scalability we focused only on those methods that require at most matrix-vector multiplication and solving linear systems directly, as intermediate steps are not allowed. The reason for this restriction is that matrix-vector multiplication is well studied in the literature and can scale very well not only on GPUs but also in distributed settings.
A well-known parallelizable version of ADMM is the Splitting conic solver (SCS) [18] and [19] which is a "matrixfree" version of ADMM which means that it involves matrix-This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3187204 vector multiplications as the most computationally demanding part, as opposed to the basic version of ADMM which is not "matrix-free". This makes SCS more scalable and better to be used for parallel and distributed computation.
Another highly scalable method is ELIPSE, proposed in [20], where the authors use a quadratic regularization term that allows applying accelerated gradient descent (AGD) and as a result, reaching O(1/k^2) convergence rate while as far as we know parallel ADMM has a convergence rate of O(1/k) [21,22]. This speed-up in ECLIPSE in turn has made solving problems up to trillions of variables a feasible task. In smaller problems, SCS outperforms ECLIPSE but as the size of the problem increases, ECLIPSE converges much faster than SCS because fewer and structured matrix-vector multiplications are used. From the privacy point of view, both SCS and ECLIPSE are similar since in both the results of multiplication of private blocks times the decision variables are to be sent to the virtual layer.
Also, from a communication overhead point of view ECLIPSE is better since fewer matrix-vector multiplications translate to a smaller number of communications. This makes ECLIPSE more scalable and more communication efficient than SCS.
From a privacy preservation point of view, although both methods are the same in that they both are in essence based on matrix-vector multiplication, having a fewer number of communications makes the proposed method less appealing to a third-party agent to try to infer about the private information of prosumers.
Another technical advantage of ECLIPSE over SCS is the ease of implementation and it can be easily programmed with low-level programming languages (like C++ or java) for Internet of things (IoT) devices.

B. Mathematical formalization of the proposed method
In this section after reformulating the problem as a smooth optimization (Eq. 14), to make the paper as self-contained as possible, a short description of the first-order method is given and the way it is implemented in the general case is presented in Algorithms 1 and 2. Then we describe in Algorithm 3 (which is the basic version of gradient descent with a fixed step-size and has (1/ ) convergence rate) how to compute the gradient for solving the regularized Lagrangian dual of the LP problem (Eq. 14). Then we describe in Algorithm 4 the accelerated gradient descent (with (1/ 2 ) convergence rate) applied to LP which is basically what has already been presented in [20]. Also, we give the full details of computing spectral bound in Algorithm 5 which is a pre-requisite to Algorithm 4. Then the main contribution of our proposed method is presented in Algorithms 6 and 7 which can be seen as the evolution of Algorithms 5 and 4 respectively, implemented in a privacy-preserving and parallel way.
To recast the problem to make it amenable to matrix-free settings like SCS or ECLIPSE, let us denote: Then (11) can be re-written as follows: which in turn can be resolved by solving the following penalized Lagrangian dual problem: where which > 0 is a parameter. Now since regularization term introduced in the Lagrangian dual makes ( ) differentiable with respect to (Lemma 3 in [20]) then the problem is amenable to firstorder methods for smooth optimization [45]. We start from the most straightforward way of smooth optimization i.e. the projected gradient descent. Basic projected gradient descent is described in Algorithm 1.
Algorithm 1: general gradient descent Step 1 Start with an initial guess 2 While ‖ ( )‖ is bigger than a threshold do: Indeed line 3 is the gradient descent and line 4 is the projection part. We can also compress by merging lines 3 and 4 as shown in Algorithm 2. Algorithm 2: general gradient descent (compressed version) Step 1 Start with an initial guess 2 While ‖ ( )‖ is bigger than a threshold do: 4 end Now we need to compute . First, note that the minimum of + 2 + ( − ) is reached when which means that the function satisfies the Lipschitz condition. Now we can use step-size 0 < ≤ 1 ‖ ‖ 2 but as in [20] we simply choose a fixed step-size = = ‖ ‖ 2 in which ‖ ‖ 2 is the maximum eigenvalue of [44]. This convexity in addition to Lipschitz continuity gives a convergence rate equal to (1/ ). In other words, if we denote the optimal value of ( ) = ≤ ≤ { + 2 + ( − )} by * then * − ( ) ≤ (1/ ).
Also, by using the two-step approach for updating (if function is convex and has Lipschitz continuity) ( 1 2 ) is guaranteed [45]. The two-step approach is called Accelerated Gradient Descent (AGD) [20,45], in which a more sophisticated way of updating +1 is used as follows: with initial parameters 1 = 1 and 1 = 0 .
Since the original ECLIPSE is proposed for applications in which the issue of privacy is not addressed, we need to figure out how to implement the different steps of ECLIPSE in a way that the privacy of prosumers is preserved. To compute ‖ ‖ 2 (the maximum eigenvalue of ) the power method described in [44] can be used as described in Algorithm 5. Since matrix A is fixed, the spectral norm ‖ ‖ 2 is calculated only once during the day-ahead optimization. Since adding an error to the forecast does not change the matrix A, but only the right-hand side b, we do not need to compute the spectral norm of A again and again during the real-time optimization.
To carry out steps 3) and 5), in a parallel and private we consider that the matrix-vector multiplication can be written as follows: On the other hand, in (11), for ≠ we have ̂= 0 hence Once is calculated we can simply do the update = (step3) as follows: By using the formulas (18)- (19) we come up with the parallelized private version for computing spectral norm as shown in Algorithm 6.

Algorithm 6: Computing spectral norm for matrix A in a parallel private fashion
Step 1 All the users independently from others start with initial configurations 1 , . . . , 2 All the users start from = −∞ 3 Each prosumer i compute and sends the resulting "vector" to the virtual layer (VL) 4 The VL computes the aggregate vector 0 = 1 1 +. . . + and broadcasts it back to all prosumers 5 Each prosumer i updates its new configuration = 0 +̂̂ (based on formula (19)) 6 Each prosumer i computes the scalar value and sends it to the VL 7 The VL computes the aggregate value ‖ ‖ = √ 1 1 +. . . + and broadcasts the resulting scalar value back to all prosumers 8 Each prosumer i updates its new configuration = ‖ ‖ 9 Each prosumer i compute =̂ and sends the scalar value to the VL 10 The VL computes = = 0 0 + 1 1 +. . . + (based on formula (18)) and then broadcasts the scalar value back to all prosumers 11 All users check whether | − | > , if yes go to step 3 otherwise the consense is reached and all of them stop. 12 All users come up with the result ‖ ‖ 2 = as the spectral norm of The details of how Algorithm 6 evolves from Algorithm 5 are given below: Steps 1) and 2) of Algorithm 6 are the same as Algorithm 5.
Steps 3), 4) and 5) of Algorithm 6 use formula (19) to privately parallelize step 3) of Algorithm 5. Steps 6), 7) and 8) of algorithm 6 are used to privately parallelize step 4) of Algorithm 5. Steps 9) and 10) of algorithm 6 use formula (18) to privately parallelize steps 5) and 6) of Algorithm 5. Finally Steps 11) and 12) of Algorithm 6 are the same as steps 7) and steps 8) of Algorithm 5. Indeed, as shown in Algorithm 6, only vectors Since the number of global constraints is three, then 1 1 , . . . , are vectors of small length so they can efficiently be communicated to VL. Also one can not infer This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3187204 from the product its components and , and the same argument applies to 0 = 1 1 +. . . + .

2)
Since (where =̂) and are scalar values, the information hidden in different components of the product is hidden.
The information flow between prosumers and the virtual layer for the calculation of the spectral norm (in a parallelized private way) is presented in Fig. 2.

FIGURE. 2. Information flow between prosumers and the virtual layer for the calculation of the spectral norm.
Now, the way we compute and in a parallel and private way (formulas (19) and (18) Therefore the step 9) in the algorithm 4 is parallelized by the following formula: The above formula implies that 'th user can update its configuration ( ) without revealing neither the information of its objective function ( ) nor the dual variables which would otherwise reveal how far she/he is from having her/his local constraints satisfied. Also, step 12 of Algorithm 4 which is the computation of the gradient ( ) = ( ) − is implemented by using formula (18)  and not the full vectors 1 , . . . , . So, the 'th prosumer only communicates the scalar value . Besides communication efficiency, this approach preserves privacy since prosumers only allow VL to know the squared norm of the dual vectors. The parallelized private AGD for LP is presented in Algorithm 7. It is easy to understand the correspondence of steps 1) and 2) of algorithm 7 to steps 1) and 2) of algorithm 4.
Steps 3) and 9) of algorithm 7 are the parallelized version of steps 3) and 9) (respectively) of algorithm 4 by using formula (21).
Step 4) of algorithm 7 is the parallelized version of steps 4) and 5) of algorithm 4 Steps 5) and (6) (Steps 10) and (11) ) of algorithm 7 are the parallelized version of step 6) (step 12)) of algorithm 4 when using formulas (22) and (23). The remaining steps of algorithm 7 are very easy to correspond to Algorithm 4. The VL only needs to communicate to prosumers the vector 0 , this in turn is of the length equal to the number of global constraints which is typically a small number, this, in turn, means that the method is communication efficient. The most demanding parts of the above algorithm from the computational point of view are the matrix-vector multiplications which can easily be conducted on CPUs or GPUs of each prosumer. Recently, due to reduced costs, lowprofile GPUs are becoming, indeed, increasingly affordable to prosumers. Indeed it is obvious that by parallelization the iteration number is still the same and the (1/ 2 ) convergence is still valid since the time for the communications between users and the VL is negligible since we need only to communicate scalar value and vectors ( ) from prosumers towards VL and back from VL to prosumers the vector 0 which are of length equal to the number of global constraints. Also, it is conceivable to apply differential privacy to the vectors ( ) by adding This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3187204 Gaussian noise which will not tangibly affect 0 = 1 1 ( )+. . . + ( ) − computed by the virtual layer. The proposed approach for privacy preservation can be categorized as data aggregation in which the different data from the meters are combined so that the adversarial agent or cyber attacker cannot identify an individual user's information. Also, this approach can be easily equipped with differential privacy.
The information flow between prosumers and the virtual layer for the parallelized private AGD for LP is presented in Fig. 3.

IV CASE STUDY AND NUMERICAL RESULTS
The baseline for the community is obtained by assuming that all batteries are charged or discharged in compliance with a typical control system used by manufacturers [11]. Considering the Italian regulation [4], [13], the ECM can present offers for up-regulation in the balancing market. The offers are subject to some constraints related to the time range (between 2 p.m. and 8 p.m.) in which the offers can be made and to the value of . , ( ) is equal to 210.58 €/MWh while , ( ) is equal to 40 €/MWh. The simulation results are obtained by using MATLAB software [46] on a PC having 24 GB of RAM and Intel (R) Core (TM) i7-4770k CPU@3.50 GHz 3.50 GHz.
Simulations have been carried out considering both the day-ahead and real-time optimizations and the scalability and privacy issues are confirmed in both cases for a different number of prosumers. As already explained, the proposed method is mainly designed for real-time optimization and the following results are obtained based on it.

A. Assessing the scalability of the real-time optimization
In this sub-section, the effectiveness of the proposed method is evaluated both in terms of accuracy and scalability. The scalability is assessed by evaluating its computation and communication time which is proportional to the number of iterations. The trend of normalized difference of the objective function and the normalized error related to the constraints are examined during different realtime optimizations. They are shown respectively in Fig. 4 and Fig. 5 for a single real-time optimization. The normalized difference of the objective function is defined as follows: where is the final value for the objective function and is the objective value at iteration number and {| |} Algorithm 7: Parallelized private AGD for LP Step 1 The virtual layer (VL) starts from = 1as well as an initial guess 0 = 0 and broadcasts = ‖ ‖ 2 as well as 0 to all prosumers 2 Each user sets = While ‖ ( )‖ is bigger than a threshold do (this is verified by the VL): 8 The VL does the update: 0 = 0 and each user sets: Each user computes ( ) = − 0 +̂+ and then does the projections ( ) = ( , ( )) and ( ) = ( , ( )) (based on formula (21)) 10 Each user computes the vectors =̂( ) − and sends the vector ( ) and scalar value to the VL 11 The VL computes 0 = 1 1 ( )+. . . + ( ) − and ‖ ( )‖ = √ 0 0 + 1 1 +. . . + (formula (23)) 12 The VL computes 0 = ( 0, 0 + 0 )

14
The VL computes = 2( −1) Since the proposed formulation needs all the constraints to be of inequality form ≤ , the normalized error related to them is defined as follows: where { 0, − } is zero if the i'th constraint is satisfied, otherwise, it is the amount of discrepancy related to the constraint. Fig. 4 shows the normalized difference between the objective function and Fig. 5 the normalized error of constraints considering 2000 prosumers for a single real-time optimization. It is worth mentioning that the warm start from the dual variables of the previous optimization allows obtaining an acceptable number of communications equal to 22. It should be mentioned that since the constraints ≤ are treated as soft ones the normalized error essentially decreases during iterations, but in some cases, it can also have a small increase to better satisfy the constraints and optimize the objective function, as shown in Fig. 3. A comparison of the proposed method with the centralized approach, considering the prosumers varying from 1000 to 5000, is shown in Table I. It is worth noting that the time necessary for each local computation in the proposed method is always lower than 1 ms and is, therefore, negligible if compared with the communication time. The computational time of the proposed method (and other parallelizable ones like SCS ADMM) is, thus, essentially almost equal to the communication time. To prove the efficiency and adequacy of the proposed method, it is, therefore, essential not only to evaluate and compare the value of the objective function and the normalized errors of constraints but also the number of communications required to converge to the optimal solution. The computational times and the number of communications are shown in Table I for the different approaches.   It is worth noting that, for the proposed method, the mean number of communications is equal to 14 and 42 communications for 1000 and 5000 prosumers, respectively. Considering a mean communication time of 10 ms for each communication, the mean communication time varies between 140 s and 420 s for each real-time optimization and is always lower than that of the centralized optimization.
A comparison in terms of the number of communications of the proposed method with SCS ADMM (both with a warm start), considering prosumers varying from 1000 to 2000, is shown in Table II. It is evident that, due to the lower number of communications, the proposed method converges much faster than SCS ADMM, as shown in Table II. As regards the achieved value of the objective function and the normalized errors of constraints, both methods achieve the same results. Also, as regards the privacy issue, it is worth noting that, even if the proposed method can keep all the variables private. The electrical powers absorbed or imported from/to the main grid should be transmitted at least every 60 seconds to the virtual layer to comply with the Italian regulation [47].

B. Evaluating the impact of the day-ahead forecast errors
To evaluate the potential impact that the errors related to the day-ahead forecast of the electrical load and the PV production can have on the real-time optimization, some simulations have been carried out considering both the day ahead and real-time optimizations. A day ahead forecast error for the load equal to 30% and the PV generation equal to 20% of the real values are assumed. The mean values in a quarter of an hour of the electrical power absorbed or imported from/to the main grid by an energy community of 5000 prosumers are shown in Fig. 6 for the day ahead and real-time optimizations. According to the results of the dayahead optimization, the ECM makes an offer to the TSO for energy flexibility in terms of Up-regulation from 2 p.m. to 6 p.m. with a constraint on the power flexibility of 2.5 MW and total energy flexibility of 10 MWh during the mentioned period. Also as shown in Fig. 6, despite the day-ahead forecast errors, the real-time optimization, by rescheduling the charging and/or discharging of the batteries and the energy exchanges between the prosumers, can precisely respect the offered flexibility, thus avoiding the payment of an economic penalty. The impact of the real-time forecast errors on the robustness of the solution obtained by the realtime optimization is evaluated by carrying out some Monte-Carlo simulations. The approach used for assessing the robustness of the real-time optimization is described and discussed in the Appendix. It is also worth noting that the adoption of exchanges in the energy community allows for reducing the costs of electrical energy for the prosumers and increasing the remuneration received for the provision of flexibility provision.

C. Evaluating the impact of energy exchanges over a day
To assess the impact that energy exchanges have on the provision of energy flexibility provided by an energy community and on the costs of electrical energy for the prosumers, a comparative study has been carried out and described in this sub-Section considering an energy community of 400 prosumers providing power flexibility for UP-regulation from 2 p.m. to 6 p.m.. According to the simulation results, shown in Table III and  Table IV, the adoption of energy exchanges in the energy community allows reducing the costs of electrical energy for the prosumers by about 7.2% with a total cost varying from around 242 €/day in the case without energy exchanges to around 225 €/day with energy exchanges. As regards the remuneration received for the provision of UP-regulation, energy exchanges determine an increase of the remuneration received by the energy community of about 3.5% with a variation from around 225 €/day without energy exchanges to around 232 €/day with energy exchanges. These economic benefits are due to the total energy exchanges of about 132.64 kWh that determine a decrease of the energy injected into the grid of about 92 kWh and a decrease of the energy imported from the grid of about 91 kWh. The mean values of the electrical power absorbed or imported from/to the main grid in a quarter of an hour by the energy community of 400 prosumers are shown in Fig. 7 considering the centralized optimization with and without energy exchanges.

V DISCUSSION
The simulation results demonstrated that, for large-scale problems involving a high number of prosumers, the proposed method for real-time optimization is superior in terms of privacy compared with the equivalent centralized optimization. To further demonstrate the efficiency and adequacy of the proposed method, it has been compared with SCS (a parallelizable and highly scalable version of ADMM) and simulation results, confirmed the higher convergence This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Although both methods can preserve the privacy of prosumers, since the proposed method allows reducing the number of communications, it is less exposed to the third party attempts to infer the private information of prosumers. Moreover, it can be easily programmed with low-level programming languages (like C++ or java) on a Raspberry platform or another similar one that represents a cheap solution for prosumers to allow the ECM to manage their batteries through the virtual layer. Being easy to be implemented it also exhibits high interoperability, which is a key property when trying to deploy the method in practice. By using the proposed method, an ECM, by carrying out the real-time optimizations with a higher frequency can be less exposed to forecast errors and sudden variations of electrical load or renewable energy production. In this way, the electrical energy-related costs of the energy community can be reduced, while the remuneration deriving from the provision of ancillary services in the ACM can be increased. The economic benefits to be shared between the prosumers and the ECM can therefore increase. The proposed method can be extended to handle different types of flexible resources, including electric vehicles and flexible loads, in an energy community. Future work will address these issues.

VI CONCLUSIONS
The issues related to the scalability and privacy of prosumers have gained much interest in recent years with regards to the realization of new energy communities. This is in part due to the rise of new emerging technologies (such as GPUs and Tensor Processing Units (TPUs) for parallel computing), as well as matured supporting mathematical tools, such as ADMM and more recently "LP solvers at scale". However, as far as we know, in the previous research, the issues of scalability and privacy of prosumers have been only partially explored, with studies mainly assessing a limited number of prosumers. Besides, the first intention behind ADMM, as well as many other parallelizable algorithms prevalent among data scientists, has not been the privacy of prosumers. Exploiting scalable and decentralized algorithms and, at the same time trying to keep the spatialtemporal information of prosumers private, is rather a more recent concept. Considering that LP problems should be solved extremely fast to manage real-time optimizations in presence of a very large number of prosumers, a new scalable and privacy-preserving method has been designed and applied to real-time distributed parallel optimization of a large-scale energy community. It can be used by an ECM to enable prosumers to provide ancillary services to the TSO by trading electrical energy with the grid in the ancillary services market. The carried-out study allows assessing the impact of energy exchanges on the flexibility services provided by an energy community.

A. Local constraints of the optimization
In this appendix, the local constraints that are considered in both the day-ahead and real-time optimizations are described. is the capacity of the battery of prosumer , ℎ / ℎ is the charging/discharging efficiency. Simultaneous charging and discharging of the energy storage system is avoided by implicit penalty in the objective function. The discharging and charging speeds of the battery are constrained by the following equations: 0 ≤ , ℎ ( ) ≤ , ℎ (31) 0 ≤ , ℎ ( ) ≤ , ℎ (32) ∀ ∈ [1, ], ∀ ∈ B. Evaluating the robustness of the real-time optimization Without losing the generality, it is assumed that the forecast errors in each quarter are normally distributed variables with a variance with a maximum percentage error of 30% for electrical load and 20% for PV production.
By comparing the proposed framework with the other robust frameworks, in the proposed formulation for real-time optimization, the uncertainty exists on the Right Hand Side (RHS) of the linear program optimization. It is worth noting that in general, while the row-wise uncertainty [48] is straightforward to deal with, tackling the issue of uncertainty is NP-hard when the uncertainty is in the RHS of a linear program [49,50], as it is the case in this paper. Also, the approach presented in [51] does not apply to this work since in practice some assumptions of the method are not respected in the considered very huge non-convex optimization problem. A practical approach based on the Monte-Carlo simulation is followed to assess the robustness of the proposed decision-making procedure and solutions similar to what has been done in [52] and [53].
Let's consider the LP optimization problem (eq. (7)): ≤ (33) ≤ ≤ where x represents the vector of decision variables and are the perturbed vectors due to forecast errors. By solving, for each sample of the Monte Carlo simulation, the following optimization problem: (34) − ≤ ≤ + , ∈ ≤ ≤ we can deduce the robustness of the solution by evaluating, for each constraint, the value that represents the allowed deviation in the decision variables, is the optimal value of the decision variable, and constant. Alternatively, if is empirically fixed to be equal to an allowed percentage of the optimal values and the optimization problem is reformulated as follows: ≤ ≤ by solving this optimization problem with i  equal to 5% of the optimal values it is possible to assess the variation of the decision variables. A Monte-Carlo simulation is carried out for 1000 iterations and 1000 prosumers. Simulation results confirmed the robustness of the solution obtained for the real-time optimization as demonstrated in Fig. 9 showing the maximum deviation of the discharging power of the battery , ℎ ( ) in Monte Carlo simulations. Three different sets of PV and demand forecast errors have been investigated and observed in Fig. 9. The conducted simulation shows that the highest deviation for , ℎ ( ) is below 4.3% for all cases. This demonstrates the robustness of the proposed approach that is also due to the limited values of the forecast errors. The proposed optimization is carried out every 15 minutes, so it is based on real-time forecasts of electrical loads and PV production. These forecasted values are, thus, sufficiently accurate and the level of forecast error is assumed to be limited. It is worth noting that in all the carried-out Monte-Carlo simulations, the values of state of charge for the battery as well as the objective function remain unchanged. It means that with some degree of perturbation of the decision variables it is still possible to achieve the same results although the input data are uncertain. This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3187204