Ksenia Syrtseva
email: k.syrtseva@edf.fr>

Welington De Oliveira

Sophie Demassey

Hugo Morais

Paul Javal

Bhargav Swaminathan

Difference-of-Convex Approach to Chance-Constrained Optimal Power Flow modelling the DSO Power Modulation Lever for Distribution Networks

Keywords:

The increasing expansion of renewable energy sources leads to the growth of uncertainty in the distribution network operation. Short-term operational planning performed by distribution system operators should evolve to address those new operating conditions, in particular to allow the efficient utilization of different flexibility levers. In this work, the use of a chance-constrained Alternating Current Optimal Power Flow (AC-OPF) is proposed to model the operational planning problem, considering the activation of several levers such as power modulation and power curtailment. The correlation between the renewable generation profiles and loads is taken into account via a joint probability constraint in the chance-constrained AC-OPF problem. The main novelty of the present manuscript is the adoption of a Difference-of-Convex approach that allows to solve the obtained optimization problem without convexification or linearization of the core OPF equations. The method starts with a reformulation of the model as a Difference-of-Convex optimization problem, and then a modified Bundle method algorithm is applied to solve it. The proposed methodology was tested in a 33 bus distribution network with 11 different values for the chance constraint satisfaction probability (safety parameter) ranging from 0.75 to 1.

Introduction 1.Context

The use of renewable energy sources (RES) is becoming increasingly important as the world seeks to transition to a more sustainable and environmentally friendly energy system. Several countries are proposing new targets to achieve carbon neutrality. For example, in Europe, several programs and initiatives have been defined such as Green Deal, Fit to 55 or RePower EU. In all documents, the targets are becoming more and more ambitious with some concrete measures such as incentives for the development of renewables and the electrification of the transport sector [START_REF] Bouckaert | From the green deal to repower EU[END_REF]. Analysing with more detail the mentioned initiatives, it is clear that the transmission and distribution systems will be neuralgic infrastructures.

Feed-in tariffs (FiT) for RES production were one of the first incentive measures implemented by governments to encourage their development [START_REF] Zhao | Feed-in tariffs, knowledge stocks and renewable energy technology innovation: The role of local government intervention[END_REF]. With these tariffs, RES production units could inject all the produced energy and receive fixed and economically attractive remuneration. Furthermore, system operators could not modulate the power injected by these units into their grid except under critical operation conditions. In practice, it implied that the system operators should be prepared to receive the contractualized capacity of production, which leads to significant investments in grid infrastructure. Depending on the circumstances, it resulted in a significant increase in cost for these producers, delays in commissioning, or even temporary limitations on the power injection during reinforcement works. Nowadays, considering the level of maturity of the RES technologies, mainly wind generation and solar photovoltaic, FiT have been replaced by other mechanisms [START_REF] Leiren | Germany: From feed-in tariffs to greater competition[END_REF]. Smart (and interruptible) connection points (SCP) and contracts have been proposed [START_REF] Anaya | Going smarter in the connection of distributed generation[END_REF] for RES and are being implemented in Europe [START_REF]Arrêté du 12 juillet 2021 d'application de l'article D. 342-23 du code de l'énergie[END_REF]. These contracts allow an increase in the hosting capacity of existing systems and thus faster integration of RES while reducing investment costs.

RES, such as solar and wind power, are inherently variable and intermittent, and therefore present significant challenges to the efficient and reliable operation of power systems. To address these challenges, system operators, mainly distribution system operators (DSOs), have to evolve from a traditional fit and forget approach in network planning and operation [START_REF] Fattaheian-Dehkordi | An incentive-based mechanism to alleviate active power congestion in a multi-agent distribution system[END_REF] to a more proactive one. Such a new approach involves optimization of multi-annual investments in human resources, infrastructure and technology, and the adoption of new methods and tools. In this context, many European DSOs are implementing a streamlined time-continuum based approach to network planning, operational planning and real-time operation, allowing a better and faster integration of new distributed energy resources (DER). This is, for example, highlighted in Enedis' (France's largest DSO) latest multi-annual network development plan [START_REF] Enedis | Network Development Plan[END_REF]. Typically, this consists of characterising large investments several years in advance, identifying and prioritising mediumterm investments, optimising works schedules in the shortterm, and operating networks in an optimized fashion. This has however, and will continue to be, highly dependent on political decisions, as is the case with above-mentioned FiT and SCP mechanisms.

New network operational planning tools specifically tailored to these needs have already been deployed and are being improved to tackle the newest problems highlighted, for example, in [START_REF] Enedis | Network Development Plan[END_REF]. These improvements include mathematical optimization based approaches to solving problems, which represent a complete shift in the operational planning paradigm for DSOs as it deals with the practical rules in use-cases and contractual obligations through mathematical optimization. DSOs like Enedis are continuing to work on deployment of newer methods that take into account other problems and opportunities like uncertainty in DER production, interaction with external actors (e.g. flexibility providers and TSOs), and integration of these new tools in their IT and OT systems.

These optimization methods rely on power flow equations in order to ensure technical constraint satisfaction. This, along with the objective of minimising costs associated with the operational solution, gives rise to an alternatingcurrent optimal power flow problem (AC-OPF). In the case of the industrial optimization solution developed by Enedis, a deterministic and linearized AC-OPF model is used [START_REF] Swaminathan | OptimGPindustrial optimisation for enedis' short-term operational planning[END_REF]. However, in a general case, an OPF problem is an optimization problem with mixed integer variables that is strongly NP-hard and nonconvex, as stated by various studies [START_REF] Bienstock | Strong NP-hardness of AC power flows feasibility[END_REF]. Several formulations and solution methods have been proposed [START_REF] Frank | Optimal power flow: a bibliographic survey i[END_REF][START_REF] Qiu | A literature survey of optimal power flow problems in the electricity market context[END_REF], including classical optimization methods that rely mainly on convex relaxations or approximations [START_REF] Low | Convex relaxation of optimal power flow-part i: Formulations and equivalence[END_REF][START_REF] Gan | Convex relaxations and linear approximation for optimal power flow in multiphase radial networks[END_REF], nondeterministic search techniques (also known as heuristic, stochastic or random search methods) [START_REF] Frank | Optimal power flow: a bibliographic survey ii. non-deterministic and hybrid methods[END_REF], and machine learning methods [START_REF] Hasan | A survey on applications of machine learning for optimal power flow[END_REF].

Optimal Power Flow (OPF) under uncertainties

As traditional OPF models assume that the system parameters are deterministic, which may not be the case in practice [START_REF] Soares | Active distribution grid management based on robust AC optimal power flow[END_REF], modelling uncertainties has become an important research topic in the field of electric power systems, particularly with the increasing penetration of RES. Solution methodologies of related optimization problems are diverse, and include meta-heuristic methods, robust and stochastic approaches. The former group has gained attention in the field of power and energy systems within the last ten years [START_REF] Chicco | Metaheuristic optimization of power and energy systems: Underlying principles and main issues of the 'rush to heuristics[END_REF], due to the ability to search for near-optimal solutions in a large solution space efficiently (an overview of methods can be found in [START_REF] Niu | A review on applications of heuristic optimization algorithms for optimal power flow in modern power systems[END_REF]). However, from a theoretical point of view, the quality of the obtained solution can be difficult to assess in practice, as no optimality certificate is provided.

Robust optimization methods rely on constructing an uncertainty set and searching for a reliable solution for any scenario realization. This approach is proposed for the management of distribution networks in [START_REF] Jabr | Adjustable robust OPF with renewable energy sources[END_REF][START_REF] Franco | AC OPF for smart distribution networks: An efficient and robust quadratic approach[END_REF][START_REF] Soares | Active distribution grid management based on robust AC optimal power flow[END_REF][START_REF] Xu | Robust transient stabilityconstrained optimal power flow with uncertain dynamic loads[END_REF] with different types of uncertain data. For instance, uncertainties are solely on renewable energy sources in [START_REF] Jabr | Adjustable robust OPF with renewable energy sources[END_REF], whereas [START_REF] Xu | Robust transient stabilityconstrained optimal power flow with uncertain dynamic loads[END_REF] deals with stochastic load composition. The authors in [START_REF] Soares | Active distribution grid management based on robust AC optimal power flow[END_REF] describe multi-period grid management applying a convex hull technique to define an uncertainty set. Solutions obtained with robust optimization methods are optimal for the worst-case scenario, and thus tend to be conservative. In contrast, stochastic optimization methods assume that the probability distribution of uncertain variables is known [START_REF] Shapiro | Lectures on Stochastic Programming: Modeling and Theory[END_REF]Chapter 1]. In particular, for the expected cost minimization framework, it allows to achieve (expected) minimal cost while preserving the required level of service for given scenarios. One widespread approach in stochastic programming is chance-constrained optimization, where one searches for a decision that minimizes costs while satisfying a set of random constraints with a prescribed probability level [START_REF] Shapiro | Lectures on Stochastic Programming: Modeling and Theory[END_REF][START_REF] Ansaripour | A chance-constrained optimization framework for transmission congestion management and frequency regulation in the presence of wind farms and energy storage systems[END_REF]. Chance-constrained optimization models are intuitive and straightforward to explain. However, they are generally difficult to solve, because they often lack essential mathematical properties such as convexity and differentiability [START_REF] Van Ackooij | A discussion of probability functions and constraints from a variational perspective[END_REF].

The use of chance (probability) constraints for energy management is discussed in detail in [START_REF] Van Ackooij | Chance Constrained Programming and Its Applications to Energy Management[END_REF], where the authors show how to ensure feasibility to the greatest extent possible, while aligning with the goals of the system operator. In the framework of OPF problems, most chance-constrained models rely on linear or convex approximations and deal with individual probability constraints, i.e., correlations between random variables are ignored. In [START_REF] Zhang | Chance constrained programming for optimal power flow under uncertainty[END_REF] an OPF problem with load uncertainties is considered, where probability constraints are represented with individual bounds on state variables. The model equations are linearized at the random vector expected value, establishing a monotonic relation between a constrained output and a random input. This fact enables to treat the probability constraints regarding the distribution of random variables. Partial linearization of power flow equations is applied in [START_REF] Roald | Chance-constrained AC optimal power flow: Reformulations and efficient algorithms[END_REF], in order to transform the individual chance constraints into deterministic ones at the forecasted operating point. One of the assumptions enabling the transformation is a relatively small forecast error on the model uncertainties. Another approximation is suggested in [10], where chance constraints are used to enforce voltage regulation with prescribed probability. The authors tackle them with convex relaxation and exploit linearization of load-flow equations. Reference [START_REF] Bienstock | Chance-constrained optimal power flow: Risk-aware network control under uncertainty[END_REF] deals with minimization of an average generation cost over random renewable power injections, while controllable generators mitigate power fluctuations. To guarantee that thermal limits are exceeded with low probability and ensure that the generation remains within imposed bounds, chance constraints are employed separately for each type of constraint. Then, the model is reduced to a deterministic convex optimization problem, more precisely, a second order cone program (SOCP). A SOCP reformulation is also used in [START_REF] Mühlpfordt | A generalized framework for chance-constrained optimal power flow[END_REF][START_REF] Bauer | Analytical uncertainty propagation for multi-period stochastic optimal power flow[END_REF] for the case of individual chance constraints. In [START_REF] Lubin | A robust approach to chance constrained optimal power flow with renewable generation[END_REF] a robust modification of the approach given in [START_REF] Bienstock | Chance-constrained optimal power flow: Risk-aware network control under uncertainty[END_REF] is presented, which addresses the uncertainty in parameters of probability distributions by restricting them to an uncertainty set. Another robust reformulation of a chance-constrained OPF modelling systems with fluctuating power sources, is obtained in [START_REF] Vrakopoulou | Probabilistic security-constrained AC optimal power flow[END_REF]. First, probabilistic bounds for the uncertainty are computed using scenario-based approach, and then a robust variant of the initial problem is solved due to these bounds. A well-known downside of individual chance constraints is that they ignore the random vector correlations, which are important statistics that should be considered when modelling real-life problems. As a strong correlation between the renewable generation profiles and loads is observed [START_REF] Verzijlbergh | Renewable energy sources and responsive demand. do we need congestion management in the distribution grid[END_REF][START_REF] Qi | Identification of critical parameters affecting voltage and angular stability considering loadrenewable generation correlations[END_REF], in the present work, the OPF problem is modelled as a joint chance-constrained optimization problem. For a given parameter 𝛼 ∈ [0, 1] and a probability measure ℙ, the problem we investigate in this work can be synthetized as follows:

𝐦𝐢𝐧 Levers { Levers activation cost } (1a)

(1c)

In this formulation, the decision variables are related to the activation of flexibility contracts (levers), which in our problem corresponds to DSO decisions on power modulation and power curtailment of distributed generation, and limitations on power supply to consumers (energy not supplied). Each grid user, who may be either a producer or a consumer, is characterized by its grid connection contract (FiT or SCP in our case). Depending on the contracts, levers activated by the DSO must satisfy specific deterministic constraints. Stochastic equations are related to technical decision feasibility: given a DSO decision and a scenario realization on power generation and loads, grid operating conditions must remain within technical limits. In practice, a decision on levers activation satisfying all the bounds and stochastic power flow equations for all the possible scenarios, may not exist. Hence, it makes sense to search for a reliable, affordable decision that satisfies the stochastic constraints with a probability level of at least 1 -𝛼. Observe that, if we take 𝛼 = 0 and a solution of problem (1) exists, then such a solution is robust over all the possible scenarios. The decision maker can thus adjust her risk aversion by appropriately setting the parameter 𝛼: small 𝛼 reflects high risk aversion.

In the proposed approach, the primary goal is to secure the required level 1 -𝛼, with 𝛼 > 0, for the system: we search for an implementable decision that is feasible with a probability of at least 1 -𝛼. As we will see shortly, the goal is obtained without any stringent assumption on the probability distribution or modelling simplifications such as linearization and convexification of core OPF equations. The deterministic case 𝛼 = 0 is not within the scope of the proposed methodology, which is applicable but not tailored for this case. In section 2, we propose a realistic mathematical model for (1) that is consequently nonconvex, nondifferentiable, and thus challenging. Nonconvexity comes from the OPF equations, whereas the requirement under probability sign in (1c) causes nondifferentiability of the problem.

Despite recent advances in theory and numerical methods for chance-constrained problems, dealing with multivariate probability functions in optimization problems generally remains a challenging task. The main difficulty arises in evaluating the probability function and computing its (sub)gradient with reasonable precision. We recall that evaluating the probability function amounts to computing numerically a multidimensional integral, a very difficult task for even moderate dimensions [START_REF] Van Ackooij | Level bundle methods for constrained convex optimization with various oracles[END_REF]. The situation becomes even more demanding for assessing (sub)gradients. For this reason, a common practice is to estimate the probability function with Monte-Carlo simulations using a finite sample of scenarios [START_REF] Pagnoncelli | Sample average approximation method for chance constrained programming: Theory and applications[END_REF]. For moderate sizes of the sample, and under the assumption that all involved functions are convex, one may model the probability constraint (1c) by using binary variables and BigM formulations (see, for instance, [START_REF] Van Ackooij | Inexact stabilized benders' decomposition approaches with application to chance-constrained problems with finite support[END_REF] and references therein). Such an approach is inappropriate for an OPF problem, because involved functions are not convex. Hence, without modelling simplifications, new mathematical approaches are required to face nonconvexity and nondifferentiability in problem [START_REF] Anaya | Going smarter in the connection of distributed generation[END_REF]. Based on the observation that the requirement under the probability in (1c) can be written as the difference of two convex functions, we apply a Difference-of-Convex (DoC) 1 optimization approach to tackle our problem. To the best of our knowledge, DoC programming has not yet been applied to the OPF class of problems.

Chance-constrained optimization and the DoC approach

A function is called DoC if it is expressible as the difference of two convex functions. As already investigated in [START_REF] Hong | Sequential convex approximations to joint chance constrained programed: A monte carlo approach[END_REF] and [START_REF] Van Ackooij | A bundle method for nonsmooth DC programming with application to chance-constrained problems[END_REF], probability functions can be approximated as accurately as one wishes by DoC functions. In Subsection 3.1, it is shown how to model (without approximation) the existence requirement in (1c) by a DoC function. Therefore, we end up with a composition of DoC functions that is itself DoC, and thus we can approximate probability constraint (1c) with a DoC constraint: ℙ [ Existence of a grid state within bounds satisfying stochastic power-flow equations

] ≥ 1 -𝛼 ≈ 𝑐 1 (𝑥) -𝑐 2 (𝑥) ≤ 0,
where 𝑐 1 , 𝑐 2 ∶ ℝ 𝑛 → ℝ are (nondifferentiable) convex functions. Given this approximation, we will show how to construct a DoC optimization model for problem (1) fitting into the following general structure:

min 𝑥∈𝑋 𝑓 (𝑥) ∶= 𝑓 1 (𝑥) -𝑓 2 (𝑥) (2a) s.t. 𝑐(𝑥) ∶= 𝑐 1 (𝑥) -𝑐 2 (𝑥) ≤ 0, ( 2b
)
where 𝑋 is a nonempty bounded polyhedron contained in ℝ 𝑛 and 𝑓 1 , 𝑓 2 ∶ ℝ 𝑛 → ℝ are convex functions as well. In the present work, we will refer to (2) as a DoC optimization problem. DoC programming forms an important subfield of nonconvex programming, as it covers a large class of nonconvex optimization problems from real-life applications. At the same time, convex analysis apparatus enables to establish optimality conditions for DoC problems and design algorithms to solve them. These facts explain the increasing interest in this field, which started in the 80s [START_REF] Thi | DC programming and DCA: thirty years of developments[END_REF]. Important facts about DoC programming such as optimality conditions and duality can be found in the tutorial paper [START_REF] De Oliveira | The ABC of DC programming[END_REF]. The survey article [START_REF] Thi | DC programming and DCA: thirty years of developments[END_REF] gives a large spectrum of examples and algorithm developments in this field. Currently, most used algorithms are based on iterative linearizations of components 𝑓 2 and 𝑐 2 [START_REF] Thi | DC programming and DCA: thirty years of developments[END_REF], on penalization technique applied to the DoC constraint [START_REF] Thi | DC programming and DCA for general DC programs[END_REF], and on improvement functions that combine constraint and objective in single level [START_REF] Van Ackooij | A bundle method for nonsmooth DC programming with application to chance-constrained problems[END_REF]. The latter is a well-known and successful strategy in the nonsmooth optimization literature [START_REF] Sagastizábal | An infeasible bundle method for nonsmooth convex constrained optimization without a penalty function or a filter[END_REF][START_REF] Apkarian | Mixed 𝐻 2 ∕𝐻 ∞ Control via Nonsmooth Optimization[END_REF][START_REF] Montonen | Bundle-based descent method for nonsmooth multiobjective DC optimization with inequality constraints[END_REF][START_REF] Van Ackooij | A bundle method for nonsmooth DC programming with application to chance-constrained problems[END_REF]. In particular, a bundle method with improvement function is proposed in [START_REF] Van Ackooij | A bundle method for nonsmooth DC programming with application to chance-constrained problems[END_REF] for dealing with DoC-constrained DoC-problems. Due to its good numerical performance reported in [START_REF] Van Ackooij | A bundle method for nonsmooth DC programming with application to chance-constrained problems[END_REF], we choose this bundle algorithm -denoted by DoC bundle method -to tackle our DoC optimization model (2) of the chance-constrained OPF problem (1).

Paper Structure and Contributions

The main contributions of the present paper are:

1. The formulation of the operational planning problem considering the activation of flexibility contracts (levers) under uncertainties, related to the probabilistic nature of nodal generation and consumption, as a joint chance-constrained OPF. The correlations among the random variables are thus taken into account. Furthermore, no model linearization or convexification is considered. 2. The design of an oracle (black-box) enabling to find a DoC decomposition of the constraint under probability sign, which imposes the OPF solution to be in the required bounds. This allows to apply DoC bundle method to the obtained chance-constrained OPF.

An illustration of the practical performance of DoC

approach on two use-cases.

The remainder of this paper is organized as follows. Section 2 is dedicated to the modelling of short-term operational planning on grid users power modulation under uncertainties as a chance-constrained OPF. A reformulation of the obtained problem as a DoC model is given in Section 3, followed by an explanation of DoC bundle method. In Section 4, numerical results for two use cases are provided, which illustrate the performance of the described method.

Chance-constrained OPF model

In what follows the chance-constrained OPF model is described by setting the decision variables, state variables and constraints. The employed notation is introduced along the way.

Decision variables and random vector

We denote by  the set of buses of the grid, and by  the set of lines. Among the end buses  , there is one slack bus and other buses with at most one connected grid user. Each grid user is either a producer or a consumer. The set of all grid users is denoted by  ⊂  , identifying a grid user 𝑖 ∈  with a corresponding bus 𝑖 ∈  . Let  and  be the subsets of producers and consumers, respectively. We wish to determine the active power 𝑝 𝑖 (positive for generation and negative for consumption) and reactive power 𝑞 𝑖 of each grid user 𝑖 ∈  , based on the following decomposition:

𝑝 𝑖 = 𝑝 𝜙 𝑖 (𝜉) -𝒑 𝜸 𝑖 -𝒑 𝝂 𝑖 , 𝑞 𝑖 = 𝑞 𝜙 𝑖 (𝜉) -𝒒 𝜸 𝑖 -𝒒 𝝂 𝑖 . ( 3
)
For buses 𝑖 ∈  ⧵  with no grid users, excluding the slack bus, we set 𝑝 𝑖 = 𝑞 𝑖 = 0. The active and reactive power demand (production and consumption), 𝑝 𝜙 𝑖 (𝜉) and 𝑞 𝜙 𝑖 (𝜉), are the uncertain parameters of our model, given as scenario realizations of a random vector 𝜉 ∈ Ξ ⊂ ℝ 𝑑 . In case of considered flexibility contracts, the decision bears on the active power modulation for SCP grid users within the bounds of their guaranteed power denoted by 𝒑 𝜸 𝑖 ; on the active power curtailment for FiT grid users and SCP grid users beyond the bounds of their guaranteed power, and on limitations on power supply to consumers (energy not supplied), denoted by 𝒑 𝝂 𝑖 . We introduce also pair variables for reactive power modulation: 𝒒 𝜸 𝑖 for SCP grid users within the bounds of their guaranteed power; 𝒒 𝝂 𝑖 for FiT grid users and SCP grid users beyond the bounds of their guaranteed power, and consumers. In fact, each reactive power variable is a function from the active, which will be defined in Subsection 2.3. This formulation was chosen as an easily modifiable for the case of other DSO levers, where the reactive power can be modulated directly. As the decision must be taken before scenario realization is known, the aforementioned decision variables do not depend directly on 𝜉. In order to distinguish different types of decision variables, we denote by  𝑆𝐶𝑃 the set of indexes for 𝒑 𝜸 and 𝒒 𝜸 , and by  𝐹 𝑖𝑇 that for 𝒑 𝝂 and 𝒒 𝝂 (note that  𝑆𝐶𝑃 ⊂  and  ⊂  𝐹 𝑖𝑇 ). To simplify the notation, if there is no need to distinguish user types, we will denote decision variables by On the contrary, the state variables |𝑉 𝑖 | and 𝛿 𝑖 representing voltage at bus 𝑖 ∈  , as well as 𝑝 𝑠𝑏 and 𝑞 𝑠𝑏 representing active and reactive power at slack bus, directly depend on 𝜉: once the scenario realization 𝜉 is known, as well as the active and reactive power 𝑝 𝑖 and 𝑞 𝑖 for all 𝑖 ∈  , state variables can be found by Gauss-Seidel, Newton-Raphson or similar methods [START_REF] Grainger | Power System Analysis[END_REF].

𝐩 ∶= ({𝒑 𝜸 𝑖 } 𝑖∈ 𝑆𝐶𝑃 , {𝒑 𝝂 𝑖 } 𝑖∈ 𝐹 𝑖𝑇 ) ∈ ℝ 𝑛 𝐪 ∶= ({𝒒 𝜸 𝑖 } 𝑖∈ 𝑆𝐶𝑃 , {𝒒 𝝂 𝑖 } 𝑖∈ 𝐹 𝑖𝑇 ) ∈ ℝ 𝑛 Ksenia Syrtseva,
Our problem thus has two main types of variables. First, the power modulation variables 𝐩 𝑖 and 𝐪 𝑖 are determined for all 𝑖 ∈  𝑆𝐶𝑃 ∪  𝐹 𝑖𝑇 before realization of the uncertain event. Then, after a scenario 𝜉 is observed, which defines the power production and consumption, 𝑝 𝜙 𝑖 (𝜉) and 𝑞 𝜙 𝑖 (𝜉) for each grid user 𝑖 ∈  , a new grid-state is determined. It is defined by the values of |𝑉 𝑖 | and 𝛿 𝑖 , 𝑖 ∈  , and 𝑝 𝑠𝑏 , 𝑞 𝑠𝑏 . The goal is to check whether an operational current transit (for a given scenario and power modulation) exists. In the next two subsections, we will detail how these two types of variables are combined in our chance-constrained OPF.

Power flow constraints under probability sign

We impose operational constraints which include bounds on voltage variables and constraints on the slack bus power variables. More specifically, we denote upper and lower bounds on voltage variables by |𝑉 𝑖 |, 𝛿 𝑖 and |𝑉 𝑖 |, 𝛿 𝑖 , respectively. Constraints on active and reactive power at slack bus 𝑝 𝑠𝑏 and 𝑞 𝑠𝑏 , are represented by a set  𝑠𝑏 defined in Demand Connection Code [START_REF]Network Code on Demand Connection DC[END_REF], Figure 1. In addition, we add thermal constraints on current transit (congestion constraints): for each line (𝑖, 𝑗) belonging to the set , we impose an upper limit (𝐼 max 𝑖,𝑗 ) 2 on a quadratic form 𝑙 from the current. All in all, we get the following set of constraints:

𝛿 𝑖 ≤ 𝛿 𝑖 ≤ 𝛿 𝑖 , ∀𝑖 ∈  (4a) |𝑉 𝑖 | ≤ |𝑉 𝑖 | ≤ |𝑉 𝑖 |, ∀𝑖 ∈  (4b) 𝑙(|𝑉 𝑖 |, |𝑉 𝑗 |, 𝛿 𝑖 , 𝛿 𝑗 ) ≤ (𝐼 max 𝑖,𝑗 ) 2 , ∀(𝑖, 𝑗) ∈  (4c) (𝑝 𝑠𝑏 , 𝑞 𝑠𝑏 ) ∈  𝑠𝑏 . ( 4d
)
In order to ease the model representation, let us define the following functions (here the sum is taken over all buses 𝑘 connected to the bus 𝑖):

L 𝑅 𝑖 (𝐩, 𝛿, |𝑉 |, 𝜉) ∶= ( 𝑝 𝜙 𝑖 (𝜉) -𝒑 𝜸 𝑖 -𝒑 𝝂 𝑖 ) + ∑ 𝑘∼𝑖 𝑌 𝑅 𝑖,𝑘 |𝑉 𝑖 ||𝑉 𝑘 | cos(𝛿 𝑖 -𝛿 𝑘 ) + ∑ 𝑘∼𝑖 𝑌 𝐼 𝑖,𝑘 |𝑉 𝑖 ||𝑉 𝑘 | sin(𝛿 𝑖 -𝛿 𝑘 ) L 𝐼 𝑖 (𝐪, 𝛿, |𝑉 |, 𝜉) ∶= ( 𝑞 𝜙 𝑖 (𝜉) -𝒒 𝜸 𝑖 -𝒒 𝝂 𝑖 ) + ∑ 𝑘∼𝑖 𝑌 𝑅 𝑖,𝑘 |𝑉 𝑖 ||𝑉 𝑘 | sin(𝛿 𝑖 -𝛿 𝑘 ) - ∑ 𝑘∼𝑖 𝑌 𝐼 𝑖,𝑘 |𝑉 𝑖 ||𝑉 𝑘 | cos(𝛿 𝑖 -𝛿 𝑘 ).
Such functions define the stochastic power flow equations of our problem.

We can now mathematically state the requirements under the probability sign in problem (1) by defining the following random set

𝑋(𝜉) ∶= ⎧ ⎪ ⎨ ⎪ ⎩ (𝐩, 𝐪) | | | | | there exist |𝑉 |, 𝛿, 𝑝 𝑠𝑏 , 𝑞 𝑠𝑏 satisfying (4), L 𝑅 𝑖 (𝐩, 𝛿, |𝑉 |, 𝜉) = 0 for all 𝑖 ∈  , L 𝐼 𝑖 (𝐪, 𝛿, |𝑉 |, 𝜉) = 0 for all 𝑖 ∈  . ⎫ ⎪ ⎬ ⎪ ⎭ (5)
Note that variables |𝑉 |, 𝛿, 𝑝 𝑠𝑏 , 𝑞 𝑠𝑏 depend directly on the random event 𝜉, whereas (𝐩, 𝐪) must be decided before realization of 𝜉. We are thus interested in finding a decision (𝐩, 𝐪) belonging to the random set 𝑋(𝜉) with probability 1 -𝛼, which is represented by a probability constraint in our optimization model.

Chance-constrained OPF formulation

We continue by defining deterministic constraints of the problem. We start with the conservative bound on generation and consumption

p𝜙 𝑖 ∶= min 𝜉∈Ξ 𝑝 𝜙 𝑖 (𝜉) for all 𝑖 ∈  p𝜙 𝑖 ∶= max 𝜉∈Ξ 𝑝 𝜙 𝑖 (𝜉) for all 𝑖 ∈ .
When the support set of the random vector 𝜉 is represented by finitely many Monte-Carlo scenarios, Ξ is a finite set, and thus computing p𝜙 𝑖 is a straightforward task. This parameter enters in our optimization problem bounding the decision variables on active power curtailment and modulation, 𝑝 𝜈 and 𝑝 𝛾 . The power curtailment must satisfy both contractual and technical bounds, for all realizations of the random vector 𝜉. Hence, we set the following constraints : This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4510304 P r e p r i n t n o t p e e r r e v i e w e d

where parameter tan 𝜙 𝑖 is given for each grid user. Meanwhile, for 𝑖 ∈  𝑆𝐶𝑃 , the bounds on power modulation are determined contractually. They can be modelled as fractions of installed or subscribed power (for producer and consumer, respectively). Due to the data used in practice, we implemented them as a fraction of p𝜙 𝑖 . All in all, we obtain the following constraints:

𝑎 - p𝜙 𝑖 ≤ 𝒑 𝜸 𝑖 ≤ 𝑎 + p𝜙 𝑖 and 𝒒 𝜸 𝑖 = tan 𝜙 𝑖 𝒑 𝜸 𝑖 , with 𝑎 -∈ [-1, 0] and 𝑎 + ∈ [0, 1].
As for the objective function, we consider the case when 𝑓 (𝐩, 𝐪) is convex, Subsection 4.3. However, our approach remains valid as long as the objective function is DoC only using a generalized version of the algorithm proposed in [START_REF] Van Ackooij | A bundle method for nonsmooth DC programming with application to chance-constrained problems[END_REF].

Based on the assumptions described in the previous paragraphs, it is now possible to present the chance-constrained model for the considered stochastic OPF problem:

min 𝐩,𝐪 𝑓 (𝐩, 𝐪) (6a) s.t. 0 ≤ 𝒑 𝝂 𝑖 ≤ p𝜙 𝑖 ∀𝑖 ∈  𝐹 𝑖𝑇 ∩  (6b) p𝜙 𝑖 ≤ 𝒑 𝝂 𝑖 ≤ 0 ∀𝑖 ∈  𝐹 𝑖𝑇 ∩  (6c) 𝒒 𝝂 𝑖 = tan 𝜙 𝑖 𝒑 𝝂 𝑖 , ∀𝑖 ∈  𝐹 𝑖𝑇 (6d) 𝑎 - p𝜙 𝑖 ≤ 𝒑 𝜸 𝑖 ≤ 𝑎 + p𝜙 𝑖 , ∀𝑖 ∈  𝑆𝐶𝑃 (6e) 𝒒 𝜸 𝑖 = tan 𝜙 𝑖 𝒑 𝜸 𝑖 , ∀𝑖 ∈  𝑆𝐶𝑃 (6f) ℙ [ (𝐩, 𝐪) ∈ 𝑋(𝜉)] ≥ 1 -𝛼. ( 6g
)

DoC reformulation and Bundle method

In this section, we first propose a DoC model for problem [START_REF] Bienstock | Chance-constrained optimal power flow: Risk-aware network control under uncertainty[END_REF], and then we show how to solve it by the bundle method of [START_REF] Van Ackooij | A bundle method for nonsmooth DC programming with application to chance-constrained problems[END_REF]. The proposed model is built upon two DoC reformulations: one exact for the mathematical requirement (𝐩, 𝐪) ∈ 𝑋(𝜉) and another approximate for the probability measure ℙ, as detailed in Subsections 3.1 and 3.2. Subsection 3.3 recalls the main details of the chosen optimization algorithm for solving the underlying problem.

DoC reformulation of the condition (𝐩, 𝐪) ∈ 𝑋(𝜉)

Let 𝜉 ∈ Ξ be fixed. For notational convenience, we will denote the decision variables (𝐩, 𝐪) of our problem by 𝑥 ∈ ℝ 2𝑛 . Our development starts by noting that the squared distance function to 𝑋(𝜉), i.e.,

𝑑 2 𝑋(𝜉) (𝑥) ∶= min 𝑧∈𝑋(𝜉) 1 2 ‖𝑧 -𝑥‖ 2 ,
yields the following useful relation:

𝑥 ∶= (𝐩, 𝐪) ∈ 𝑋(𝜉) ⟺ 𝑑 2 𝑋(𝜉) (𝑥) = 0, 𝑥 ∈ ℝ 2𝑛 .
The squared distance function is not convex as the set 𝑋(𝜉) given in ( 5) is not convex. However, it is a DoC function as the following development shows:

𝑑 2 𝑋(𝜉) (𝑥) = min 𝑧∈𝑋(𝜉) { 1 2 ‖𝑥‖ 2 -⟨𝑧, 𝑥⟩ + 1 2 ‖𝑧‖ 2 } = 1 2 ‖𝑥‖ 2 -max 𝑧∈𝑋(𝜉) { ⟨𝑧, 𝑥⟩ - 1 2 ‖𝑧‖ 2 } = 𝑔 1 (𝑥, 𝜉) -𝑔 2 (𝑥, 𝜉).
Indeed, for any arbitrary 𝜉 fixed, the convexity of both functions

𝑔 1 (𝑥, 𝜉) ∶= 1 2 ‖𝑥‖ 2 (7a) 𝑔 2 (𝑥, 𝜉) ∶= max 𝑧∈𝑋(𝜉) { ⟨𝑧, 𝑥⟩ - 1 2 ‖𝑧‖ 2 } (7b)
with respect to 𝑥 directly results from the definition of convexity. Moreover, observe that 𝑔 1 -𝑔 2 is non-negative.

As a result, we have the following DoC reformulation for checking whether (𝐩, 𝐪) satisfies the constraints in [START_REF] Bauer | Analytical uncertainty propagation for multi-period stochastic optimal power flow[END_REF]:

𝑥 = (𝐩, 𝐪) ∈ 𝑋(𝜉) ⟺ 𝑔 1 (𝑥, 𝜉) -𝑔 2 (𝑥, 𝜉) ≤ 0.
While 𝑔 1 is a simple quadratic function, 𝑔 2 (⋅, 𝜉) is more involved: it is the optimal value of a nonconvex optimization problem and is thus nondifferentiable. Its subdifferential at a given point 𝑥 is given by

𝜕𝑔 2 (𝑥, 𝜉) ∶= { 𝑠 2 (𝜉) ∈ ℝ 2𝑛 ∶ 𝑔 2 (𝑧, 𝜉) ≥ 𝑔 2 (𝑥, 𝜉) + ⟨𝑠 2 (𝜉), 𝑧 -𝑥⟩ ∀ 𝑧 ∈ ℝ 2𝑛 } = arg max 𝑧∈𝑋(𝜉) { ⟨𝑧, 𝑥⟩ - 1 2 ‖𝑧‖ 2 } = arg min 𝑧∈𝑋(𝜉) 1 2 ‖𝑧 -𝑥‖ 2 .
If the set 𝑋(𝜉) were convex (and closed), then the above projection problem would have had a unique solution and, thus, 𝑔 2 would have been differentiable. Note that evaluating 𝑔 2 at a given point 𝑥 and computing one of its subgradient (an element of the subdifferential) amounts to projecting 𝑥 onto the set 𝑋(𝜉) given in [START_REF] Bauer | Analytical uncertainty propagation for multi-period stochastic optimal power flow[END_REF]. This task can be accomplished (at least approximately) by OPF tools because this projection problem is indeed an OPF with a quadratic objective function. A numerical procedure (oracle) for accessing the DoC function 𝑔 1 (𝑥, 𝜉) -𝑔 2 (𝑥, 𝜉) and its first-order information is as follows. Here, we use the notation (𝐩 𝑘 , 𝐪 𝑘 ) to highlight that this point is fixed, conceivably given by an algorithm at its iteration 𝑘. The most time-consuming task in Oracle 1 is Step 1. To compute ( p, q), we may employ standard solvers from the OPF literature and leverage computational burden by first checking whether (𝐩 𝑘 , 𝐪 𝑘 ) satisfies the constraints in [START_REF] Bauer | Analytical uncertainty propagation for multi-period stochastic optimal power flow[END_REF]. The latter amounts to solving a system of load-flow equations and verifying if the computed solution (𝛿, |𝑉 |, 𝑝 𝑠𝑏 , 𝑞 𝑠𝑏 ) satisfies the bounds in (4). If it is the case, then (𝐩 𝑘 , 𝐪 𝑘 ) automatically solves the OPF of Step 1: there is no need for calling an OPF solver, but only a (simpler) load-flow algorithm. However, if the computed solution of the load-flow equations does not satisfy (4), then we say that (𝐩 𝑘 , 𝐪 𝑘 ) is infeasible for the future event 𝜉. In this case, an OPF solver must be applied to compute a point ( p, q) that is feasible for the scenario 𝜉 and as close as possible to (𝐩 𝑘 , 𝐪 𝑘 ). Such a step thus amounts to answering the following question: given that (𝐩 𝑘 , 𝐪 𝑘 ) is infeasible for scenario 𝜉, what is the smallest necessary perturbation on (𝐩 𝑘 , 𝐪 𝑘 ) to render it feasible?

The answer is ( p, q) -(𝐩 𝑘 , 𝐪 𝑘 ), which is nothing but the opposite direction to the subgradient 𝑠 𝑘 1 (𝜉) -𝑠 𝑘 2 (𝜉). This fact explicitly reveals the practical role of the subgradients computed by Oracle 1: it guides the optimization process to seek for a better candidate solution (at the next iteration), and the use of an efficient OPF solver allows to improve numerical performance. However, caution is necessary: in this analysis, the oracle only sees the given individual scenario 𝜉. It is thus necessary to account for all the scenarios and the probability measure in a higher-level oracle. That is the goal in the following subsection.

DoC reformulation of probability constraint

This subsection aims at presenting a DoC approximation for the probability constraint (6g). From the previous subsection we have

ℙ[(𝐩, 𝐪) ∈ 𝑋(𝜉)] ≡ ℙ[𝑔 1 (𝐩, 𝐪, 𝜉) -𝑔 2 (𝐩, 𝐪, 𝜉) ≤ 0],
with 𝑔 1 , 𝑔 2 two convex functions given in [START_REF] Bienstock | Strong NP-hardness of AC power flows feasibility[END_REF]. Hence,

ℙ[(𝐩, 𝐪) ∈ 𝑋(𝜉)] ≥ 1 -𝛼 is equivalent to ℙ[𝑔 1 (𝐩, 𝐪, 𝜉) -𝑔 2 (𝐩, 𝐪, 𝜉) > 0] ≤ 𝛼.
Next, we follow the lead of [START_REF] Van Ackooij | A bundle method for nonsmooth DC programming with application to chance-constrained problems[END_REF] to approximate the probability measure by a DoC function. To this end, let 𝑣(𝜉) = 𝑔 1 (𝐩, 𝐪, 𝜉) -𝑔 2 (𝐩, 𝐪, 𝜉) be the random variable of interest, 𝔼[⋅] the expected value operator w.r.t. ℙ, and let 𝟏 (0,∞) (⋅) denote the characteristic function of the segment (0, ∞), that equals to 1 if 𝑣 > 0, and 0 if 𝑣 ≤ 0. Recall that the following useful equivalence:

ℙ[𝑣(𝜉) > 0] = 𝔼[𝟏 (0,∞) (𝑣(𝜉))].
The main source of difficulties is that 𝟏 (0,∞) (⋅) is not convex and, even worse, it is discontinuous at 0. As in [START_REF] Hong | Sequential convex approximations to joint chance constrained programed: A monte carlo approach[END_REF] and [START_REF] Van Ackooij | A bundle method for nonsmooth DC programming with application to chance-constrained problems[END_REF],

we now approximate the characteristic function by a DoC one. Given a small parameter 𝑡 > 0, the discontinuous characteristic function can be approximated by the continuous one

𝜁 𝑡 (𝑣) ∶= ⎧ ⎪ ⎨ ⎪ ⎩ 0, if 𝑣 ≤ 0 𝑣 𝑡 , if 0 < 𝑣 ≤ 𝑡 1, if 𝑡 < 𝑣. (8)
Observe that lim 𝑡↓0 𝜁 𝑡 (𝑣) = 𝟏 (0,∞) (𝑣) and 𝜁 𝑡 (𝑣) has the following DoC decomposition:

𝜁 𝑡 (𝑣) = max { 𝑣 𝑡 , 0 } -max { 0, 𝑣 -𝑡 𝑡 } .
Hence, the following expected value

𝔼[𝜁 𝑡 (𝑔 1 (𝐩, 𝐪, 𝜉) -𝑔 2 (𝐩, 𝐪, 𝜉))]
is an approximation of the probability

ℙ[𝑔 1 (𝐩, 𝐪, 𝜉) -𝑔 2 (𝐩, 𝐪, 𝜉) > 0].
Such approximation is as good as one wishes: the smaller is the parameter 𝑡 > 0, the closer 𝜁 𝑡 (⋅) is to 𝟏 (0,∞) (⋅). Furthermore, the composition of DoC functions under the expectation is itself a DoC function: by their sample average approximations. The justification of such an approach is well documented in the literature (e.g., [START_REF] Pagnoncelli | Sample average approximation method for chance constrained programming: Theory and applications[END_REF]), and specialized to the the DoC setting in [START_REF] Hong | Sequential convex approximations to joint chance constrained programed: A monte carlo approach[END_REF] and [START_REF] Javal | Integrating uncertainties in short-term operational planning[END_REF]Subsection 7.7.2]. Hence, we can approximate the probability constraint (6g) with the following DoC constraint

𝜁 𝑡 (
𝑐 1 (𝐩, 𝐪) -𝑐 2 (𝐩, 𝐪) ≤ 0, where ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝑐 1 (𝐩, 𝐪) ∶= 1 𝑁 𝑁 ∑ 𝑗=1 max{𝑔 1 (𝐩,
s.t. 𝑐 1 (𝐩, 𝐪) -𝑐 2 (𝐩, 𝐪) ≤ 0, (10a)
where

𝑋 ∶= {(𝐩, 𝐪) ∈ ℝ 2𝑛 ∶ (6b) -(6f)} (10b)
is a polyhedral set. As discussed in Subsection 1.3, this optimization problem fits the structure of ( 2), and we can apply the DoC bundle method of [START_REF] Van Ackooij | A bundle method for nonsmooth DC programming with application to chance-constrained problems[END_REF] to tackle it.

DoC Bundle method

We start by highlighting that oracles for the functions in problem (10) are readily available. Indeed, 𝑓 is a convex function and thus simple. Furthermore, an oracle for 𝑐 1 and 𝑐 2 , providing their values and first-order information, is readily implementable thanks to Oracle 1 and assumption that we have finitely many 𝑁 scenarios to represent the future random events: given 𝑥 ∶= (𝐩, 𝐪) ∈ 𝑋, an oracle provides (𝑐 1 (𝑥), 𝑠 1 ∈ 𝜕𝑐 1 (𝑥)) and (𝑐 2 (𝑥), 𝑠 2 ∈ 𝜕𝑐 2 (𝑥)). Note, however, that such an oracle is not a straightforward one: it requires calling Oracle 1 𝑁 times for the same 𝑥 = (𝐩, 𝐪) given. In other words, 𝑁 deterministic OPF (or load-flow problem) must be solved by Oracle 1 to compute (via [START_REF]Network Code on Demand Connection DC[END_REF] the function values and a pair of subgradients for 𝑐 1 and 𝑐 2 .

Given these oracles, we can go further and briefly present the solving methodology, which has many numerical advantages: it does not require a feasible initial point, does not need penalty parameters, and numerical experience suggests that the approach is likely to escape bad-quality critical points. The interested reader is referred to [START_REF] Van Ackooij | A bundle method for nonsmooth DC programming with application to chance-constrained problems[END_REF] for a detailed presentation of DoC bundle method, as well as its mathematical properties. In that paper, the algorithm is given for a general case of DoC objective function, while the version presented below is adapted for convex objective. We start with the improvement function definition, an essential tool for presenting the algorithm: 𝐻 τ (𝑥) = max{𝑓 (𝑥) -τ𝑓 , 𝑐 1 (𝑥) -𝑐 2 (𝑥) -τ𝑐 }. [START_REF] De Oliveira | The ABC of DC programming[END_REF] In this definition, τ ∈ ℝ 2 is a parameter. The best choice for its value is τ = (𝑓 * , 0), with 𝑓 * the optimal value of (10): in this case, any solution of the problem min 𝑥∈𝑋 𝐻 τ (𝑥) solves (10). However, as 𝑓 * is unknown we update τ iteratively as follows: given x ∈ 𝑋, a candidate point to solve (10) produced by the algorithm, we set τ ∶= ( τ𝑓 , τ𝑐 ) = (𝑓 ( x)+𝜌 max{𝑐( x), 0}, 𝜎 max{𝑐( x), 0}), [START_REF]Arrêté du 12 juillet 2021 d'application de l'article D. 342-23 du code de l'énergie[END_REF] where 𝜌 > 0 and 𝜎 ∈ [0, 1) can be freely chosen. The DoC bundle algorithm can be depicted as in Algorithm 1.

Algorithm 1 DoC bundle method

1: Given 𝑥 0 ∈ 𝑋, 𝛼 ∈ [0, 1], choose 𝜇 > 0 and 𝜅 ∈ (0, 1). Set x ← 𝑥 0 and τ as in (12) 2: Compute (𝑓 (𝑥 0 ), 𝑠 0 𝑓 ∈ 𝜕𝑓 (𝑥 0 )) and (𝑐 𝑖 (𝑥 0 ), 𝑠 0

𝑖 ∈ 𝜕𝑐 𝑖 (𝑥 0 )), 𝑖 = 1, 2 3: for 𝑘 = 0, 1, 2 … do ⊳ Trial point 4:
Let 𝑥 𝑘+1 be the 𝑥-part solution of the quadratic program

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ min 𝑥∈𝑋,𝑟∈ℝ 4 𝑟 4 -⟨𝑠 𝑘 2 , 𝑥⟩ + 𝜇 2 ‖𝑥 -x‖ 2 s.t. 𝑓 (𝑥 𝑗 ) + ⟨𝑠 𝑗 𝑓 , 𝑥 -𝑥 𝑗 ⟩ ≤ 𝑟 1 , 𝑗 = 0, … , 𝑘 𝑐 1 (𝑥 𝑗 ) + ⟨𝑠 𝑗 1 , 𝑥 -𝑥 𝑗 ⟩ ≤ 𝑟 2 , 𝑗 = 0, … , 𝑘 𝑐 2 (𝑥 𝑗 ) + ⟨𝑠 𝑗 2 , 𝑥 -𝑥 𝑗 ⟩ ≤ 𝑟 3 , 𝑗 = 0, … , 𝑘 𝑟 1 + 𝑟 3 -τ𝑓 ≤ 𝑟 4 𝑟 2 -τ𝑐 ≤ 𝑟 4 ⊳ Stopping test 5: if ‖𝑥 𝑘+1 -x‖ ≤ 𝚃𝚘𝚕 then 6:
Stop and return x 7:

end if ⊳ Oracles call 8: Compute (𝑓 (𝑥 𝑘+1 ), 𝑠 𝑘+1 𝑓 ∈ 𝜕𝑓 (𝑥 𝑘+1 ))
and

(𝑐 𝑖 (𝑥 𝑘+1 ), 𝑠 𝑘+1 𝑖 ∈ 𝜕𝑐 𝑖 (𝑥 𝑘+1 )), 𝑖 = 1, 2 ⊳ Descent test 9: if 𝐻 τ (𝑥 𝑘+1 ) ≤ 𝐻 τ ( x) -𝜅 𝜇 2 ‖𝑥 𝑘+1 -x‖ 2 then 10:
Set x ← 𝑥 𝑘+1 and update τ as in [START_REF]Arrêté du 12 juillet 2021 d'application de l'article D. 342-23 du code de l'énergie[END_REF].

11:

end if

12: end for

We highlight that Algorithm 1 is a simplified version of the one presented in [START_REF] Van Ackooij | A bundle method for nonsmooth DC programming with application to chance-constrained problems[END_REF]. For instance, the prox-parameter 𝜇 > 0 can be updated iteratively, the number of constraints in the quadratic (master) program can be kept bounded, and the objective function can be DoC. The convergence analysis given in [48, Section 4]

Results

Network

The current research study utilizes a medium voltage distribution network, inspired by [START_REF] Baran | Network reconfiguration in distribution systems for loss reduction and load balancing[END_REF][START_REF] Sousa | Intelligent energy resource management considering vehicle-to-grid: A simulated annealing approach[END_REF], consisting of 33 buses, accommodating 31 loads and with a total peak consumption of approximately 16.6 MW. The network incorporates three distributed generation (DG) units, comprising two biomass plants (one with SCP contract connected in bus 12 with an installed capacity of 6.97 MW and another with FiT contract connected in bus 29 with an installed capacity of 3.86 MW), and one wind farm with a FiT contract connected in bus 32 with an installed capacity of 0.47 MW. The network is connected to the high-voltage network in bus 1 which is considered the slack bus in the present study. The consumption and generation values were defined to have the network operating near by their technical limits creating the need for flexibility activation. The schematic representation of the network is illustrated in Figure 3.

Scenario generation

In order to test the proposed DoC approach on the test network described in Subsection 4.1, we use Enedis Open Data [START_REF] Enedis | Enedis open data[END_REF] on July 27, 2020, to construct load and generation profiles for 𝑁 = 1000 scenarios. Three types of grid users are considered: biomass generation, wind generation and consumption (data based on small and medium enterprises). For each grid user 𝑖 ∈  , we denote the provided data by ( p𝜙 𝑖 , q𝜙 𝑖 ). Next, we attribute an individual variance 𝜎 2 𝑖 : 0.0248 for biomass generation type, 0.01 for wind generation type, and 4 different values in the range between 6.01 ⋅ 10 -5 and 0.01 for consumers. Consider covariance matrix Σ with Σ 𝑖,𝑗 = 𝑟 𝑖,𝑗 𝜎 𝑖 ⋅ 𝜎 𝑗 , where 𝑟 𝑖,𝑗 = 1 if the types of users 𝑖 and 𝑗 coincide, and 0 otherwise; and normalization matrix 𝐴 with 𝐴 𝑖,𝑗 = 0.5 if the user 𝑖 = 𝑗 is a generator, 𝐴 𝑖,𝑗 = 5000 if the user 𝑖 = 𝑗 is a consumer, and 0 otherwise. Applying the procedure described in [START_REF] Higham | Computing a nearest symmetric positive semidefinite matrix[END_REF], we calculate the nearest symmetric positive semi-definite matrix Σ to 𝐴𝜎. Finally, we generate vectors 𝑝 𝜙 and 𝑞 𝜙 following multivariate Gaussian distributions with the means p𝜙 and q𝜙 , respectively, and covariance Σ.

Parameters of DoC approach

Here and in what follows 𝑃 𝑏𝑎𝑠𝑒 = 1 MW, 𝑉 𝑏𝑎𝑠𝑒 = 12.66 kV and 𝐼 𝑏𝑎𝑠𝑒 = 78.99A.

The upper and lower bounds in (4b) are defined in the network design manual of a French DSO [13, Annex 1.3] specifically related to connection contracts for grid users, and represent ±5% of nominal voltage for MV network. The upper and low bounds on 𝛿 in (4a) are set to ± 𝜋 2 . For all (𝑖, 𝑗) ∈ , thermal constraints (4c) are reduced to |𝐼 𝑖,𝑗 | 2 ≤ (𝐼 max 𝑖,𝑗 ) 2 . In order to simplify computations of Oracle 1, we use a convexification of the feasible set  𝑠𝑏 in (4d), Figure 1. Such a convexification is defined by

 = { (𝑝 𝑠𝑏 , 𝑞 𝑠𝑏 ) ∈ ℝ 2 ∶ 𝑝 min 𝑠𝑏 ≤ 𝑝 𝑠𝑏 ≤ 𝑝 max 𝑠𝑏 , 𝑞 min 𝑠𝑏 ≤ 𝑞 𝑠𝑏 ≤ 𝑞 max 𝑠𝑏 , 𝑞 𝑠𝑏 ≥ -0.48𝑝 max 𝑠𝑏 -𝑝 min 𝑠𝑏 + 0.25𝑝 max 𝑠𝑏 𝑝 𝑠𝑏 + 0.48𝑝 max 𝑠𝑏 𝑝 min 𝑠𝑏 -𝑝 min 𝑠𝑏 + 0.25𝑝 max 𝑠𝑏 } .
and corresponds to the union of green and blue sets on Figure 4. The objective function is convex and has the following structure (all the coefficients are non-negative):

𝑓 (𝐩, 𝐪) = 𝑓 1 (𝐩) + 𝑓 2 (𝐩), (13a)
where and

𝑓 1 (𝐩) = ∑ 𝑖∈ 𝑆𝐶𝑃 𝐶 𝑖 𝑆𝐶𝑃 |𝒑 𝜸 𝑖 | + ∑ 𝑖∈ 𝐹 𝑖𝑇 ∩ 𝐶 𝑖 𝐹 𝑖𝑇 ,𝑔 |𝒑 𝝂 𝑖 | + ∑ 𝑖∈ 𝐹 𝑖𝑇 ∩ 𝐶 𝑖 𝐹 𝑖𝑇 ,𝑙 |𝒑 𝝂 𝑖 | (13b)
𝑓 2 (𝐩) = ∑ 𝑖∈ 𝑆𝐶𝑃 𝑐 𝑖 𝑆𝐶𝑃 (𝒑 𝜸 𝑖 ) 2 + ∑ 𝑖∈ 𝐹 𝑖𝑇 ∩ 𝑐 𝑖 𝐹 𝑖𝑇 ,𝑔 (𝒑 𝝂 𝑖 ) 2 + ∑ 𝑖∈ 𝑐 𝑖 𝐹 𝑖𝑇 ,𝑙 (𝒑 𝝂 𝑖 ) 2 . ( 13c
)
The values of coefficients in (13b) and (13c) are given in Table 1. Our motivation for this choice is as follows. First, the component 𝑓 1 (𝐩) reflects the cost of levers activation (since coefficients in Table 1 are set as dimensionless quantities, this cost is expressed in pu). As power modulation of a FiT producer is more expensive compared to that of a SCP grid user, the inequality2 𝐶 𝑖 𝐹 𝑖𝑇 ,𝑔 ≫ 𝐶 𝑗 𝑆𝐶𝑃 should be respected for 𝑖, 𝑗 ∈ . Moreover, as limitation on power supply to a consumer is more expensive than power curtailment of producer, the inequality 𝐶 𝑖 𝐹 𝑖𝑇 ,𝑙 ≫ 𝐶 𝑗 𝐹 𝑖𝑇 ,𝑔 should be satisfied for 𝑖, 𝑗 ∈  𝐹 𝑖𝑇 . Thus, the following relation 𝐶 𝑖 𝐹 𝑖𝑇 ,𝑙 ≫ 𝐶 𝑗 𝐹 𝑖𝑇 ,𝑔 ≫ 𝐶 𝑘 𝑆𝐶𝑃 is respected in (13b). Further, as penalty for the energy not supplied is defined by a regulator, we set 𝐶 𝑖 𝐹 𝑖𝑇 ,𝑙 = 𝐶 𝑗 𝐹 𝑖𝑇 ,𝑙 for all consumers 𝑖, 𝑗 ∈ . Since the contracts may differ depending on generation technology, we fix different coefficients for biomass generation and wind generation. Meanwhile, the quadratic terms are introduced as a penalty that encourages fairness in power modulation among the same types of users, as the minimum in (13c) is attained at a point where activation of levers is equal among grid users with equal coefficients. Thus, we set 𝑐 𝑖 𝐹 𝑖𝑇 ,𝑙 = 𝑐 𝑗 𝐹 𝑖𝑇 ,𝑙 for 𝑖, 𝑗 ∈  𝐹 𝑖𝑇 , i.e. we impose the coefficients to be equal for FiT consumers. We also assume that the quadratic coefficients are equal among FiT producers of the same generation type. Unless otherwise specified, the approximation parameter 𝑡 from ( 8) is set to 10 -5 . The choice of parameters in DoC Bundle method [START_REF] Van Ackooij | A bundle method for nonsmooth DC programming with application to chance-constrained problems[END_REF] is as follows: 𝜌 = 10 7 , 𝜎 = 0.5, 𝜅 = 0.9, 𝛿 𝑇 𝑜𝑙 = 10 -7 , 𝜇 min = 10 -6 and 𝜇 max = 10 6 . To test the performance of the presented DoC approach for the chance-constrained problem (6), we set 11 values of safety parameters 1 -𝛼 ranging from 0.75 to 1 with a step size of 0.025 (the case 𝛼 = 0 is tested as an extreme one, the algorithm is not designed for deterministic framework). We consider two cases: one where only voltage constraints are detected (buses 9, 10, 11, 15, 16 and 17) in initial state of the grid, i.e. without levers activation, and another with an additional congestion constraint. For the second case, we set an upper limit on current for the line connecting buses 2 and 19. The constraint value is chosen in such a way that it is violated for 533 scenarios in initial state of the grid (without levers activation).

Case 1: Voltage constraints

For each value of safety parameter, the algorithm manages to find a feasible critical point with average execution time of 1665 seconds (in a personal laptop) ranging from 1059 (1 -𝛼 = 0.875) to 2279 (1 -𝛼 = 0.775) seconds. In order to check validation of the chance constraint, we compare targeted value of safety parameter with the ratio of scenarios satisfying load-flow equations. In initial state of the grid, without levers activation, the latter is equal to 0.545. After the optimization is performed, it gets close to targeted safety parameter with a tendency to approach it as the parameter increases, Figure 5. It remains below the targeted value due to approximation of probability function used for DoC formulation (see Subsection 3.2). Varying parameter 𝑡 that participates in DoC approximation of probability constraint (see eq. ( 8)), 𝑡 = 10 -4 , 10 -5 and 10 -6 , for a fixed safety parameter value 1 -𝛼 = 0.9, the ratio of scenarios satisfying load-flow equations becomes 0.839, 0.876 and 0.893, respectively. This illustrates, as expected, an increase in accuracy as 𝑡 goes to zero. The volume of active power modulation grows linearly up to 1 -𝛼 = 0.925 and accelerates afterwards, as illustrated in Figure 6. The same tendency is valid for the volume supplied by FiT grid users, whose share in total power curtailment remains within 69 -77%. Meanwhile, a SCP grid user supplies the remaining part of the volume within the bounds of her guaranteed power for all values of safety parameter. Power modulation cost, represented by (13b) in the objective function, follows the same upward trend as power modulation volume. It is consistent with the price formation, which is in linear dependence on the latter.

Therefore, a slight decrease in safety parameter implies a reduction in total cost of power modulation. This result aligns with our expectations from probabilistic modelling of power flow constraints. Computing maximum and average amplitudes of constraint violation over all scenarios where constraint violations are detected, we compare them with average volume of power modulation supplied by one grid user. Without levers activation, the maximum amplitude of constraint violation is 0.0339 pu, whereas the average is 0.0057 pu. Figure 7 reveals a notable decrease in maximum amplitude of constraint violation for 1 -𝛼 > 0.825 up to 0.0102 pu (𝛼 = 0), and a significant growth in average volume of power modulation while increasing the system reliability. This confirms that the higher volume of levers is activated, the less amplitude of constraint violation is. To ensure the system reliability with 1 -𝛼 ≥ 0.9, average power modulation volume exceeds maximum amplitude of constraint violation. The gap between two figures drastically increases as 1-𝛼 approaches 1 due to stochastic character of our problem: if we did not have to deal with uncertainties and the scenario realization was known, covering the maximum amplitude of constraint violation would have been sufficient.

Meanwhile, the average amplitude of constraint violation fluctuates near the value 0.0045 pu with a slight tendency to decrease. These upward and downward trends in amplitude of constraint violation lead to the following conclusions. More reliable but costly solutions cover risky scenarios which allow higher amplitude of constraint violation. At the same time, they tolerate those where constraint violation amplitude is closer to average. This result is based on simulations, and does not follow from chance-constrained formulation as it says nothing about amplitude of constraint violation. Moreover, this conclusion contradicts mathematical intuition that we would have had for a convex problem: maximum amplitude of constraint violation grows, as risky scenarios are more costly to cover.

Case 2: Voltage and congestion constraints

The algorithm finds a critical point with average execution time of 1947 seconds ranging from 1353 (1 -𝛼 = 0.85) to 2850 (1 -𝛼 = 0.825) seconds. However, for deterministic case 𝛼 = 0, the algorithm does not manage to find a feasible solution. The ratio of scenarios satisfying load-flow equations is equal to 0.266 for the initial state of the grid. As Figure 8 shows, the latter approaches the targeted value of safety parameter once levers activation is optimized. However, for 1-𝛼 > 0.95, the difference between targeted and obtained value of safety parameter increases. This can be due to the fact that a robust solution may not exist. The difference between targeted and obtained value of safety parameter represents 4.6% of targeted value for 𝛼 = 0, whereas, initially, it was only 3.4% for 1 -𝛼 = 0.75. The gap for 𝛼 = 0 is significant due to infeasibility of obtained solution. This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4510304 P r e p r i n t n o t p e e r r e v i e w e d bounds of her guaranteed power. However, compared to the previous case, the growth acceleration of power modulation cost and volume is less pronounced as safety parameter goes to 1. Moreover, these values become almost stable at 1 -𝛼 = 0.975 and 𝛼 = 0, which is consistent with the fact that the corresponding number of covered scenarios is very close (953 for 1-𝛼 = 0.975 and 954 for 𝛼 = 0). Thus, except the deterministic case with 𝛼 = 0 when obtained solution is not feasible, a small decrease in safety parameter enables to reduce the total cost of levers activation. As in the previous case, we compute average and maximum amplitudes of voltage constraint violation. These values are calculated over all scenarios where voltage constraint violation is detected (scenarios where only congestion constraint is violated, are not included). Without levers activation, they remain the same as in the case without congestion constraint. Meanwhile, maximum and average amplitudes of congestion constraint violation calculated over all scenarios where congestion constraint violation is detected, are 0.3851 pu and 1.5155 pu, respectively. As the order of magnitude of congestion constraint violation is greater than that of voltage constraint violation, the corresponding values are plotted separately, Figure 10 and Figure 11. All in all, they illustrate the same trends as in the previous case, namely, a steady decrease in maximum amplitude of constraint violation, both for voltage and congestion constraints. At the same time, Figure 10 reveals a significant growth in average volume of power modulation while increasing the system reliability. Thus, the conclusion that an increase in volume of levers activation reduces the amplitude of constraint violation, remains valid.

Comparing curves of average power modulation volume on Figure 7 and Figure 10, we observe that it is higher for the case with congestion constraint. Moreover, average power modulation volume always exceeds maximum amplitude of voltage constraint violation in the latter case. This is due to stochastic character of our problem, but also to additional amplitude of congestion constraint violation that should be covered by power modulation. Meanwhile, we see the opposite trend for the curves of maximum amplitude of constraint violation (without congestion constraint). It can be explained by the curtailment of grid users downstream of the congestion constraint for all values of safety parameter, which is aimed to cover an essential part of scenarios where congestion constraint is detected. This curtailment weakens voltage constraints at other buses of the grid, and reduces amplitude of their violation for not covered scenarios. Furthermore, the ratio of scenarios where congestion constraint is violated among all not covered scenarios, remains within 47 -59%, and tends to slightly increase as safety parameter goes to 1. In other words, more reliable but costly solutions are more prone to cover scenarios without congestion constraints.

Comparing activation of levers for grid users downstream the line connecting buses 2 and 19 (congestion constraint in Case 2) in two considered use cases, we observe that there is no power curtailment of those grid users for all values of safety parameter 1 -𝛼, 𝛼 > 0 in Case 1. This result is coherent with the expectations, as they do not contribute to resolution of existing voltage constraints. Meanwhile, they are curtailed in Case 2 for all values of safety parameter 1 -𝛼, 𝛼 > 0. As the voltages constraints are the same in both use cases, we conclude that their curtailment resolve This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4510304 P r e p r i n t n o t p e e r r e v i e w e d only the congestion constraint. Moreover, we also observe that the power curtailment is not equal among those users, although all of them are consumers. Consequently, the soft constraint imposed by quadratic terms (13c) in the objective function, is not sufficient to force an equal power curtailment inside the same group of grid users. This is due to nonconvex nature of the optimization model.

Conclusion

This paper presented a novel formulation of the distribution network operational planning problem as a chanceconstrained AC-OPF modelled using a DoC approach. This approach takes into account uncertainties related to nodal production and consumption, and proposes solutions to the chance-constrained problem through the activation of flexibility (power modulation) from SCP and FiT grid users. The solution methodology consists in model reformulation as a DoC problem with the use of a special numerical procedure (oracle), and the subsequent use of DoC Bundle method. The main advantages of the employed approach are:

1. The joint probability constraint in the obtained chanceconstrained OPF model allows to take into account correlation between RES generation profiles and load profiles. 2. The use of DoC approach enables to solve the optimization problem avoiding linearization and convexification of power flow equations. 3. The method allows to consider wide range of objective functions and constraints, as long as they can be represented as DoC functions. In particular, it can be applied for modelling other DSO levers, e.g. reactive power modulation. 4. Different OPF solvers/strategies can be employed in the oracle that computes function values and subgradients. Integration of an efficient solver improves the numerical performance of DoC Bundle method.

The perspectives of this work include application of a generalize version of DoC Bundle method presented in [START_REF] Syrtseva | Minimizing the difference of convex and weakly convex functions via bundle method[END_REF] for solving DoC reformulation of the chance-constrained OPF. From a theoretical point of view, this algorithm provides a solution satisfying stronger optimality condition, but has potentially lower computational efficiency as no explicit DoC decomposition of involved functions is used. Another perspective consists in adjusting the objective function and deterministic constraints of the model for the practical goals of the DSO, as the solution methodology gives a large spectre of options. However, the use of binary variable can be required for modelling some of DSO rules, for instance, an accurate requirement of the fairness in power modulation among the same types of grid users. For these cases, a modification of the solution methodology is necessary, as DoC Bundle method is adapted only for continuous variables.

[10] Dall'Anese, E., Baker, K., Summers, T., 2017. Chance-constrained AC optimal power flow for distribution systems with renewables. This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4510304 P r e p r i n t n o t p e e r r e v i e w e d

[

  s.t. Activated levers satisfy contractual constraints (1b) ℙ Existence of a grid state within bounds satisfying stochastic power-flow equations ] ≥ 1 -𝛼.

  Paul Javal, Welington de Oliveira, Sophie Demassey, Hugo Morais, Bhargav Swaminathan: Preprint submitted to Elsevier with 𝑛 = | 𝑆𝐶𝑃 | + | 𝐹 𝑖𝑇 |.

Figure 1 :

 1 Figure 1: Feasible 𝑃 -𝑄 diagram for slack-bus following Demand Connection Code.

0

  ≤ 𝒑 𝝂 𝑖 ≤ p𝜙 𝑖 and 𝒒 𝝂 𝑖 = tan 𝜙 𝑖 𝒑 𝝂 𝑖 , 𝑖 ∈  𝐹 𝑖𝑇 ∩  p𝜙 𝑖 ≤ 𝒑 𝝂 𝑖 ≤ 0 and 𝒒 𝝂 𝑖 = tan 𝜙 𝑖 𝒑 𝝂 𝑖 , 𝑖 ∈  𝐹 𝑖𝑇 ∩ . Ksenia Syrtseva, Paul Javal, Welington de Oliveira, Sophie Demassey, Hugo Morais, Bhargav Swaminathan: Preprint submitted to Elsevier Page 5 of 15

Oracle 1 .

 1 Black-box for the DoC function 𝑔 1 (𝑥, 𝜉)-𝑔 2 (𝑥, 𝜉).

Figure 2 :

 2 Figure 2: Function 𝜁 𝑡 (𝑣) for 𝑡 = 10 -4 , 𝑡 = 10 -5 and 𝑡 = 10 -6 .

Figure 3 :

 3 Figure 3: Medium voltage distribution network, 32 buses.

Figure 4 :

 4 Figure 4: Convexification of the feasible set for slack-bus constraint.

Figure 5 :

 5 Figure 5: Comparison of targeted safety parameter 1 -𝛼 with the obtained one (for t=10 -5 ).

Figure 6 :

 6 Figure 6: Cost and volume of power modulation.

Figure 7 :

 7 Figure 7: Volume of power modulation and amplitudes of constraint violation.

Figure 8 :

 8 Figure 8: Comparison of targeted safety parameter 1 -𝛼 with the obtained one.

Figure 9 :

 9 Figure 9: Cost and volume of power modulation.

Figure 10 :

 10 Figure 10: Volume of power modulation and amplitudes of constraint violation (excluding congestion constraint violation).

Figure 11 :

 11 Figure 11: Amplitudes of congestion constraint violation.

Ksenia Syrtseva ,

 Syrtseva Paul Javal, Welington de Oliveira, Sophie Demassey, Hugo Morais, Bhargav Swaminathan: Preprint submitted to Elsevier Page 12 of 15

  DoC Approach for Chance-Constrained OPF modeling the DSO Power Modulation Lever for Distribution Networks

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4510304 P r e p r i n t n o t p e e r r e v i e w e d

1 :

 1 Given (𝐩 𝑘 , 𝐪 𝑘 ) ∈ ℝ 2𝑛 and event/scenario 𝜉 ∈ Ξ, let ( p, q) be a solution of the OPF DoC Approach for Chance-Constrained OPF modeling the DSO Power Modulation Lever for Distribution Networks

	min (𝐩,𝐪)∈𝑋(𝜉)	1 2	‖(𝐩, 𝐪) -(𝐩 1 2 ‖𝐩 𝑘 ‖ 2 + 1 2 ‖𝐪 𝑘	‖ 2 and 𝑠 𝑘 1 (𝜉) ←	𝐪 𝑘 ( 𝐩 𝑘	)
	𝑠 𝑘 2 (𝜉) ←	(	p q)

𝑘 , 𝐪 𝑘 )‖ 2 , with 𝑋(𝜉) given in

[START_REF] Bauer | Analytical uncertainty propagation for multi-period stochastic optimal power flow[END_REF]

2: Set 𝑔 1 (𝐩 𝑘 , 𝐪 𝑘 , 𝜉) ← 3: Set 𝑔 2 (𝐩 𝑘 , 𝐪 𝑘 , 𝜉) ← ⟨𝐩 𝑘 , p⟩ -1 2 ‖ p‖ 2 + ⟨𝐪 𝑘 , q⟩ -1 2 ‖ q‖ 2 and 4: Return the first order information (𝑔 1 (𝐩 𝑘 , 𝐪 𝑘 , 𝜉), 𝑠 𝑘 1 (𝜉)) and (𝑔 2 (𝐩 𝑘 , 𝐪 𝑘 , 𝜉), 𝑠 𝑘 2 (𝜉)). Ksenia Syrtseva, Paul Javal, Welington de Oliveira, Sophie Demassey, Hugo Morais, Bhargav Swaminathan: Preprint submitted to Elsevier Page 6 of 15

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4510304 P r e p r i n t n o t p e e r r e v i e w e d

  𝐪, 𝜉 𝑗 ), 𝑔 2 (𝐩, 𝐪, 𝜉 𝑗 )}

					+𝑡(1 -𝛼),
	𝑐 2 (𝐩, 𝐪) ∶=	1 𝑁	𝑁 𝑗=1 ∑	max{𝑔 1 (𝐩, 𝐪, 𝜉 𝑗 ), 𝑔 2 (𝐩, 𝐪, 𝜉 𝑗 ) + 𝑡}
					(9)
	are convex functions. As a result, we have our DoC optimiza-
	tion model for the chance-constrained problem (6):
	{	min

(𝐩,𝐪)∈𝑋 𝑓 (𝐩, 𝐪)

  DoC Approach for Chance-Constrained OPF modeling the DSO Power Modulation Lever for Distribution Networks conditions, Algorithm 1 provides a critical point (i.e., a point satisfying necessary optimality conditions) of problem[START_REF] Ansaripour | A chance-constrained optimization framework for transmission congestion management and frequency regulation in the presence of wind farms and energy storage systems[END_REF].

ensures that, under certain Ksenia Syrtseva, Paul Javal, Welington de Oliveira, Sophie Demassey, Hugo Morais, Bhargav Swaminathan: Preprint submitted to Elsevier Page 8 of 15

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4510304 P r e p r i n t n o t p e e r r e v i e w e d

Table 1

 1 DoC Approach for Chance-Constrained OPF modeling the DSO Power Modulation Lever for Distribution Networks Coefficients in the objective function[START_REF] Enedis | Principes d'étude et de développement du réseau pour le raccordement des clients consommateurs et producteurs BT[END_REF] for each grid user (GU).

	GU	Contract	Type	Coeff. in	Coeff. in
				𝑓 1 (𝐩)	𝑓 2 (𝐩)
	12	SCP	Biomass	4.2 ⋅ 10 -5 0
			generation
	12, 29 FiT	Biomass	4.2 ⋅ 10 -3 0.01
			generation
	32	FiT	Wind	0.02	0.1
			generation
	Others FiT	Consumption 1	1

Ksenia Syrtseva, Paul Javal, Welington de Oliveira, Sophie Demassey, Hugo Morais, Bhargav Swaminathan: Preprint submitted to Elsevier Page 9 of 15

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4510304 P r e p r i n t n o t p e e r r e v i e w e d

  IEEE Transactions on Power Systems 32, 3427-3438. doi:10.1109/ TPWRS.2017.2656080.

Ksenia Syrtseva, Paul Javal, Welington de Oliveira, Sophie Demassey, Hugo Morais, Bhargav Swaminathan: Preprint submitted to Elsevier Page 13 of 15

The standard abbreviation for Difference-of-Convex is DC. To avoid confusion with "Direct Current" in OPF problems, we refer to the former as DoC.Ksenia Syrtseva, Paul Javal, Welington de Oliveira, Sophie Demassey, Hugo Morais, Bhargav Swaminathan: Preprint submitted to Elsevier

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4510304 P r e p r i n t n o t p e e r r e v i e w e d DoC Approach for Chance-Constrained OPF modeling the DSO Power Modulation Lever for Distribution Networks

The relation "𝑎 ≫ 𝑏" means that 𝑎 is much greater than 𝑏.

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4510304

6. Acknowledgements W. Oliveira acknowledges financial support from the Gaspard-Monge Program for Optimization and Operations Research (PGMO) project "Scalable Optimization for Learning and Energy Management." H. Morais was supported by Portuguese national funds through FCT, Fundação para a Ciência e a Tecnologia, under the project UIDB/50021/2020. K. Syrtseva acknowledges financial support from EDF R&D for the PhD thesis work.

CRediT authorship contribution statement


domain: shs sde sde.es

text:
Nicolas Bernigaud

Frédérique Bertoncello

Marie-Jeanne Ouriachi

Bondeau Alberte

Joel Guiot

Marie-Jeanne Ouriachi

Agent

Frédérique Azur

Marie-Jeanne Bertoncello

Ouriachi

Laure

M.-J Ouriachi

A Bondeau

Agent-Based Modelling, Archaeology and Paleoenvironment

Keywords:

what is this?

▪Software simulating interactions between agents and their environment ▪Interactions produce complex phenomena that we are trying to reproduce (or which are unexpected)

42e

  Rencontres Internationales d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 3 General elements of an agent-based model A typical agent (Figures after Macal and North 2014) Agent-based modeling (ABM) or Multi-Agent System (MAS):

  d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 5 ABM & Archaeology: kick off from USA in the 1990' I-General principles, historiographical overview •Growing Artificial Societies, Social Science from the Bottom Up (Joshua M. Epstein, Robert Axtell) •Dynamics in Human and Primate Societies, Agent-Based Modeling of Social and Spatial Processes (Timothy A. Kohler & George J. Gumerman) •Population growth and collapse in a multiagent model of the Kayenta Anasazi in Long House Valley (Axtell et al. 2002) •The Model-Based Archaeology of Socionatural Systems (Kohler & van der Leeuw 2007) al. 2002) 42e Rencontres Internationales d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 6 Artificial Anasazi (AA) (Axtell et al. 2002) https://ccl.northwestern.edu/netlogo/models/ArtificialAnasazi ENKIMDU (Wilkinson et al. 2005)

Agent▪▪➢

  principles, historiographical overview ▪ Complexity of human-environment dynamics (Human Ecodynamics) ▪ Understanding interactions between societies and environment as factor of evolutions/transformations ▪ Identify the internal and external causes of these changes: social, cultural, demographic and environmental factors (climate, erosion, vegetation, etc.) 42e Rencontres Internationales d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 7 principles, historiographical overview •2018 « Human History and Digital Future » Computer Applications and Quantitative Methods in Archaeology (CAA) ▪2020 # 415 « From Abacus to Calculus. Computational Approaches to Roman Economy » (I. Romanowska, Ph. Verhagen) •2021 # 367 « Agent-Based Modelling of socioecological Systems in Archaeology Towards a New Research Community » (D. Daems) European Association of Archaeologists (EAA) ▪ NASSA (Network for Agent-based modelling in Socioecological Systems in Archaeology) (Coord. Dries Daems) -Collect and compile ABM code elements -Organize them as an openly accessible modelling library 42e Rencontres Internationales d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 8 Mining with Agents (Kowarik et al. 2008) • How mine of salt were exploited in Hallstatt in Bronze Age ? Patagonian Hunter-gatherers (Barceló et al. 2010) • Emergence of ethnicity & cultural differentiation "Celtic Oppida" (Cimler et al. 2012, ) • Hypotheses related to agricultural practices around Celtic oppida, Carrying capacity, population growth MERCURY (Brughmans & Poblome 2015) • Market Economy and Roman Ceramics Redistribution ROMFARM (Joyce 2019, Verhagen 2022) • Modelling Agricultural Strategies in the Dutch Roman Limes I-General principles, historiographical overview Archaeological ABMs in Europe: some examples ▪ Varied periods: Prehistory, Bronze Age, Roman period,… ▪ Varied topics: Salt Mining, Agriculture, Social phenomenon, Roman Economy ➢ Not really a socionatural perspective for most models 42e Rencontres Internationales d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 9 I-General principles, historiographical overview ABM in France: geographers open the path ▪ Still few models! ANR OBRESOC (2009-2011) • Trajectoire du Néolithique Rubané • Coord. J.-P. Bocquet-Appel ANR TransMonDyn (2011-2014) • Modéliser les grandes transitions de l'évolution du peuplement dans l'Ancien et le Nouveau Monde • Coord. L. Sanders Project A*Midex RDMed (2018-2020) • Resilience and Adaptation to Droughts and Extreme Climate Events in the Mediterranean • Coord. J. Guiot SimPop (Guérin-Pace et al. 1993, Pumain et al. 1995) •Settlement pattern evolution •Genesis of urban phenomenon in southern France OBRESOC/ BEAN (Bocquet et al. 2015, Dubouloz et al. 2017, Zanotti 2016) •Spread of agriculture/LBK culture in Europe in early Neolithic •[ANR OBRESOC] AmphorABM (Crabtree 2016) •Modeling the Trade of Wine in the Bronze to Iron Age Transition in Southern France •[ANR TransMonDyn] ModelAnSet (Bertoncello et al. 2018) •Settlement pattern and dynamics during the Roman period in South-Eastern France (impact of social & environmental factors) SIMFEODAL (Cura 2020, Tannier et al. 2020) •Concentration of settlement in Middle Age •[ANR TransMonDyn] ROMCLIM/FER AGRI (Bernigaud et al. 2022) • Impact of Climate Change on gaulish and roman societies in southern Gaul •[RDMed project] 42e Rencontres Internationales d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 ▪ ModelAnSet simulates the processes of creation, growth, abandonment or maintenance of the rural settlements ▪ Agents: • 3 types of landowners [FARMERS] [BIG LANDOWNERS], [ARISTOCRATS] with different economic power, social status and behavior • 2 types of agricultural settlement [FARM] [VILLA] ▪ Model dynamics: • The profit made by each settlement depends on its environmental context, which is influenced by the climatic conditions, and on the transportation cost to the nearest marketplace • According to their economic power and to the profit they derive from each of their settlements, landowners decide to maintain, enlarge, improve or abandon it, or to create a new farm/villa Explore the respective roles of social & environmental factors in the settlement pattern and dynamics during the Roman period in the roman city of Forum Julii (Fréjus, Var) Instanciation of the model: the territory of the Roman colony of Forum Iulii (Fréjus, Var, France) (Bertoncello et al. d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 14 According to their social status (farmer, big landowner, aristocrat), landowners have different behaviour (spatial, symbolic) and perceptions (of land fertility, of the macro-economic context)  The repeated landowners' decision-making produces a changing macro-level settlement pattern, in terms of number, type and location of the settlements Needed improvements: • Introduce "cognitive agents" (vs "reactive agents" acting according to fixed rules) who can diversify their actions according to their beliefs (about themselves and about their environment) • Refine the calculation of the settlement's profit • More accurate simulation of the climate impacts on the agricultural yield II-Case studies (Bertoncello et al. 2018) ▪ ROMCLIM simulates the evolution of potential profits made by virtual agricultural exploitations (wine, olive and wheat production) in southern Gaul (Gallia Narbonensis), according Climate Change ▪ Agents: [WINERIES] [OLIVE GROVES] [CEREAL FARMS] [MAJOR TOWNS] (programmed according to agronomic treatises) ▪ Potential yields (wine, olive, cereals) are calculated by an emulation the agro-ecosystemic model LUND-POTSDAM-JENA-MANAGED-LAND (LPJML) ▪ Paleoclimatic data needed by LPJmL are processed for each century between the VI th c BC and the VII th PC from an actual climate reconstruction (ALADIN) modified with paleoclimatic anomalies (from Guiot & Kaniewski 2015) 15 42e Rencontres Internationales d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 Screenshot of the interface of the ROMCLIM Model II-Case studies Did the Roman Climate Optimum (RCO) boost the economy of the Roman Empire? /kg/km * Distance to the nearest capital Benefits Surplus price -costs Harvest price -costs Action Disappear if annual harvest can't feed the family Disappear if benefits < 0 Main characteristics of the agents: 42e Rencontres Internationales d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 17 II-Case studies ➢ Strong positive impact of RCO on the potential agricultural profits (peak in the II nd c. BC), slight decrease during Late Antiquity (IV th -VII th cd'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 18 II-Case studies Map of potential profit for wine Map of potential profit for olive oil ▪ For winemaking archaeological data for wine roughly match with the area of highest potential profits (Antibes-Fréjus area, middle et lower Rhone valley, languedocian plains) ▪ For olive oil, main archaeological data are concentrated in Languedoc and mainly in the Var region (thanks to J.-P. Brun works), where the potential appears medium. The area with highest potential profits in the languedocian hinterland doesn't match with archaeological data, but researches are not very intensive in this area. Difficulty to validate the outputs of the model with archaeological data due to their spatial heterogeneity -> but model have a stimulating predictive value for archaeological research Archaeological testimony 42e Rencontres Internationales d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 19 d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 ▪ Impact of Climate Change on potential cereal production and demographical growth are simulated ▪ Cereals yields vary according to paleoclimatic data (LPJmL emulation) II-Case studies Interface of FER-AGRI Model (Bernigaud et al. in press) Density of settlement in Northern Gaul in Iron Age and climatic fluctuations (Malrain et al. 2013) What was the impact of Climate Change on the Iron Age demography in Gaul between the VI th c. BC and the I st c. BC ? Temperatures Precipitations Population grows each year according to a demographic rate ➢ Agricultural workers (workforce) make a harvest per year for 6 centuries ➢ When the harvest and stocks are insufficient in one cell to feed everybody, population disappears Numbers of sites per regions 42e Rencontres Internationales d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 21 FER-AGRI Model: outputs, results and interpretation

heuristic tool for testing hypothesis (or a "tool to think")

	Scientifical
	questions:
	How, how many,
	Why, When, etc.
	➢Understanding the social and environmental
	causes of changes seen by archaeologists in
	archaeological records
	➢Test in ABM theorical or data driven
	hypothesis by simulations
	➢Confronting outputs of the model with ▪More than fifteen computer programs are dedicated to archaeological records for the ABMs (NetLogo, Gama, Repast, Swarm, Mason, etc.) validating/invalidating hypotheses ▪Numerous fields of application: Agriculture, ➢If all hypotheses are invalided, try to make a Anthropology, Biomedical Research, Crime Analysis, Ecology,… new hypothesis!

I-General principles, historiographical overview 42e Rencontres Internationales d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 4 I-General principles, historiographical overview Why using ABM in archaeology ? ABM: a

  Assessments and Conclusions: ABM: The End of the Beginning? ▪ Powerful tool for explore hypothesis & rethinking historical narratives ▪ Use of ABMs grows in Europe, but probably still in a learning & experimental phase in France ▪ Development of ABMs in France still shy ▪ Various difficulties entertain methods improvements ▪ Heuristic virtues: test hypotheses & models ▪ Integrative power: stimulates collaborations between communities of archaeologists, environmentalists and computer scientists (cf OBRESOC/RDMed projects) ▪ Environmental change still few integrated in most models 23 42e Rencontres Internationales d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022

		III-Assessments & perspectives III-Assessments & perspectives
	Global assessment	Global perspectives
	Strengths:	Perspectives:
	▪ social modeling for gallo-roman society [ModelAnset] ▪ climate-agricultural yields reconstructions	Simulation ▪Development for research in landscape archeology Number of archaeological sites ▪ Hybridizing the strengths of (test hypothesis for the causes of fluctuations ModelAnSet & ROMCLIM/FER-between cultivated/wooded areas) AGRI in a socio-environmental ▪Development for research in the socionatural
	[ROMCLIM/FER-AGRI]	perspective in the long term & transitions between perspective chronoculturals periods
		▪But necessity of strengthen collaboration with
		geoscientifics & paleoenvironmentalists
		(archaeobotany, archaeozoology, geoarcheology,..)
		Archaeodemographical modelling
		(Isoardi 2008)

II-Case studies

▪ Noticeable RCO effect on cereals surplus, but not really on demographical growth ▪ Cooling of more than -2° C have a negative impact on cereal yields and cause a demographical collapse (system resilient under this threshold) 42e Rencontres Internationales d'Archéologie et d'Histoire de Nice Côte d'Azur, 12-14 octobre 2022 22

-> Assessment of our models

Epstein & Axtell 1996Kohler & Gumerman 2000