Alexander Engelmann
email: alexander.engelmann@ieee.org

Sungho Shin
email: sshin@anl.gov

François Pacaud
email: francois.pacaud@minesparis.psl.eu

Victor M Zavala
email: victor.zavala@wisc.edu

Scalable Primal Decomposition Schemes for Large-Scale Infrastructure Networks

The real-time operation of large-scale infrastructure networks requires scalable optimization capabilities. Decomposition schemes can help achieve scalability; classical decomposition approaches such as the alternating direction method of multipliers (ADMM) and distributed Newtons schemes, however, often either suffer from slow convergence or might require high degrees of communication. In this work, we present new primal decomposition schemes for solving large-scale, strongly convex QPs. These approaches have global convergence guarantees and require limited communication. We benchmark their performance against the off-the-shelf interior-point method Ipopt and against ADMM on infrastructure networks that contain up to 300,000 decision variables and constraints. Overall, we find that the proposed approaches solve problems as fast as Ipopt but with reduced communication. Moreover, we find that the proposed schemes achieve higher accuracy than ADMM approaches.

Introduction

The operation of infrastructure networks such as power systems, district heating grids or gas networks is challenging. In many cases, these networks are large and composed of many complex subsystems such as lower-level networks or buildings. Operation is often based on numerical optimization due to its flexibility and recent advances in solver development, which allows to solve large-scale problems quickly and to a high accuracy. For large networks, however, a centralized solution is often not desirable since, a), the problem becomes computationally challenging, even with state-of-the-art solvers; b), information collection in a central entity should be avoided due to confidentiality and privacy concerns, and, c), the responsibility for operation and updates in modeling should stay mainly in the subsystems.

One line of research addresses the above challenges via aggregation. Here, the idea is to simplify the subproblems by projecting the constraint set on the coupling variables of the infrastructure network. Examples for this can be found for power systems [START_REF] Capitanescu | TSO-DSO interaction: Active distribution network power chart for TSO ancillary services provision[END_REF][START_REF] Kalantar-Neyestanaki | Characterizing the Reserve Provision Capability Area of Active Distribution Networks: A Linear Robust Optimization Method[END_REF]. A drawback of this approach is a loss of optimality. Moreover, aggregation is often not straightforward, feasibility is hard to guarantee and disaggregation requires solving additional local optimization problems.

A second line of research is based on distributed optimization. Prominent approaches are duality-based first-order algorithms such as Lagrangian dual decomposition and the Alternating Direction Method of Multipliers (ADMM) [START_REF] Everett | Generalized Lagrange multiplier method for solving problems of optimum allocation of resources[END_REF][START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF][START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF]. Application examples range from the operation of power systems [START_REF] Erseghe | Distributed Optimal Power Flow Using ADMM[END_REF][START_REF] Kim | A comparison of distributed optimal power flow algorithms[END_REF], over gas networks [START_REF] Shin | Graph-Based Modeling and Decomposition of Energy Infrastructures[END_REF], district heating systems [START_REF] Huang | Coordinated dispatch of electric power and district heating networks: A decentralized solution using optimality condition decomposition[END_REF][START_REF] Cao | Decentralized Operation of Interdependent Power Distribution Network and District Heating Network: A Market-Driven Approach[END_REF], to water networks [START_REF] Coulbeck | A hierarchical approach to optimized control of water distribution systems: Part I decomposition[END_REF]. With their at most linear rate of convergence, these approaches often require many iterations to converge even for a modest solution quality. This is often prohibitive for real-time implementation.

Distributed second-order methods exhibit faster convergence. Here, classical approaches aim at decomposing the block-structure of the Karush-Kuhn-Tucker (KKT) system within interior-point algorithms [START_REF] Chiang | Structured nonconvex optimization of large-scale energy systems using PIPS-NLP[END_REF][START_REF] Zavala | Interior-point decomposition approaches for parallel solution of large-scale nonlinear parameter estimation problems[END_REF] or sequential quadratic programming [START_REF] Varvarezos | Multiperiod design optimization with SQP decomposition[END_REF]. Alternative second-order methods based on augmented Lagrangians can be found in [START_REF] Engelmann | Toward Distributed OPF Using ALADIN[END_REF][START_REF] Houska | An Augmented Lagrangian Based Algorithm for Distributed NonConvex Optimization[END_REF]. As a drawback, these approaches typically require an expensive central coordination, although it is possible to partially alleviate the computational load by decentralizing the Newton-steps [START_REF] Engelmann | Decomposition of Nonconvex Optimization via Bi-Level Distributed ALADIN[END_REF][START_REF] Engelmann | An essentially decentralized interior point method for control[END_REF][START_REF] Stomberg | Decentralized non-convex optimization via bi-level SQP and ADMM[END_REF]. Moreover, ensuring global convergence is difficult in a distributed setting, since this requires monitoring a merit function centrally.

Primal decomposition approaches [START_REF] Geoffrion | Primal Resource-Directive Approaches for Optimizing Nonlinear Decomposable Systems[END_REF][START_REF] Demiguel | A local convergence analysis of bilevel decomposition algorithms[END_REF][START_REF] Demiguel | On Decomposition Methods for a Class of Partially Separable Nonlinear Programs[END_REF] come with the advantage of achieving a high degree of feasibility in a small number of iterations. Here, the idea is to partition the problem in a hierarchical fashion, i.e. to construct lower-level problems coordinated by one upper-level problem, where the upper-level problem considers the lower-level problems by their optimal value functions. Primal decomposition has been very successful in solving large-scale problems from chemical engineering [START_REF] Zavala | A fast moving horizon estimation algorithm based on nonlinear programming sensitivity[END_REF][START_REF] Yoshio | A Nested Schur decomposition approach for multiperiod optimization of chemical processes[END_REF] and some of the largest Quadratic Programs (QPs) and Nonlinear Programs (NLPs) from power systems [START_REF] Tu | A Two-Stage Decomposition Approach for AC Optimal Power Flow[END_REF][START_REF] Petra | Solving realistic security-constrained optimal power flow problems[END_REF][START_REF] Curtis | A Decomposition Algorithm for Large-Scale Security-Constrained AC Optimal Power Flow[END_REF]. Moreover, these approaches allow to use specialized, domain-specific solvers to solve the subproblems and the master problem efficiently [START_REF] Demiguel | A local convergence analysis of bilevel decomposition algorithms[END_REF].

In this work, we propose two primal decomposition schemes for solving large-scale strongly convex QPs, with global convergence guarantees. Both methods rely respectively on augmented Lagrangians and exact ℓ1-penalties for ensuring feasibility in the subproblems. Whereas similar ℓ1-penalty based approaches have been proposed in previous works [START_REF] Demiguel | A local convergence analysis of bilevel decomposition algorithms[END_REF], the augmented-Lagrangian framework is new to the best of our knowledge. We show that the augmented Lagrangian formulation exhibits improved performance compared to the ℓ1 formulation. Moreover, we demonstrate that the algorithms are faster than off-the-shelf interior-point solvers. We benchmark our algorithms against a distributed ADMM and the nonlinear solver Ipopt with sparse and parallel linear algebra. As benchmarks, we consider the operation of HVAC systems in a city district with a variable number of buildings and with up to 300, 000 decision variables and inequality constraints.

Notation

Given A ∈ R m×n , [A] j denotes the jth row of A and nr(A) . = m corresponds to the number of rows of A. The Lagrange multiplier λ ∈ R ng associated to the constraint function g : R nx → R ng is written as g(x) = 0 | λ. Given a vector v ∈ R n , D = diag(v) ∈ R n×n is a matrix with the elements of v on the main diagonal. For a tuple of matrices (A, B), blkdiag(A, B) denotes their block-diagonal concatenation.

Problem Formulation

Many infrastructure network problems can be formulated as strongly convex QPs over a set of subsystems S = {1, . . . , S},

min {x i } i∈S ,y i∈S 1 2 x i y ⊤ H xx i H xy i H xy⊤ i H yy i x i y + h x i h y i ⊤ x i y (1a) subject to A x i A y i x ⊤ i y ⊤ ⊤ -b i = 0, i ∈ S, (1b)
B x i B y i x ⊤ i y ⊤ ⊤ -d i ≤ 0, i ∈ S, (1c)
A y y -b y = 0, B y y -d y ≤ 0. (1d)
Here, the global decision variable vector

x = [x 1 , . . . , x S ] ⊤ is composed of local decision variables x i ∈ R n xi
, where each x i belongs to one subsystem i ∈ S. The decision variables y ∈ R ny are "global" in the sense that they belong to the interconnecting infrastructure network, described by the constraints (1d). Each coefficient matrix/vector in the objective (1a) and the constraints (1b), (1c) belongs to one i ∈ S.

Primal Decomposition Schemes

In contrast to duality-based techniques such as ADMM or dual decomposition, primal decomposition decomposes entirely in the primal space, i.e. no dual variables are updated in the solution process. The main idea here is to replace the subproblems in (1) by their optimal value functions. Specifically, one reformulates (1) as min

y i∈S ϕ i (y), s.t. A y y -b y = 0, B y y -d y ≤ 0, (2)
where for all i ∈ S, the value function ϕ i is defined as

ϕ i (y) . = min x i 1 2 x i y ⊤ H xx i H xy i H xy⊤ i H yy i x i y + h x i h y i ⊤ x i y (3a) subject to A x i A y i x ⊤ i y ⊤ ⊤ -b i = 0, i ∈ S, (3b)
B x i B y i x ⊤ i y ⊤ ⊤ -d i ≤ 0, i ∈ S. (3c)
The key idea is to apply standard algorithms for solving (2) by optimizing only with respect to the coupling variables y. Doing so can lead to enhanced robustness, as the complexity of the subproblems is not exposed to the algorithm solving [START_REF] Kalantar-Neyestanaki | Characterizing the Reserve Provision Capability Area of Active Distribution Networks: A Linear Robust Optimization Method[END_REF]. Algorithms for solving (2) typically require first-order and possibly second-order derivatives of all {ϕ i } i∈S . Since all {ϕ i } i∈S are non-smooth because of the inequality constraints, one typically relies on smooth reformulations. Inspired by interior-point methods [START_REF] Demiguel | A local convergence analysis of bilevel decomposition algorithms[END_REF], we introduce log-barrier functions and slack variables s i ∈ R n ii , which approximate (3) by

Φ δ i (y) . = min x i ,s i 1 2 x i y ⊤ H xx i H xy i H xy⊤ i H yy i x i y + h x i h y i ⊤ x i y -δ1 ⊤ ln(s) (4a) subject to A x i A y i x ⊤ i y ⊤ ⊤ -b i = 0, i ∈ S, (4b)
B x i B y i x ⊤ i y ⊤ ⊤ -d i + s i = 0, i ∈ S, (4c)
where δ ∈ R + is a barrier parameter, 1 . = [1, . . . , 1] ⊤ , and the ln(•) is evaluated component-wise. Note that lim δ→0 Φ δ i (y) = ϕ i (y), and that Φ δ i is smooth1 . A basic primal decomposition strategy with smoothing is summarized in Algorithm 1.

Algorithm 1: A basic primal decomposition scheme.

Initialize y 0 , δ 0 . while not terminated do 1) Solve (2) for ϕ i ≡ Φ δ i with a NLP solver; in case the NLP solver calls

(∇ y Φ δ i , ∇ 2 y Φ δ i ), compute them locally for all i ∈ S. 3) Decrease δ. end Return y k , {x k i } i∈S .

Dealing With Infeasibility

An issue in Algorithm 1 is that the subproblems (3) may be infeasible for a given y.

One way of circumventing this is to introduce auxiliary variables z i ∈ R ny and to use relaxation techniques either based on augmented Lagrangians or on exact ℓ1-penalties. Consider a set of auxiliary variables {z i } i∈S and introduce additional constraints z i = y for all i ∈ S. Then, we reformulate (4) equivalently as2

Φ δ i (y) = min x i ,s i ,z i 1 2 x i y ⊤ H xx i H xy i H xy⊤ i H yy i x i y + h x i h y i ⊤ x i y -δ1 ⊤ ln(s) (5a) subject to A x i A y i x ⊤ i z ⊤ i ⊤ -b i = 0, i ∈ S, (5b)
B x i B y i x ⊤ i z ⊤ i ⊤ -d i + s i = 0, i ∈ S, (5c)
z i = y, i ∈ S, (5d)
which can still be infeasible, but paves the way for augmented Lagrangian and exact ℓ1 relaxations.

Augmented Lagrangian Relaxation

A simple way of making [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF] feasible for all y is to relax (5d) via a quadratic penalty [START_REF] Demiguel | A local convergence analysis of bilevel decomposition algorithms[END_REF]. However, in this case, large penalty parameters might lead to numerical difficulties and feasibility can in general not be guaranteed for a finite penalty parameter. Hence, we use an Augmented Lagrangian (AL) approach to solve (2) for late outer iterations with a constant barrier parameter δ. Assigning the Lagrange multiplier λ i to (5d), we relax (5d) in an Augmented Lagrangian fashion by adding the terms λ k⊤ i (y

-z i ) + ρ 2 ∥y -z i ∥ 2 2 to the objective: Φ δ,ρ i (y, λ k i ) . = min x i ,s i ,z i 1 2   x i y z i   ⊤   H xx i H xy i 0 H xy⊤ i H yy i + ρI -ρI 0 -ρI ρI     x i y z i   +   h x i h y i + λ k i -λ k i   ⊤   x i y z i   -δ1 ⊤ ln(s i ) subject to A x i A y i x ⊤ i z ⊤ i ⊤ -b i = 0, | γ i B x i B y i x ⊤ i z ⊤ i ⊤ + s i -d i = 0, | µ i . (6)
Here, (γ i , µ i ) are Lagrange multipliers corresponding to the constraints in the same line.

Note that the solution of (5) can be forced to be arbitrarily close to that of (6) by letting ρ → ∞, if (5) is feasible for y. In addition, if one has a good Lagrange multiplier estimate λ k ≈ λ ⋆ , ( 5) and ( 6) are equivalent for a finite ρ < ∞ [28, Sec. 3.2.1], [START_REF] Nocedal | Numerical Optimization[END_REF]Thm. 17.5]. We will exploit this fact in the following. Primal decomposition based on the augmented Lagrangian works as follows: In phase 1, the barrier parameter δ and the penalty parameter ρ are increased/decreased simultaneously. In phase 2, when δ/ρ are sufficiently small/large, both are held constant and a standard augmented Lagrangian algorithm is applied to the resulting optimization problem to obtain feasibility in (5d). We use the standard first-order update rule from augmented Lagrangian algorithms [29, Chap. 17.3]

λ k+1 i = λ k i + ρ(y k -z k i ). ( 7
)
The resulting scheme is summarized in Algorithm 2.

Algorithm 2: AL-based primal decomposition.

Initialize y 0 , δ 0 , ρ 0 ; λ i = 0, i ∈ S. while phase 1 do 1) Solve (2) for ϕ i ≡ Φ δ,ρ i with a NLP solver; in case the NLP solver calls (∇ y Φ δ,ρ i , ∇ 2 yy Φ δ,ρ i ), compute them locally for all i ∈ S. 3) Decrease δ, increase ρ. end while phase 2 do 1) Solve (2) as in phase 1.

2) Update λ k i according to [START_REF] Kim | A comparison of distributed optimal power flow algorithms[END_REF]. end Return y k , {x k i } i∈S .

ℓ1-penalty relaxation

A second variant to ensure feasibility is to relax (5d) via an ℓ1-penalty function. This has the advantage that it is exact also for a finite penalty parameter λ ∈ R + without the need for Lagrange-multiplier estimation. By doing so, the objective becomes nonsmooth. However, the non-smoothness can be eliminated by using an elastic relaxation [29, p.535]: the ℓ1-penalty min y,z i ∥y -z i ∥ 1 is reformulated by introducing two nonnegative auxiliary variables v i , w i ∈ R ny

+ as min y,z i ,v i ,w i v i + w i subject to y -z i = v i -w i . The corresponding reformulation of (5) reads Φ δ, λ i (y) . = min x i ,s i ,z i v i ,w i 1 2 x i z i ⊤ H xx i H xy i H xy⊤ i H yy i x i z i + h x i h y i ⊤ x i z i + λ1 ⊤ (v i + w i ) -δ(1 ⊤ ln(s i ) + 1 ⊤ ln(v i ) + 1 ⊤ ln(w i )) (8) subject to A x i A y i x ⊤ i z ⊤ i ⊤ -b i = 0, | γ i , B x i B y i x ⊤ i z ⊤ i ⊤ + s i -d i = 0, | µ i , y -z i -v i + w i = 0, | χ i ,
where the bounds (v i , w i ) ≥ 0 are replaced by log-barrier functions. If one chooses λ < ∞ large enough, ( 5) and ( 8) are equivalent [START_REF] Nocedal | Numerical Optimization[END_REF]Thm 17.3].

Primal decomposition based on the ℓ1-penalty solves (2) using Φ δ, λ i with a fixed λ larger than a certain threshold and decreases the barrier parameter δ during the iterations. The overall algorithm is summarized in Algorithm 3.

Algorithm 3: ℓ1-based primal decomposition.

Initialize y 0 , δ 0 , λ large enough. while not terminated do 1) Solve ( 2) for ϕ i ≡ Φ δ, λ i with a NLP solver; in case the NLP solver calls (∇ y ϕ i , ∇ 2 yy ϕ i ), compute them locally for all i ∈ S.

3) Decrease δ. end Return y k , {x k i } i∈S .

Computing sensitivities

Next, we show how to compute ∇ y Φ δ i and ∇ 2 yy Φ δ i under standard regularity assumptions based on the implicit function theorem [START_REF] Demiguel | On Decomposition Methods for a Class of Partially Separable Nonlinear Programs[END_REF]. Reformulate (4) by

Φ δ i (y) = min x i ,s i f δ i (x i , s i ; y) (9)
subject to

g i (x i ; y) = 0 | γ i , h i (x i ; y) + s i = 0 | µ i ,
where f δ i is defined by (4a), and g i and h i are defined by (4b), (4c). Define the Lagrangian to [START_REF] Huang | Coordinated dispatch of electric power and district heating networks: A decentralized solution using optimality condition decomposition[END_REF],

L δ i (x i , s i , γ i , µ i ; y) . = f δ i (x i , s i ; y) + γ ⊤ i g i (x i ; y) + µ ⊤ i (h i (x i ; y) + s i ).
Assume that ( 9) is feasible for a given y and that regularity conditions such as LICQ, SOSC and SCC hold, cf. [START_REF] Demiguel | On Decomposition Methods for a Class of Partially Separable Nonlinear Programs[END_REF]. Then, the KKT conditions to (9) form an implicit function in form of

F δ i (x ⋆ i , s ⋆ i , γ ⋆ i , µ ⋆ i ; y) = 0
, where the superscript (•) ⋆ indicates a KKT stationary point. Thus, by the implicit function theorem, there exist a neighborhood around y for which there exists functions p ⋆ i (y)

. = (x ⋆ i (y), s ⋆ i (y), γ ⋆ i (y), µ ⋆ i (y)) such that F δ i (p ⋆ i (y); y) = 0. Hence, we can rewrite (9) as Φ δ i (y) = f i (x ⋆ i (y), s ⋆ i (y); y) = L i (p ⋆ i (y); y) since p ⋆ i (y) is feasible.
Applying the total derivative and the chain rule yields

∇ y Φ δ i (y) =∇ y L δ i (p ⋆ i (y); y)+∇ p ⋆ i L δ i (p ⋆ i (y); y)∇ y p ⋆ i (y).
By the KKT conditions, we have that ∇ p ⋆ i L δ i (p ⋆ i (y); y) = 0 and thus

∇ y Φ δ i (y) =∇ y L δ i (p ⋆ i (y); y). (10)
Again by the total derivative, the Hessian can be computed by

∇ 2 yy Φ δ i (y) = ∇ 2 yy L i (p ⋆ i (y); y) + ∇ 2 yp ⋆ i L i (p ⋆ i (y); y) ∇ y p ⋆ i (y). ( 11
)
It remains to derive an expression for ∇ y p ⋆ i (y). The KKT conditions of (9) read

F δ i (x ⋆ i , s ⋆ i , γ ⋆ i , µ ⋆ i ; y) =     ∇ x i f i (x ⋆ i , y) + ∇ x i g i (x ⋆ i , y)γ ⋆ i + ∇ x i h i (x ⋆ i , y)µ ⋆ i -δS ⋆-1 i 1 + µ ⋆ i g i (x ⋆ i , y) h i (x ⋆ i , y) + s ⋆ i     ! = 0,
where

S ⋆ i = diag(s ⋆ i )
. By the total differential and the chain rule we have

∇ y F δ i (p ⋆ i (y), y)+ ∇ p ⋆ i F δ i (p ⋆ i (y), y)∇ y p ⋆ i (y) = 0.
Hence, we can compute the Jacobian ∇ y p ⋆ i (y) by solving the system of linear equations

∇ p ⋆ i F δ i (p ⋆ i (y), y) ∇ y p ⋆ i (y) = -∇ y F δ i (p ⋆ i (y), y). (12)
Observe that ( 12) is a system of linear equations with multiple right-hand sides. In summary, we can compute ∇ 2 yy Φ δ i (y) locally for each i ∈ S by combining ( 11) and ( 12). The corresponding formulas for the gradient and the Hessian of Φ δ,ρ i and Φ δ, λ i from ( 6) and ( 8), i.e. of the AL relaxation and the ℓ1 relaxation ( 9) are given in Appendix A.

Solving the Master Problem and Globalization

An important question is how to solve the master problem (2) for different variants of ϕ i . In general, this can be done by any sensitivity-based NLP solver. We proceed by showing how to obtain a simple globalized version of Algorithm 1 based on a linesearch scheme; here, the idea is to show global convergence for the relaxed problem (2) with ϕ i ∈ {Φ δ,ρ i , Φ δ, λ i } for fixed penalty and barrier parameters. This leads to converge of a solution to the original problem (1) by standard results from penalty and barrier methods [START_REF] Nocedal | Numerical Optimization[END_REF]Thms. 17.1,17.6].

Define the objective of (2), ψ(y) . = i∈S ϕ i (y), as a global merit function, where

ϕ i ∈ {Φ δ,ρ i , Φ δ, λ i }.
The basic idea is to employ a Sequential Quadratic Programming (SQP) scheme, where we ensure a sufficient decrease in ψ at each step via the Armijo condition. The overall algorithm is summarized in Algorithm 4. Similar to the general primal decomposition scheme from Algorithm 1, the master problem solver evaluates the sensitivities (∇ y ϕ i , ∇ 2 yy ϕ i ) in step (i), in order to construct a quadratic approximation of (2) in step (ii). Solving this approximation yields a search direction ∆y. The stepsize α is updated with a backtracking line-search with the Armijo condition as termination criterion.

Global Convergence

We now establish global convergence 3 of Algorithm 4 to a minimizer of the relaxed problems [START_REF] Kalantar-Neyestanaki | Characterizing the Reserve Provision Capability Area of Active Distribution Networks: A Linear Robust Optimization Method[END_REF]

for ϕ i ∈ {Φ δ,ρ i , Φ δ, λ i }.
Assume that the following regularity assumptions Algorithm 4: A simple master problem solver.

Initialize y 0 , ϵ, δ, ζ ∈ (0, 1), ρ or λ, ϕ i ∈ {Φ δ,ρ i , Φ δ, λ i }. while ∥∆y∥ > 0 do (i) compute (∇ y ϕ i , ∇
. = [p ⋆ i (y ⋆ )] i∈S . Assumption 1 (Regularity). Assume that for all i ∈ S a) ∆x ⋆ i ∆y ⋆ ⊤ H xx i H xy i H xy⊤ i H yy i ∆x ⋆ i ∆y ⋆ > 0, for all ∆x ⋆ i ∆y ⋆ ̸ = 0 with A x i A y i ∆x ⋆ i ∆y ⋆ = 0; b) A x i A y i B x i B y i has full row rank; c) [µ ⋆ i ] j + [B x i B y i ][x ⋆⊤ i y ⋆⊤ ] ⊤ j ̸ = 0, ∀j = 1, . . . , nr(B x i ).
Line-search methods require that the search direction ∆y is a descent direction, i.e. ∆y i∈S ∇ i ϕ δ i (y) < 0. This can be ensured by showing that i∈S ∇ 2 yy ϕ δ i ≻ 0 for all variants of ϕ i , which we do with the next lemma.

Lemma 1 (Positive definite Hessians). Let Assumption 1 hold and assume that

(s i , µ i , v i , w i ) > 0. Then, a), the Hessian ∇ 2 yy Φ δ,ρ i is positive definite for all ρ > 0. Moreover, b), ∇ 2 yy Φ δ, λ i is positive definite if λ is larger than all multipliers associated to the elastic constraints λ > max j |[χ ⋆ i ] j |.
The proof of Lemma 1 is given in Appendix B. Now we are able to show global convergence of Algorithm 4 to the solution of problem (2) with ϕ i ∈ {Φ δ,ρ i , Φ δ, λ i } for fixed penalty and barrier parameters.

Theorem 1 (Convergence of Algorithm 4 for fixed δ, ρ, λ). Consider Algorithm 4 with either fixed (δ, ρ, λ) if ϕ i = Φ δ,ρ i or with fixed (δ, λ) if ϕ i = Φ δ, λ i . Let Assumption 1 and the conditions of Lemma 1 hold. Then, Algorithm 4 converges to the global minimizer of problem [START_REF] Kalantar-Neyestanaki | Characterizing the Reserve Provision Capability Area of Active Distribution Networks: A Linear Robust Optimization Method[END_REF]

with ϕ i ∈ {Φ δ,ρ i , Φ δ, λ i }.
Proof. The unconstrained minimizer to [START_REF] Zavala | Interior-point decomposition approaches for parallel solution of large-scale nonlinear parameter estimation problems[END_REF], ∆y = -i∈S ∇ 2 yy ϕ -1 i (y) i∈S ∇ y ϕ i (y), is a descent direction for the merit function ψ(•) since i∈S ∇ y ϕ i (y) ⊤ ∆y < 0 by the positive definiteness of the Hessians from Lemma 1. Observe that ∆y = 0 is feasible for (13) by feasibility of y. This shows that either ∆y = 0, or ∆y is a descent direction. Hence, by [START_REF] Nocedal | Numerical Optimization[END_REF]Lem 3.1], there exists an α ∈ (0, 1] such that ψ(y) -ψ(y + α∆y) ≥ -σα∇ y ψ(y) ⊤ ∆y is satisfied and thus the inner while loop is well defined. Moreover, by the convergence of line-search methods [29, Thm 3.2], Algorithm 4 will either converge to a stationary point of ψ returning ∆y = 0. Alternatively, if the unconstrained search direction is blocked by the constraints in [START_REF] Zavala | Interior-point decomposition approaches for parallel solution of large-scale nonlinear parameter estimation problems[END_REF], ∆y = 0 is returned since ∆y = 0 is feasible.

We now show that y is optimal for (2) if ∆y = 0. The KKT conditions associated to (13) read

         i∈S ∇ 2 yy ϕ i (y)∆y + ∇ y ϕ i (y) ⊤ + A y⊤ γ y + B y⊤ µ y = 0 A y (y + ∆y) -b y = 0, B y (y + ∆y) -b y ≤ 0, (B y (y + ∆y) -b y ) ⊤ µ = 0, µ ≥ 0,
which are precisely the KKT conditions for (2) if ∆y = 0. Since ( 2) is convex and Assumption 1 holds, the assertion follows.

Combining Theorem 1 with the convergence results for the ℓ1-penalty method [29, Them 17.3] or the augmented Lagrangian method [30, Prop 2.7] implies convergence of Algorithm 4 to the minimizer of the original problem (1) for sufficiently large penalties.

Remark 1 (Satisfying the assumptions of Theorem 1). Observe that the assumptions for Theorem 1 are standard regularity assumptions from nonlinear programming [22, Ass. 1-4]. Moreover, (s i , µ i , v i , w i ) > 0 is always ensured when using interior-point solvers for solving (5) and (6) even in the case of early termination.

Implementation Aspects Solving Local Subproblems

The evaluation of the sensitivities of ϕ i ∈ {Φ δ,ρ i , Φ δ, λ i } requires solving local optimization problems [START_REF] Erseghe | Distributed Optimal Power Flow Using ADMM[END_REF] or [START_REF] Shin | Graph-Based Modeling and Decomposition of Energy Infrastructures[END_REF] for fixed δ, ρ, λ. Observe that this can be done using specialized and optimized interior-point solvers, if they allow termination once a certain barrier δ is reached. Moreover, interior-point solvers factorize the KKT matrices ∇ p i F δ i (cf. ( 23), ( 26)) at each inner iteration and these factorizations can be re-used for Hessian computation via [START_REF] Chiang | Structured nonconvex optimization of large-scale energy systems using PIPS-NLP[END_REF]. Here we provide two variants: our own interior-point QP solver In early iterations, it is typically not necessary to solve the local problems to a high precision, since the barrier parameter δ are still large and the penalty parameters (ρ, λ) are still small. Hence, we solve the subproblems to an accuracy measured in the violation of the optimality conditions

i∈S v k i ≤ v v k 1 v k |S| . . .
∥F δ i (p k i , y k )∥ ∞ and terminate if ∥F δ i (p k i , y k )∥ ∞ < min(δ, 1/ρ) or ∥F δ i (p k i , y k )∥ ∞ < min(δ, 1/ λ)
. This is inspired by the termination of inexact interiorpoint methods [START_REF] Byrd | On the Local Behavior of an Interior Point Method for Nonlinear Programming[END_REF]. Warm-starting the local solves with the solution of the previous iteration reduces computation time significantly.

Numerical Linear Algebra

Efficient numerical linear algebra is crucial for performance. The most expensive steps in terms of memory and CPU time are the factorization of the local KKT matrices and the backsolves in [START_REF] Chiang | Structured nonconvex optimization of large-scale energy systems using PIPS-NLP[END_REF]. Since we only consider QPs, large parts of the KKT matrices in [START_REF] Chiang | Structured nonconvex optimization of large-scale energy systems using PIPS-NLP[END_REF] are constant over the iterations, which can be exploited for pre-computation. Here, we make heavy use of the Schur complement. Details on how to achieve this for the AL formulation are given in Appendix A.

Updating Penalties

We use simple update rules for the penalty parameters. They are

δ k+1 = 0.2δ k , λk+1 = 2 λk , ρ k+1 = 3ρ k ,
with initial values δ 0 = 0.1, ρ 0 = 10 3 , and λ0 = 100. After a fixed amount of 8 iterations, phase 2 starts and the above values stay constant.

Numerical Case Studies

We consider an optimal control problem for a city district with a scalable number of commercial buildings connected via a electricity grid with limited capacity. The building data is from [START_REF] Rawlings | Economic MPC and real-time decision making with application to large-scale HVAC energy systems[END_REF]. We neglect the waterside HVAC system and assume that the buildings are equipped with heat pumps with a constant coefficient of performance.

Building Model

The evolution of the temperature of the mth zone in the ith building reads

C m,i (T k+1 m,i -T k m,i ) = -H m,i (T k m,i -T k a ) - n∈Z i G mn,i (T k m,i -T k n,i ) -Qck m,i + Qdk m,i , (14)
where at time step k, T k m,i is the temperature of zone m and T k a the ambient temperature, C i,m is the thermal capacity, H m,i and G mn,i are heat transfer coefficients with the ambient and between two zones. Moreover, ( Qck m,i , Qdk m,i ) are the controllable/uncontrollable heat influxes from the heat pump and from sources of disturbance such as solar irradiation and occupancy. Eq. ( 14) can be written in compact form as

C i (T k+1 i -T k i ) = -H i (T k i -T k a ) -G i T k i -Qck i + Qdk i ,
where

T ⊤ i . = [T 1,i , . . . T |Z i |,i ], Qck ⊤ i . = [
T k+1 i . = z k+1 i = I -C -1 i (H i + G i ) T k i -C -1 i Qck i + C -1 i [H i I][T k⊤ a Qdk⊤ i ] ⊤ (15) . = A i x k + B i u k i + E i d k i . (16)
Stacking the above over N time steps yields zi = 0 0

I N -1 ⊗ A i 0 zi + 0 I N -1 ⊗B i ūi + 0 I N -1 ⊗ E i di +z f i . = Āi zi + Bi ūi + Ēi di + zf i ,
where z⊤

i . = [z 0⊤ i , . . . , z N ⊤ i ], ū⊤ i . = [u 1⊤ i , . . . , u N ⊤ i ], d⊤ i . = [d 1⊤ i , . . . , d N ⊤ i ], zf⊤ i .
= [z f i , 0 ⊤ , . . . , 0 ⊤ ] and z f i are the initial temperatures. Define the total energy consumption of building i ∈ S at time step k by v k i = 1 ⊤ u k i , and v⊤

i . = [v 0⊤ i , . . . , v N ⊤ i ].
Then, the above is equivalent to

(I -Āi ) -Bi 0 z⊤ i ū⊤ i v⊤ i ⊤ = Ēi di + zf i . (17)
The grid coupling between all subsystems i ∈ S induces an upper-bounded energy supply writing as a global constraint:

1 ⊤ v k i ≤ v for all times k. Moreover, we have local comfort constraints T 1 ≤ z k i ≤ 1 T .

Optimal Control Problem

Assume that the goal of each building is to minimize its cost of energy consumption over N time steps respecting all constraints. The cost function to buy the power form the utility is given by

f k i (u k i ) . = 0.5 c k (u k i ) 2 +g k u k i
, where g k is a linear and c k > 0 is a (small) quadratic cost coefficient. This yields a discrete-time Optimal Control Problem (OCP) for building i,

ϕ i (v i ) . = min zi ,ū i ,v i 1 2   zi ūi vi   ⊤   0 0 0 0 cI 0 0 0 0     xi ūi vi   +   0 1 ⊤ ⊗ g 0   ⊤   zi ūi vi   subject to (17), I 0 0 -I 0 0   zi ūi v i   ≤ T 1 -T 1 , vi = I ⊗ 1 ⊤ ūi .
The overall OCP -including global grid constraints -reads min

v 1 ,...,v |S| i∈S ϕ i (v i ), subject to (1 ⊤ ⊗ I) [v ⊤ 1 , . . . , v ⊤ |S| ] ⊤ ≤ v 1. (18)
To obtain a problem in form of (1), define

x ⊤ i = [x ⊤ i , ū⊤ i , v⊤ i ], i ∈ S, y ⊤ = [v ⊤ 1 , . . . , v ⊤ |S| ], H x i = blkdiag(0, cI, 0), h x i = [0, 1 ⊤ ⊗g, 0] ⊤ , h y i = 0, H xy i = 0, B y i = 0, d ⊤ i = [ T 1 ⊤ , -T 1 ⊤ ], A x i = (I -Āi ) -Bi 0 0 I ⊗ 1 ⊤ 0 , A y i = 0 -e ⊤ i ⊗I , B x i = I 0 0 -I 0 0
,

where e i is the ith unit vector, b ⊤ i = ( Ēi di + zf i ) ⊤ 0 , A y = 0, b y = 0, B y = 1 ⊤ ⊗ I, and d y = v1.

Technical Setup

We benchmark our algorithms against ADMM (as one of the most popular algorithms for decomposition) and against Ipopt (as one of the most prominent centralized NLP solvers). 4 The particular variant of ADMM is described in Appendix C. We rely on OSQP v0.6.2 [START_REF] Stellato | OSQP: An operator splitting solver for quadratic programs[END_REF] for solving subproblems and the coordination problem in ADMM. In primal decomposition, we rely on our own interior-point solver for the subproblems and on Algorithm 4 for coordination, where we solve (13) via Ipopt.

We perform all simulations on a shared-memory virtual machine with 30 cores and 100GiB memory. The underlying hardware is exclusively used for the case studies. All algorithms are parallelized via Julia multi-threading. Ipopt v3.14.4 runs partially in parallel via the parallel multi-frontal linear algebra solver MUMPS v5.4.1.

Numerical Results

We compare the numerical performance of all algorithms on OCP [START_REF] Engelmann | An essentially decentralized interior point method for control[END_REF] for |S| ∈ {30, 180, 300} buildings. Table 2 shows the corresponding number of local/global decision variables n x /n y , the number of local equality/inequality constraints n e /n i , and the number of global equality/inequality constraints n ey /n iy . We employ ADMM from Appendix C with penalty parameters ρ ∈ {10 1 , 10 2 , 10 3 , 10 4 }. 2 illustrates the numerical performance of both primal decomposition variants and ADMM. Figures 34show the AL formulation only, since the ℓ1 formulation runs out of memory for these cases. The constraint violations for the equality constraints (1b),

[A x A y ][x ⊤ y ⊤ ] ⊤ -b ∞ , for the inequality constraints (1c), max([B x B y ][x ⊤ y ⊤ ] ⊤ - b)
, and the value of the cost function f ([x ⊤ , y ⊤ ] ⊤ ) from (1a) are displayed, where the xaxis shows the iteration count. One can observe that the primal decomposition schemes achieve a high degree of feasibility in less than 10 iterations for all cases. Moreover, the optimality gap (f ([x ⊤ , y ⊤ ] ⊤ )-f ([x ⋆⊤ , y ⋆⊤ ] ⊤ ))/f ([x ⋆⊤ , y ⋆⊤ ] ⊤ ) is below 0.01% in less than 10 iterations for both primal decomposition variants and for all |S|, where [x ⋆⊤ , y ⋆⊤ ] ⊤ is computed via Ipopt. For ADMM, infeasibility and the optimality gap stay large independently of the choice of the tuning parameter ρ.

Remark 2 (Scaling of the ℓ1 formulation). The reason for the poor scaling of the ℓ1formulation is two-fold: First, the relaxation (8) introduces 2n y additional slack variables and inequality constraints. Hence, the KKT system in the subproblems defined via [START_REF] Petra | Solving realistic security-constrained optimal power flow problems[END_REF] has a larger size than the KKT system we get with the AL formulation [START_REF] Zavala | A fast moving horizon estimation algorithm based on nonlinear programming sensitivity[END_REF]. Moreover, the additional inequality constraints potentially lead to smaller stepsizes due to the fraction-to-boundary rule [START_REF] Nocedal | Numerical Optimization[END_REF]Eq 19.9]. Hence, more iterations in the subproblems are required compared to the AL formulation.

Discussion of Algorithmic Properties

Next, we discuss algorithm properties in view of the desirable properties from Section 1.

Computation Times

Comparing computation times between algorithms is difficult, since the numerical performance strongly depends on the implementation. Nonetheless, we would like to provide some timing information to underline the potential of primal decomposition methods. Table 2 shows the computation times and the number of iterations for all algorithms.

Ipopt is terminated at optimality with default settings and the computation time for the primal decomposition schemes are evaluated once an optimality gap of 10 -4 and a maximum constraint violation for equality/inequality constraints of 10 -5 is reached. ADMM is terminated after a fixed amount of 100 iterations. One can observe that the AL formulation and Ipopt with parallel linear algebra have similar computation times independently of |S|. Although fast computation is not our primary focus here, this indicates the potential of primal decomposition for large-scale optimization.

Ipopt requires more than 100 iterations in all cases, which is unusual for second-order interior-point methods. Both ADMM and the open-source solver OSQP are not able to solve these problems to a sufficient accuracy with a reasonable number of iterations (4,000 for OSQP). This indicates that these problems are rather difficult, which might be due to the large number of inequality constraints coming from the temperature bounds. On the other hand, the primal decomposition schemes require far less iterations. A possible explanation for the behavior of Ipopt is that interior-point methods typically [START_REF] Nocedal | Numerical Optimization[END_REF]Eq 19.9]. This can lead to a slow progress since in this case only small steps are taken. In primal decomposition schemes, each subproblem chooses "its own" stepsize independently which might be one explanation for faster progress.

Internal timings for the AL formulation and different sizes |S| are shown in Table 3. Here, the time spent in the coordination problem [START_REF] Zavala | Interior-point decomposition approaches for parallel solution of large-scale nonlinear parameter estimation problems[END_REF] and in the local solvers stays relatively constant for varying |S|. The time spent for sensitivity computation, however, increases significantly with |S|. One explanation for that is that the complexity in the backsolves for computing the Hessian via [START_REF] Chiang | Structured nonconvex optimization of large-scale energy systems using PIPS-NLP[END_REF] increases with O(n 3 y ). In case the backsolves can be parallelized, for example if the subproblems themselves use multithreading in a cluster environment, computation time can be reduced.

Communication Requirements

Both primal decomposition schemes require the communication of the gradient and the Hessian of the optimal value functions in each outer iteration. In the line-search, only the current iterate has to be communicated to the subsystems and the resulting optimal value of the local problem is communicated back. Hence, the required amount of communication in one iteration is higher than in ADMM, where only the current iterates of the coupling variable y has to be communicated. However, the difference in total communication-i.e. the number of bits to be communicated until a certain solution accuracy is reached-is still often smaller since the primal decomposition schemes require a significantly smaller number of iterations for a given target solution quality.

Feasibility and Optimality

Reaching feasibility fast is often crucial in the context of infrastructure systems to ensure system stability. Here, primal decomposition can shine. For all case studies, one can observe a high degree of primal feasibility in 10-20 iterations. ADMM requires far more iterations to reach a sufficient degree of feasibility, which is a known limitation of ADMM [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF].

Conclusion and Outlook

We have presented two primal decomposition schemes to solve large-scale QPs for the operation of infrastructure networks. The developed methods are proven to converge globally to the optimal solution. Numerical experiments have demonstrated their potential for solving large-scale QPs in a small number of iterations to a high degree of feasibility and optimality, which distinguishes them from classical distributed methods such as ADMM. Moreover, we have shown that primal decomposition based on augmented Lagrangians has numerical benefits compared to the classical ℓ1-formulation.

Future work will further improve implementation aspects of the developed primal decompositions schemes. Here, sparse backsolves or quasi-Newton Hessian approximations have the potential to greatly reduce the computation time spent in Hessian computation.

A. Sensitivities and Precomputation for Augmented Lagrangians

Observe that for computing ∇ y Φ δ i and ∇ 2 yy Φ δ i in ( 10) and ( 11), the partial derivatives of the implicit function F δ i and L i are required. Next, we derive these quantities for the two relaxed local problems ( 6) and [START_REF] Shin | Graph-Based Modeling and Decomposition of Energy Infrastructures[END_REF].

For [START_REF] Erseghe | Distributed Optimal Power Flow Using ADMM[END_REF], the Lagrangian (omitting arguments) reads

L δ,ρ i . = 1 2   x i y z i   ⊤   H xx i H xy i 0 H xy⊤ i H yy i + ρI -ρI 0 -ρI ρI     x i y z i   +   h x i h y i + λ k i -λ k i   ⊤   x i y z i   -δ1 ⊤ ln(s i ) + γ ⊤ i A x i A y i x ⊤ i z ⊤ i ⊤ -b i + µ ⊤ i ( B x i B y i x ⊤ i z ⊤ i ⊤ + s i -d i ).
Hence, the local KKT conditions read

T δ,ρ i (q ⋆ i , y) . =        H xx i x ⋆ i + H xy i y + h x i + A x⊤ i γ ⋆ i + B x⊤ i µ ⋆ i ρ(z ⋆ i -y) -λ k i + A y⊤ i γ ⋆ i + B y⊤ i µ ⋆ i -(S ⋆ i ) -1 δ1 + µ ⋆ i A x i A y i x ⋆⊤ i z ⋆⊤ i ⊤ -b i B x i B y i x ⋆⊤ i z ⋆⊤ i ⊤ + s ⋆ i -d i        = 0,
where

q ⊤ i . = x ⊤ i , z ⊤ i , s ⊤ i , γ ⊤ i , µ ⊤ i . Moreover, ∇ q i T δ,ρ i (q i , y) =       H xx i 0 0 A x⊤ i B x⊤ i 0 ρI 0 A y⊤ i B y⊤ i 0 0 S -1 i M i 0 I A x i A y i 0 0 0 B x i B y i I 0 0       , (19)
∇ y T δ,ρ i (q i , y) = H xy⊤ i -ρI 0 0 0 ⊤ , (20)
where M i = diag(µ i ). Moreover, by [START_REF] Cao | Decentralized Operation of Interdependent Power Distribution Network and District Heating Network: A Market-Driven Approach[END_REF],

∇ y Φ δ,ρ i (y) = ∇ y L i δ,ρ (q ⋆ i (y); y) = (H yy i + ρI)y + H xy i x ⋆ i + h y i + λ k i -ρz ⋆ i .
Furthermore, by [START_REF] Coulbeck | A hierarchical approach to optimized control of water distribution systems: Part I decomposition[END_REF],

∇ 2 yy Φ δ,ρ i (y)=H yy i + ρI + [H xy⊤ i -ρI 0 0 0] ∇ y q ⋆ i (y), (21)
where q ⋆ i (y) can be computed by the system of linear equations

∇ q i T δ,ρ i (q ⋆ i , y) ∇ y q ⋆ i (y) = -∇ y T δ,ρ i (q ⋆ i , y) (22)
Precomputation for Hessian Evaluation

Next, we show how to precompute matrices to make [START_REF] Demiguel | On Decomposition Methods for a Class of Partially Separable Nonlinear Programs[END_REF]

      H xx i 0 0 A x⊤ i B x⊤ i 0 ρI 0 A y⊤ i B y⊤ i 0 0 S ⋆-1 i M ⋆ i 0 I A x i A y i 0 0 0 B x i B y i I 0 0             ∇ y x ⋆ i ∇ y z ⋆ i ∇ y s ⋆ i ∇ y γ ⋆ i ∇ y µ ⋆ i       =       -H xy i ρI 0 0 0       . ( 23
)
By the third block-row, we have that

∇ y s ⋆ i = -M ⋆-1 i S ⋆ i ∇ y µ ⋆ i . This yields     H xx i 0 A x⊤ i B x⊤ i 0 ρI A y⊤ i B y⊤ i A x i A y i 0 0 B x i B y i 0 -M ⋆-1 i S ⋆ i         ∇ y x ⋆ i ∇ y z ⋆ i ∇ y γ ⋆ i ∇ y µ ⋆ i     =     -H xy i ρI 0 0     .
Since H xx i is invertible, we have

∇ y x ⋆ i ∇ y z ⋆ i = H xx i 0 0 ρI -1 . =P -1 i -H xy i ρI - A x i A y i B x i B y i ⊤ . =K ⊤ i ∇ y γ ⋆ i ∇ y µ ⋆ i . ( 24
)
Employing the Schur complement with respect to the first two block rows yields

K i P -1 i K ⊤ i + 0 0 0 M ⋆-1 i S ⋆ i . =W i ∇ y γ ⋆ i ∇ y µ ⋆ i = -A x i H xx-1 i H xy i + A y i B x i H xx-1 i H xy i + B y i . =R i .
Observe that H xx-1 i , K i and R i can be precomputed. Moreover, the above system of linear equations has significantly less decision variables compared to [START_REF] Demiguel | On Decomposition Methods for a Class of Partially Separable Nonlinear Programs[END_REF] and is in addition positive definite under Assumption 1. This allows to use the Cholesky or Bunch-Kaufmann (LDL) factorization instead of LU.

Precomputation for the Local Solvers

Similar to the above, we can use precomputation for speeding up the solution of the local optimization problems (6) in case interior-point methods are used. Here, we need to compute Newton steps ∇ q i T δ,ρ i (q i , y) ∆q i = -T δ,ρ i (q i , y), i.e.

    H xx i 0 A x⊤ i B x⊤ i 0 ρI A y⊤ i B y⊤ i A x i A y i 0 0 B x i B y i 0 -M -1 i S i         ∆x i ∆z i ∆γ i ∆µ i     = (25) -      H xx i x i + H xy i y + h x i + A x⊤ i γ i + B x⊤ i µ i ρ(z i -y) -λ k i + A y⊤ i γ i + B y⊤ i µ i A x i A y i x ⊤ i z ⊤ i ⊤ -b i B x i B y i x ⊤ i z ⊤ i ⊤ + s i -d i + M -1 i (δ1 -S i µ i )      . =     g 1 g 2 h 1 h 2     ,
where we have eliminate the third block-row via ∆s i = -M -1 i (S i ∆µ i -δ1 + S i µ i ). Solving for the first two block-rows yields

∆x i ∆z i = P -1 i g 1 g 2 -K ⊤ i ∆γ i ∆µ i ,
where P -1 i and K i are from [START_REF] Yoshio | A Nested Schur decomposition approach for multiperiod optimization of chemical processes[END_REF]. Again, employing the Schur-complement with respect to the first two block rows yields

K i P -1 i K ⊤ i + W i ∆γ i ∆µ i = K i P -1 i g 1 g 2 - h 1 h 2 .
Observe that the above is again a small system of linear equations with positive definite coefficient matrix, which allows for using Cholesky factorization and precomputed H xx-1 i and K i .

Sensitivities for the ℓ1 Formulation

The Lagrangian to (8) reads

L i = 1 2 x i z i ⊤ H xx i H xy i H xy⊤ i H yy i x i z i + h x i h y i ⊤ x i z i -δ(1 ⊤ ln(s i ) + 1 ⊤ ln(v i ) + 1 ⊤ ln(w i )) + γ ⊤ i A x i A y i x ⊤ i z ⊤ i ⊤ -b i + µ ⊤ i ( B x i B y i x ⊤ i z ⊤ i ⊤ + s i -d i ) + λ1 ⊤ (v i + w i ) + χ ⊤ i (y -z i -v i + w i ).
Hence, the KKT conditions require

T δ, λ i (u ⋆ i , y) . =              H xx i x ⋆ i + H xy i z ⋆ i + h x i + A x⊤ i γ ⋆ i + B x⊤ i µ ⋆ i H yy i z ⋆ i + H xy⊤ i x ⋆ i + h y i -χ ⋆ i + A y⊤ i γ ⋆ i + B y⊤ i µ ⋆ i -S -1 i δ1 + µ ⋆ i -δV -1 i 1 + ( λ1 -χ ⋆ i ) -δW -1 i 1 + ( λ1 + χ ⋆ i ) A x i A y i x ⋆⊤ i z ⋆⊤ i ⊤ -b i B x i B y i x ⋆⊤ i z ⋆⊤ i ⊤ + s ⋆ i -d i y -z ⋆ i -v ⋆ i + w ⋆ i              ! = 0, where u ⊤ i . = x ⊤ i , z ⊤ i , s ⊤ i , v ⊤ i , w ⊤ i , γ ⊤ i , µ ⊤ i , χ ⊤ i . Thus, ∇ u i T δ, λ i (u i , y) =              H xx i H xy i 0 0 0 A x⊤ i B x⊤ i 0 H xy⊤ i H yy i 0 0 0 A y⊤ i B y⊤ i -I 0 0 S -1 i M i 0 0 0 I 0 0 0 0 V -1 i ( λI -X i ) 0 0 0 -I 0 0 0 0 W -1 i ( λI + X i ) 0 0 I A x i A y i 0 0 0 0 0 0 B x i B y i I 0 0 0 0 0 0 -I 0 -I I 0 0 0              , (26)
where

V i = diag(v i ), W i = diag(w i ), and X i = diag(χ i ). Moreover, ∇ y T δ, λ i (u i , y) = [0 0 0 0 0 0 0 I] ⊤ . (27)
Furthermore, ∇ y Φ δ, λ i (y) = ∇ y L i = χ i , ∇ yy L i = 0, and ∇ yu ⋆ i L i = [0 0 0 0 0 0 0 I]. Thus, by [START_REF] Coulbeck | A hierarchical approach to optimized control of water distribution systems: Part I decomposition[END_REF],

∇ yy Φ δ, λ i (y) = ∇ y χ ⋆ i (y). (28)

B. Proof of Lemma 1

First, we will show that Z

⊤ (ZCZ ⊤ ) -1 Z = Z ⊤ ZC -1 Z ⊤ Z for a regular, symmetric C ∈ R n×n , Z ∈ R m×n with m < n. Consider a re-ordered eigendecomposition C = QΛQ ⊤ and partition Q = [Q 1 Q 2 ], Λ i = blkdiag(Λ 1 , Λ 2 ) such that Q 2 is a nullspace-basis of Z, i.e. ZQ 2 = 0. Hence, we have ZCZ ⊤ = Z[Q 1 Q 2 ] blkdiag(Λ 1 , Λ 2 )[Q 1 Q 2 ] ⊤ Z ⊤ i = ZQ 1 Λ 1 Q 1 Z ⊤ since ZQ 2 = 0. Thus, Z ⊤ (ZCZ ⊤ ) -1 Z =Z ⊤ ZQ 1 Λ -1 1 Q ⊤ 1 Z ⊤ Z. Again, since ZQ 2 = 0, we can expand this expression to Z ⊤ Z[Q 1 Q 2 ] blkdiag (Λ -1 1 , Λ -1 2 )[Q 1 Q 2 ] ⊤ Z ⊤ Z = Z ⊤ ZC -1 Z ⊤ Z.
Proof of a): By (21), we need ∇ y q ⋆ i (y) for computing ∇ yy Φ δ,ρ i , where ∇ y q ⋆ i (y) is defined by [START_REF] Demiguel | On Decomposition Methods for a Class of Partially Separable Nonlinear Programs[END_REF]. Define C

i . = blkdiag(H xx i , ρI, S -1 i M i )), D i . = A x i A y i 0 B x i B y i I
, and E i . =

[H xy⊤ i -ρI 0] ⊤ . Consider [START_REF] Geoffrion | Primal Resource-Directive Approaches for Optimizing Nonlinear Decomposable Systems[END_REF] and parametrize (∇ y x ⋆⊤ i , ∇ y z ⋆⊤ i , ∇ y s ⋆⊤ i ) ⊤ . = Z i P i , where Z i is a nullspace matrix to D i , i.e., the columns of Z i form an orthogonal basis of the nullspace of D i and P i ∈ R (n xi +ny-nr(A x i ))×ny is an auxiliary matrix. Using the above parametrization in [START_REF] Demiguel | On Decomposition Methods for a Class of Partially Separable Nonlinear Programs[END_REF] and multiplying with Z ⊤ i yields Z ⊤ i C i Z i P i = -Z ⊤ i E i by Z ⊤ i D ⊤ i = 0. Since s i , µ i > 0 and Assumption 1 holds, we have C i ≻ 0 and thus Z ⊤ i C i Z i is invertible by full rank of Z i . Hence, by [START_REF] Demiguel | A local convergence analysis of bilevel decomposition algorithms[END_REF] and the above derivation, ∇ yy Φ δ,ρ i (y) =

H yy i +ρI -E ⊤ i Z i (Z ⊤ i C i Z i ) -1 Z ⊤ i E i = H yy i +ρI -E ⊤ i Z i Z ⊤ i C -1 i Z i Z ⊤ i E i . Notice that Z i Z ⊤
i is a diagonal matrix with rank(Z i ) ones and dim(C i ) -rank(Z i ) zeros. Hence, since C i is positive definite, it suffices to show that ∇ yy Φ δ,ρ i (y) ≻ 0 for the worst case, i.e. Z i Z ⊤ i = I (no constraints). In this case we have ∇ yy Φ δ,ρ i (y) = H yy i -H xy⊤ i (H xx i ) -1 H xy i ≻ 0 by the definition of E i , C i , by Assumption 1 a) and the Schur-complement Lemma [36, A.14].

Proof of b): By [START_REF] Bertsekas | Nonlinear Programming[END_REF], we need to show that ∇ y χ ⋆ i (y) ≻ 0, which can be computed by the system of linear equations ( 26), [START_REF] Curtis | A Decomposition Algorithm for Large-Scale Security-Constrained AC Optimal Power Flow[END_REF]. Define F i = blkdiag i G i (∇ y γ ⋆⊤ i , ∇ y µ ⋆⊤ i , ∇ y χ ⋆⊤ i ). Thus, G ⊤ i F -1 i G i (∇ y γ ⋆⊤ i , ∇ y µ ⋆⊤ i , ∇ y χ ⋆⊤ i ) = [0 0 I] ⊤ . Since F -1 i ≻ 0 and by full rank of G i from Assumption 1,

H xx i H xy i H xy⊤ i H yy i , S -1 i M i , V -1 i ( λ- X i ), W -1 i ( λ+X i
G ⊤ i F -1 i G i ≻ 0 and thus (∇ y γ ⋆⊤ i , ∇ y µ ⋆⊤ i , ∇ y χ ⋆⊤ i ) = (G ⊤ i F -1 i G i ) -1 [0 0 I] ⊤ . Since (G ⊤ i F -1 i G i ) -1
≻ 0, all leading principle minors of this matrix must be positive definite by Sylvester's criterion [START_REF] Horn | Matrix Analysis[END_REF]Col 7.1.5]. By variable reordering, the assertion follows.

C. Solution of (1) via ADMM

We derive a distributed ADMM version for (1) as a baseline for numerical comparison. Consider (1), introduce auxiliary variables z i ∈ R ny and consensus constraints y = z 1 = • • • = z S | λ 1 , . . . , λ S . This yields min x,y,z i∈S 1 2

x i z i ⊤ H xx i H xy i H xy⊤ i H yy i x i z i + h x i h y i ⊤ x i z i (29a) subject to A x i A y i x ⊤ i z ⊤ i ⊤ -b i = 0, i ∈ S, (29b)
B x i B y i x ⊤ i z ⊤ i ⊤ -d i ≤ 0, i ∈ S, (29c)
A y y -b y = 0, B y y -d y ≤ 0 (29d)

z i = y, i ∈ S. (29e)
The augmented Lagrangian with respect to y -z i = 0 reads

L ρ = i∈S ϕ x i (x i , z i ) + ϕ y (y) + λ k⊤ i (y -z i ) + ρ 2 ∥y -z i ∥ 2 ,
where ϕ x i are defined by (29a)-(29c) and ϕ y is the indicator function for (29d). Minimizing L ρ w.r.t. (x i , z i ) for fixed (y k , λ k i ) yields for all i ∈ S (x k+1 i , z k+1 i ) = arg min

x i ,z i 1 2 Finally, the Lagrange multiplier update reads

x i z i ⊤ H xx i H xy i H xy⊤ i H yy i + ρI x i z i + h x i h y i -λ k i -ρy k ⊤ x i z i subject to A x i A y i x ⊤ i z ⊤ i ⊤ -b i = 0, B x i B y i x ⊤ i z ⊤ i ⊤ -d i ≤ 0. ( 30
λ k+1 i = λ k i + ρ(y k+1 -z k+1 i ), i ∈ S. (32)
The update rules ( 30)-( 32) define the ADMM iterations. Note that ( 30) and ( 32) can be executed locally for all i ∈ S, whereas (31) defines the global coordination step.

Figure 1 :

 1 Figure 1: Buildings connected via a network with limited capacity.

Figure 2 :

 2 Figure 2: Convergence for |S| = 30 buildings.

Figure 3 :

 3 Figure 3: Convergence for |S| = 180 buildings.

Figure 4 :

 4 Figure 4: Convergence for |S| = 300 buildings.

Figure

  Figure2illustrates the numerical performance of both primal decomposition variants and ADMM. Figures3-4show the AL formulation only, since the ℓ1 formulation runs out of memory for these cases. The constraint violations for the equality constraints (1b),[A x A y ][x ⊤ y ⊤ ] ⊤ -b ∞ , for the inequality constraints (1c), max([B x B y ][x ⊤ y ⊤ ] ⊤b), and the value of the cost function f ([x ⊤ , y ⊤ ] ⊤ ) from (1a) are displayed, where the xaxis shows the iteration count. One can observe that the primal decomposition schemes achieve a high degree of feasibility in less than 10 iterations for all cases. Moreover, the optimality gap (f([x ⊤ , y ⊤ ] ⊤ )-f ([x ⋆⊤ , y ⋆⊤ ] ⊤ ))/f ([x ⋆⊤ , y ⋆⊤ ] ⊤) is below 0.01% in less than 10 iterations for both primal decomposition variants and for all |S|, where [x ⋆⊤ , y ⋆⊤ ] ⊤ is computed via Ipopt. For ADMM, infeasibility and the optimality gap stay large independently of the choice of the tuning parameter ρ.

  0 -I I   . By Assumption 1, s i , v i , w i , µ i > 0, and λ > max j |[χ i ] j |, we have that F i ≻ 0. Hence, ∇ y x ⋆⊤ i , ∇ y z ⋆⊤ i , ∇ y s ⋆⊤ i , ∇ y v ⋆⊤ i , ∇ y w ⋆⊤ i ⊤ = -F -1

  )Minimising L ρ w.r.t. y for fixed (x k+1 i to A y y -b y = 0, B y y -d y ≤ 0.

  A y (y +∆y) -b y =0, B y (y +∆y) -b y ≤0.

		2 yy ϕ i ) locally for all i ∈ S;
	(ii) solve the coordination problem
	min ∆y	i∈S	1 2	∆y∇ 2 yy ϕ i (y)∆y ⊤ + ∇ y ϕ i (y) ⊤ ∆y	(13)
	s.t. α = 1
	while ψ(y) -ψ(y + α∆y) ≥-σα∇ y ψ(y) ⊤ ∆y do
	α ← ζα
	end
	y ← y + α∆y
	end
	Return y, {x

i } i∈S .

hold at the optimal solution at p ⋆

Table 1 :

 1 #decision variables and #constraints for |S| ∈ {30, 180, 300} buildings.

	|S|	n x	n y	n e	n i n ey	n iy
	30	28,200	690	15,090	28,800	0	1,403
	180 169,200 4,140	90,540 172,800	0	8,303
	300 282,000 6,900 150,900 288,000	0 13,823

Table 2 :

 2 Timing and #iter for |S| ∈ {30, 180, 300} buildings, 30 cores.

	Ipopt	ADMM	ADMM

Table 3 :

 3 Internal timing (%) for the AL formulation, 30 cores. |S| sensitivity eval. local sol. coord.

	line search other

Under standard regularity assumptions [22, A1-C1].

Observe that we have replaced y by zi in the constraints here but not in the objective. Exchanging y in the objective is possible but might lead to a different numerical behavior.

We use the definition of global convergence in the context of NLPs, i.e. the convergence to a KKT point from an arbitrary initialization[START_REF] Nocedal | Numerical Optimization[END_REF] Chap. 3]. However, since problems (1) and (2) are strongly convex and we assume regularity in Assumption 1, every KKT point is also a global minimizer.

The ADMM-based QP solver OSQP solver did not converge for the problems presented here.