Advances in Energy System Optimization

This paper investigates the role of feasible initial guesses and large consensus-violation penalization in distributed optimization for Optimal Power Flow (OPF) problems. Specifically, we discuss the behavior of the Alternating Direction of Multipliers Method (ADMM). We show that in case of large consensus-violation penalization ADMM might exhibit slow progress. We support this observation by an analysis of the algorithmic properties of ADMM. Furthermore, we illustrate our findings considering the IEEE 57 bus system and we draw upon a comparison of ADMM and the Augmented Lagrangian Alternating Direction Inexact Newton (ALADIN) method.

optimization is considered to be helpful as it allows splitting large OPF problems into several smaller subproblems; thus reducing complexity and avoiding the exchange of full grid models between subsystems. We refer to [2] for a recent overview of distributed optimization and control approaches in power systems.
One frequently discussed convex distributed optimization method is the Alternating Direction of Multipliers Method (ADMM) [3], which is also applied in context of AC OPF [1,4,5]. ADMM often yields promising results even for large-scale power systems [4]. However, ADMM sometimes requires a specific partitioning technique and/or a feasible initialization in combination with high consensusviolation penalization parameters to converge [1]. The requirement of feasible initialization seems quite limiting as it requires solving a centralized inequalityconstrained power flow problem requiring full topology and load information leading to approximately the same complexity as full OPF. In previous works [6][7][8] we suggested applying the Augmented Lagrangian Alternating Direction Inexact Newton (ALADIN) method to stochastic and deterministic OPF problems ranging from 5 to 300 buses. In case a certain line-search is applied [9], ALADIN provides global convergence guarantees to local minimizers for non-convex problems without the need of feasible initialization. The results in [6] underpin that ALADIN is able to outperform ADMM in many cases. This comes at cost of a higher per-step information exchange compared with ADMM and a more complicated coordination step, cf. [6].
In this paper we investigate the interplay of feasible initialization with high penalization for consensus violation in ADMM for distributed AC OPF. We illustrate our findings on the IEEE 57-bus system. Furthermore, we compare the convergence behavior of ADMM to ALADIN not suffering from the practical need for feasible initialization [9]. Finally, we provide theoretical results supporting our numerical observations for the convergence behavior of ADMM. The paper is organized as follows: In Sect. 2 we briefly recap ADMM and ALADIN including their convergence properties. Section 3 shows numerical results for the IEEE 57-bus system focusing on the influence of the penalization parameter ρ on the convergence behavior of ADMM. Section 4 presents an analysis of the interplay between high penalization and a feasible initialization.

ALADIN and ADMM
For distributed optimization, OPF problems are often formulated in affinely coupled separable form min x∈R nx i∈R where the decision vector is divided into sub-vectors x = [x 1 , . . . , x |R| ] ∈ R n x , R is the index set of subsystems usually representing geographical areas of a power system and local nonlinear constraint sets X i := {x i ∈ R n xi | h i (x i ) ≤ 0}.

Parallelizable
Step: Solve for all i ∈ R locally 2. Update dual variables: 3. Throughout this work we assume that f i and h i are twice continuously differentiable and that all X i are compact. Note that the objective functions f i : R n xi → R and nonlinear inequality constraints h i : R n xi → R n hi only depend on x i and that coupling between them takes place in the affine consensus constraint i∈R A i x i = 0 only. There are several ways of formulating OPF problems in form of (1) differing in the coupling variables and the type of the power flow equations (polar or rectangular), cf. [4,6,10].
Here, we are interested in solving Problem (1) via ADMM and ALADIN summarized in Algorithms 1 and 2 respectively. 1,2 At first glance it is apparent that ADMM and ALADIN share several features. For example, in Step (1) of both algorithms, local augmented Lagrangians subject to local nonlinear inequality constraints h i are minimized in parallel. 3 Observe that while ADMM maintains multiple local Lagrange multipliers λ i , ALADIN considers one global Lagrange multiplier vector λ. In Step (2), ALADIN computes sensitivities B i , g i and C i (which often can directly be obtained from the local numerical solver without additional computation) whereas ADMM updates the multiplier vectors λ i . 1 We remark that there exist a variety of variants of ADMM, cf. [3,11]. Here, we choose the formulation from [9] in order to highlight similarities between ADMM and ALADIN. 2 Note that, due to space limitations, we describe the full-step variant of ALADIN here. To obtain convergence guarantees from a remote starting point, a globalization strategy is necessary, cf. [9]. 3 For notational simplicity, we only consider nonlinear inequality constraints here. Nonlinear equality constraints g i can be incorporated via a reformulation in terms of two inequality constraints, i.e. 0 ≤ g i (x i ) ≤ 0.

Parallelizable
Step: Solve for all i ∈ R locally 2. Compute sensitivities: Compute Hessian approximations B k i , gradients g k i and Jacobians of the active constraints C k i , cf. [9]. 3. Consensus Step: Solve the coordination problem min Δx,s i∈R 4. Update variables: z k+1 ← x k + Δx k , λ k+1 ← λ QP . Similarity to ADMM: Remove C k i in (5), set B k i = ρA i A i , g k i = A i λ k and Σ i = A i A i , set μ = ∞ (i.e. s = 0) and use dual ascent step (3) for updating λ k in (4), cf. [9].

In
Step (3), both algorithms communicate certain information to a central entity which then solves a (usually centralized) coordination quadratic program. However, ALADIN and ADMM differ in the amount of exchanged information: Whereas ADMM only communicates the local primal and dual variables x i and λ i , ALADIN additionally communicates sensitivities. This is a considerable amount of extra information compared with ADMM. However, there exist methods to reduce the amount of exchanged information and bounds on the information exchange are given in [6]. Another important difference is the computational complexity of the coordination step. In many cases, the coordination step in ADMM can be reduced to an averaging step based on neighborhood communication only [3], whereas in ALADIN the coordination step involves the centralized solution of an equality constrained quadratic program.
In the last step, ADMM updates the primal variables z i , while ALADIN additionally updates the dual variables λ. Differences of ALADIN and ADMM also show up in the convergence speed and their theoretical convergence guarantees: Whereas ALADIN guarantees global convergence and quadratic local convergence for non-convex problems if a certain globalization strategy is applied [9], few results exist for ADMM in the non-convex setting. Recent works [12,13] investigate the convergence of ADMM for special classes of non-convex problems; however, to the best of our knowledge OPF problems do not belong to these classes.

Numerical Results
Next, we investigate the behavior of ADMM for large ρ and a feasible initialization to illustrate performance differences between ADMM and ALADIN. We consider power flow equations in polar form with coupling in active/reactive power and voltage angle and magnitude at the boundary between two neighbored regions [6]. We consider the IEEE 57-bus system with data from MATPOWER and partitioning as in [6] as numerical test case. Figures 1 and 2 show the convergence behavior of ADMM (with and without feasible initialization (f.)) and ALADIN for several convergence indicators and two different penalty parameters ρ = 10 4 and ρ = 10 6 . 4 Therein, the left-handed plot depicts the consensus gap Ax k ∞ representing the maximum mismatch of coupling variables (active/reactive power and voltage magnitude/angle) at borders between two neighbored regions. The second plot shows the objective function value f k and the third plot presents the distance to the minimizer x k − x ∞ over the iteration index k. The right-handed figure shows the nonlinear constraint violation g(z k ) ∞ after the consensus steps (4) and (5) of Algorithms 1 and 2 respectively representing the maximum violation of the power flow equations. 5 In case of small ρ = 10 4 , both ADMM variants behave similar and converge slowly towards the optimal solution with slow decrease in consensus violation, nonlinear constraint violation and objective function value. If we increase ρ to ρ = 10 6 with results shown in Fig. 2, the consensus violation Ax k ∞ gets smaller in ADMM with feasible initialization. The reason is that a large ρ forces x k i being close to z k i leading to small Ax k ∞ as we have Az k = 0 from the consensus step. But, at the same time, this also leads to a slower progress in optimality f k compared to ρ = 10 4 , cf. the second plot in Figs. 1 and 2.
On the other hand, this statement does not hold for ADMM with infeasible initialization (blue lines in Figs. 1 and 2) as the constraints in the local step  Fig. 3 Convergence behavior of ADMM with infeasible initialization, ADMM with feasible initialization (f.) for ρ = 10 12 and ALADIN and the consensus step of ADMM enforce an alternating projection between the consensus constraint and the local nonlinear constraints. The progress in the nonlinear constraint violation g(z k ) ∞ supports this statement. In its extreme, this behavior can be observed when using ρ = 10 12 depicted in Fig. 3. There, ADMM with feasible initialization produces very small consensus violations and nonlinear constraint violations at cost of almost no progress in terms of optimality.
Here the crucial observation is that in case of feasible initialization and large penalization parameter ρ ADMM produces almost feasible iterates at cost of slow progress in the objective function values. From this, one is tempted to conclude that also for infeasible initializations ADMM will likely converge, cf. Figs. 1, 2, and 3. This conclusion is supported by the rather small 57-bus test system. However, it deserves to be noted that this conclusion is in general not valid, cf. [1].
For ALADIN, we use ρ = 10 6 and μ = 10 7 . Comparing the results of ADMM with ALADIN, ALADIN shows superior quadratic convergence also in case of infeasible initialization. This is inline with the known convergence properties of ALADIN [9]. However, the fast convergence comes at the cost of increased communication overhead per step, cf. [6] for a more detailed discussion. Note that ALADIN involves a more complex coordination step, which is not straightforward to solve via neighborhood communication. Furthermore, tuning of ρ and μ can be difficult.
In power systems, usually feasibility is preferred over optimality to ensure a stable system operation. Hence, ADMM with feasible initialization and large ρ can in principle be used for OPF as in this case violation of power flow equations and generator bounds are expected to be small and one could at least expect certain progress towards an optimal solution. Following this reasoning several papers consider A(x k − z k ) ∞ < as termination criterion [1,4]. However, as shown in the example above, if ρ is large enough, and ADMM is initialized at a feasible point, this termination criterion can always be satisfied in just one iteration if ρ is chosen sufficiently large. Consequently, it is unclear how to ensure a certain degree of optimality. An additional question with respect to ADMM is how to obtain a feasible initialization. To compute such an initial guess, one has to solve a constrained nonlinear least squares problem solving the power flow equations subject to box constraints. This is itself a problem of almost the same complexity as the full OPF problem. Hence one would again require a computationally powerful central entity with full topology and parameter information. Arguably this jeopardizes the initial motivation for using distributed optimization methods.

Analysis of ADMM with Feasible Initialization for ρ → ∞
The results above indicate that large penalization parameters ρ in combination with feasible initialization might lead to pre-mature convergence to a suboptimal solution. Arguably, this behavior of ADMM might be straight-forward to see from an optimization perspective. However, to the best of our knowledge a mathematically rigorous analysis, which is very relevant for OPF, is not available in the literature.

Proposition 1 (Feasibility and ρ → ∞ imply x k+1
i − x k i ∈ null(A i )) Consider the application of ADMM (Algorithm 1) to Problem (1). Suppose that, for all k ∈ N, the local problems (2) have unique regular minimizers x k i . 6 Fork ∈ N, let λk i be bounded and, for all i ∈ R, zk i ∈ X i . Then, the ADMM iterates satisfy lim ρ→∞ x k i (ρ) − xk i ∈ null(A i ), ∀k >k.
Proof The proof is divided into four steps. Steps 1-3 establish technical properties used to derive the above assertion in Step 4.
Step 1. At iterationk the local steps of ADMM are xk i (ρ) = argmin Now, by assumption all f i s are twice continuously differentiable (hence bounded on X i ), λk i is bounded and all zk i ∈ X i . Thus, for all i ∈ R, lim ρ→∞ xk i (ρ) = zk i + vk i with vk i ∈ null(A i ).
Step 2. The first-order stationarity condition of (6) can be written as where γk i is the multiplier associated to h i . Multiplying the multiplier update formula (3) with A i from the left we obtain A i λ k+1 Combined with (7)  Step 3. Next, we show by contradiction that Δxk i ∈ null(A i ) for all i ∈ R and ρ → ∞. Recall the coordination step (4b) in ADMM given by min Δx i∈R Observe that any Δxk i ∈ null(A i ) is a feasible point to (8) as i∈R A i xk i = 0.
Consider a feasible candidate solution Δx i / ∈ null(A i ) for which i∈R A i (xk i + Δx i ) = 0. Clearly, λk +1 i A i Δx i (ρ) will be bounded. Hence for a sufficiently large value of ρ, the objective of (8) will be positive. However, for any Δx i ∈ null(A i ) the objective of (8) is zero, which contradicts optimality of the candidate solution Δx i / ∈ null(A i ). Hence, choosing ρ sufficiently large ensures that any minimizer of (8) lies in null(A i ).
Corollary 1 (Deterministic code, feasibility, ρ → ∞ implies x k+1 i = x k i ) Assuming that the local subproblems in ADMM are solved deterministically; i.e. same problem data yields the same solution. Then under the conditions of Proposition 1 and for ρ → ∞, once ADMM generates a feasible point xk i to Problem (1), or whenever it is initialized at a feasible point, it will stay at this point for all subsequent k >k.
The above corollary explains the behavior of ADMM for large ρ in combination with feasible initialization often used in power systems [1,4]. Despite feasible iterates are desirable from a power systems point of view, the findings above imply that high values of ρ limit progress in terms of minimizing the objective. ALADIN for ρ → ∞) Note that for ρ → ∞, ALADIN behaves different than ADMM. While the local problems in ALADIN behave similar to ADMM, the coordination step in ALADIN is equivalent to a sequential quadratic programming step. This helps avoiding premature convergence and it ensures decrease of f in the coordination step [9].

Conclusions
This method-oriented work investigated the interplay of penalization of consensus violation and feasible initialization in ADMM. We found that-despite often working reasonably with a good choice of ρ and infeasible initialization-in case of feasible initialization combined with large values of ρ ADMM typically stays feasible yet it may stall at a suboptimal solution. We provided analytical results supporting this observation. However, computing a feasible initialization is itself a problem of almost the same complexity as the full OPF problem; in some sense partially jeopardizing the advantages of distributed optimization methods. Thus distributed methods providing rigorous convergence guarantees while allowing for infeasible initialization are of interest. One such alternative method is ALADIN [9] exhibiting convergence properties at cost of an enlarged communication overhead and a more complex coordination step [6]. 12 The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Introduction
There is a world-wide trend towards increasing the share of renewable energy sources for electricity generation. Since the yield of renewable energy sources like e.g. solar and wind power is time-varying and depends on geographic location there is an increasing demand for large-scale power transfer over long distances. In recent years there have been more and more high-voltage direct current (HVDC) links for power transmission commissioned and more are yet to be planned and installed. In [1] an extensive overview of the advantages of DC technology for power transmission is given. For instance, in Germany there are several embedded HVDC M. Giuntoli · S. Schmitt ( ) ABB Corporate Research Center Germany, Ladenburg, Germany e-mail: marco.giuntoli@de.abb.com; susanne.schmitt@de.abb.com links (i.e. they connect synchronous AC nodes) planned in order to transfer wind power from the north to the high load regions in the south, see [2].
Optimal power flow is a widely discussed optimization problem to determine the optimal operation of controllable equipment in a power grid, see e.g. [3,4]. For transmission grid operation the so-called N − 1 criterion as defined in [5] is a requirement. It means that for any single outage of a grid component (e.g. a generator, a line or a power transformer) the remaining grid must be able to operate safely within its limits. This leads to security-constrained optimal power flow problems (SC-OPF). For this class of problems an extensive overview is provided in [6]. SC-OPF problems can be formulated either in the preventive way, i.e. there is a single set of control variables that has to be feasible in all considered contingency cases or in the corrective way, i.e. for each contingency case there may be a different set of control variables, and their values can be reached by short-term control actions when a contingency case is taking place, see [7].
Optimal power flow formulations for hybrid AC/DC grids including security constraints are presented e.g. in [1,8] and [9]. In [10] a corrective SC-OPF methodology is employed to investigate the capabilities of embedded HVDC systems to compensate for contingency cases in the AC system and thus to reduce redispatch by generators.
In this work we use SC-OPF to study the effects of the N − 1 criterion with emphasis on contingency cases within the embedded HVDC systems. To our knowledge the effects of N − 1 cases inside the DC systems have not yet been studied extensively in the literature. As in [10] we focus on the capabilities of the grid components and in particular the embedded HVDC systems to compensate for contingency cases and to avoid generator redispatch.
This work is structured as follows: in Sect. 2 we present the mathematical model for SC-OPF for hybrid AC/DC grids. Then, in Sect. 3 we describe the case study grid and scenario used as well as the implementation of the approach (Sect. 4). The case study results are presented and discussed in Sect. 5, and we give conclusions and suggestions for future work in the final Sect. 6.

The Mathematical Model
The mathematical model of the SC-OPF is used to perform a steady-state optimization for hybrid AC/DC grids, where the two types of grids are coupled via a set of power converters. The model is formulated such that just one single optimization problem can be used to solve the AC and DC grids simultaneously; this means in particular that no iterative algorithm is needed to reach the final solution.
The model is written to address one-phase electric systems, i.e. it can evaluate a mono-phase system or a multiphase system but in the latter case the system must be balanced and symmetric to consider only the positive sequence (Fortescue transformation).
All the variables and equations are expressed in complex polar form and they are differentiable and continuous. This means that discrete variables (e.g. limited number of tap changer positions or discontinuous active power range for the generators) are not considered. The set of variables are split into two parts: control (or independent) variables u representing power setpoints (active and reactive) for controllable devices as well as tap changer and phase changer positions for power transformers. State (or dependent) set of variables x consisting of voltage magnitude and phase at AC nodes, voltage at DC nodes as well as power flows over AC and DC branches.
Regarding the N − 1 criterium, the model is able to take into account the outage of a branch (overhead line, power transformer or cable) or a generic equipment inside the system (power generator, power converter). Both the preventive and the corrective N − 1 methodology can be used. Furthermore, there is no particular restriction about the grid topology. This means that it is possible to solve meshed AC and DC (e.g. multiterminal DC) synchronous and asynchronous grids.
The mathematical model can be formulated as the following nonlinear optimization problem: where N c is the number of contingency cases considered (c = 0 corresponds to the base case where no contingencies have happened). The set of equality constraints g c (u c , x c ) describes the power balance for each node and equipment, while the set of inequality constraints h c (u c , x c ) describes the physical limits of each component. The deviation of the control variable in case of corrective N − 1 method can be accordingly limited with the upper Δu c and lowerΔu c bounds (e.g. active power setpoint deviation between the contingency cases and the base case).
The objective function f may be chosen as the generator dispatch cost in the system, system losses, voltage magnitude or angle variability, reactive power generated or another power system objective. It is also possible to provide flexible weights for the different contingency cases under consideration. In this work the objective function is defined as the economic dispatch cost in quadratic form: Here N g is the number of generators in the grid, P g is the active power generated by the generator g, while α g , β g , γ g are cost coefficients for the individual generators. Δp + g,c and Δp − g,c are the active power deviation (positive and negative) for each contingency and they are calculated with respect to the base case and weighed with the relative linear cost coefficients c + g and c − g . The grid model is composed of the following items: nodes (AC and DC), power loads, branches and active equipments (generators, converters and reactive power sources like e.g. static VAR compensators (SVC)). Two power balance equations are given for the AC nodes and one for the DC nodes. In the subsequent paragraphs an overview over the models of the grid components is provided.
The following equation describes the active power balance for an AC node: where the subscripts n and c denote the n-th node and the c-th contingency case. The term p(v n,c , θ n,c ) is the active power flow into the n-th AC node where v and θ are the voltage magnitude and phase; the term p g,c is the power supplied by the g-th active source; the term p h,ac,c is the power supplied by the h-th converter; the term p d is the power demand. In the same way it is possible to describe the reactive power balance for the AC nodes: For the DC nodes, only one power balance equation is needed: since in this case no power generators and loads are present. For both AC and DC nodes the voltage magnitude is constrained to its lower and upper limits. The branches are modeled as generic passive components able to transfer current from one node to another according to Kirchhoff's circuit laws. The double bi-pole model with PI scheme has been used. This model can be applied to overhead lines, cables, power transformers (with and without tap changer and/or phase shifter) and to DC overhead lines and cables. Without loss of generality we assume the tap changer and phase shifter to be installed on the from-node of the branch and constrained with their limits. Both can be modelled with a pure gain δ and a complex rotation e jθ . For a generic branch b, contingency c and defining i as from-node and j as to-node, it is possible to describe the relationship between the current and voltage phasors as follows: Exploiting the symmetry of the model, the matrix component with subscript ii is equal to the component with jj and, in the same way, the component with ij is equal to the component with j i. Those values arise from the longitudinal and transversal susceptance and conductance parameters. The equations that describe a double bi-pole model are used to model the active and reactive power flow through the branches, and thus those flows are used inside the nodes power balance equations (3), (4), (5). Obviously, the power flow (active power or apparent power or current module) can be constrained at both sides of the branch with an upper limit. The power generators are generic active components able to inject active and/or reactive power into the system. The generator capability is defined as a domain created by the intersection of a rectangle (active and reactive upper and lower bounds) and a circle (apparent power bound). Additional linear bounds between active and reactive power can be defined in order to customize the capability domain.
In a similar way, by neglecting the active power, voltage controlling reactive power sources can be used inside the model: in this case, a linear relationship between voltage magnitude and reactive power output is used.
The power converters are able to transfer power in both directions between AC and DC nodes. The converter model takes into account the power balance between the AC and DC ports with the losses due to its power electronics. Moreover, the converters may have capability domains described in the same way as for generators. The converter power balance equation is written as follows: where the subscripts loss denotes the power losses inside the h-th converter and the c-th contingency. The AC side current is determined by the following equation that describes the power balance in the AC converter side: Finally, the power losses are calculated by a formula that creates a quadratic relationship between the current i h,ac,c and losses p h,loss,c (see e.g. [11] and [12]) with parameters α h , β h and γ h that may be individual for each converter: This model is suited for steady-state analyses, so that transient behaviour, e.g. fault and stability analysis cannot be investigated with it.

Case Study
The case study in this work is based on the IEEE118 test grid (see [13]).
The generator costs have been modified in order to encourage power transfer from the (low cost) western zone to the (high cost) eastern zone (see Fig. 1). In particular, the generation cost in the western zone have been set to zero to emulate renewable generation, while the cost coefficients in the western zone have not been modified.
The focus of our analysis lies on one additional power corridor between the nodes 8 and 65 (highlighted in red in Fig. 1). There are three options for this considered corridor: an AC corridor, a monopolar HVDC corridor as well as a bipolar HVDC corridor with identical total power ratings of 500 MVA and total length of 500 km. The electrical parameters of the AC corridor as well as the resistance of the DC lines are taken from [14]. Due to the huge length of the corridor the AC line has been split in three line segments. The HVDC converters exhibit a loss model with a quadratic relationship to the AC current (see [15]).
For all three configurations the effects of security constraints on the corridor are analysed. For both the AC corridor and the monopolar HVDC corridor the considered contingency implies the outage of the whole corridor, while for the bipolar configuration an N − 1 contingency case involved the outage of one of the two poles.
In order to discourage redispatch of generators in the contingency case, we require the generator active and reactive power output to be identical in the base case and the contingency case, while HVDC converter setpoints may differ between base case and contingency case. This means that the generators follow a preventive security strategy while HVDC converters follow a corrective security strategy. The same behaviour can be achieved by imposing very high redispatch cost on the generators. The possibly changing losses after an outage are balanced by adapted voltage profiles. In addition, the setpoints for transformers are fixed.

Implementation
A software tool for solving the security-constrained optimal power flow problem as described in Sect. 2 has been developed in C++ using Microsoft Visual Studio TM . Inside this software the open-source library IPOPT (see [16]) is used as solver for non-linear programming based on an interior point method. To speed up the optimization process and to increase accuracy, we derived the analytic forms of Jacobian and Hessian matrices of the objective function and the constraints and pass them directly to the solver instead of using numerical differentiation.
The input data as well as the optimization results including all control and state variables are given as text files. In the results discussion below (Sect. 5) we summarize these results in the form of the key performance indicators dispatch Illustration if the IEEE118 test grid, the two generator cost zones and the additional corridor under investigation [17] cost, power flow over the considered corridor, active power losses and voltage angle deviation.

Results
The result of the SC-OPF is analysed in terms of dispatch cost (Fig. 2), power flow over the considered corridor (Fig. 3), active power losses in the system (Fig. 4) and voltage angle difference along the corridor (Fig. 5).
The first observation is that in the cases without the N − 1 criterion, all the values for the monopolar and the bipolar HVDC corridors coincide, which is obvious since they have identical parameters and rating. When looking at the dispatch cost in Fig. 2 it can be observed that the cost is lowest for both the HVDC variants without N − 1 criterion. With N − 1 criterion both HVDC variants have lower dispatch cost than the AC variant, while the bipolar configuration has the lowest dispatch cost. This can be explained by the higher efficiency of the HVDC systems compared to the AC corridor (lower losses over long distances) and the higher redundancy of the bipolar configuration compared to the monopolar one. This way more power can be transferred from the low cost generation zone to the high cost generation zone which enables a generation dispatch at lower total cost. This explanation is confirmed by the power flow over the corridor (Fig. 3), where the bipolar HVDC configuration is able to transfer the most power from the low cost generation zone to the high cost generation zone. The monopolar HVDC corridor is transferring significantly less power in the N − 1 case since it has to prepare for the outage of the full corridor. The AC corridor is not even able to maintain the power flow direction in the N −1 case. The significant deviations in the power flows between the AC and monopolar HVDC configuration can be explained by voltage limits and differences in the generator dispatch.
The active power losses in the system can be seen in Fig. 4. There we did not observe a clear pattern, which can be explained by the fact that in all the scenarios the generator dispatch and thus the power flows through the whole system are very different, and the losses were not the objective of the optimization performed.
When looking at the angle deviations in Fig. 5 is can be observed that the AC corridor exhibits a large angle deviation in the case without security criterion, while all the HVDC configurations are able to maintain smaller angle deviations which indicates a higher system stability.
In Table 1 the voltage magnitudes and reactive power injected at the corridor terminal nodes 8 and 65 are listed for the N − 1 scenario before and after the contingencies. It can be observed that the bipolar HVDC configuration is able to

Conclusion and Future Work
A mathematical model and a software for security-constrained optimal power flow for generic hybrid AC/DC grids has been developed. In this work it has been used to study the effects of the N − 1 criterion on embedded HVDC systems. The above results show that for the bipolar HVDC configuration the power plant dispatch with lowest cost could be established when the N −1 criterion is employed. This study has been performed on the IEEE118 grid with a single corridor under consideration. The results indicate the benefits of employing bipolar HVDC systems in order to increase the social welfare.
In order to draw conclusions for real transmission systems, additional studies on real grid configurations with possibly several embedded HVDC links would have to be performed.
Additionally, an investigation including contingency cases both in the AC and DC parts of the power system would be interesting to investigate the capabilities of embedded HVDC to compensate for general N − 1 cases and thus to reduce redispatch cost for grid operators.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Introduction
The increasing penetration of Renewable Energy Sources (RES) leads to a power system operation closer to its operational limits. Enhanced methods and tools will be crucial for a secure, but also efficient and cheap planning and operation of the power grid. A Transmission System Operator (TSO) will have larger assets of controllable devices at hand, such as Voltage Source Converter (VSC-) based HVDC-systems and energy storage systems. Combined with the need to assure N-1 security, this requires a Security-Constrained Dynamic Optimal Power Flow (SC-D-OPF), which incorporates both N-1 security and multiple time steps into the optimization. First research in that area was done in [1], where not only power but also necessary energy reserves can be determined to ensure N-1 security. In [2], successive linear algorithms and approximated power flows are used to solve large-scale SC-D-OPF problems in a European context. The authors of [3] propose an SC-D-OPF model including uncertainty of wind power or equipment availability, which is solved in a two-stage stochastic program. In [4], SC-D-OPF is used as inner iteration in a hierarchical approach to optimize a central deployment signal which is sent to the units able to provide a reserve. Furthermore, [5][6][7] include energy storage systems in an SC-D-OPF calculation. A further extension is shown in [8], where SC-D-OPF is solved in a hybrid AC-DC grid with energy storage.
Due to the increasing complexity of the power system and an operation closer to network limitations, central coordination in large scale networks comes with major computational burdens. Furthermore, privacy can be a concern for each transmission system operator (TSO) controlling a certain region. Subsequently, the interest in distributed optimization, also referred to as multi-area optimization has substantially grown in recent years. An early overview of distributed OPF algorithms can be found in [9] and the most recent developments are examined in detail in [10]. The Alternating Direction of Multipliers Method (ADMM) [11], has become a popular branch to tackle the non-convex AC OPF problem [12][13][14]. Each area uses variable duplicates from neighboring areas and is solved to optimality. Penalty terms, which are calculated from the coupling variable deviation and Lagrangian multipliers, are added to the objective function to push two neighbors towards a common boundary solution. ADMM was applied to a hybrid AC-DC grid in [15]. Multiarea optimization is applied to D-OPF with storage in [16] using OCD, but to the best of our knowledge, no attempts have been made toward distributing SC-D-OPF.

Security-Constrained Dynamic Optimal Power Flow (SC-D-OPF)
The optimization horizon is defined by a set of time steps T = {1, . . . , T } with T the total number of considered time steps. Furthermore, we define a scenario set C = {0, . . . , C}. Here, scenario c = 0 defines the base case, i.e. the original fault-free network. A total of C possible contingencies is added by including modified versions of the base case to the scenario set. A modification could be the outage of an AC line, a DC line or a converter. With x t,c the optimization variables of time step t and scenario c, the full set of optimization variables is collected in The full problem has the form (2): It minimizes the sum of weighted costs f over all scenarios. Factor p t,c defines a probability of each contingency to enable a risk-based OPF. Constraints h base describe the basic constraints of a conventional OPF, which must be fulfilled during each scenario (t, c). That includes AC and DC node power balance; thermal AC and DC branch limits; AC and DC voltage limits; quadratic loss function of AC-DC converters; and operational limits of generators, converters or sheddable loads. Time-dependent constraints, such as storage energy limits or generator ramping limits, are collected in h dyn . Those constraints only apply to the base case (c = 0), since set points after an outage are coupled solely to the respective base case of the identical time step. This N-1 coupling is modeled by the constraints h N-1 and guarantees N-1 security. Here, the coupling between each contingency and the base case is explicitly modeled. That is, set point deviation limits from before to after the occurrence of an outage can be defined. Note that for the sake of readability, we omit equality constraints which can be interpreted as reformulated inequalities as well. We refer to [8] for more details on the modeling.

Distributed Formulation
For the sake of readability, we write (2) in the form (3): Let R non-overlapping regions be defined in R = {1, . . . , R}, that is, each node belongs to exactly one region. An equivalent problem as (3), but in separable form, can be written as (4): Here, x k ∈ R l k , h k : R l k → R m k and f k : R l k → R are optimization variables, non-linear constraints and objective function, respectively, in a region k. Matrix A k ∈ R n×l k maps x k onto the full set of n coupling constraints and enforces consensus between regions. Thus, (4b) forms independent local SC-D-OPF constraints and (4c) form the coupling constraints which recuperate a feasible overall power flow.

Network Decomposition in AC-DC Grids
To allow for a separable formulation, the network is first partitioned and then decomposed. Network partitioning can be crucial for distributed algorithms to achieve good performance. We believe that network partitions are inherently given by structural responsibilities. For example, a region or control area could represent one TSO, multiple TSOs or a whole country. Thus, the number and dimension of AC regions are historically known. However, responsibilities in overlaying DC networks are yet to be defined and various options are thinkable [15]. We choose a Shared-DC approach. Here, we define R = R AC hybrid AC-DC regions, where each AC region may also contain DC nodes. Thus, each existing TSO gains control over converters in its own control area. Let N k identify all nodes in region k. Herewith, N AC k collects all AC nodes and N DC k all DC nodes in region k.

Decoupling of AC or DC Tie Line
In principle, the decoupling of AC and DC tie line is identical, except for extended consensus constraints due to complex voltage and power in the AC case (Fig. 1).
The original line is cut into two halves; auxiliary nodes (m, n) and auxiliary generators (a, b) are added at both open ends. Thus, node sets are augmented to To guarantee a feasible power flow, the voltage Region A Region B Fig. 1 Decoupling model of a tie line (AC or DC) between nodes i and j . Two auxiliary nodes (m and n) and two auxiliary generators (a and b) are added at the middle. Reactive power source Q G is only added for an AC tie line must be equal at nodes m and n. Furthermore, the generators must produce the same amount of power of the opposite sign. In the case of an AC tie line (i ∈ N AC A , j ∈ N AC B ), this leads to boundary conditions including complex voltage and both active and reactive power. For time step and scenario (t, c), we have: In the case of a DC tie line (i ∈ N DC A , j ∈ N DC B ), only real voltage and active power must meet the constraints. For time step and scenario (t, c), we have: Equations (5) and (6) are the only boundary constraints.

Building Consensus Constraint Matrix A
Since x t,c k uses the augmented node sets of region k, the auxiliary generator variables are included. Thus, consensus constraints (5)-(6) between region A and B can be written toÃ t ∈ T c ∈ C In a more compact form, (8) can be written as

Implemented ADMM Algorithm
The general idea of ADMM, see also [11], is the following. Augmented regional OPFs are solved and the deviation of boundary variables from a fixed auxiliary variable z, which is information stemming from neighboring regions, is penalized. The regions then exchange information and z is updated. The update is re-distributed to the local agents (areas) for a new OPF calculation until consensus between regions is achieved. It is constructed an augmented Lagrangian of the form with penalty parameter ρ ∈ R. Dual variables of the consensus constraints are denoted with λ k ∈ R n×1 , and W ∈ R n×n is a positive definite, diagonal weighting matrix, 1 where each entry is related to one coupling constraint. We refer to W (S) if the entry is related to a power variable, and to W (V ) if the entry is related to a voltage variable. In ADMM literature, W is the identity matrix. The main steps during one iteration of the solving process are . The first step (11a) minimizes a non-linear problem, where constraint region enforces local constraints (4b). Since z is fixed, (11a) is in fact a series of R independent problems which can be calculated in parallel. The second step minimizes a coupled quadratic problem, where enforces consensus of auxiliary variables z. In fact, (11b) calculates the average value between two consensus variables of neighboring regions [11,13] and the problem can be reduced to Wind PV G3 G2 G1 The weighting factors ρ and W can be neglected in (14), if the same weight is assigned to two coupled variables, which is a reasonable choice. Once a region has gathered necessary neighbor information, this step can be calculated locally as well [12]. In the third step (11c), dual variables are updated based on a weighted distance between x and z. Again, with given local x k and z k , each region can calculate λ k independently (Fig. 3).
In ADMM, an update rule on ρ can be useful to enforce consensus. A ρ k ∈ R is assigned to each region, which can be increased depending on the local residual, see (18b). If the residual has not decreased sufficiently compared to the previous iteration (indicator 0 < Θ ∈ R < 1), the penalty is increased by a constant factor of τ ∈ R > 1. The penalty parameter must be chosen carefully since it is widely known to be crucial for good convergence behavior [17]. An overview of the implemented ADMM is given in Algorithm 1.

Test Scenario
We use the illustrative AC-DC test system from [15], see Tables 1 and 2 for line and generator parameters, respectively. The generator at Node 5 is interpreted as a wind park and a PV park is connected to Node 3. RES upper power limit and load are variable over time. Forecast power profiles are taken from the Belgian transmission system operator Elia, which provides detailed RES generation and load data on its website. Adapted profiles for usage in the 5-bus test system are shown in Figs. 6 and 7. We assume quadratic generator cost functions as in [18] (see Table 3) and one Algorithm 1 ADMM 1: Initialization: Weighting matrix W , tolerance ; for all k ∈ R: initial guesses z k , penalty parameters ρ k = ρ, dual variables λ k = 0, local solutions x k = ∞, local residues Γ k = ∞. 2: while ||Ax|| ∞ > and ||x − z|| ∞ > do 3: Solve for all k ∈ R the decoupled NLPs Solve the coupled averaging step minimize z k∈R 6: Calculate local residues and penalty parameter updates for all k ∈ R Update Γ k ← Γ + k for all k ∈ R 8: end while cost coefficient for all reactive power injections. RES power injection is prioritized by assigning a low fuel price. Two energy storage systems with characteristics from Table 4 are connected to Node 3 and 5, respectively. Line flow limits of 400 MVA for Line 1-2 and 240 MVA for Line 4-5, respectively, are used. Furthermore, the capacity of each VSC is set to 200 MVA. Allowed voltage ranges lie between 0.9 and 1.1 pu on the AC side, and between 0.95 and 1.05 pu on the DC side.

Multi-area SC-D-OPF
To show applicability of the distributed approach to N-1 secure and dynamic OPF problems, a small example is presented in the following. An SC-D-OPF dispatch optimization is performed for the scenario set C = {0, Line 1-2 } and T = {1, . . . , 4}. Storages and VSCs are allowed curative control, that is, operating set points are allowed to move freely to a new set point after an outage. Contrarily, conventional generators are operated preventively, that is, the power set point must be valid both before and after the outage. However, to cope with changing system losses after an outage, a deviation of 1 % of installed capacity is allowed.
General convergence behavior is shown in Fig. 4. Consensus error and objective suboptimality show similar behavior compared to conventional OPF cases. Consensus error (||Ax|| ∞ ) which depicts cross-border feasibility, falls below the threshold of 10 −4 after 344 iterations. The deviation from centrally computed costs f * is denoted withf = |1 − f/f * |. The objective value error falls below 0.01 % when convergence is achieved.
The progress of storage variables is shown in Fig. 5 for all time steps after the outage. Dashed lines denote results from the centralized solution. Storage systems provide positive (ESS 2 ) and negative (ESS 1 ) power reserve after the outage. The amount of reserve is increased with time for a growing network stress level, but energy levels remain within limits. The distributed solution is represented by solid lines. Those oscillate around and eventually approach the target levels.

Conclusion
This paper presents an approach to calculate SC-D-OPF in a distributed manner. That is, each control area computes a local N-1 secure and dynamic optimal dispatch which respects boundary constraints from neighboring areas for each time step and topology. After iteratively exchanging information and updating the local solution, a consensus is found which allows for a feasible cross-border power flow for all possible scenarios. Additionally, the result is close to the central optimal solution. Optimal curative storage or VSC set points can be determined for a contingency scenario which is in a foreign area. In this work, we apply the method to a small test system -the next steps are to increase system size, time horizon and number of contingencies.    The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Introduction
While today's power grid infrastructure has been designed for centralised power generation in conventional power plants, the ever increasing share of renewable energy sources (RES) leads to an increasingly uncertain, volatile and decentralised generation. In order to ensure a reliable grid operation and to maintain today's security of supply, it is necessary to develop methods to simulate accurately a power grid at high temporal resolutions. This can be done by efficient methods for power grid optimisation, including an accurate consideration of the non-linear and nonconvex AC power flow constraints. And these are to be improved.
The task to find the optimal operating state of a given power grid, also known as Optimal Power Flow (OPF), is formalized as the minimisation of a cost function with respect to a vector of continuous optimisation variables like the generated active power and the generated reactive power. Dynamic Optimal Power Flow (DOPF) considers several OPF problems, each corresponding to one specific point in time. In this case, the power grid's power demand is time-dependent and additional constraints must be taken into account. These constraints couple some of the optimisation variables corresponding to different time steps. Among others, these couplings are introduced by energy storage facilities and ramping constraints for conventional power plants [20]. Since DOPF is a large-scale non-linear and nonconvex optimisation problem, solving it in parallel is of crucial importance.
For this kind of problems Primal Dual Interior Point Methods (PDIPM) have proven to be among the most powerful algorithms. The main computational effort, when PDIPM is applied, lies in solving linear systems of equations [14]. And this a task suitable to be carried out in parallel on multi-core CPUs.
There exist a few works on solving linear systems, that arise from PDIPM applied to DOPF, in parallel. One approach is based on parallel matrix factorization by means of Schur complement techniques [6,20]. Other strategies use Benders Decomposition to decompose the complete optimisation problem into smaller ones [1].
Our contribution is the use of a Schwarz Domain Decomposition Method as preconditioner for Krylov-type iterative methods for solving these linear systems of equations in parallel. The original Schwarz method was formulated in 1870 as theoretical tool for proving existence of elliptic partial differential equations (PDEs) on complicated domains [17]. Later on, modifications of it have been used as standalone iterative methods for solving PDEs and have become a standard technique for preconditioning Krylov-type methods in the context of PDEs [18]. We apply this technique to DOPF by decomposing the time domain of interest into several smaller subdomains, which allows the use of multiple processors for computation.

Dynamic Optimal Power Flow
We will now formulate the DOPF problem in a mathematically rigorous way. This formulation has already been stated in more detail by several authors, e.g. in [12,20] and [16]. We would like to mention that we consider a continuous formulation of the DOPF problem, i.e. one without any discrete decision variables. A non-linear mixed-integer formulation of an OPF is considered, for example, in [19].

Formulation of the DOPF
In order to formulate the DOPF problem for a given time period of interest [0, T ], we consider a uniform partition 0 The power grid consisting of N B nodes denoted by B := {1, . . . , N B } and N E transmission lines denoted by E ⊂ B × B, is described by an admittance matrix The admittance matrix is symmetric, i.e. Y = Y T , and furthermore Y lk = Y kl = 0 if and only if there is a branch connecting node k and l, i.e., (k, l) ∈ E and (l, k) ∈ E. Therefore, Y is a sparse matrix for most real world power grids.
The complex voltage at node k ∈ B for time step t ∈ T is given as the complex number with real part E t k ∈ R and imaginary part F t k ∈ R. We define E t := [E t k ] k∈B and F t := [F t k ] k∈B . Moreover, we define the following sets, which number: The resulting vector of variables corresponding to active power and energy for time step t is given by with P t ∈ R N P . In a similar way, the vector of reactive power injections is defined by Every variable related to active and reactive power can be assigned to a specific node of the power grid by means of the power injection and extraction matrices C P ∈ {−1, 0, 1} N B ×N P and C Q ∈ {−1, 0, 1} N B ×N Q [21, p.23]. Since the states of charge (SOC) of storage facilities do not directly influence the power flow equations, the corresponding rows in C P are equal to the zero vector.

Constraints
At first, we consider the equality constraints. We begin with the power flow equations: Denoting the active and reactive power load for time step t by L t p and L t q , respectively, and assigning the loads to the nodes with the help of the matrix C L ∈ {0, 1} N B ×N L , the AC power flow equations are given by [22] where which describe the sum of power flows at node k. Note that we dropped the index t to improve readability. In the following, we abbreviate the equations (5) and (6) with the help of the following "symbolic" equation: Besides the power flow equations (7), one has to take into account an equality constraint for the slack bus, given by F t 1 = 0 for all t, and for generators k, l ∈ DC modelling a DC line kl: Here, line losses are taken into account by means of the factor ν DC ∈ (0, 1). Moreover, an additional set of equality constraints has to be imposed for modelling the state of charge, SOC t i , of storage facilities. These constraints are given by with factors ν l , ν g ∈ (0, 1) describing the efficiency of the loading process P t sl,i and generation process P t sg,i , respectively. Secondly, we take a look at the inequality constraints. Due to technical restrictions arising from power plants, renewable energy sources, storage facilities and DC transmission lines, it's necessary to impose lower and upper bounds for active and reactive power injections: In our application the only bounds, which are time dependent, correspond to the power injected by renewable energy sources, i.e. P t r . Also the node voltages and the active power flow over transmission lines, denoted by p kl , must be constrained: Finally, we impose ramping constraints for conventional power plants to bound the rate of change for injected power between consecutive time steps by α is the ramping factor in percentage per unit time.

Cost Function
The cost function f accounts for expenditure and income related to active power processes for the complete time interval [0, T ] and is composed of contributions from each time step t ∈ T [5]: with f t (P t ) = i∈C (a C,i,2 (P t g,i ) 2 + a C,i,1 P t g,i ) + i∈R (a R,i,2 (P t r,i ) 2 + a R,i,1 P t r,i ) and coefficients a * , * , * ∈ R.

Optimisation Problem
With the definitions of the previous sections at hand, the DOPF problem can be stated in an abstract way: where the vector of optimisation variables is given by Moreover, we used and will use the following definitions: n x,t := 2N B + N P + N Q and n x := N T n x,t and, as of now, let n λ,t be defined as the number of equality constraints corresponding to time step t, i.e. λ t ∈ R n λ,t with λ t = (λ t I , λ t S ) denoting the Lagrangian multipliers corresponding to g t I and g t S , respectively. Analogously, μ t = (μ t I , μ t R ) ∈ R n μ,t denotes the Lagrangian multipliers for the inequality constraints h t I and h t R . Consequently, we define n λ := N T n λ,t and n μ := N T n μ,t . Note that the ramping constraints h R given by (13) and the storage constraints g S given by (9) establish the only couplings between the variables corresponding to different time steps. Without them, a solution to (16) could be obtained by solving N T independent OPF problems. The optimisation problem (16) is non-linear due to the AC equations and due to the voltage and line flow inequality constraints. 1 Moreover, the latter ones and the AC equations are the sources of non-convexity in (16).
To simplify the notation in the following sections, we abbreviate the optimisation problem (16) to with quadratic functions

Primal Dual Interior Point Method for DOPF
In this section, we briefly describe the application of the Primal Dual Interior Point Method (PDIPM) to a DOPF problem. The description is in line with the implementation found in the MATLAB package MATPOWER [21].

The PDIPM Algorithm
For the general optimisation problem (18), the corresponding Lagrangian function is Assuming additional constraint qualifications, e.g., the linear independence constraint qualification (LICQ) [14], for every local minimum x * of (18), there exist corresponding Lagrangian multipliers λ * , μ * such that (x * , λ * , μ * ) are in accord with the Karush-Kuhn-Tucker (KKT) conditions [5]: where s ∈ R n μ is a vector containing slack variables and S ∈ R n μ ×n μ denotes a diagonal matrix whose diagonal elements are S kk = s k . The PDIPM is an algorithm to find solutions to Eq. (21). Mathematically it solves a slightly modified version of this equation by applying to it Newton's method. The modification consists in adding a "term", which decreases with iteration index k. This term keeps the slack variables s away from zero and, thus, the inequality constraints inactive. This means that solutions (x (k) , s (k) , λ (k) , μ (k) ) to this modified version of Eq. (21) are inside the feasible region and not on its surface, hence, the method is called interior point method. More rigorously formulated, we find for a sequence of barrier parameters For details we refer the inclined reader to [5,8] and [21].
Remark A very basic PDIPM is applied, i.e. there is no globalisation strategy.
Since we focus on parallel linear algebra techniques for DOPF problems, a locally convergent algorithm is sufficient to our purpose.

KKT Matrix in PDIPM
The application of Newton's method to Eq. (22) yields iterations in which the following linear systems of equations must be solved: Henceforth we omit the iteration index k. Assuming that s, μ > 0, the Newton system (23) can be transformed into a reduced, symmetric saddle point form with diagonal matrix Σ = diag μ 1 s 1 , . . . , μ n μ s n μ and The KKT matrices A, arising in every PDIPM iteration k, may become more and more ill-conditioned when a KKT point is approached. This is caused by the sequence of scaling matrices Σ (k) . If strict complementary holds at the KKT point [14]. 2 For this reason, many IPM software packages (like IPOPT [7]) use direct methods such as LDL T -factorizations to solve the appearing linear systems. These methods are less sensitive to ill-conditioned matrices.
In contrast, iterative linear solvers like GMRES [15] are very sensitive to illconditioned matrices. A good preconditioner is needed to lower the matrices' condition numbers. In exchange, they offer a higher potential of parallelization and allow the use of inexact Interior Point methods, see Sect. 3.3.
The distribution of the spectrum of a matrix K, defined as is an indicator of the convergence behaviour of GMRES applied to Kv = b. By rule of thumb, a spectrum which is clustered away from 0 is beneficial to the rate of convergence of GMRES. [11,Th. 4.62] For general non-linear optimisation problems with the corresponding KKT matrix the so-called constraint preconditioner is an example of a preconditioner. It is given byK withH being an approximation to H , that is easy to factorize, e.g.,H = diag(H ) [14]. One can show that σ (K −1 K) = {1} ∪σ with |σ | = n x − n λ . Therefore, at most n x −n λ eigenvalues ofK −1 K are not equal to 1. The distribution of these remaining eigenvalues depends on how wellH approximates H on the nullspace of B [10]. 3 In Sect. 4, we describe how to take advantage of the special structure given by (16) to construct a parallel preconditioner. Furthermore, we will see that there is a close relation between our proposed preconditioner and the constraint preconditioner defined above, see Sect. 4.4.

Inexact PDIPM
In contrast to direct methods, iterative solvers allow solving linear systems with prescribed accuracy. One can exploit this fact by means of inexact Interior Point methods. Within these methods, the linear systems corresponding to the first PDIPM iterations are solved with a rather low accuracy and therefore a low number of iterations. As the PDIPM iterate approaches the exact solution, the accuracy of the linear solver is increased. We use the approach from [3] for determining the accuracy in each PDIPM iteration. Here, (23) has to be solved with residualr (k) , i.e, Therefore, "γ k +η k can be viewed as forcing sequence of inexact Newton methods" [3]. In our numerical experiments, we setη k = 0.1 for all k and computed the relative tolerance for the iterative solver as follows .
We stopped to decrease the relative tolerance when it reached a value smaller than 10 −10 .

Termination Criteria for PDIPM
We use the following criteria for monitoring the progress of PDIPM [21]: with x − denoting the previous iterate.

Schwarz Domain Decomposition Method for DOPF
The key idea behind Schwarz domain decomposition methods applied to DOPF problems is the decomposition of the set of time steps T into several smaller subsets. Due to this decomposition, it is possible to define corresponding submatrices of the global KKT matrix defined in (24) that are smaller in size and therefore faster to factorize. Furthermore, the factorization of these submatrices can be done in parallel. After computing subsolutions corresponding to these submatrices, one can reconstruct an approximate solution to the global system (24), which is used as preconditioner within a Krylov-type iterative solver. In Sect. 4.1 we introduce the mathematics, which describe the additive Schwarz Method (ASM) in the context of DOPF problems. In Sect. 4.2 we briefly describe some issues related to its implementation. Section 4.3 is intended to provide a deeper insight into the subproblems that have to be solved when the ASM is applied. Finally, we describe the relation between the ASM and the above mentioned constraint preconditioner in Sect. 4.4.

Mathematical Formulation of the Additive Schwarz Method
For the application of the ASM as a preconditioner for the KKT matrix A(x, s, λ, μ) ∈ R (n x +n λ )×(n x +n λ ) defined by (24), it is needed to decompose the set of time steps T = {1, . . . , N T } into q non-overlapping subdomains: Afterward, each subdomainT l is expanded by additional s ol time steps on both ends, yielding an overlapping decomposition of T (see Fig. 1): with Henceforward, all expressions involving˜refer to the non-overlapping subdomains T l . When restricting the optimisation variables to the components, which belong to T l , we write The combination of x [l] and λ [l] yields the following definition: This vector contains the components of the solution vector of (24) belonging to the time steps in T l . Moreover, let R l ∈ {0, 1} n l ×n be a restriction matrix corresponding to the subset T l defined by its action: The matrices R x l and R λ l are composed of either zero or identity blocks. As an example, consider the case of l = 1 and l = q: Note that the size of the zero matrices varies with each restriction matrix and each l. Analogously, we defineR l as restriction operator corresponding to the nonoverlapping subdomainT l .
With these definitions at hand, it is possible to extract submatrices A l from the KKT matrix A by multiplying it with R l : In Sect. 4.3, we explore the structure of these submatrices.
On the important assumption, that all submatrices A l are non-singular, we are able to formulate the multiplicative Schwarz Method (MSM) for approximately solving AΔ R = b. We do this in the form of an algorithm [18]: In every iteration l, the current error with respect to the exact solution Δ R , given by By restricting r l−1 to subdomain T l and solving with A −1 l , one obtains an approximationê to the exact error R l e l−1 on subdomain T l . Prolongation with R T l yields a correction δ l = R T lê l to the current approximation Δ l−1 R . Thus, the multiplicative Schwarz method can be seen as "defect correction algorithm". 4 Unfortunately, the the MSM is a sequential algorithm. In order to parallelize it, one omits residual updates in each iteration, yielding the ASM [18]: which can be written as The right preconditioned linear system for solving AΔ R = b is then given by Since M ASM is symmetric but in general not positive definite, we use the Generalized Minimal Residual (GMRES) method for solving the non-symmetric system (37). In general, the ASM computes a less accurate approximation to the solution of AΔ R = b and therefore needs an increased number of GMRES iterations compared to the MSM given by Algorithm 1. However, the ASM offers a higher potential of parallelization, which results usually in better parallel scaling. 4 In contrast to "classical defect correction" algorithms, there is no converging sequence

Implementation
In our implementation, we use one processor per subdomain T l and distribute A such that every processor storesR l A in its local memory.R l A contains the nonoverlapping portion of rows of A l = R l AR T l . To set up M ASM once, every processor first has to form its local submatrix A l , i.e., in this step the overlapping part of A l has to be communicated to processor l by processor l − 1 and l + 1. Afterward, each processor computes an LU -factorization of A l , i.e., A l = L l U l . This step doesn't involve any inter-processor communication and can be done in parallel.
The employment of the ASM as a preconditioner inside the GMRES method, requires to compute M ASM v (k) in each iteration k. v (k) denotes the k-th basis vector of the Krylov subspace K k (AM ASM , b). On the assumption that the basis vector's data layout is the same as the matrix's one, i.e. every process stores v (k) is done by one forward substitution with L l and a backward substitution with U , respectively. This step does not involve any communication. After this, the local solution is prolonged back to obtain the global result of the matrix vector product M ASM v (k) . This again requires communicating the overlapping part of v (k) l to process l − 1 and l + 1. Finally, it's necessary to multiply A with M ASM v (k) . Due to A's data layout and its off-diagonal elements, corresponding to intertemporal couplings, additional communication is required.
One can further improve the performance of ASM by using the so-called restricted version of ASM [4], which is given by For this preconditioner, just the non-overlapping part of the local solution v (k) l is prolonged instead of the entire (overlapping) vector. Experiments show a beneficial behaviour in terms of GMRES iterations compared to standard ASM [4]. Furthermore, prolongation byR l doesn't involve any communication.
To further stabilise the LU -factorization, a small regularisation parameter = 10 −8 is added to the KKT matrix A, i.e. A + I .

Structure of Local KKT Matrices
In this section, we investigate the structure of the local submatrices A l defined in Eq. (36). A good way to get an idea of their structure is to permute the KKT matrix A in such a way that all rows and columns belonging to one time step are grouped together. This yields a block diagonal matrix with some off-diagonal elements. Every block represents one time step and the off-diagonal elements arise due to the time-dependent constraints. Extracting the time steps belonging to the subdomain T l results in a permuted version of the submatrix A l . 6 Its structure is exactly the same as the one of the original KKT matrix A. If the off-diagonal elements are taken into account as "boundary conditions", it is possible to interpret the submatrix A l as representing an optimisation problem in its own, i.e. the original dynamical OPF problem is reduced to the subdomain T l . Figure 2 tries to demonstrate this. As a consequence the submatrices A l tend to be ill-conditioned.

Relationship Between the ASM and the Constraint Preconditioner
To see a relation between the additive Schwarz method and the constraint preconditioner, as introduced in Sect. 3.2, it is necessary to consider a dynamical OPF problem without time-dependent equality constraints. In our case, we had to drop the storage facilities g t S to do this. Furthermore, it is important to consider a nonoverlapping decomposition of the time domain, i.e. s ol = 0. In this case, the corresponding ASM preconditioner, reduces to a block Jacobi method with omitted couplings between variables belonging tot + l andt − l+1 for each subdomain l. Since there are only time-dependent inequality constraints left, the coupling between the variables belonging to different time steps has to arise in the (1, 1)-block of the KKT matrix, namely in the (∇h) T Σ(∇h) part. This means that multiplyingM ASM with b is like solving a modified KKT system, or equivalently, H denotes the KKT matrix's modified (1, 1)-block. Thus,M ASM has the form of a constraint preconditioner. At this point, it also becomes clear that the ASM can only be interpreted as a constraint preconditioner, if there are no time-dependent equality constraints. If there existed time-dependent equality constraints, the (2, 1)-and (1, 2)-block of M −1 ASM would have changed too. As pointed out in Sect. 3.2, 1 is an eigenvalue ofM ASM A of multiplicity 2n λ and the remaining n x −n λ eigenvalues are solutions of a generalized eigenvalue problem of the following form [10]: For a low number of subdomains q, one can expectH to be a good approximation to H , leading to a clustered spectrum σ (M ASM A) close to one. However, the more subdomains, the more neglected couplings and the less accurate doesH approximate H . This results in a more scattered eigenvalue distribution ofM ASM A, but still a large number of eigenvalues are equal to 1.

Numerical Experiments
In this section, we present some results for our proposed method applied to a DOPF problem. In its first part, we discuss the parallel speedup of our solver. In the second one, we investigate the way in which the KKT matrix's spectrum σ (A) changes with the iterations of the interior point method. The test grid, which we used for our experiments, represents a possible variant of the German Transmission Grid in the year 2023 and consists of 1215 nodes, 2247 AC lines, 8 DC lines, 181 conventional power plants, 355 renewable energy sources and 67 storage facilities. For more details, we refer the reader to [12]. For solving (16), we use the PDIPM algorithm mips which is written in MATLAB code and part of MATPOWER [21]. In this algorithm, we replace the standard MATLAB backslash operator \ for solving linear systems by our own linear solver. Our solver applies a right preconditioned GMRES method. We use the ASM as the preconditioner. For computing the LU -factorizations of the local systems A l , we use the software package Intel R MKL PARDISO. Our solver is written in C++ and makes use of the KSPGMRES and PCASM methods provided by PETSc [2], which is compiled in Release mode with the Intel compiler 17.0. The computation of the eigenvalues of the KKT matrices and the preconditioned KKT matrices have been done with SLEPc [9] and MATLAB. Inter-process communication is realised by means of the Message Passing Interface (MPI). The numerical experiments were carried out on the BwForCluster MLS&WISO Production at Heidelberg University with two Intel Xeon E5-2630v3 processors i.e. 16 CPU cores at 2.4 GHz per node, and 64 GB working memory.

Parallel Speedup
In our numerical experiments, we concatenated the 48 h load profile, which was part of the described test grid, to obtain a test case comprising 48 days with a temporal resolution of 1 h, i.e. N T = 960 time steps. This causes the linear system (24), which has to be solved in every PDIPM iteration, to have the dimension n x + n λ ≈ 6.168 · 10 6 . Furthermore, we set α = 0.4, which led to ramping constraints limiting the change of power generation of conventional power plants by 40%/h of their respective maximum power generation.
We set f eas = grad = cost = comp = 10 −5 as PDIPM termination criteria and limited the number of PDIPM iterations to 300. The number of GMRES iterations was restricted to 1500. The GMRES method was allowed to restart after 300 iterations. Its relative residual tolerance rtol was determined as described in Sect. 3.3.
The ASM method was not restricted and the overlap s ol was set to 3 for all tests. The parallel speedup on q processors for the first k PDIPM iterations is defined as s q (k) := T 1 (k) T q (k) . We computed the sequential reference time T 1 (k) by solving the linear systems (24) using a LU -factorization and adding up the times needed to solve them up to the k-th PDIPM iteration. MKL PARDISO's sequential LUfactorization algorithm was applied.
Analogously, T q (k) was obtained by adding up the times needed to solve the same linear systems, but this time iteratively, i.e. the GMRES method and the ASM were used. In Fig. 3 we plotted the results of our speedup computations for different numbers of processors. The dashed graph represents the speedup obtained when the KKT systems (24) are solved directly with MKL PARDISO's parallel implementation of the LU -factorization using 16 processors. For this amount of processors we observed the best speedup.
Besides, we scaled the y-axis logarithmically to ease drawing a comparison with the LU -solver and to emphasise the large speedups during the first iterations. 7 We observed for our proposed iterative solver, no matter the number of processors q, a clear speedup in comparison with the LU -solver. The parallelization, due to the ASM, shows its power during the first PDIPM iterations. Unfortunately, this initial speedup decreases rapidly. After 40 iterations we end up with a total speedup of five. This decrease is accompanied by an increase in GMRES iterations. To demonstrate this, we plotted in Fig. 4 the number of GMRES iterations needed to solve the KKT system at the k-th PDIPM iteration. The first thing to note is, that the number of GMRES iterations does not only change with PDIPM iterations, but also with the number of processors used to solve the KKT systems. It's likely that the increase in GMRES iterations with PDIPM iterations is caused by a change in the spectrum of the original KKT matrix A. The increase with the number of processors can be explained as follows: The more proccesors we use, the smaller the subproblems, described by A l , become and, hence, the more intertemporal couplings are neglected, which reduces the "approximity" of M ASM to A −1 . And the closer M ASM is to A −1 , the more the spectrum of the preconditioned operator σ (AM ASM ) is clustered around 1.
Before we start to actual investigate the spectrum of the matrices in the next part of this section, we would like to summarise our findings. Even though we observed a rapid decrease in the speedup, it was possible to obtain a moderate speedup with only 16 processors and in comparison to a standard parallel LU -factorization also using 16 processors, our solver was more than a factor 2 faster.

Eigenvalue Distribution
Since the GMRES method converges faster if the eigenvalues of matrix A are clustered around 1, i.e. it needs less iterations, we would like to investigate how much the ASM, applied as preconditioner, improves the eigenvalues' distribution. Therefore, we calculated the eigenvalues with the largest and smallest magnitude of the KKT matrix A and the preconditioned KKT matrix AM ASM , respectively. We Seeing that it takes very long to compute, or that it is even not possible to compute the eigenvalues for the complete N T = 960 time steps, we decided to calculate them for N T = 96 time steps and for each PDIPM iteration. We set the overlap to s ol = 1 and used q = 32 subdomains. The results are presented in Fig. 5. The red lines depict the development of λ A max and λ A min of the KKT matrix, respectively. Thus, they form the limits for all possible magnitudes of eigenvalues corresponding to A. The latter is indicated by the red shaded area. Analogously, the blue lines and the blue area represent the spectrum of the preconditioned KKT matrix.
The spectrum of the preconditioned KKT matrix, is clustered around 1 in every interior point method iteration. Looking at the magnified part of the graph shows

Conclusion
In this work we proposed a way of solving linear systems arising from Dynamic Optimal Power Flow (DOPF) problems in parallel by means of an overlapping Schwarz domain decomposition method, namely the additive Schwarz method. It was shown how to apply this method in the context of DOPF problems and that the obtained submatrices correspond to localised formulations of a DOPF problem with additional boundary conditions. Numerical tests on one large-scale DOPF problem showed that the combination of the Generalized Minimal Residual (GMRES) method and the Additive Schwarz Method (ASM) is able to solve DOPF problems. Furthermore, we were able to show that this combination yields good parallel speedup for some of the PDIPM iterations. And, in a small example, we demonstrated that the ASM is a good preconditioner in the sense, that it clusters the eigenvalues around 1.
Unfortunately, the proposed solver had problems with some of the KKT systems arising in the large-scale DOPF problem. Therefore, it is necessary to improve it and this will be one of our main goals in future work. 19 Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Introduction
Renewable power sources have an ever increasing share of all power sources. Though renewable energy has been developed in recent years with great success, its intermittent and unpredictable nature raises the difficulty to balance the energy production and consumption [6,18]. A frequent proposal is to use gas turbine plants to compensate for sudden drops in power of renewable sources because these plants are relatively flexible in comparison to coal or nuclear plants. A welcome advantage of this approach is the possibility to run gas plants with fuel produced from renewable electricity via power-to-gas plants, thereby reducing or even negating carbon emissions of the gas plants. For a review of power-to-gas capability see [6].
It is desirable to have a joint optimal control framework for power and gas sector of the energy system to model this compensation. So far only steady-state flow in the gas network has been considered [3,18,20], which may be too coarse for several applications. Therefore, we focus on an optimal control strategy for the instationary gas network model [1] coupled to a power grid [2] via compressor stations, see for example [8,15,17]. The mathematical foundation for the gas-to-power coupling has been recently introduced in [5], where conditions for the well-posedness have been derived and proved. This work is the first attempt to model this interaction and yields an understanding of the underlying equations. The next step will be realworld scenarios.

Optimal Control Problem
The gas dynamics within each pipeline of the considered gas networks are modeled by the isentropic Euler equations, supplemented with suitable coupling and boundary conditions. For the power grid, we apply the well-known powerflow equations. The coupling between gas and power networks at gas-driven power plants is modeled by (algebraic) demand-dependent gas consumption terms. To react on the demand-dependent influences on the gas network, controllable devices as compressor stations are considered within the gas network. The aim is to fulfill given state restrictions like pressure bounds whereas at the same time the entire fuel gas or power consumption of the compressor stations is to be minimized.
The network under consideration is similar to the one presented in [5] and is depicted in Fig. 1. In addition, there is now a compressor with a time-dependent control u(t), which is used to satisfy pressure bounds. The optimal control problem is to minimize compressor costs while satisfying power demand and gas dynamics: Mathematically, this is an instationary nonlinear optimization problem constrained by partial differential equations, see [9] for an overview. To solve the problem, we make use of a first-discretize-then-optimize approach and apply the interior point solver IPOPT [16]. The necessary gradient information for IPOPT, i.e., gradients with respect to all controllable devices, is efficiently computed via adjoint equations. Here, the underlying systems can be solved time-step-wise (backwards in time), where additionally the sparsity structure is exploited.
S0 S17 S4 S8 S20 S25 Fig. 1 Gas network connected to a power grid. Red nodes are PQ/demand nodes, green nodes are generators (PV nodes) and the blue node is the slack bus (also a generator, with gas consumption of the form ε(P ) = a 0 + a 1 P + a 2 P 2 ). The circle symbol indicates a compressor station We remark that the cost function is given in Sect. 2.3, while the bounds on the pressure are introduced as box constraints within the numerical optimization procedure in Sect. 3. Within the Sects. 2.1, 2.2, 2.3, 2.4, and 2.5, we now describe the constraints of the optimal control problem in detail and focus on the technical details.

The Isentropic Euler Equations
The gas network is modeled by the isentropic Euler equations (see [1,4]), which govern gas flow in each pipeline between nodes, where ρ is the density, q = ρv is the flow, p(ρ) = κρ γ is the pressure function. In our example we use γ = 1 and κ = c 2 = (340 m s ) 2 , that is, the isothermal Euler equations with speed of sound c. The Euler equations must be met for 0 ≤ x ≤ l e and t ≥ 0, where l e is the length of the pipe. Furthermore, S is a source term, where λ is the solution to the Prandtl-Colebrook formula, The Reynolds number Re is given by The roughness k e and diameter d e of the pipes are all the same in our example and given by k e = 5 · 10 −4 m, d e = 6 · 10 −1 m.

Coupling at Gas Nodes
At the nodes we use the usual Kirchhoff-type coupling conditions: The pressure is the same near the node in all pipes connected to it and the flows must add up to zero (where the sign for inflow is positive, the sign for outflow negative), ingoing pipes The example gas network we are using is a small part of the GasLib-40 network [10]. In Table 1 the only remaining parameter, the length l e of each pipe, is gathered. Note that no length is given for the arc connecting S0 and S17, because only a compressor is situated between these nodes.

Compressor Stations
To compensate for pressure losses in the gas network, we consider compressor stations, which are also modelled as arcs. Those arcs have (time-dependent) inand outgoing pressure (p in , p out ) and flux values (q in , q out ). In general, the two separate flux values allow the modelling of fuel gas consumption of the compressor station, whereas we will consider an external power supply for the compressor and therefore have q in = q out . The power consumption is modelled as a quadratic function of the power required for the compression process. Denoting this as function P (p in (t), q in (t), p out (t), q out (t)), our objective function in the optimal control problem is of the form Note that the power consumption does not influence the network dynamics and is therefore only of interest for the optimization procedure. For our investigations below it is sufficient to know that the consumption and therewith the costs increase if the ratio p out /p in increases. For the details of the power consumption model, we refer to [17,Section 3.2.3]. The influence of the compressor station on the network dynamics is modelled by the control u(t) of the pressure difference:

Power Model
For the power grid we use the AC powerflow equations (see [7] for an introduction), where P k , Q k are real and reactive power at node k,|V k | is the voltage amplitude, φ k is the phase (and φ k,j = φ k − φ j ) and B kj , G kj are parameters of the transmission lines between nodes k and j or of the node k for B kk and G kk . Each node is either the slack bus (in our case N1; V and φ given), a generator bus (N2 and N3; V and P given) or a load bus (N4 through N9; P and Q given). All in all for N nodes we get 2N equations for 2N variables.
The considered power grid is taken from the example "case9" of the MAT-POWER Matlab programming suite [19]. A per-unit system is used, whose base power and voltage are 100 MW and 345 kV respectively. The corresponding node and transmission line parameters are gathered in Table 2, these are the entries of the nodal admittance matrix (see [7]).

Coupling
The last ingredient is a model for converting gas to power at a gas power plant. In our example, it will be situated between the nodes S4 and N1 and convert a gas flow ε into a real power output P according to The flow ε must be diverted from the pipeline network, hence the coupling condition at node S4 must be changed to The details of simulation of such a combined network and the treatment of all arising mathematical issues can be found in [5]. We now showcase a concrete example.

Problem Setup
As already noted, we consider a small part of the GasLib-40 network from [2] consisting of 7 pipelines with a total length of 152 km. This network is extended by a compressor station and additionally connected to a power grid with 9 nodes by a gas-to-power generator. For this coupled gas-power network, we simulate a sudden increase in power demand within the power grid and study its effect on the gas network. The considered compressor station is supposed to compensate part of the pressure losses in the gas network such that a given pressure bound is satisfied all the time, while power consumption of the compressor is minimized.
To complete the problem description, the following initial and boundary conditions are given: -P (t) and Q(t) at load nodes, -P (t) and V (t) at generator nodes, -V (t) and φ(t) at the slack bus, p(t) at S5, q(t) at S25, p(x, 0), q(x, 0) for all pipelines. Node  More precisely, the initial conditions for the power network are given in Table 3. These remain constant over time except for the power demand at node N5, which changes linearly between t = 1 h and t = 1.5 h from 0.9 p.u. to 1.8 p.u. for the real power and from 0.3 p.u. to 0.6 p.u. for the reactive power, see also Fig. 2.
For the gas network the incoming pressure at S5 is fixed at 60bar, the outflow at S25 is fixed at q = 100 m 3 s · ρ 0 A e , where ρ 0 = 0.785 kg m 3 and A e = π d 2 e 4 . The fuel consumption parameters we use in equation (5) are given by a 0 = 2, a 1 = 5, a 2 = 10. Since the data for the considered gas and power networks are taken from different sources, the parameters of the gas-to-power generator are chosen in such a way that a significant influence is caused.
Further, the pressure at S25 is supposed to satisfy a lower pressure bound of 41bar, i.e., for all times t, where we consider a time horizon of T = 12 h.
As we will see below, the compressor station between nodes S0 and S17 will have to run at a certain time, i.e., u(t) > 0, to keep this pressure bound. In general, a solution to the described optimal control problem consists of the control u(t) and the entire network state for all times t, i.e., for all nodes in the power grid, p(x, t), q(x, t) for all pipelines in the gas network, fulfilling the given model equations and the pressure constraint.

Discretization and Optimization Schemes
To solve the described optimal control problem, we follow a first-discretizethen-optimize approach. The model equations of the power grid only require a discretization in time, which means that the given boundary conditions and the powerflow equations (4) hold for discrete times t j = j Δt with Δt = 15 min and j ∈ {0, . . . , 48} in our scenario. The discretization of the isentropic Euler equations within the pipelines of the gas network additionally requires a spatial grid (here with grid sizes Δx e ≈ 1 km) and an appropriate discretization scheme. Here we apply an implicit box scheme [14], which allows the application of large time steps as Δt = 15 min for the considered spatial grid. Considering the isentropic Euler equations as a system of balance laws of the form the applied scheme reads The numerical approximation of the balance law is thought in the following sense: Together with the algebraic equations modelling the compressor station and the coupling and boundary conditions, the discretization process results in a system of nonlinear equations for all state variables of the coupled gas-power network. For simulation purposes, the entire discretized system is solved with Newton's method. Note that the system can be solved time-step per time-step and that we exploit the sparsity structure of the underlying Jacobian matrices.
So far we can only compute the state of the considered gas-power network (for a given compressor control u(t)) and evaluate quantities of interest like the power consumption of the compressor station or the pressure constraint within the time horizon of the simulation. For the given time discretization, the compressor costs (3) are approximated by the trapezoidal rule and formally contained in J (u, y(u)) below. Next, we want to solve the (discretized) optimal control problem, i.e., to find control values u(t j ) such that the pressure constraint is satisfied, while the power consumption of the compressor is minimized. For this purpose, we apply the interior point optimization code IPOPT [16] to the (reduced) discretized optimal control problem. In addition to our simulation procedure to evaluate quantities of interest for a given control, IPOPT further requires gradient information of those quantities with respect to the control. Based on the considered discretization, such information can be efficiently computed by an adjoint approach, which we briefly describe in the following. Therefore, we consider the discretized optimal control problem in the following abstract form: The vector u contains all control variables (here the compressor control u(t j ) for all times t j = j Δt ∈ [0, 24]) and y contains all state variables of the coupled gaspower network of all time steps t j . The mapping u → y(u) is implicitly given by our simulation procedure. Thus, IPOPT does not have to care about the model equations formally summarized in E(u, y), but only about the further constraints g(u, y) and minimizing the objective function J (u, y). Accordingly, we need to provide total derivatives of J and g with respect to u.
In the following, we consider the computation of these derivatives only for the objective function, since the procedure is identical for the constraints, and we follow the description given in [12]. First of all, the chain rule yields While the partial derivatives of J with respect to u and y can be directly computed, the derivatives of the states y with respect to the control u are only implicitly given. Differentiating the model equations E(u, y(u)) = 0 yields and therewith (formally) Even though the partial derivatives on the right-hand-side of (8) can be directly computed, one would have to solve Nu systems of linear equations here. Instead of that, we insert (8) into (7) and get With the so-called adjoint state ξ as the solution of the adjoint equation we finally have It is the fact that (9) is a single linear system and has a special structure, which can be easily exploited (see for instance [11,13]), which makes the computation of derivatives via the presented adjoint approach very efficient. Nevertheless note that for given control variables u one still has to solve the model equations to get y(u), before one may compute the gradient information via (9) and (10).

Results
We first discuss the simulation with inactive compressor. In the course of the simulation, due to the increase in power demand at node N5, the power demand at the slack bus rises as well and leads to increased fuel demand at node S4. This increases the inflow into the gas network, as can be seen in Fig. 3. Also the pressure in the final node S25 decreases and violates the pressure bound after approximately 4 h, see Fig. 4. After the optimization procedure, the compressor station compensates part of the pressure losses in the gas network such that the pressure bound is satisfied all the time. Since the power consumption of the compressor station is minimized within the optimization, the pressure constraint is active after roughly 4 h (see again Fig. 4), i.e., the compressor station applies as little as possible energy.

Conclusion
The proposed optimization model allows to predict pressure transgressions within a coupled gas-to-power framework. Simulation and optimization tasks are efficiently solved by exploiting the underlying nonlinear problem structure while keeping the full transient regime. This makes it possible to track bounds much more accurately than with a steady state model, thereby achieving lower costs. 17 The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Utilising Distributed Flexibilities in the European Transmission Grid
Manuel Ruppert, Viktor Slednev, Rafael Finck, Armin Ardone, and Wolf Fichtner Abstract In this paper, we investigate the effect of distributed flexibilities on the operation of the transmission grid. The flexibilities considered are heat pumps, electric vehicles, battery energy storage systems and flexible renewable generation. For this purpose, we develop a two-stage approach of first determining an optimal electricity market solution considering the optimal dispatch of each generation element and flexibility. In the second step we determine the required dispatch adjustments due to transmission grid constraints and investigate the effect of integrating battery energy storage systems into the adjustable generators to solve congestions. In our case study, we investigate the central European transmission grid for a scenario based on the Distributed Generation scenario of the Ten-Year Network Development Plan for the year 2030. Integrating distributed flexibilities leads to a strong increase in the security of supply, while the overall effect on the generation adjustment is small. A comparison of the results for an AC and DC formulation shows that both approaches differ significantly in individual cases.
Keywords AC optimal power flow · Battery energy storage systems · DC optimal power flow · Distributed flexibilities · Distributed generation · Flexible demand · Transmission grid

Introduction
Today's change in the electricity system is driven by the decarbonisation initiative developed to meet emission reduction targets defined in international agreements in order to limit the global temperature increase. To reduce the carbon intensity of the electricity generation, an increasing amount of electricity generation from renewable sources (RES-E) is being installed throughout Europe. The two most  [1]. On the demand side, decarbonisation efforts have begun to lead to an increased electrification in the heat and transportation sectors. In some countries, rapid market growth of such elements which can be categorised as distributed flexibilities can be already be observed today. Future energy scenarios hint towards a further increase of these developments until 2030. With potentially high concurrencies of individual, decentralised demand patterns, uncontrolled or market-driven operation of these flexibilities will increase local as well as aggregated consumptions peaks in the power sector.
These changes pose significant challenges to the operation of the European transmission grid. First, balancing supply and demand will require additional flexibilities in a system dominated by distributed, intermittent RES-E and new sources of electrical demand. Second, RES-E at a specific point in time will be largely determined by the predominant meteorological conditions with high spatial concentration, increasing average distance between generation and consumption and thus increase the stress on the grid infrastructure. Lastly, market requirements to create a cost-minimal dispatch might contradict those necessary to enable a congestion-free transmission grid.
In the past, extensive research has been performed on different flexibility types concerning the technologies, their modelling and their contribution to renewable integration. A comprehensive overview on the current developments in electrical energy storage technologies and their potential for application in power system operation is given in [2] and [3]. Current reviews of different technological options can be found for power-to-heat applications [4], residential photovoltaic-battery energy storage (PV-BESS) [5], heat pumps (HP) in smart grids [6], and battery electric vehicles (BEV) [7]. Most literature focuses on the local balancing of demand and supply, either on a single household level or a community. This neglects the implications on the grid level or focuses on the regional effects in micro grids while including only a small set of technological options.
There are numerous models, which take into account flexibilities and simplified transmission grid constraints. In most cases, these models have been developed to analyse electricity markets, e.g. for investment decision support or operation decision support. In most studies liberalised markets in the United States or Europe are analysed and grid constraints are represented in a linearised way or by import/export constraints between market zones. Hence, most techno-economic models only take these constraints into account.
Few models, which consider flexibilities in a grid context take into account the full AC formulation of transmission grid constraints. Besides the group of studies and models, which consider a linearised form of grid constraints, some models can integrate full AC constraints. However, they are limited to either a very short time horizon or their focus lies on transient technical aspects like fault level detection, transient stability, harmonic analysis, reliability and power flow. Another group of models takes into account many flexibilities and AC restrictions but only on a distribution grid level, thus limiting system size and thereby complexity significantly [8]. A comprehensive overview on modelling approaches and their consideration of the grid can be found in [9].
Another area where transmission constraints are explicitly considered is the co-optimisation of transmission and generation capacities. The complexity of the resulting problem forces modellers to apply linearised transmission grid constraints or reduce the number of modelled flexibilities [10]. Furthermore, such models are often only validated on small scale or test grids, making it difficult to assess their suitability for large real-world power grids [11]. In most cases, the focus lies on investment decisions rather than on utilising operational flexibility in future transmission systems. [12] provide an overview of requirements and approaches to address this problem.
In this paper, we present a modelling framework to investigate the effect of distributed flexibilities on the transmission grid operation and investigate the benefit of utilising distributed flexibilities to reduce grid congestion while taking into account the full AC model of the transmission grid.

Modelling the Interaction of Market and Grid Operation
In order to model today's interaction between the electricity market and transmission grid operation, we developed a two-step approach to model market and grid operation. In the first step, we use a linear optimisation approach to determine the minimal-cost, copperplate-based results of the interconnected electricity market. In the second step, we determine the required dispatch adjustments due to congestion in the transmission grid using a multi-period nonlinear optimisation model. In the following section, the general modelling of the electricity market is described, followed by the integration of distributed flexibilities and the application of the method for modelling the transmission grid operation.

Electricity Market
To model the interconnected electricity market consisting of multiple market areas, we use a linear optimisation approach with the objective function of minimal total annual system cost under the assumptions of perfect foresight and perfect competition. 1 The system cost consist of the aggregated variable cost of electricity production C g for the electricity generation p g,t of generator g in time step t and the cost of load shedding C LS for the amount of load shedding p LS d,t of demand d in time step t. min P g,t ,P d,t g∈G,t ∈T In every market area m, the electricity demand P d,t in every time step must be balanced out against generation, load shedding (LS) and interconnector flows into the market area p c,t .
The bounds of each generator are determined by the minimum capacity P min g,t and maximum capacity P max g,t of the respective generating unit. While the boundaries of conventional generation are determined by the technical parameters of the generator, the boundaries of RES-E are determined by the available generation due to the intermittent nature. Available interconnection flows between market areas are restricted by the available exchange capacity P min c , P max c . In order to ensure feasibility of the problem, LS up to total demand can be performed at each time step.

Modelling Flexibilities
In the following, we apply a set of general equations to model various types of time-coupled flexibilities, ranging from different power storage technologies to consumption and renewable flexibilities. Ignoring self-discharging losses, the energy level v s,t of a storage s ∈ S at time step t ∈ T = {1, . . . NT } equals its previous level, the external power inflow ζ in s,t or outflow ζ out s,t and the sum of charged and discharged power provided by connected generators p g,t ∈ G S ⊆ G: v s,t = v s,t−1 + g∈G inS p g,t * η g − g∈G outS p g,t /η g + ζ in s,t −ζ out s,t ∀s ∈ S, t ∈ T For a simpler notation, the set of storage-connected generators G S is split such that G outS denotes those generators supplying power to the electricity system by discharging a storage and G inS vice versa (charging of a storage). The charging and discharging efficiency is denoted by η g . Furthermore, all variables might be restricted to lower and upper bounds and the starting condition for the storage is either a coupling of the first and last time steps or the definition of an initial storage level:

Storage Technologies
Modelling the flexibility provided by actual storage technologies S S ⊆ S based on the equations shown above is straightforward and could easily be applied to hydro storages and batteries. Pure hydro pump storages and large-scale battery systems allow a rather simple definition of variable bounds. Storage volumes and generator capacities are non-negative and time independent, so that the equations for the energy level and power output can be reformulated.
While large-scale battery systems allow to ignore external power inflow and outflow, which do not result from indirectly connected generators or consumption, this is only true for non-cascading pump storage systems. However, we will ignore the case of cascading systems in the following and refer to the description of its modelling [13]. Including a time series of natural inflow and the restriction of a potential minimal outflow, allows for modelling mixed-hydro pump storages and/or pure seasonal storages, with an empty set G inS in the last case. 2 Concerning small scale batteries, such as BESS combined with a PV system, the mathematical modelling is identical to large-scale battery systems, as long as no further restrictions apply to the battery utilisation. In practical application, however, regulators might apply specific restrictions for the operation of such systems. In Germany, for example, a widely used support scheme for PV-BESS is limiting the grid feed-in up to 50% of the PV capacity. 3 Given a demand P d,t and a maximum grid feed-in P lim g , these restriction might be easily implemented:

Demand Flexibilities
In general, consumption processes with the ability to store energy can be directly modelled as single storages. Under the assumption that a certain share of a specific load might be shifted within a certain time range, the storage restrictions, however, might be applied to a broader range of demand flexibilities or for modelling the flexible potential of aggregated loads. Considering our target to model the impact of demand flexibilities at high and extra high voltage levels, we focus on modelling the flexibility for large-scale consumption process or of aggregated loads as storages S L ⊆ S. When looking at the load shifting potential of an EV fleet, aggregation allows to neglect technical restrictions of single units and to derive the load shifting potential from the specific load profile [14]. This finding is supported by [15] where an aggregated modelling of BEV as flexible loads led to almost the same results in a unit commitment problem as their individual representation, after the BEVs where clustered accordingly. Therefore, we model the demand flexibility of utility scale HP and the aggregation of numerous BEV and household heat pumps (HH-HP) by defining certain share α of a specific load P d,t as flexible By restricting the storage outflow to the flexible demand profile We define that the flexible demand either directly translates to a system load or is stored. In case that load shedding is not allowed, Z min s,t equals Z max s,t . Assuming that the energy of a flexible demand profile has to be consumed within a time interval T T m , the bounds of the storage volume are defined by the integral of the profile over T T m and equal zeros in the last time step.

Renewable Flexibilities
Given a generation profile, a flexible operation of RES technologies, which depend on a storable resource, such as hydro-or bioenergy, might be modelled analogously to the case of flexible demand. This might be useful in cases where the actual storage parameters are unknown or for modelling virtual power plants, comprising multiple small units. The general idea is to model the flexibility by enabling a shift of the initial generation profile within a certain time interval. If a certain share α of the renewable generation is based on stored energy, we first split the generation variable into a flexible and inflexible part: We define that the flexible renewable generation is either directly feed-in or stored.
In case that dumping the flexible renewable generation is not allowed, Z min s,t equals Z max s,t . Assuming that the energy of a flexible renewable generation profile has to be fed in within a time interval T T m , the bounds of the storage volume are defined by the integral of the profile over T T m and equal zero in the last time step.
Concerning the bounds for the storage discharge, the minimal or maximal values of the profile time series within the certain time interval T T m or the whole time horizon can be used.
∀g ∈ G outS s , t ∈ T In order to model the specific flexible renewable technologies, only the definition of the time intervals and α f lex , α + , α − has to be adjusted.

Grid Operation
Based on the dispatch result for each generation unit and time step from the market model, we determine the necessary dispatch adjustment due to grid constraints in the second step using a nonlinear optimisation approach. As a conventional cost minimising objective function under consideration of grid constraints would lead to a nodal pricing result, which would discard the dispatch determined by the market, we establish the minimum amount of dispatch adjustment in the entire network as the optimisation objective and evaluation figure to assess the efficiency in the transmission grid. Therefore, the generalised objective function can be formulated as follows: Positive and negative deviation of generation units from the base solution P g,t are considered equivalently towards the objective, while grid-induced load shedding P LS d,t is unidirectional by definition. We solve the resulting optimisation problem with a multi-period nonlinear formulation of the ACOPF utilising an augmented Lagrangian formulation for coupling of time-dependent storage units as described in [16]. In order to eliminate the possibility of results representing significantly suboptimal solutions caused by converging into local optima due to the nonlinear nature of the problem, we benchmark our results with a DCOPF formulation, as described in [17].

Case Study
In the following, we apply the presented modelling approach on a scenario derived from the scenario "Distributed Generation" (DG) of the European Network of Transmission System Operators for Electricity (ENTSO-E) Ten-Year Network Development Plan (TYNDP) 2016. In our study, we restrict the scope of both electricity market and grid to the central European countries and their transmission networks, thus not considering the effect of flows from further interconnected countries. The considered area contains the countries France, Belgium, Netherlands, Switzerland, Austria, Czech Republic, Poland, Germany and Luxembourg. These countries form individual market areas, with the exception of Germany and Luxembourg, which form a common market area.
In addition to the scenario data from the DG scenario, the closely related high renewables scenario C of the german network development plan 2030 was used for Germany due to the fact that this source includes data with a higher level of detail. Based on the projected national RES capacity development of the scenario an optimal allocation planning model as described in [18] was run in the first step for regionalisation. In detail, a cost optimal expansion planning with a high spatial resolution until 2050 is performed, considering national and regional lower and upper bounds for the development of single renewable technologies. Due to the longer optimization horizon, exceeding the forecast horizon of the underlying scenarios, a minimal RES-E share on the demand is defined in the later years, such that an 80% renewable share is achieved in Europe until 2050.
In the second step, the capacity adequacy is analysed based on a Monte-Carlo simulation of the generation availability. Based on an integrated European approach, accounting for a cross-border balancing, and the availability of a flexible demand shift (BEV and heap pumps in households and district heating) and flexible renewables (hydro storage and biomass) we computed the needed dispatchable capacity. 4 Afterwards we utilised results obtained from a generation expansion planning problem of the European power system in order to meet the calculated capacity demand for a secure operation of the system [19].
In the current scenario, we assumed that existing power plants are decommissioned at the end of their technical lifetime. For coal-fired power plants, a premature decommissioning pathway until 2040 is assumed and individually adjusted for each unit with regard to its technical parameters and national regulations, aligned as much as possible with the political guidelines on coal phase out known today. The portfolio of hydro power plants with storage or pump-storage capability is based on today's generation and announced construction or expansion projects under the assumption, that generators reaching the end of their lifetime are replaced with the same parameters.
In the presented case study, 555 existing or currently projected and 145 new built thermal power plants and 110 hydro power plants are operational in the scenario in 2030. The resulting thermal, hydro pump and hydro turbine capacity and the location of their connection to the transmission grid are shown in Fig. 1.
The case study year with input data for transmission grid as well as generation and demand is chosen as 2030 and the simulations are performed with an hourly time granularity. For the transmission grid part, a weekly optimization horizon of transmission grid without rolling horizon was chosen. For this, storage states of each unit at start and end time were fixed to the state given by the previously determined state after market-based dispatch.

Transmission Grid Model
The transmission grid model contains the interconnected AC network of the same countries considered on the market side, including overhead lines and cables of voltage levels above 200 kV. For the studied area, this includes the nominal voltage levels of 220 and 380 kV. Additionally, high voltage DC (HVDC) lines are included in the dataset. HVDC lines connecting offshore wind generation are considered at the point of RES-E calculation and allocation and subsequently reduced from the dataset. All AC and HVDC lines are connected to busbars, which are assigned to georeferenced substations. The number of busbars at each substation is limited to the number of voltage levels present at the respective substation. In addition to the transformation between the extra high voltage levels, the transformation to the high voltage levels between 60 and 150 kV is considered and used for the connection of small-scale generation and demand as described in [18].
The data include the present state of the transmission grid as well as projected expansion measures in terms of deconstruction, replacement and construction of substations, busbars, lines and transformers until the year 2030. The technical data for each grid element was derived from publicly available sources whenever possible and otherwise approximated based on available information from predominant comparable equipment in the same geographical and organisational area. Among the sources used for completing the dataset are the static grid models of the transmission grid operators of the Central Western Europe Region (CWE), as well as the grid dataset of the TYNDP 2016 and Open Street Map (OSM). The future expansion of the national grid networks was obtained from various sources such as national network development plans of the respective grid operators and the projects listed in the TYNDP 2016. The model represents a snapshot of the projected grid topology in 2030 and includes a total of 3010 busbars, distributed over 1513 substations. The busbars are connected by 4103 AC lines and transformers, as well as six HVDC lines, with the exception of the interconnector between Belgium and Germany being located in the German market area. Furthermore, 237 reactive power compensation elements are located in the entire system, with either positive or negative reactive power provision. The snapshot of the transmission grid topology for 2030 is shown in Fig. 2.
The interconnection capacities between market areas are determined based on the available thermal capacities of interconnecting lines from the transmission grid dataset, assuming that 40% of the total thermal capacity of interconnecting AC lines between the market areas and the full thermal capacity of interconnection HVDC lines is being made available for market operation. Generators from nonintermittent sources are assumed to have reactive power provision capabilities of cosϕ = 0.925 w.r.t the installed generator capacity. Furthermore, RES-E with a nodal distribution of installed capacity as shown in Fig. 3 is assumed to be able to

Parameters for Flexibility Modelling
Considering the modelling of new consumers, their flexibility and their regionalisation the following assumptions were made:

Heat Pumps in Households
The modelling of electricity consumed due to heat demand in households from HPs follows the approach of the Munich City Utilities (SMW) for calculating the standard load profile of HPs within their grid, based on the yearly household load and the temperature. Assuming a full flexibility (α f lex = 1, α + = 1, α − = 0) within a 3 h horizon ( 00 : 00 − 03 : 00) , 03 : 00 − 06 : 00) , . . . .), the regionalisation follows the same approach as that of the household demand, restricted to single and two family buildings.

Heat Pumps in District Heating
We apply the same approach for modelling the generation profile and its flexibility as it is the case for HH-HP. Due to the availability of a heating grid and the possibility to store larger amounts of energy, a full flexibility within a 24 h horizon 00 : 00 − 24 : 00) is assumed. The regionalisation is in general based on the capacity of combined heat and power (CHP) units within a country. In Germany furthermore the data of the German District Heating Association (AGFW) for the actual district heating grids demand and location where combined with the CHP database for a more detailed analysis.

BEV PBEV
Profiles are based on [14,20] for a direct charging at arrival and adjusted based on Table 1.
The regionalisation corresponds to that of the household demand in all European countries except of Germany, where the vehicle registration numbers at NUTS 3 level were used as a distribution key instead of household number and electricity demand. Concerning the flexibility, we assumed the potential to shift the full demand (α f lex = 1, α + = 1, α − = 0) within a 12 h horizon either during the day or the night ( 06 : 00 − 18 : 00) , 18 : 00 : −06 : 00).

PV-BESS
In the current analysis, decentral BESS installed with PV on household level are modelled as simple battery storages and differ only in their sizing and regionalization from classical battery storages. Based on a linear relation between household demand and self-optimised PV-BESS capacity, obtained from a mixed integer optimization of PV-BESS sizing of households [21] in northern Germany and southern Germany, an upper bound for the regional distribution of PV-BESS is calculated. By adjusting the BESS capacity according to the actual installed rooftop PV smaller 15kWp in single and two family buildings, the key for the final PV-BESS distribution is computed.

RES-E Flexibilities
For the flexibility of hydro storage and biomass generation we assumed α f lex = 0.5, α + = inf , α − = 0, meaning that half of the generation profile is assumed to be flexible and restricted to the maximum value of the total time series of the flexible part. Furthermore, the potential to shift the generation within 24 h was assumed in case of biomass and of 168 h in case of hydro storage.

Scenarios
In the following, we will present the results obtained from the simulations performed with the input dataset described above. The results are presented for four scenarios, which differ due to the amount of utilisation of the flexibilities both on the market and the transmission grid side. In this, we divide the available flexibilities into

The Impact of Distributed Flexibilities
The yearly aggregated results are shown in Fig. 4. Transmission lines which reach their continuous thermal current limit during the time horizon are colored  Fig. 4. The highest number of fully loaded hours is found on the HVDC lines, implying that the start and end locations of these lines are generally well suited for relieving the stress on the AC grid. In the AC grid, areas with significant congestion can be found in Southern France, the border of France and Belgium, Northern Germany and the border area of Poland, the Czech Republic and Austria. Generally, the required decrease in generation at certain locations is reached by renewable curtailment rather than power decrease of thermal power plants. A trend of required power reduction in the southern part of the network can be observed over the entire year, while power increase is more focussed in the northern and eastern areas. When interpreting the results, possible discrepancies between the generation and demand scenario on the one hand and the grid model on the other hand have to be considered. For example, the grid expansion considered based on the TYNDP in France is significantly lower than the evolution of the grid in Germany in the time frame between today and 2030. Thus, the increase of both RES-E and new electricity demand types has an higher impact on a system that does not undergo a significant transformation. The impact of the utilisation of distributed flexibilities on the security of supply can be seen in the reduction of the required load shedding in Table 2. While both the Base and the NoFlex scenario require load shedding in order to achieve a feasible solution in more than 100 single hours during the year, this value is greatly decreased in the Flex scenario and even more in the DynFlex scenario. The main reason for this is the reduction in peak demand and generation due to the market-based dispatch of flexibilities as described in Sect. 2. Furthermore, the additional control capabilities -albeit with a small capacity over a longer time horizon -lead to less congestion events, which cannot be prevented without curtailing demand.
In the Base scenario, the average hourly positive dispatch adjustment requirement in the entire transmission grid area considered is 8.672 GW, with peaks reaching adjustments of up to 18.5 GW. As a permanent additional power provision is required in order to balance the transmission losses over lines and transformers, the positive adjustment does not decrease to zero even in cases without any violation of voltage or thermal limits. In the NoFlex scenario, the hourly requirement reduces by 48 MW while the peak increases to 20.1 GW. The hourly average decreases further to 8.54 GW in the Flex scenario and 8.497 GW in the DynFlex scenario, with the peak decreasing to 17.7 and 17.6 GW. The average negative dispatch adjustment required is not distorted by additional grid losses, over the scenarios it decreases by 3.6%, 2.1% and 1.0% from the NoFlex to the DynFlex scenario. Similarly to the positive adjustments, the negative yearly peak values are reaching the highest absolute values for the NoFlex scenario and the lowest for the DynFlex scenario.
Overall, the results show that as expected more flexibilities available both in the market and transmission grid lead to a lower amount of required redispatch. Yet, the different utilisations of both central and distributed flexibilities have a lower impact on the grid operation than expected. Among flexibilities the inclusion of largescale hydro in the Flex scenario has the largest impact, as observed in the yearly values of negative dispatch adjustment. Operating distributed flexibilities according to central market signals leads to a small improvement in terms of transmission grid congestion. A similar observation can be made for the inclusion of BESS in the adjustable generation on the grid side. The comparably small difference between the scenarios might be explained by the discrepancies in the transmission grid model and the generation and demand data for 2030, with structural congestion accounting for a large part of the observed overloading events. However, the resulting security of supply is largely increased by utilising distributed flexibilities with load shedding mostly rendered unnecessary, as can be seen in Table 2. In order to improve the accuracy of findings for the total amount, more grid expansion in the input data or a better coordination of renewable expansion and demand increase on the one hand and the transmission grid on the other hand might be needed.

AC vs. DC Formulation
The comparison of hourly results of the DynFlex scenario with the results of an identical simulation performed with the DC restrictions in the load flow equations is shown in Fig. 5. For the sake of an equivalent comparison, grid losses are subtracted from the positive dispatch adjustment of the AC solution, and the demand of the DC model was increased by today's average national transmission grid losses for all market areas considered. Otherwise, both model formulations, restrictions and input data were kept consistent. The DC problem was solved using the commercial solver Gurobi. Even though the set of DC constraints form a relaxed formulation of the problem, the AC solution has a lower requirement for adjustments after correcting the AC results for the required adjustments to account for grid losses. This is due to the additional positive generation increase needed to account for grid losses in the AC case, which are part of the optimisation problem in the formulation chosen in this paper and thus contribute to relieving existing congestions. As a result they lower the additionally needed adjustment after the subtraction of the losses, leading to lower requirements, even though the linearised DC formulation is commonly referred to as a relaxation of the AC formulation. In the section up to 5 GW, the effect of voltage limits can be seen, which leads to higher AC adjustments as voltage is constant in the DC formulation. Overall, individual deviation in single hours can be significant while the entire trend and structure are similar. The correlation between the positive adjustment value series is 91.1%. This shows that the DC approach is suitable as an estimator for the AC equations for the investigation of general trends, while individual results and peaks might differ significantly due to the simplifications undertaken. For these cases, the AC formulation leads to better results.

Conclusion
In this paper, we presented a two-step approach to investigate the influence of distributed flexibilities on the operation of the transmission grid in the central European area. For this, we have developed a market optimisation formulation, which enables us to obtain dispatch results which include different kind of distributed flexible generation and demand sources. Subsequently, we analysed the effect of four scenarios with different implementation rates of flexibilities in market and grid operation based on a case study using data for the year 2030. Furthermore, we compared the results of our multi-period AC formulation to determine grid congestion with a linearised DC formulation.
We conclude based on the required adjustments over the simulation horizon of 8760 h, that a scenario like the DG scenario of the TYNDP, with high RES-E increase and the additional introduction of demand-side flexibilities into the electricity system is generally feasible for the anticipated transmission grid of the year 2030. Locally concentrated, larger adjustments can be explained by the gap between known grid expansion in selected countries and the ambitious targets of the chosen scenario. Also, distributed flexibilities show to have a strong effect on solving grid congestions that lead to load shedding, which is necessary in the Base scenario of our study in more than 200 h in order to achieve a feasible solution. This is almost entirely resolved when optimising distributed flexibilities with the market and enabling the adjustment of BESS during grid operation.
Furthermore, the results of the four scenarios investigated show that the impact of increased participation of both central and distributed flexibilities in the market and grid operation has a positive but comparably low influence on the overall required adjustment to operate the grid within its technical limits in our chosen case study. The largest change on the average adjustments required occurs when adding central storage units to the set of generators and consumers available for adjustment in the grid simulation. The comparison of the AC and DC results showed that the required adjustment can be lower in the AC case, if the provision of grid losses is included into the AC formulation and generation increase is performed at locations which are suitable for lowering stress on the grid. Additionally, the individual results of each hour showed a generally well reproduced trend using the DC relaxation, however single hours can differentiate significantly in both cases.
In the current state a simplified approach of modelling the flexible operation of shiftable demands and renewables based on their profile was integrated in the market model. In the future we are planning to integrate the decentral flexibilities modelling into the grid model und validate it against a more detailed representation of individual units. As demand and supply uncertainty have a major impact on the utilisation of storages and load shifting potentials, we are furthermore planning to extend our modelling to a stochastic optimization. Finally, the extension of the current modelling approach to a unit commitment problem would allow us to fully evaluate the contribution of flexibilities for the future power system. 16 The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Introduction
Buildings in the commercial and residential sectors account for about 40% of the total energy consumption in European Union countries [1,2]. Heating, Ventilation and Air Conditioning (HVAC) systems represent one of the most significant loads contributing to electrical energy consumption. HVAC systems are particularly suited to be controlled by making the most of: (i) the thermal inertia of spaces to be conditioned (leading to some time dissociation between the energy consumption and the energy service provided), and (ii) the flexibility of consumers' preferences regarding the provision of energy services in face of dynamic time-of-use tariffs (displaying some willingness to bear the indoor temperature somewhat farther from the reference settings for a limited period of time). Since time-differentiated tariff schemes are expected to become relevant commercial offers in electricity retail markets in smart grids, demand response actions should be developed to achieve optimal operation of HVAC systems. Optimal HVAC control can be beneficial for consumers (by reducing the energy bill without jeopardizing comfort), retailers (by managing buying and selling prices) and grid operators (by contributing to release the strain in distribution networks and enhancing the utilization of renewable sources). Consumer engagement in demand response programs can be made operational by means of automated energy management systems endowed with optimization algorithms and the capability to control appliances. Therefore, adequate and tractable models for optimal thermostat programming should be developed as operational tools (to be embedded with sensors and control) for demand response programs.
The potential of thermostatic loads for demand response actions has been exploited in the literature, due to their significant consumption and the possibility of being controlled. For this purpose, mathematical models have been proposed, in particular mixed integer linear programming (MILP) formulations, i.e. with integer and continuous variables. The authors in [3] consider heat pumps for heating residential buildings with a floor heating system using a linear state space model in an Economic Model Predictive Control framework. Thermal models are developed in [4] for a smart house for determining the value of thermal inertia in demand response dynamic pricing using a MILP formulation. The authors in [5] present a user-centric demand response control for scheduling the electric space heating load under price and load uncertainty to minimize a weighted-sum of the expected payment, loss of comfort, and financial risk of a customer while considering the end-user's preferences. The household thermal behavior is modeled by means of a two-capacity building model. An HVAC system is considered in [6] that is controlled (in combination with other type of loads, PV generation and storage) by a home energy management system, thus enabling residential consumers to participate in demand response programs. A price-based demand response strategy for multi-zone office buildings to optimize the energy cost of HVAC units and the thermal discomfort of occupants is formulated in [7] as a MILP model. The authors in [8] develop an approach based on a partial-differential equation model of thermal diffusion to determine the thermostat settings to minimize the electricity bill for a consumer with energy time-of-use and power prices, in which the optimal thermostat programming for HVAC is formulated as a constrained dynamic optimization problem. The authors in [9] consider the optimization of the investment and operation planning of a decentralized energy system, which is subject to different sources of uncertainties, encompassing photovoltaic generators and load flexibility using heat pumps in combination with thermal energy storage units for space heating and domestic hot water, which is tackled by two-stage stochastic programming.
This paper presents novel comprehensive MILP models in which the physical characteristics of a thermostatically-controlled air conditioning system are considered. These models have the common aim of determining the optimal thermostat operation using different forms of thermostat control. It will be discussed that the combinatorial nature of these models makes it difficult to solve them by a state-of-the-art exact solver. This fact impairs embedding realistic mathematical models in automated building energy management systems when a fine-grain time discretization of the planning period is considered.
The study is made in the perspective of the building owner or tenant, i.e. the one who aims to minimize the energy bill. These models can also be useful to optimize demand response strategies that can be of interest for the grid operator (possibly through the mediation of an aggregator of end-users' flexibility that uses this asset to participate in ancillary service provision markets). Thus, this paper adopts an Operational Research methodology focusing on the development of mathematical programming models.
The paper begins by developing a simple thermodynamic model to derive a timediscretized equation expressing the instantaneous indoor temperature (at instant t i ) as a function of the indoor temperature, the external temperature and the air conditioning system operation at the preceding instant (t i−1 ). MILP models are then presented accounting for the hysteresis behavior of the thermostat between minimum and maximum temperatures defining the dead-band around a reference temperature (set-point). It is shown that, if just modelling this behavior, the model is solved in an acceptable computation time. However, if a minimum ON or OFF period is considered (also to avoid excessive switching), offering increased operation flexibility, the computation effort becomes impracticable.
In Sect. 2 the building thermal model is developed. MILP models of a thermostatically-controlled air conditioning system with different features are described in Sect. 3. Section 4 reports extensive computational experiments of the different MILP models. Conclusions are drawn in Sect. 5.

Development of the Building Thermal Model
A space, a building unit, is considered for being heated/cooled by means of an air conditioning (AC) system. Assuming the building unit as a (homogeneous) control volume, whose instantaneous indoor temperature is uniform and denoted by θ in (t), the overall thermal energy balance on a heat rate basis (at some generic instant t) is [10]: for AC in heating mode, and (2) for AC in cooling mode. In these expressions: q ext (t) is the instantaneous external load [kW]; q g (t) is the instantaneous internal load (rate of heat generation, by occupants, lighting, appliances) [kW]; C = ρc p V is the overall thermal capacity [kJ/ • C], in which ρ and c p are, respectively, the mass density [kg/m 3 ] and the specific heat at constant pressure [kJ/(kg. • C)] of indoor air, if no other thermal inertia of the building is considered, or some weighted values for indoors, and V is the volume of the building unit [m 3 ]; q AC (t) is the instantaneous air conditioning heating/cooling load [kW].
Considering a finite-difference approach over a time-step Δt to integrate equation (1) in a fully explicit time discretization -all variables at the "previous" instant t i−1 , except for temperature θ in (t i ) at the "present" instant of Δt -it can be assumed that the rate of the energy storage, i.e. the right-hand side in (1) is given by U is the (weighted-average) overall heat transfer coefficient of the building unit envelope [kW/(m 2 · • C)] and A is the surface area of the envelope [m 2 ]. Therefore, U.A represents the overall thermal conductance of the unit envelope [kW/ • C].
Introducing the expressions above in (1), equations (3), (4), (5) are obtained: In a dynamic simulation of buildings, the internal thermal load q g [kW] is usually defined according to the operation profile (e.g., daily/weekly profile of occupation, of lighting, etc.), i.e. all internal loads (thermal power) associated with the type of utilization of the building. Without loss of generality for the application of the model, and for the sake of simplification in this context, the following assumptions are considered: 1. the internal heat loads are neglected: where COP is the AC coefficient of performance, P demand and P AC nom are the power demand and the AC nominal power, respectively, and s t i−1 is the AC control variable  (5) can be rewritten as: Therefore, since COP , P AC nom , Δ t and C are always positive quantities, the sign in equation (10) is positive as this equation has been derived from (1) for heating mode. Otherwise, γ = − Δ t C COP (negative) for cooling mode. The coefficients α, β and γ derive from the building characteristics (area, envelope, etc.) and the COP of the AC. Considering a building unit with a 225 m 2 floor area and 3 m height, and the properties of the air inside (density of 1.19 kg/m 3 and specific heat of 1.007 kJ/(kg. • C)), the air thermal capacity is C ≈ 810 kJ/ • C. Taking U -values of 0.35 and 0.30 W/(m 2 . • C) for the building facades and roof, respectively, a value of U.A ≈ 0.129 kW/ • C is obtained. Assuming a COP = 2.5, typical values of parameters α, β and γ are displayed in Table 1 for different discretization intervals of the planning period.

Mathematical Models of a Thermostatically-Controlled AC System
This section presents different MILP models aimed at determining the optimal operation of the AC within a planning period to minimize energy costs, considering distinct forms of control. Although other types of loads could be considered (such as shiftable and interruptible appliances), the goal herein is to focus on the mathematical modelling of the AC operation. Shiftable loads include dishwashers, laundry machines and clothes driers, i.e. appliances for which the operation cycle cannot be interrupted once initiated. The supply of interruptible loads, which include electric water heaters and the battery of electric vehicles, can be interrupted within preferred time slots provided a certain amount of energy is supplied. In general, shiftable and interruptible loads require MILP models, which can be solved very fast by state-of-the-art solvers (e.g. Cplex), although requiring many binary variables. Examples of such models are the ones presented in [11] and [12]. However, the MILP modelling of the thermostat behavior is more complex and imposes a higher computational effort to obtain the optimal solution. In particular, MILP models are presented accounting for the hysteresis behavior of the thermostat between minimum and maximum temperatures defining the dead-band around a reference temperature (set-point) established by the user according to his comfort preferences.
The planning period consists of T time intervals t = 1, . . . , T of a given duration; this discretization can be, for instance, 15-, 5-or 1-min. The length of the time interval is denoted by Δt, in hours. For instance, if a 1-min discretization is used then Δt = 1/60 h. Figure 1 displays the behavior of the thermostat of an AC in heating mode. The hysteresis operation of the thermostat prevents excessive switching when the indoor temperature is around the set-point, which may be established as the midpoint within the dead-band defined by the comfort range of indoor temperature [θ min ,θ max ] specified by the user. These comfort bounds may depend on t, i.e. may be different during the planning period ([θ min t , θ max t ], t = 1, . . . , T ); for the sake of clarity, we will consider the comfort bounds constant throughout the planning period in the formulations presented below. Hysteresis determines the AC control variable: s t = 1, i.e. the AC continues ON as the indoor temperature θ in t increases above θ min until it reaches θ max (a in Fig. 1); s t = 0, i.e. the AC continues OFF as the indoor temperature decreases below θ max until it reaches θ min (b in Fig. 1).

Modelling the Thermostat Behavior
Model 1A forces s t = 1 when the indoor temperature (θ in t ): -is below the lower bound of the comfort band (θ in t < θ min ) OR Fig. 1 Behavior of an AC thermostat for heating mode -is within the comfort band (θ min ≤ θ in t ≤ θ max ) AND in the previous instant the AC was also ON (s t −1 = 1) Model 1A (heating mode): s.t.: The binary variables s t control the thermostat switching: The auxiliary binary variables y t and w t establish the consistency conditions associated with thermostat switching to the ON position (s t = 1): M is a positive large number, which is used in (13) and (14) to enforce the logical conditions. M must be larger than any temperature value; for instance, M = 100 (or any higher value) is suitable as temperatures are in • C. The constraints of Model 1A ensure that: 1. if θ in t < θ min , then y t = 1 and w t = 1 by (13) and (14); hence, by y t +w t −s t ≤ 1 (16), s t = 1; 2. if θ min ≤ θ in t < θ max and s t −1 = 1, then y t = 1 and w t = 1 by (13) and (15); hence, by y t + w t − s t ≤ 1 (16), s t = 1; 3. in other cases, s t is free to be 0 or 1. The objective function pushes these variables to 0, in order to turn the AC to the OFF state and minimize cost.
Additional input information to be specified include: -The temperatures defining the thermostat dead-band, i.e. the minimum temperature θ min ( • C) below which the thermostat should be ON (s t = 1) and the maximum temperature θ max ( • C) above which the thermostat should be OFF (s t = 0). -The initial indoor temperature θ in 0 ( • C) and the initial status (0, 1) of the thermostat s 0 . This model assumes that the "natural state" of the control variable s t is 0 (since it has a positive coefficient in the objective function to be minimized). The objective function forces the binary variables to 0 whenever the constraints do not impose these variables to be 1. However, in this model it may happen that, when indoor temperature drops below the upper bound temperature (θ max ) after having been above it (and hence s t = 0), the control variable switches from 0 to 1 within the thermostat comfort band (note that there are no constraints forcing s t = 0 because this would be the normal state of the variable). This may be beneficial for the objective function by avoiding switching on at some later time intervals when the electricity cost is higher. Therefore, Model 1A does not replicate exactly the thermostat hysteresis behavior. This model guarantees the maintenance of the ON state within the comfort dead-band (a in Fig. 1) but not the maintenance of the OFF state (b in Fig. 1). This freedom to switch ON to benefit the cost objective function is the reason why this model takes so much computation time. A model forcing the control variables to 0 or 1 according to the thermal balance equation and thermostat hysteresis, i.e. a rule-based system modelling the thermostat switching, is almost instantaneously solved. This is implemented in Model 1B. In addition to Model 1A, Model 1B (also for heating mode) explicitly forces s t = 0 when the indoor temperature (θ in t ): -is above the upper bound of the comfort band (θ in t > θ max ) OR -is within the comfort band (θ min ≤ θ in t ≤ θ max ) AND in the previous instant the AC was also OFF (s t −1 = 0) (b in Fig. 1) Model 1B (heating mode): s t variables are rigidly determined by the conditions that establish θ in t and the thermostat operation (Fig. 1). In this case, there is a single feasible solution (the one that complies with the thermostat hysteresis switching rules), so obtaining the solution to this model is almost instantaneous.

Modelling Minimum ON/OFF Duration Periods
A new model (Model 2) has been developed to offer the possibility of controlling the AC to make the most of time-differentiated tariffs, i.e. being ON during lower electricity price periods when it is not strictly necessary to heat the building space so that thermal comfort is achieved even limiting the time ON when electricity prices are higher. In face of this freedom given to the model, it is necessary to inhibit excessive switching, which is accomplished by imposing a new set of constraints. In the comfort band, the AC control variables s t are determined to minimize costs in face of dynamic tariffs, but once there is a switch from 1 to 0 or from 0 to 1, then the new state ON/OFF should be kept for at least d time intervals (these minimum ON/OFF periods could be different for each state). These new constraints establish: This model is also quite demanding from the computation time point of view due to offering the possibility of switching from 1 to 0 or from 0 to 1 in the comfort band. Model 2 (heating mode): s.t.:

Modelling an Inverter AC
Inverter technology AC appliances have growing acceptance because of increased efficiency, extended life and elimination of abrupt load and temperature variations. In this type of appliances, a variable-frequency drive adjusts the compressor and the cooling/heating output. The previous models can be extended to accommodate inverter units, for instance considering they can be operated at different levels with respect to the nominal power. This can be modelled by introducing additional binary decision variables: For instance, for R = 5 with the power levels 20-40-60-80-100% of the nominal power, the AC operation is modelled in Model 3. In this model, θ min and θ max should be interpreted as absolute bounds the user is willing to endure. Note that the sets of constraints in the previous models could also have been used here (for simplification purposes they were omitted). p AC t is the power required to operate the AC at time t of the planning period, i.e. the AC is either OFF (p AC t = 0) or supplied at one of the specified power levels.

Model 3 (heating mode):
be combined with the thermostat modelling of Models 1 and 2 and/or with the cost vs. comfort trade-off.
In the next section, illustrative results of the exact solution (whenever possible with a given computation budget) of these models using a MIP solver are presented drawing attention to the computational effort.

Data
A planning period of 24 h with a time discretization of 1 min is adopted, i.e. 1440 time steps, which imposes a heavier computational burden. The thermostat deadband is defined by θ min = 20 • C and θ max = 24 • C. The initial indoor temperature θ in 0 = 18 • C in all models, except in model 3 in which it is θ in 0 = 20 • C. The initial thermostat status s 0 = 0. The external temperatures θ ext t ( • C) for a winter day (from a meteorological station in Coimbra, central Portugal, on January 1 st , 2012) are available at http://www.uc.pt/en/org/inescc/publications/files/temp_out_ winter.xlsx.
The nominal power of the AC, P AC nom = 1.5 kW. The time-differentiated energy prices c t (e/kWh), considering 10 sub-periods, are given in Table 2. Experiments have been made on an Intel(R) Core(TM) i5, 3.33 GHz computer with GAMS ver. 24.0.2 and Cplex ver. 12.5.0.0.

Results
A computational budget of 7200 s or relative MIP gap of 0.5% was established. The behavior of thermostat switching as well as the evolution of indoor temperature are displayed in Figs. 2 and 3 for Model 1A and Model 1B, respectively. The cost of the corresponding optimal solutions for Model 1A and Model 1B, as well as the dimension of the models and the Cplex time to obtain these solutions are given in Table 3 Table 4. Considering that the AC can operate at power levels 20-40-60-80-100% of P AC nom in Model 3, the cost of the optimal solution as well as the dimension of the model and the Cplex time to obtain these solutions are given in Table 5. The corresponding solution of the AC operation and indoor temperature variation is displayed in Fig. 6.  The frequency of switching in Model 1A is higher than in Model 1B due to the freedom given by that model to minimize cost. In Model 1A, the AC is ON for a longer time in pricing sub-period P 3 to account for the increase in price in subsequent sub-periods (Fig. 2). Note that the solution presented in Fig. 2 and Table 3 for Model 1A still presents a MIP gap of 43.47% after 2 h of computation, due to the combinatorial difficulties of this model. Model 1B, as expected, is solved to optimality instantly.
The temperature peaks in Figs. 4, 5 and 6 are due to the capability of the models to accommodate for comfort by making the most of lower price sub-periods. The MIP gaps in Model 2 (Table 4) are still significant after 2 h of computation.
The optimal solution to Model 3 is better than the solutions obtained for the other models due to the possibility of the inverter AC to work at different power levels.

Conclusion
This paper presented a set of mixed integer linear programming models in order to optimize the operation of a thermostatically-controlled air conditioning system (AC), aiming at minimizing the energy costs. These models capture different physical control characteristics: -Model 1A models the hysteresis operation of the thermostat (in heating mode), preventing from excessive switching when the indoor temperature is around a setpoint. This model gives, however, some freedom to thermostat switching when the indoor temperature is within the dead-band, enabling to change from OFF to ON if this benefits the cost objective function. -Model 1B removes the flexibility given by the previous model, forcing the AC to keep its state when the indoor temperature is within the dead-band. This model implements, in a MILP formulation, a rule-based system of logical implications.
The feasible region is reduced to one solution and the model can be solved very quickly, unlike Model 1A which is computationally very demanding. -Model 2 gives some flexibility while inhibiting excessive switching, by considering that, when the indoor temperature is within the dead-band, the thermostat must keep its state (ON or OFF) for at least a predefined number of time intervals. This model is also quite computationally demanding. -Model 3 introduces a new feature that can be used to extend the previous models: it models inverter units, considering that the AC can operate at different levels of nominal power.
Modelling cost reduction at the expense of accepting worsening the indoor temperature by some degrees was also introduced in this paper. It has been shown that the combinatorial nature of these models imposes a severe computation burden, which in some cases impairs obtaining the optimal solutions in a practical, acceptable computational time by a state-of-the-art solver. Therefore, this work lays the foundations for understanding those computational difficulties and gives insights for the development of other approaches, namely dedicated meta-heuristics customized for the features of the models, in order to offer good solutions with a suitable computational effort. The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Introduction
Renewable energies are a resource that strain the grid through intermittent availability and difficulty in prediction. Demand side management addresses the issue by reversing the paradigm of grid operation and controlling power consumers instead of only power generators.
Within demand side management there is a need for robust algorithms that face the unpredictability aspect of renewable energies, which makes real-time algorithms with no forecast information favorable.
Finding solutions that improve the residual power problem is only the first step. Said problem is the question of how unused renewable supply should be handled. Using the Greedy algorithm significant improvements can be accomplished. However, the results indicate that distress would be unevenly spread in the consumer population. If distress is distributed unevenly, compensation would also be distributed unevenly, which means that the market is not design towards fairness. This paper introduces Weighted Fair Queuing (WFQ) as a scheduling algorithm for deferrable loads in smart grids. Three variant implementations are presented that feature different concepts of fairness: serving more and simultaneously smaller loads, fairly distributing wait times, and treating any number of groups equally.

Model
The experimental environment is designed in AnyLogic and uses an interface to control power consuming processes, as presented in [1]. The interface is suitable for any process that is capable of prolonging periods of its activity or inactivity as can be seen in Fig. 1. This requires no intimate knowledge of the process: when an algorithmic decision is made to alter the behavior of an individual load, the request is passed on to the agent representing the load. In accordance with its internal conditions it then accepts or rejects. This is to replicate the sovereignty of loads, as in actuality algorithms are informed of the current states and consider them accordingly.
Every process is divided into 4 states: activity and inactivity, and for each of those one deferrable and one non-deferrable. Processes can alter their power demand from step to step or follow a mathematical function-anything that describes the real behavior. The model captures distinct periods in which switching can be deferred.  The deferral and operation times are a matter of survey in the companies (interviews, data sheets and surveillance). The total load shape is based on real loads that are replicated with the simulation. Figures 2 and 3 are two examples of the power demand and the repeated activation of loads.
The analyzed algorithms make no use of forecasting of any kind; they operate based on information of the "now". The power supply is processed data-point-bydata-point, left to right.

Population
The model population consists of approximately 300 individual loads. These encompass beverage production, metal working, glass finishing, plastics production, textile treatment, sewage treatment, electric vehicles and compressors for both pressure applications and freezing. The switching behavior of a load for the purposes of the simulation is detailed. The realistic accuracy was verified by a test implementation in [2]. The composition of the population, i.e. the proportions of types of loads are elected to reflect a midsized city by approximation, based on statistics found in [3,4] and [5]. The overall peak power demand across the whole population is 25 MW, with single members requiring power between 1 and 350 kW. The average population power demand is 5.8 kW (cf. Fig. 7). State changes are deferrable by 15% to 50% and depend on the current and future state (arrows in Fig. 1) of the individual load. At any given time an average of approximately 20% of loads is deferrable.

Objective Function
Considering the negative impact of renewable energies, optimizing the load population toward a renewable energy supply profile seems adequate. As the volatility of the renewables is the straining property, immediate consumption by the deferrable loads would lead to the minimization of the residual power function and thereby stabilization. Residual power (R) is the amount of power leftover from renewable supply after subtracting the demand. A sample from historic supply data is picked, that consists of wind and solar energy over 24 h with a resolution of one data point per 3 s. This deviates from the resolution of simulated demand, which is continuous. Figure 4 shows the sample. At this resolution the objective function incorporates no useful gradient information.
However, selecting a random supply function is not sufficient: First the data points of the supply function must be scaled according to the demand of the population of deferrable loads. The first simulation experiment is a simple observation of how much power the consumer population requires when there is no optimization at all. This is the so called "natural run" that allows any and all requests to switch on or off. This unimpeded 24 h simulation returns a power demand function (P (t) Demand ) In Eq. (2) the integrals of the natural run and the supply functions which are W Demand and W Supply are equated. This yields a result for c, the scaling factor. The result of this approach is a function of power supply (P (t) Supply ) that can satisfy the demand in terms of energy (kWh), but in terms of the progression of power over time (kW) any algorithm must manage transient deficits.

W Demand
This scenario seems manageable concerning not only algorithms, but is also realistic with regard to resource allocation. The latter meaning that purely renewable supply seems unrealistic and would probably feature a safety factor to account for fluctuations in demand. At the same time, however, in an economic sense a rough balance between supply and demand is sensible, because of investment and operative costs. At the least, equal supply and demand appears to be a good initial metric.
To facilitate the comparison the summation from equation (3) is used, approximating the integral of the residual power for discrete time steps: This condenses the result of one 24 h simulation into one manageable value.

Greedy
The Greedy strategy (Algorithm 1) is a heuristic search algorithm. It is incomplete, but fast and very simple. As [6] shows it generally yields good results, though it is not capable of systematically finding global optima. I chose this algorithm as a starting point for its popularity, non-specificity and low computational requirements. I have previously outlined the use of the Greedy algorithm in [1] and [2]. The algorithm sorts the available candidates according to size and attempts to activate as many as possible, as can be seen in Algorithm 1.

Weighted Fair Queuing
Weighted Fair Queuing (WFQ) is a sophisticated algorithm designed by Demers et al. (cf. [7]), building on the work of Nagle [8]. WFQ manages congestion in gateway queuing of datagram networks. The problem occurs when participants in a network attempt to send data packets, but some entities always send larger packets than others. WFQ replaces the pragmatic FIFO (first in, first out) paradigm by categorically distinguishing datagram sources by packet size and permanently assigning them to queues with predetermined weights. Weight, packet size and a queues history determine which packet gets serviced next. Identifiable parallels to the scheduling of deferrable loads lead to the application of WFQ in a novel way. The deferrable loads also feature a well-defined requirement for resources which is their power demand in Watts. This parameter is equatable to the packet size, as is the network bandwidth to the usable residual power. The significant reason that makes a more sophisticated algorithm desirable is the unfair dissemination of resource with the Greedy strategy. By sorting according to size, Greedy categorically prefers large loads, which disadvantages smaller loads. deactivate; 8 e n d 9 end 10 end WFQ (2) employs a heuristic called virtual finish time 1 F (Eq. 4), which is calculated using the packet size, or in this case the wattage of a load l and the weight w of each queue i. Every load is usually permanently assigned to a queue, as are weights to queues.
This yields a numeric value for each load that wants to schedule for activation. For each queue a so-called session virtual finish time SF i is calculated (Eq. 5).
Within every queue, all elements are queued on a FIFO basis. The element of the queue i to be scheduled next is F i . All previously successfully activated elements are added with the summation term F n . At every point in time the algorithm selects the queue with the lowest SF i value. The significant difference in comparison to the original WFQ implementation is the explicit check (Algorithm 2, ln. 10) whether the current residual power is sufficient for the regarded load. In the network realm decreasing bandwidth will not cause a packet to be rejected service-it will, however, take increasing amounts of times to submit. This "flow" property is not given in the case of deferrable loads, therefore "choking" cannot be applied, and switching decisions are discrete.

Fairness
In datagram networks fairness is the property of allotting a certain resource to any participant according to a predefined key. Such a key maybe simple and assign equal fractions to everyone. In the case of WFQ, however the distribution scheme accounts for a priori observations in terms of the typical resource requirements of each participant. With WFQ bandwidth cannot be manipulated which means that the amount of sheer data transmitted will at best remain the same (though more packets cause more gaps). In fact the objective is to service a greater number of participants.

Fairness of Wattage
This exact thought is transferred to demand side management applied to a load's wattage. In contrast to the original network application participants do not forfeit after multiple trials in vain. It is implausible to assume that a load be denied access altogether. Loads forcing themselves to power on is unfavorable, but in line with the conventional grid operation paradigm of supplementation in such a case. As the key influenceable property is the length of part of the "on" and "off" episodes of each load, disproportionately large loads compensate with their "on" time. This means that larger loads are operated for shorter periods of time by the algorithm. The objective of employing this queuing scheme is to serve more participants in the same time interval. Here fairness is targeting decongestion.

Fairness of Distress
"Wait time" is a second, additional fairness measure worthy of investigation. Wattage is a strictly objective parameter. If we consider the period that a process can wait to be served, this requires very intimate knowledge of the process to determine the validity or honesty of time parameters as declared by a machine operator. A company might easily misrepresent this information. For the simulation experiment at hand the processes were personally surveyed and documented during live observations. Especially as there was no motivation presented to falsify data, the time parameters can considered to be sufficiently accurate.
In this regard fairness would disregard the wattage of a load and consider the time it has been waiting to get served. Fairness in this sense is the equal distribution of waiting-distress among the entirety of the load population.
For this variant fixed queue assignments are replaced with dynamic ones. The longer a process waits, the further it moves into queues that make it more likely to be picked next. To this end standardized wait time is favorable, i.e. the fraction of time waited in proportion to the maximum amount of time the load can wait (Eq. 6). Without this measure loads that provide less flexibility are irrationally rewarded.
Fairness Group-by-Group As for a third measure of fairness the load population is subdivided into regions of equal size, representative of quarters or districts. In this configuration each queue represents one region, and regions compete for the available resources. This option is designed to grant each group the same right to resources. Any attribute can be used to queue in this manner. This generic queuing scheme targets an equal distribution of deferral times and power to independent groups, but not within the groups. As the virtual finish time paradigm is still applied, smaller loads are preferred.

Weights and Queues
There is no definitive answer on how weights should be chosen. Existing procedures such as that outlined in [9] are not applicable because the prevailing problem of node congestion in networks is not a relevant phenomenon. Therefore a mode of selecting weights should rely on the available parameters. Each of the fairness measures suggested above requires its own weighing scheme.

Weight by Power (WFQ-Power)
Considering that the resource "power" is engrossed correlating with wattage, a load's wattage should inversely determine its likelihood of being elected. As a number of queues to which all loads are assigned, 5 was elected as it is typical. The fifth of the population with the highest wattage is assigned to the according quintile, and so forth. The weight for each queue is determined by computing its average wattage. Following equation (4), fairness is established when the virtual finish times of all individual loads are equal. This means solving for w and equating the average wattages of all queues. Table 1 shows which weights are assigned to queues, and the respective average wattage within the queue. For example: At 250 kW load # 1 is a member of the queue weighted at 0.09, i.e. the group of highest wattages. Compared to bottom quintile member load # 2 of 20 kW and its weight of 1. Load # 1 has a virtual finish time 138 times higher. The mismatch in wattage reflects proportionately in the weights.

Weight by Wait Time (WFQ-Time)
Queuing by standardized wait time is carried out dynamically. In contrast to wattage-queuing there is no a priori assignment. With increasing wait time, a load is assigned to queues with higher weights that equate to higher likelihoods of getting picked. Once again there are five queues, each designated for a 20% interval of wait time completed. Table 2 shows how weights are governed for each queue. Decreasing weights augment urgency.

By Region (WFQ-Geo)
In the case of grouping by region equal weights are assigned to all queues. As the fictitious regions comprise equal portions of every category of the overall population, and given the equal weights, this queuing scheme acts comparable to a round-robin. This is true for all types of WFQ with equal weights (cf. [8]).

Model Limitations
The proposed model requires highly automated processes. The underlying work plans are machine accessible so that the ensuing steps and their duration can be projected. Semi-automation or increased human interaction would hinder this procedure. Even though they were not encountered during survey, processes that can continuously adapt the level of their power input can only be replicated as a series of discrete steps. On the contrary this model is a response to the lack of continuously adaptable processes. Yet this is a limitation.
In addition the size of the consumer population cannot be arbitrarily increased. To this end an arbitration mechanism appears feasible. Multiple optimization algorithms would work on a fraction of a larger problem which in turn is consolidated by optimization layers higher up in the hierarchy. This approach could also help to reduce the amount of datagrams that are sent between the loads and the optimizer.

Results
In each simulation there is a population of approximately 300 manageable loads that must be scheduled over 24 h following a randomly selected supply signal. There are five experiments: 1. Natural run: All loads operate without any algorithmic interference. This is a reference trial. The algorithm in experiment 3, with weights defined by power demand, aims to serve more participants, especially smaller ones. Experiment 4 entails queuing by wait time, which is designed to equally distribute the wait time percentages over the entire population. Queuing by region, experiment 5, is designed to allocate the same resources to arbitrarily defined subsets of the population.
In order to judge the fulfillment of these secondary objectives additional parameters must be analyzed. Cycles counts all activations which is the sum of every individual power-on that was granted-in the datagram sense this is analogous to serviced packages. The standardized average wait time is indicative of how long every activation was deferred in percent of the maximum. Deferral length is the sum of all deferrals as time. The variance of this value is representative of how evenly deferrals are distributed among the population. Table 3 summarizes the results of experiments 1-4. As outlined in Sect. 2.3 the significant data is the power consumption function that every algorithm causes. Figure 5 is the consumption profile caused by the Greedy strategy. Consumption curves for all algorithms including the natural run (experiment 1) can be found in the appendix (Figs. 7,8,9,10,and 11). The visual differences in the residual curves are exiguous due to the algorithmic similarity. The key difference is the residual energy. Systematic differences do not express discernibly in the graphs. Showing which load   is activated at any given time offers no insights from algorithm to algorithm, as no pattern emerges. The Greedy algorithm performs best by reducing the residual by 54.8%. WFQ-Time improves the deferral distribution by a factor of 2.08, while maintaining a very similar deferral length. Table 4 shows the results of all queues of experiment 5, each representing one region. As this algorithm is specified not to improve the global population, but to serve queues equally, its results cannot be directly compared to the other experiments. The only exception is the residual. Queue deferral length is on average 11% higher compared to WFQ-Power. Queue by queue comparison shows that results group close together for all parameters, but worse in general.  Fig. 6 The WFQ algorithms can compete with the versatile Greedy algorithm in terms of the residual optimization (blue). However Greedy focuses 80% of distress on only 69% of the population (green). In addition WFQ-Power rejects especially few loads (red) Figure 6 summarizes the advantages of WFQ: it is an improvement on the Greedy algorithm which is fast and versatile. At the cost of slightly worse results with regard to the optimization of the power residual, WFQ algorithms distribute distress more equally, and serve more loads all else being equal.

Improvement
As to be expected the natural run that features no algorithm or adaptation to the objective function causes a high residual. The Greedy algorithm delivers the greatest improvement in terms of residual. At the same time this algorithm unevenly distributes the wait time distress across the population which is expressed in the variance of the average deferral time. WFQ-Time significantly narrows the variance of deferral, following its secondary objective. The WFQ-Power variant is capable of placing more individual activations which can be seen in its cycle count. Keeping in mind that the objective function cannot be fully attained, as indicated in Sect. 2.3, the algorithm still reduces the gap in serviced packages by approximately 60% in comparison to Greedy.
The results from the algorithm WFQ-Geo in experiment 5 are inferior in general as it causes longer wait times and the variance thereof indicates that wait times are heterogeneously distributed within queues. The objective of equal queue treatment is achieved at the cost of worse results all in all.

Distress
Distress distribution is a desirable trait in demand side management that is not trivially deducible. Wait time is a parameter that cannot be easily subverted, if for example the maximum deferral time is repaid as an incentive. Weighted Fair Queuing is capable of delivering comparable result to a Greedy strategy in the primary objective of residual power improvement. In addition it can be adapted to pursue additional objectives. I presented three possible measures of fairness: equality between multiple queues or regions (WFQ-Geo), the even distribution of wait times (WFQ-Time), and serving more participants (WFQ-Power)-the latter being most similar to the original algorithm. This shows that WFQ can be used to introduce fairness to the selection of deferrable loads.

Methodological Contribution & Policy Advice
Despite its complexity Weighted Fair Queuing is a staple of network congestion management. It is highly functioning and ubiquitous because of its performance and its ability to accomplish fairness. The application to energy distribution problems is novel. At the same time it introduces the problem of fairness into the matter, which was disregarded because of the focus on residual power optimization, monetary compensation or welfare gain.
The optimization problem presented and solved in this paper is significant in its size regarding multiple dimensions. Firstly the time resolution exceeds the commonly quoted 15 min intervals that stem from the tertiary grid balancing realm. Grid stabilization however can only be achieved when the shrinking number of spinning masses can be replaced. The realization that compensatory measure must move below the 30-seconds-threshold is at the core of this research. Secondly, the population is substantial with 300 members and 25 MW peak. Each load is replicated as an agent with internal processes, variables and decision making. Furthermore their states are not estimated or stochastically approximated-decision making and transmission is acute, meaning that statuses are inquired, evaluated, decided on and requests dispatched. Thirdly, no simplifying assumptions are applied to the target function (power supply) such as smoothing. The complexity of the problem is embraced to the extent that the resolution of the objective function (supply) is the restricting factor. The argument is in favor of preparedness, as methodologically the question answered here is that demand side management features the necessary scalability and agility.
The suggested algorithms are implementations of centralized algorithms. As the technology is still at an early stage centralization is acceptable. Especially as the methodological exploration of this topic-which this research is part of-is still ongoing. For future applications, however, this is not advisable, as a single entity will possess all information on all clients (loads). The controlling entity would be vulnerable to attacks . The presented model, is acting on a request basis which means that a load can reject any suggestion to change its state. Designing intrinsically safe systems increases complexity, but is quintessential policy advice.

Practical Implications
For companies owning deferrable loads the practical implications are low, especially if processes are highly automated. The test implementation presented in [1] required no manual interaction. While loads were in interference the optimization paused and resumed automatically. The more a deferrable load is interwoven into the process of a company, the more difficulty can arise from automated unexpected deferrals. Deferrable loads are an effort that aim to improve the usage of renewable loads. Some resulting distress is to be expected. Although not insinuated here, market design should compensate for this.
The suggested interface is at the core of the research presented, as it allows for the load to negotiate deferral and reject it. From a load's perspective WFQ-Time is the most favorable, as it guarantees that no load is overburdened with deferral. The Greedy algorithm is unfavorable, as some loads are asked to defer significantly more often than others (unfairness).

Algorithm Comparison
There is a multitude of optimization algorithms that can be used. The advantage of the above presented algorithms is their simplicity. They can be categorized as local search algorithms, which means that they use a type of heuristic, like SF i in Eq. (5). As there are no simplifying assumptions, and every load is designed to be addressed with a direct request, the problem presents itself to be massive and without useful gradient information, which excludes classic optimization algorithms.
A suitable alternative would for example be "Simulated Annealing" (cf. [10]). It is an algorithm that begins by amply exploring the solution space with random solutions. A solution in this case is an allowable sequence of activation and pause for all loads. The algorithm continues by rejecting or accepting results, based on the metaheuristic property ("temperature"). In principle this algorithm converges on the global optimum solution, which local search algorithms are not systematically capable of. The significant disadvantage is considerably higher run time. A preliminary test implementation shows that, after 14 h Simulated Annealing delivers results, which are comparable to the Greedy strategy which takes approximately 15 min on the same desktop computer (2.5 GHz Intel Core i7 with 16 GB 1600 MHz DDR3 memory).

Related Work
Gellings laid out the fundamentals of demand side management in [11].
The objective function and its definition are based on [12] and [13].
Weighted Fair Queuing and the approach of Generalized Processor Sharing are based on [7,8] and [14]. These original sources were used for the adaptation to demand side management.
Most related work is based on [15] where an optimal deferrable load control problem is defined. The focus is placed on market parameters which is why price bounds are employed. [16] focuses on the problem laid out in this publication. To solve this problem a decentralized algorithm is employed, which communicates the power residual to all participants. This is in contrast to centralized solution efforts, but capable of scheduling multiple instances of the same load. [17] utilizes a similar approach but only for singular placement of loads in 24 h.
[18] uses a stochastic model which defers one activation of each load and no reactivation. The underlying idea is very strongly connected to the paradigm of grid operation. This means that the time scales of balancing power in the European grid are adopted and all activity takes places in time steps of 15 min, and deferral times are at least 4 h. The significant difference is the lack of understanding and anticipating load behavior. The publication focuses on the simulation of wind power only and the cost of demand side management.
[19] has a similar approach in that only one activation instance of every participant is considered. The multitude of loads is approached as a combinatorial problem. Here the smallest time step is one hour. In neither of these loads can extend their operative time.
Formulated as a constraint problem regarding the availability of the demand side, [20] manages a day-ahead scenario by interacting with energy markets. By the introduction of welfare they move away from monetary quantification. They outline the difference in time scale between markets (larger steps) and the load behavior (smaller steps).

Introduction
A ZEN is a neighborhood that has a net zero emission of CO 2 over its lifetime. Many aspects are embedded in the idea of ZEN. Energy efficiency, materials, users behavior, energy system integration are all aspects that need to be accounted for in this concept. In addition, different parts of the life cycle can be included but in this paper we only consider the operation phase and no embedded emissions. Two types of action exist to make neighborhoods more sustainable. One is to act on the demand, via better insulation, user behavior or other efficiency measures. The other is to act on the supply and have a local energy system minimizing the CO 2 emissions. There is consequently a need for a way of designing the energy system of such neighborhoods. The questions to be answered are, which technologies are needed to satisfy the demand of heat and electricity of a neighborhood, and how much of it should be installed so that it is as inexpensive as possible. The problem is then to minimize the cost of investment and operation in the energy system of a neighborhood so that it fulfills the ZEN criteria. This paper presents an optimization model to solve such problems with a focus on operations research methodology.

State of the Art and Contribution
The ZEN concept is specific to this particular project, however similar topics have been studied in different settings either at the neighborhood level, the city level or the building level, for example during the research center on Zero Emission Building. In this context, K B Lindberg studied the investment in Zero Carbon Buildings [2] and Zero Energy Buildings [3] which are variations around the concept of Zero Emission Buildings. In both papers an optimization based approach is used to study the impact of different constraints on the resulting design. The second one [3] in particular uses binary variables to have a more realistic representation of the operation part (part load limitation and import/export). In [4], Gabrielli et al. tackle the problem of investment and operation of a neighborhood system and show an approach allowing to model the system complexity while keeping a low number of binary variables. It also constrains the total CO 2 emissions. It uses design days and proposes two methods for allowing to model seasonal storages while keeping the model complexity and reducing the run time. In [5], Hawkes and Leach look at the design and unit commitment of generators and storage in a microgrid context using 12 representative days per season in a linear program. It is particular in that it defines how much the microgrid would be required to operate islanded from the main grid and include this in the optimization. It also discusses the problematic of market models within microgrids. In [6], Weber and Shah present a mixed integer linear programming tool to invest and operate a district with a focus on cost, carbon emission and resilience of supply. A specificity of this tool is that it also designs the layout of the heat distribution network taking into account the needs of the buildings and the layout of each area. It uses the example of a town in the United Kingdom for its case study. In [7], Mehleri et al. study the optimal design of distributed energy generation in the case of small neighborhoods and test the proposed solution on a Greek case. Emphasis is put on the different layouts of the decentralized heating network. In [8], Schwarz et al. present a model to optimize the investment and the energy system of a residential quarter, using a two stage stochastic MILP. It emphasizes on how it tackles the stochasticity of the problem in the different stages, from raw data to the input of the optimization, and on the computational performance and scalability of the proposed method. In [9], he also studies the impact of different grid tariffs on the design of the system and on the self-consumption of the PV production. In [10], Li et al. separate the investment and the operation into a master and a follower problem. The master problem uses a genetic algorithm to find the optimal investment while a MILP is used to find the operation in the follower problem. In [11], Wang et al. also use a genetic algorithm, but at the building level and using a multi objective approach focused on environmental considerations. A life cycle analysis methodology as well as exergy consumption are used to assess the design alternatives. In [12], Mashayekh et al. uses a MILP for sizing and placement of distributed generation using a MILP approach including linearized AC-power flow equations. In [13], Yang et al. also use a MILP approach for the placement and sizing problem but consider discrete investment in technologies at the district scale. These papers give us an overview of different methods for optimal investment in the energy system of neighborhoods or buildings, but none apply the ZEN concept and the influence of tight requirements on the CO 2 emissions on the modelling and on the results has not been demonstrated.
In this paper, the focus is put on getting a fast yet precise solution that can take long term trends, such as cost reduction of technologies or climate. To this end, the proposed model uses a full year representation, ensuring a correct representation of seasonal storage of heat and electricity, and allows to divide the lifetime of the neighborhood into several periods, each represented by one year. It is also different by using the Zero Emission framework on a neighborhood level as a guide for the emission reduction constraint. This adds an integral constraints coupling each timestep and increasing the complexity of the problem. The use of binary variables is limited to the minimum.

ZENIT Model Description
ZENIT stands for Zero Emission Neighborhoods Investment Tool. It is a linear optimization program written in Python and using Gurobi as a solver. It minimizes the cost of investing and operating the energy system of a ZEN using periods, with a representative year in each period. Different technologies are available, both for heat and for electricity. It is most suited for greenfield investment planning but can also take into account an existing energy system. The objective function is presented below: The objective is to minimize the cost of investing in the energy system as well as its operation cost.
The operation phase can be separated in different periods during the lifetime of the neighborhood, and one year with hourly time-steps is used for each period. In addition to technologies producing heat or electricity, there is also the possibility to invest in a heating grid represented by the binary b hg that also gives access to another set of technologies that would be inappropriate at the building level. In the equation above, the E represent discount factors either global for the whole study (3) or for each period (2). They are calculated in the following way: The calculation assumes that reinvestment in this technology is made for the whole lifetime of the neighborhood, and is discounted to year 0. The salvage value is also accounted for. The formula used is: In the objective function, y exp t,p represent the total export from the neighborhood. It is simply the sum of all exports from the neighborhood: The most important constraint, and what makes the specificity of the "Zero Emission" concept, is the CO 2 balance constraint. It is a net zero emission constraint of CO 2 over a year. It takes into account the emissions from the used fuels and electricity with the corresponding CO 2 factors for the emission part and the exports of electricity for the compensation part. In this study the same factor is used for imports and for exports of electricity. This constraint is expressed below, ∀p: In the particular ZEN framework of this study, the idea behind the compensation is that the electricity exported to the national grid from on-site renewable sources allows to reduce the national production, and thus to prevent some emissions from happening. The corresponding savings, the compensation, stand on the right-hand side of the equation. In the ZEN framework, this constraint is set as an annual constraint. It can however also be used for shorter periods of time.
Other necessary constraints are the different electricity and heat balances which guarantee that the different loads are served at all times. The electricity balance is represented graphically in Fig. 1. The corresponding equations are also written below. The electricity balance is particular because, we want to keep track of the origin of the electricity sent to the battery. It is managed by representing each battery as a combination of two other batteries: one is linked to the on-site production technologies, while the other is connected to the grid. It allows to keep track of the self-consumption and to differentiate between the origin of the energy for the CO 2 balance. Node I (8) represents the main electric balance equation while II (9) and V (10) are only related to the on-site production of electricity. Node II (9) describes that the electricity produced on-site is either sold to the grid, used directly or stored, while node V (10) states that at a given time step what is stored in the batteries is equal to what is in excess from the on-site production. Electricity Electricity balance V: ∀t, p g g ch t,p,g = est y pb_ch t,p,est Heat also has its own balance, that guarantees that the demand of each building is met: Note that the demand is not divided between domestic hot water (DHW) and space heating (SH). The batteries are represented, as mentioned earlier and as seen on Fig. 1, as two entities: one on the on-site production side and the other on the grid side. This means that we have two "virtual" batteries with their own set of constraints as well as constraints linking the two.
The first constraint is a "reservoir" type of constraint and it represents the energy stored in the battery at each time-step: ∀t ∈ T * , p, est Equations (14), (16) and (17) link both batteries. They make sure the sum of the stored energy in the "virtual" batteries is less than the installed capacity, and making sure the rate of charge and discharge of the battery is not violated.
The storage level at the beginning and the end of the periods should be equal.
The heat storage technologies also have the same kind of equations as the batteries, for example: ∀t ∈ T * , p, hst Equations (14) to (18) In order to not add additional variables, the mutual exclusivity of import and export is not explicitly stated. It is still met however due to the price difference associated with importing and exporting electricity.
In addition to the above equations, different constraints are used to represent the different technologies included. The maximum investment possible is limited for each technology. ∀i: The amount of heat or electricity produced is also limited by the installed capacity: ∀q, t, p : q q,t,p ≤ x q (22) ∀g, t, p : g g,t,p ≤ x g The amount of fuel used depends on the amount of energy provided and on the efficiency of the technology: respectively ∀γ ∈ F ∩ Q, p, t and ∀γ ∈ E ∩ Q, p, t For CHPs technologies, the Heat to Power ratio is used to set the production of electricity based on the production of heat. ∀t, p g CH P ,t,p = q CH P ,t,p α CH P For the heat pumps, the electricity consumption is based on the coefficient of performance (COP).
∀hp, b, t, p The heat pumps are treated differently from the other technologies because they are not aggregated for the whole neighborhood but are separated for each building. This is because the COP depends on the temperature to supply, which is different in passive buildings and in older buildings and which is also different for DHW and for SH, and dependent on the temperature of the source. The source is either the ground or the ambient air depending on the type of heat pump. The COP is then calculated using a second order polynomial regression of manufacturers data [3] and the temperature of the source and of the outside timeseries. The possibility to invest in insulation to reduce the demand and improve the COP of heat pumps is not considered. The global COP is calculated as the weighted average of the COP for DHW and SH.
The solar technologies, solar thermal and PV, also have their own set of specific constraints. ∀t, p: The hourly efficiency of the PV system is calculated based on [14], and accounts for the outside temperature and the irradiance. This irradiance on a tilted surface is derived from the irradiance on a horizontal plane that is most often available from measurements sites by using the geometrical properties of the system: azimuth and elevation of the sun and tilt angle and orientation of the panels.
The irradiance on the horizontal plane data comes from ground measurements from a station close to the studied neighborhood which can for example be obtained from Agrometeorology Norway. 1 The elevation and azimuth of the sun is retrieved from an online tool. 2 This calculation takes into account the tilt of the solar panel and its orientation. Several assumptions were necessary to use this formula. Indeed, the solar irradiance is made up of a direct and a diffuse part and only the direct part of the irradiance is affected by the tilt and orientation. However there is no good source of irradiance data that provides a distinct measurement for direct and diffuse parts in Norway as far as the authors know. Thus we make assumptions that allow us to use the complete irradiance in the formula. We assume that most of the irradiance is direct during the day and that most is diffuse when the sun is below a certain elevation or certain azimuths. This assumption gives a good representation of the morning irradiances while still accounting for the tilt and orientation of the panel during the day. On the other hand, this representation overestimates the irradiance during cloudy days, when it is mostly indirect irradiance. Obtaining direct and diffuse irradiance data would solve this problem.

Implementation
The model presented in the previous section has been implemented in the case of campus Evenstad, which is a pilot project in the ZEN research center [15]. This implementation of the model and the parameters used are presented in this section. Campus Evenstad is a university college located in southern Norway and is made up of around 12 buildings for a total of about 10,000 m 2 . Most of the buildings were built between 1960 and 1990 but others stand out. In particular two small buildings were built in the nineteenth century and the campus also features two recent buildings with passive standards. The campus was already a pilot project in the previous ZEB center and one of those buildings was built as a Zero Emission Building. In addition, on the heating side a 100 kW CHP plant (40 kW electric) and a 350 kW Bio Boiler both using wood chips were installed along with 100 m 2 of solar collectors, 10,000 L of storage tank, 11,600 L of buffer tank and a heating grid. On the electric side, the same CHP is contributing to the on-site generation as well as a 60 kW photovoltaic system. A battery system is already planned to be built accounting for between 200 and 300 kWh. Based on this we assume in the study an existing capacity of 250 kWh. We keep those technology in the energy system of the neighborhood for one part of the study. In addition, the heating grid is kept in all cases (b hg = 1). The technologies included in the study are listed in Table 1 along with the appropriate parameters.
Two main sources for the parameters and cost of the technologies are used as references for the study. Most of the technologies' data is based on a report made by the Danish TSO energinet and the Danish Energy Agency [16] on technology data for energy plants. The other source includes the technology data sheets made by IEA ETSAP [17] and is used in particular for the gas and the biomass CHP. The cost of PV is based on a report from IRENA [18]. The two efficiencies reported for the CHP plants correspond to the thermal and electrical efficiency, noted by a subscript ( th for thermal and el for electrical). Note that: at the neighborhood level, only ground source heat pump is considered (Table 2) and that PV is only considered at the building level but the roof area limit to the size of the PV is not implemented.
The heat storage values are based on a data sheet by ETSAP [17] while the values used for the batteries are based on a report from IRENA [19].  The values in Table 3 come from different sources. The cost of biomass comes from EA Energy Analyses [20], the cost of gas is based on the cost of gas for non household consumers in Sweden 3 (we assume similar costs in Norway). For the technologies in Table 1, the O&M costs, expressed as a percentage of the investment costs, are respectively: 1, 1.3, 1, 1.3, 2, 0.8, 2.3, 4, 5.5, 1, 1 and 5. For the storage technologies in Table 2, the operating cost is 0. The CO 2 factors of gas and electricity for Norway are based on a report from Adapt Consulting [21] and the CO 2 factor for biomass is based on [22].
The electricity prices for Norway are based on the hourly spot prices for the Oslo region in 2017 from Nordpool. 4 On top of the spot prices, a small retailer fee and the grid charges are added. 5 The prices are rather constant with a fair amount of peaks in the winter and some dips in the summer. This cost structure is close to the actual structure of the electricity price seen by consumers. We assume hourly billing due to its relevance to prosumers and its emergence in Norway.
The irradiance on the horizontal plane and temperatures are obtained and used in the calculations as described in the previous section. The ground station used to retrieve data is Fåvang, situated 50 km to the west of Evenstad. The electric and heat load profiles for the campus are derived from [23]. The load profiles are based on the result of the statistical approach used in these papers and the ground floor area of each type of building on the campus. In addition, the domestic hot water (DHW) and Space Heating (SH) are derived from the heat load based on profiles from a passive building in Finland where both are known [24].
The problems are solved on a Windows 10 laptop with a dual-core CPU (i7-7600U) at 2.8 GHz and 16 GB of RAM. Each case typically has about 450,000 rows, 600,000 columns and 2,400,000 non-zeros. They are solved using the barrier method in Gurobi in about 150 s each.

Results
The optimization was run several times with different conditions. It was run with a yearly CO 2 balance with and without including the energy system that already exists at Evenstad. When the pre-existing energy system is included, the preexisting amounts of heat storage, PV, solar thermal and biomass heating (CHP and boilers) represent the minimum possible investments in those technologies for the optimization. The energy systems resulting from those optimizations are presented on Fig. 2.
Both cases are interesting. Indeed the case with the pre-existing technologies included in the optimization allows to know in which technology to invest to move towards being a ZEN for the campus Evenstad while the case that does not include the pre-existing technologies allows to see how it would look like if it was built today from the ground up using the optimization model presented here and the given ZEN restrictions.
A first observation from Fig. 2 is that the technologies already installed (heat storage ST, biomass boiler BB, CHP, battery) do not get additional investments, except for PV which gets a lot of additional investments to meet the ZEN criteria. In addition to the large investment in PV the only additional investment for Evenstad appears to be a heat pump. In the case without any pre-installed technologies the system is quite different. There is still a need for investment in PV, though it is slightly lower and the optimization does not chose to invest in a battery. On the heating part the chosen design uses heat pumps and electric boiler in addition to a heat storage smaller than already installed in Evenstad. The results highlight the predominance of PV in the results. This shows that the other possible designs are not cost competitive. Alternative designs, for example relying on biomass CHP, could be incentivized to obtain a better mix of technology. The amount of the incentive could be explored by a sensitivity analysis using this model, but this remains as future work.
On Fig. 3, the self consumption and the total demand of electricity is presented while on Fig. 4 it is the imports (red) and exports (blue) of electricity that are presented. Both figures show a week for the case of the yearly balance and including pre-existing technologies. In the summer the neighborhood produces electricity in excess and needs to send it to the grid. The battery, that is part of the pre-existing technologies, is used but is not large enough to allow for relying on self produced electricity during the night. It is also not large enough to limit the amount of electricity sent to the grid. Figure 4a illustrates this: the exports during the days have highs peaks that represent around four times the night imports in terms of peak power. This has implications on the sizing of the connection to the grid and is especially important in the context of the introduction of new tariffs based on peak power in Norway. Indeed, the introduction of smart-meters enables the use of more complex grid tariff structures. Such tariffs would promote avoiding large peaks in consumption. This may be beneficial to highly flexible neighborhoods such as ZENs and might promote investment in batteries. Investigating the impact of grid tariffs on the design of ZENs remains as future work. Outside of the ZEN context, a positive impact of certain grid tariff designs has been shown on self-consumption and peak electricity import [9]. In the winter, some of the electricity is still self consumed due to the CHP that is part of the pre-existing technologies. This self consumption stays limited and no electricity is exported.
Ultimately, all resulting designs require huge investment in PV to attain the status of ZEN. In those systems, which rely heavily on electricity, heat pumps and electric boilers appear to be the preferred heating solution.

Limitations
This study has several limitations, on the methodology and on the case study. For the case study, assumptions were necessary due to the lack of data, in particular for the loads or the insolation (diffuse and direct). For the methodology, the will to limit the use of binary variables meant leaving out constraints such as part load limitations which would be needed to have a better representation of some technologies. In addition, using an hourly resolution leads to an underestimation of the storages and possibly of the heating technologies size. There is a trade off between the solving time and the precision of the results and the resolution needs to be chosen accordingly. Additionally, being deterministic, the model leaves out several uncertainties. Those uncertainties concern the evolution of the price of the technologies, the electricity price or the price of other fuels and the climate conditions. Those can be partially addressed by specifying additional periods in the model. The short-term uncertainties are not included either and induce an overly optimistic operation of the system. Despite those limitations it provides insights in the design methodology that can be used to design the energy system of a ZEN. The choice of CO 2 factors for electricity is also greatly impacting the results and this should be studied in more detail in future work.

Conclusion
This paper presented in detail the ZENIT model for investment in Zero Emission Neighborhoods as well as its implementation and the results on a realistic case study of campus Evenstad in Norway, with a focus on methods from the field of operations research. The model is formulated as a MILP, using as few binaries as possible. The Zero Emission constraint complexifies the problematic of designing the energy system of a neighborhood and the long term trends can be accounted for by defining periods. For Evenstad, the results suggest that additional investments, mainly in PV, are necessary in order to attain the status of ZEN. Investments happen at both levels but mainly at the building level. When the technologies already installed at Evenstad are not included, they are not invested in (except for heat storage). The optimal choice in order to become Zero Emission for Evenstad in the current ZEN framework thus appears to be a massive investment in PV and a heating system fueled by electricity. Further work includes disaggregating the heat part of the model and a more detailed operation part in the optimization.
There are key takeaways for policy makers in this study, in particular for Norway due to the setting of the case study. The results suggest that the Zero Emission constraint used in this study is sufficient to get PV investment without any additional incentive. However, under the CO 2 factor assumptions used in this study, huge investment in PV are made which would be problematic in case of a largescale application of the concept of ZEN. This suggests the need for incentives in alternative technologies such as biomass CHPs in case the concept of ZEN becomes more common. The methodology presented in this paper can be used to assess such policies and their potential effect on investments in ZEN. Other policies such as the grid tariff structure can also be studied with this model. Finally, the hourly limitation on electricity export from prosumers has recently been replaced in Norway by a tariff on exported electricity. The results of this paper suggest that without this change in policy, ZEN would become even more expensive due to the necessity of large batteries to make the exports more constant. We thus recommend continuing on the path of facilitating the development of the number of prosumers for example with the implementation of capacity grid tariffs. For countries other than Norway, similar methodology can be used to assess the cost and design of ZEN. Further policy recommendations cannot be drawn from this study due to the specificity of the Norwegian electricity mix, that is reflected in its electricity CO 2 factor.

t (T )
Timestep The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Introduction
Grid storage systems on megawatt scale play a vital role for the integration of renewable energies into electricity markets and grids. Several investigations focus on the development of optimized battery operation strategies [1][2][3][4][5]. For several reasons, existing grid storage systems usually consist of multiple batteries and inverters. For reasons of economies of scale, hardware manufacturers offer components in limited number of size classes, allowing lower production costs. Furthermore, the use of modular components supports the systems' scalability. In addition, the size of modular systems can be changed over their lifetime. Moreover, modular systems avoid single points of failure, which leads to higher fault tolerance. Although the modularity of existing grid storage systems is well known, most of the analyses describe the storage system as a single battery combined with a single inverter [1][2][3][4][5][6][7][8][9]. Few studies are known that analyze the modularity of grid-connected battery systems and the related effect on system efficiency during operation and the influence of these energy losses on the operating strategy of the system [10,11].
The consideration of this architectural aspect in the model-based analyses of battery operation provides a degree of freedom in optimizing the overall yield of a grid storage system. Figure 1 shows an example for a grid storage system with a high modularity. In this setup, three inverters and batteries are connected to one point of common coupling.
Storage systems operated together with a renewable energy source are most likely not charged and discharged at their nominal power. Figure 2 shows the frequency of the operational inverter power over the course of a year for a large battery that is operated together with a wind farm with the purpose of supporting market integration of wind energy [14].
During 79% of the time, the system is in standby mode. During 5% of the time the system is operated at less than 50% rated power, while during 14% the power rating is between 80% and 100%.
This suggests the importance of energy conversion efficiency not only in fullload operation, but also in standby mode and in partial-load operation. The inverter efficiency, however, is a nonlinear function in terms of power flow (Fig. 3 with  In partial-load operation, the relative inverter losses are significantly higher in comparison to relative losses at full-load operation. Because of the possibility to use several inverter units, the efficiency is higher in high modularity systems. At low power flows, using one smaller inverter unit can limit the losses. Figure 3 shows the effect of number of inverters on the total efficiency curve. Here, we assume the power rating of the total storage system to be constant. Therefore, the power rating of the inverter and the constant losses scale with 1 N . The impact on the efficiency is evident up to a number of inverters of 4. Constant losses become less relevant with the increase from 4 to 8 inverters. Figure 3 provides an overview of the influence on the efficiency of modular battery systems. For operating a system with higher modularity, with discrete battery-inverter pairs, the operating strategy is more complex. Without an optimization strategy, which combines knowledge about power transfer losses, the SOC of the batteries and requests from the grid, suboptimal operating states can be realized, which increases losses. Hence, the operating strategy has a significant influence on the yield of a grid storage system. In this study we avoid choosing an individual optimization strategy by formulating an optimization problem. This optimization problem is formulated in terms of power flows and minimizes the power losses of the system. So, changes in the operating strategy, caused by increasing the number of inverter-battery pairs, are identified by solving the optimization problem.
The target function optimizing the operating strategy of the system is purely technical driven, i.e. we optimize in terms of power losses. Since most of the business models for grid storage systems rely on the power in and outflows of the system for economic evaluation, an optimization in terms of losses is always of benefit for the system operator.
As there is no direct connection between the batteries on the DC-side the energy management needs to take the state of charge of each battery into account. This increases the complexity of the energy management strategy.
This study addresses the questions to which extent the segmentation of the system into multiple battery-inverter subsystems can reduce operating losses and aims at quantifying this effect. Therefore, we develop a mixed-integer linear program that represents the operational strategy of an energy management system of such a modular system as it could be applied in real applications. The mixed-integer linear program determines the system operation for meeting a given schedule for the whole system at the point of common coupling by minimizing inverter losses. As a result, the required energy to be fed into or stored from the grid is allocated over the different battery-inverter-subsystems. The resulting SOC of the batteries is investigated and ideas for further investigations are derived.
Furthermore, the overall system topology has an influence on the power losses and the operating strategy. This study reduces the complexity by looking for identical battery-inverter pairs, which are coupled on the AC-side. A coupling on the DC-side has not been investigated in this study. This paper is structured as follows. Section 2 describes the mathematical model and presents the data chosen for the analysis. Section 3 presents the results and compares the efficiencies of systems with different degrees of modularity. Section 4 draws conclusions for the design of grid-connected battery systems.

Power Flow Model
The grid storage system is described as a set of nodes, representing the discrete time steps t with duration of the time steps d, with transfer of energy from one node to another [12,13]. In this model, each pair of inverter and battery represents a single storage node. The inverter is only described by its influence on the power transfer from the battery into the grid and vice versa. The SOC represents the battery state of charge and interlinks the periods between each other, adding energy flows, both for energy charged and energy discharged during the previous period. N storage systems S i are connected to one point of common coupling. G t represents an additional connection to the grid. It serves as a backup when power requests to the storage cannot be met and shall avoid infeasibility of the model.
Hereinafter, a power flow from node A to node B is described as A B and the efficiency of this transfer is described as η AB (A B ). Per definition, transfer losses are assigned to the sink. The underlying power model is generic and initially independent of the technical implementation of the storage system. Different technologies might add different boundary conditions and parameters to the power flow. All parameters underlying the assumed technical realization are described in 2.1.3 to 2.1.5.

Consideration of Losses
The power transfer losses can be divided into constant, linear, and quadratic terms. This results in the common representation of the inverter efficiency (1): In case of battery storage systems, constant losses a A B have their origin in auxiliary power, e.g. for the battery management system, active cooling and thermal control of the battery cells, and other subsystems. These losses are independent from the status of the inverter. Some systems can reduce the constant losses during standby operation by switching off subsystems. We consequently assume that the systems considered in our analysis can be turned off, avoiding standby losses during operation, as it is shown to be possible in Munzke et al. [14].
Linear losses b A B are proportional to the rated power. They are mainly heat losses. Furthermore, we consider battery efficiency. This is reflected by losses that are proportional to the charging and discharging power.
Quadratic losses c A B represent losses caused by nonlinear saturation effects in the inductance. In this work, the quadratic terms are not taken into account because non-linear losses cannot be clearly observed in all inverter systems [14].
Self-discharge of the battery cells, here represented as power loss of the storage path η SS , is usually very low and thus neglected [14][15][16]. Moreover, we do not consider the power supply of the battery management system (BMS).
We furthermore do not take into account losses which are independent of the degree of modularity, such as the power consumption of sensors or the climate control. Moreover, we do not consider transformer losses, neglecting the effect of additional power flows, introduced for ensuring model feasibility, on transformer losses.
Inverter loss data is obtained based on a curve fitting approach of the efficiency curve of SMA's 250 kVA inverter "Sunny Central" [17]. The battery loss parameters are approximated based on data measured in Munzke et al. [14]. Table 1 shows the exogenously given parameters used in the model. Table 2 lists the decision variables as well as their lower and upper bounds. Power values are always given in Watt. The penalty parameter p and two decision variables for additional, unplanned power flows to and from the grid (G t L and G t C ) are introduced for reducing deviations from the predetermined schedule and at the same time ensuring the model solvability. The same penalty parameter is applied both for charging from and discharging into the grid.

Target Function
The target function minimizes all losses and penalizes additional power flows (2). The solution is obtained by carrying out one optimization for every hour, i.e. four 15-min time steps at a time.
The target function determines the behaviour of the optimized power management strategy. In a modular system, the given power demand can be met by a subset of inverters. Given the example of two inverters in Fig. 3, only one inverter is used if it can provide the requested power, thus limiting constant losses. When power demand exceeds the single inverter's power rating, the second subsystem is used. In this case, there is no incentive to operate the modules at different power. Therefore, if charge and discharge power requests are sufficiently high, we expect a nearly equal distribution on the SOC.

Constraints
Two binary variables are introduced for charging and discharging (γ G S i , γ S i G ∈ {0; 1}) in order to ensure that the storage system is not charged and discharged at the same time (3).
The power flow at the point of common coupling is defined by Eq. (4): Equation (5) describes the energy flow for each storage system. It is solved for each hour independently. Therefore, the SOC t =0 is set as an external parameter. For the subsequent iterations, the SOC-value of the first time step is set as parameter according to the solution of the optimization of the previous optimization.
In addition, Eqs. (6) and (7) assure that battery charging or discharging is only realized if the binary variable is equal to 1.
Equations (8) and (9) set the boundaries for the power demand at the point of common coupling. These boundaries are calculated outside of and prior to the optimization.

Implementation
The model is implemented in MATLAB, using CPLEX as a solver. The solution is obtained on an hourly basis in 15-min resolution with a MIP gap of 0.5%.

Data
The battery schedule was calculated by Ried et al. [4]. It results from a planned operation of a 50 MW/100 MWh battery which is connected to a 50 MW wind farm. The system participates in the German day-ahead and tertiary control reserve markets. The battery schedule is a time series in 15-min resolution and obtained by a mixed-integer linear program maximizing the contribution margin of the system. For this study, the schedule is scaled to a 2 MW/2 MWh battery system and used as the power demand at the point of common coupling G t . Since it does not account for detailed losses, the battery used for this study is scaled 30% larger. All assumptions are given in Table 3. Based on this data, the model is applied to six systems with varying degrees of modularity (Table 4).

Results and Discussion
In this section, we compare the yearly losses and resulting battery operation for the different systems, and discuss potential implications for the system layout. Figure 4 shows the losses of systems with varying degrees of modularity relative to the losses of the system with N = 1. In accordance with the assumption shown in Fig. 3, the losses scale with the number of inverter-battery subsystems. By increasing the number of inverters to 32, the operating losses can be reduced by 38%. The gradual decrease of operating losses declines as the modularity increases. This effect can be explained if the distribution of charge and discharge power over the course of one year for one and 32 inverters are compared. In a single inverter system, 50-70% of the charge and discharge power, neglecting standby mode, lies between 70% and 80% of the rated inverter power (Fig. 5). While 91% of the energy is charged at a rated power above 70%, 88% of the energy is discharged above 70% power rating. Power below 50% is requested during 29-31% of the time, corresponding to 5-6% of the energy flow. This is the mode of operation in which the individual inverter is less efficient and losses increase.

Losses
As the number of inverters increases, the maximum power of each inverter decreases. Figure 6 shows the resulting distribution of charge and discharge power. As expected, the most likely power is the maximum power of the inverter, which corresponds to 1 N of the maximum power of the single inverter setup. Only very few events occur, where low power is requested.

Battery Operation
A difference between the systems with one and 32 inverters is a different SOC-level during the course of the year. Figure 7 shows the distribution of the state of charge over the course of a year for systems consisting of one and 32 sub-systems. In the non-modular system, the battery is operated at an average 13% SOC, while in the high-modularity system the average SOC is 27%. While 86% of the SOC-values in the non-modular system range between 10-20%, 57% of the SOC-values are in this range in the system consisting of 32 inverters and batteries. Due to lower operating losses, less energy is wasted and higher SOC-values occur more often.
Due to a larger number of batteries, there are relatively more downtimes of the batteries in high-modularity systems. While the system consisting of one inverter stands still during 79% of all periods, the sub-systems of the high-modularity systems are in standby-mode during 87% of the time. This suggests a higher redundancy of the system consisting of 32 batteries, which could be beneficial in case of failure, especially when power output must be guaranteed.

Conclusions
In summary, we have performed a study on the operating losses of grid-connected battery storage systems on megawatt-scale consisting of several battery-inverter subsystems. Therefore, we applied a mixed-integer linear program determining the optimal operation of all sub-systems for a given schedule for the point of common coupling of a case study where a battery is operated together with a wind farm.
The analysis shows that by increasing the system modularity, the frequency of operating points at low power decreases, resulting in a reduction of operating losses. The higher degree of freedom in modular systems thus allows following a given schedule at lower losses. By increasing the modularity from one to 32 sub-systems, operating losses can be reduced by almost 40%. The effect of modularity on efficiency is most evident for systems consisting of few sub-systems. With increasing modularity, the gradual decrease of operating losses declines. The resulting state-of-charge of the batteries shows a similar operation pattern for the investigated systems. The avoidance of losses in high-modularity systems, however, leads to a shift towards higher SOC levels.
We believe that this work motivates a modular layout of grid-connected battery systems. A higher efficiency should translate into lower operating cost. Moreover, the higher redundancy of modular battery systems can be advantageous particularly in applications that must ensure the system's availability. The related investment cost is another aspect to consider when determining the system design. On the one hand, very large inverter and battery systems might translate into lower investment cost due to less balance of system components. A higher degree of modularity could thus be associated with higher costs. The opposite could also be the case, if the modularity enables the usage of standardized components produced in higher volumes, overcompensating the higher costs for the peripheral components.
Further analyses could address the effect of considering monetary aspects within the target function. Moreover, the penalty parameter ensures model feasibility, but does not fully avoid additional power flows. Both the additional power flows and the choice of the penalty parameter could be further investigated. In order to penalize any deviation from the predetermined schedule, the same penalty parameter was applied both for charging from and for discharging into the grid. However, especially deviations leading to battery discharging might be less critical and sometimes even desirable from energy system perspective, possibly leading to positive revenues. This could be the case if the system participates in additional electricity markets. Consequently, the value of the penalty parameter could be different for positive and negative deviations and should be subject to further research.
Maybe the most interesting potential for further research lies in the operation of the different batteries, the effect on battery ageing and implications on technology choice and system design. Our results show higher average SOC-levels of the batteries in high-modularity systems. The average SOC has an impact on the calendar ageing, so that higher SOC-levels could lead to reduced battery lifetime, at least for most technologies. Another aspect worth investigating is the effect of higher C-rates on battery degradation. Because of allocating the same requested power to fewer batteries with lower capacities, the C-rates might be higher in high-modularity systems. The modular system architecture might thus be more suitable for some battery technologies than for others. The selection of one or more appropriate battery technologies could positively influence battery lifetime.
Lastly, the battery management system determining the operational strategy could be adapted in order to take quadratic losses into account or vary operating patterns, SOC levels, and power requests for different storages. Moreover, asymmetric system topologies would allow for stressing the batteries at different depths of discharge and C-rates, even when running the same schedule. This might be an argument for hybrid storage systems consisting of different battery technologies with different characteristics and furthermore improve cycle ageing of the batteries, another important factor for battery degradation.