# A Reverse Stackelberg Model for Demand Response in Local Energy Markets

Juan Sepúlveda, Luce Brotcorne, Hélène Le Cadre

<figure><img src="image_1.png" /></figure>

## To cite this version:

Juan Sepúlveda, Luce Brotcorne, Hélène Le Cadre. A Reverse Stackelberg Model for Demand Response in Local Energy Markets. 2024. hal-04525404

HAL Id: hal-04525404

https://hal.science/hal-04525404

Preprint submitted on 28 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# A Reverse Stackelberg Model for Demand Response in Local Energy Markets

Juan Sepúlveda, Luce Brotcorne, Hélène Le Cadre

<figure><img src="image_1.png" /></figure>

## To cite this version:

Juan Sepúlveda, Luce Brotcorne, Hélène Le Cadre. A Reverse Stackelberg Model for Demand Response in Local Energy Markets. 2024. hal-04525404v2

HAL Id: hal-04525404

https://hal.science/hal-04525404v2

Preprint submitted on 4 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Reverse Stackelberg Model for Demand Response in Local Energy Markets *

Juan Sepúlvedaª, Luce Brotcorneª and Hélène Le Cadreª

ªUniv. Lille, Inria, CNRS, Centrale Lille, UMR 9189 CRISTAL, Lille, F-59000, France

## ARTICLE INFO

### Keywords:

- Time-of-use tariff

- Pricing

- Reverse Stackelberg Game

- Bilevel optimization

- Demand response

- Local electricity markets

- Reactive power

## Abstract

In an era where renewable energy resources are increasingly integrated into our power systems, and consumer-centric approaches gain traction, local energy markets emerge as a pivotal mechanism for empowering prosumers. This paper presents a novel bilevel optimization model that uniquely blends the dynamics of peer-to-peer energy markets with the physical realities of power distribution networks. The innovation steams from introducing a tariff design approach based on affine functions to shape prosumer behavior towards operationally efficient and secure energy exchanges. This is critical as previous market designs often overlooked the physical constraints of power flows, leading to potential risks in voltage regulation and economic efficiency. The lower level of the model encapsulates the interactions among prosumers in a generalized Nash equilibrium problem (GNEP), modeling active and reactive power injections of prosumers. The upper level, representing the role of the distribution system operator, strategically computes tariffs to steer the market to an operationally efficient equilibrium. The paper relies on the classical Nikaido–Isoda (NI) reformulation to characterize the GNEP, a key aspect in leveraging a proof of strong stability of the lower-level solution. Computational experiments on various IEEE test feeder instances reveal the model’s capacity to efficiently align prosumer behavior with operational objectives, utilizing only the tariff information, thereby simplifying the decision-making process in complex distribution systems.

## 1. Introduction

The energy transition represents a path toward reshaping the global energy sector from its fossil-fuel foundations to a carbon-neutral state by the latter half of this century. This journey demands transformative shifts within power systems, particularly at the distribution level. One central aspect involves the widespread adoption of distributed energy resources such as renewable generators, storage systems, adaptable loads, and electric vehicles. Simultaneously, future regulations are expected to pivot towards a consumer-centric model, leveraging local electricity markets and empowering the prosumers to actively manage their energy assets.

While the potential benefits of this transition are substantial, it also presents notable challenges. The uncertain and fast-changing nature of renewable resources, coupled with the autonomy of prosumers, may compromise the efficient and secure operation of the distribution network. To address this concern, the distribution system operators (DSOs) could potentially implement strategies to steer the behavior of prosumers. We believe these strategies must go beyond the control of energy usage patterns and extend to reactive power injections as well as participation in local electricity markets, therefore extending the notion of demand-side management.

The objective of the work is to design an incentive policy in the form of a tariff function to promote prosumer behavior that is both operationally safe and economically efficient. The tariff applied to each prosumer is a function of its net active and

reactive power output, thereby being a price signal that acts on local variables of the system. The tariff design problem considers two classes of decision-makers: (i) the DSO whose interest is to achieve efficient and secure operation of the distribution network, and (ii) the prosumers, selfishly participating in a peer-to-peer market, agnostic of the global system operation. Naturally, we adopt a Stackelberg game setting in which the DSO is the leader (upper level) and the prosumers are the followers (lower level). The main research questions revolve around evaluating the performance of the proposed tariff design specially when compared against conventional constant tariffs.

Recent studies have increasingly explored tariff design as a method for demand-side management, capitalizing on the responsiveness of prosumers to price signals. These approaches are often based on Stackelberg games and bilevel optimization to reflect the complex interplay between the different classes of decision makers involved. Notably, Schittekatte et al. (2018) pioneered the use of a best-response algorithm to achieve a Stackelberg equilibrium between the DSO and the prosumers, a methodology that was expanded by Hoarau and Perez (2019) to encompass electromobility considerations. Pricing techniques based on bilevel optimization have been proposed in the works Askeland et al. (2020, 2021) and Anjos et al. (2021). In an innovative approach, Aussel et al. (2020) developed a customized solution concept for a single-leader multi-follower trilevel problem, considering energy price and quantity as prosumers variables. Similarly, the authors in Tushar et al. (2014) examined a single-leader multi-follower Stackelberg game where the prosumers were modeled by a generalized Nash equilibrium problem. While Maharjan et al. (2016) applied an iterative approach for a multi-leader multi-follower Stackelberg game.

The works mentioned above reveal gaps in incorporating the impact of prosumers' actions on physical variables such

<footnote>* This work emanates from research supported by the ERA-NET Cofund grant under the CHIST-ERA IV Joint Call on Novel Computational Approaches for Environmental Sustainability (CES) project "Supporting Energy Communities-Operational Research and Energy Analytics" (SEC-OREA) and ANR.</footnote>

<footnote>ORCID(s):</footnote>
# A Reverse Stackelberg Model for Demand Response in Local Energy Markets *

Juan Sepúlvedaª, Luce Brotcorneª and Hélène Le Cadreª

ªUniv. Lille, Inria, CNRS, Centrale Lille, UMR 9189 CRISTAL, Lille, F-59000, France

## ARTICLE INFO

**Keywords:**

Time-of-use tariff

Pricing

Reverse Stackelberg Game

Bilevel optimization

Demand response

Local electricity markets

Reactive power

## Abstract

In an era where renewable energy resources are increasingly integrated into our power systems, and consumer-centric approaches gain traction, local energy markets emerge as a pivotal mechanism for empowering prosumers. This paper presents a novel bilevel optimization model that uniquely blends the dynamics of peer-to-peer energy markets with the physical realities of power distribution networks. The innovation steams from introducing a tariff design approach based on affine functions to shape prosumer behavior towards operationally efficient and secure energy exchanges. This is critical as previous market designs often overlooked the physical constraints of power flows, leading to potential risks in voltage regulation and economic efficiency. The lower level of the model encapsulates the interactions among prosumers in a generalized Nash equilibrium problem (GNEP), modeling active and reactive power injections of prosumers. The upper level, representing the role of the distribution system operator, strategically computes tariffs to steer the market to an operationally efficient equilibrium. The paper relies on the classical Nikaido-Isoda (NI) reformulation to characterize the GNEP, a key aspect in leveraging a proof of strong stability of the lower-level solution. Computational experiments on various IEEE test feeder instances reveal the models capacity to efficiently align prosumer behavior with operational objectives, utilizing only the tariff information, thereby simplifying the decision-making process in complex distribution systems.

## 1. Introduction

The energy transition represents a path toward reshaping the global energy sector from its fossil-fuel foundations to a carbon-neutral state by the latter half of this century. This journey demands transformative shifts within power systems, particularly at the distribution level. One central aspect involves the widespread adoption of distributed energy resources such as renewable generators, storage systems, adaptable loads, and electric vehicles. Simultaneously, future regulations are expected to pivot towards a consumer-centric model, leveraging local electricity markets and empowering the prosumers to actively manage their energy assets.

While the potential benefits of this transition are substantial, it also presents notable challenges. The uncertain and fast-changing nature of renewable resources, coupled with the autonomy of prosumers, may compromise the efficient and secure operation of the distribution network. To address this concern, the distribution system operators (DSOs) could potentially implement strategies to steer the behavior of prosumers. We believe these strategies must go beyond the control of energy usage patterns and extend to reactive power injections as well as participation in local electricity markets, therefore extending the notion of demand-side management.

The objective of the work is to design an incentive policy in the form of a tariff function to promote prosumer behavior that is both operationally safe and economically efficient. The tariff applied to each prosumer is a function of its net active and reactive power output, thereby being a price signal that acts on local variables of the system. The tariff design problem considers two classes of decision-makers: (i) the DSO whose interest is to achieve efficient and secure operation of the distribution network, and (ii) the prosumers, selfishly participating in a peer-to-peer market, agnostic of the global system operation. Naturally, we adopt a Stackelberg game setting in which the DSO is the leader (upper level) and the prosumers are the followers (lower level). The main research questions revolve around evaluating the performance of the proposed tariff design specially when compared against conventional constant tariffs.

<footnote>* This work emanates from research supported by the ERA-NET Cofund grant under the CHIST-ERA IV Joint Call on Novel Computational Approaches for Environmental Sustainability (CES) project Supporting Energy Communities-Operational Research and Energy Analytics (SEC-OREA) and ANR.</footnote>

<footnote>ORCID(s):</footnote>
as voltage levels, congestion, and reactive power injections. The integration of distributed energy resources (DERs) and the decentralization of power grids have brought to the forefront the critical role of reactive power support Dall'Anese et al. (2015). DERs are increasingly recognized for their potential to provide reactive power, a capability that is essential for the efficient operation of networks leveraging these resources. Recent literature, including studies like Potter et al. (2023), Bozionek et al. (2022), and Wolgast et al. (2022) underscores the growing emphasis on creating reactive power markets and tariffs designed to address reactive power contributions. We highlight works like Pediaditis et al. (2021), which aim to incorporate physical aspects such as flow constraints and reactive power injection as part of the tariff design problem while considering a peer-to-peer market of energy, illustrating a move towards more comprehensive solutions.

Current literature in tariff design remains in the scope of classical approaches such as Labbé et al. (1998), where the tariff is synthesized as a constant coefficient of a linear expression in the objective function of the prosumers, thereby introducing a linear cost expression. While recognizing state-of-the-art information systems and local energy management technologies Jin et al. (2017). These advancements allow for a reevaluation of conventional assumptions about the intelligibility of tariffs, as the ones in Labbé et al. (1998) and Pediaditis et al. (2021).

As a result, we move beyond constant tariffs based on discrete levels, opting for tariff coefficients that are seen by the prosumers as affine functions of their decision variables. Despite this increased complexity, the focus on simplicity is concretized by ensuring that the decision-making model addressed by prosumers is convex and computationally tractable. This approach finds its origins in works like Zheng and Basar (1982) and is based on the control-theoretic view of incentives presented in Ho et al. (1980). There, under a Stackelberg game setting, the leader calculates optimal cost terms on the objective function of the follower by means of announcing its own reaction function in advance—in our case, a tariff function—thereby called reverse Stackelberg game Groot et al. (2016). Although the finer distinctions within the definitions of the game—theoretic framework are discussed in greater detail in Ho et al. (1981), in the end, the computation of solutions boils down to solving a mathematical program with equilibrium constraints (MPECs).

The MPECs are a class of problems, of which bilevel optimization is a special case, that allows for modeling two layers of decision-making, in which an equilibrium problem (the lower-level problem) is nested into an optimization problem (the upper-level problem); see Luo et al. (1996). They have been applied in tariff design and taxation through modeling single-leader multiple-follower games Hart et al. (2021). The formulation of these programs involves the integration of complementarity conditions into optimization models, enabling the detailed analysis of single-leader multi-follower interactions in competitive markets Raghunathan and Biegler (2003); Flegel and Kanzow (2005); Ye (2005).

Different notions of equilibrium might be chosen to represent the lower-level problem of a MPEC. We focus on the generalized Nash equilibrium problem (GNEP) for its ability to represent a peer-to-peer energy market equilibrium, as done in Le Cadre et al. (2020). This solution concept extends that of the classical Nash equilibrium problem, representing the equilibrium of multiple optimization problems coupled through their objective functions and feasibility sets; detailed reviews can be found in Facchinei and Kanzow (2007) and Fischer et al. (2014). Jointly-convex GNEPs appear as an important family of problems for which there are well-known and useful reformulations: (i) the KKT reformulation, in the case standard constraint qualifications are met, and (ii) the Nikaido–Isoda reformulation; both reformulations are detailed in Facchinei and Kanzow (2007). The Nikaido–Isoda reformulation links jointly convex GNEPs and traditional optimization problems, allowing the analysis of GNEPs at the lower-level of MPECs as parametric optimization problems. Consequently, this allows the application of the notion of stability in parametric optimization, as developed in Dempe (2002), whose importance in an application context has been highlighted in Caruso et al. (2020).

Our paper pioneers the application of reverse Stackelberg games to demand response in power distribution systems, exploring the dynamics of peer-to-peer markets alongside voltage and reactive power management. We address this through a bilevel programming framework with second-order conic constraints that uniquely overcomes the common issue of non-uniqueness in lower-level solutions. Notably, our research is the first within the domain of electricity market demand response to offer a proof of stability for the solutions. Our main contributions are:

1. The development of an inclusive demand-response model that integrates a peer-to-peer market, battery operations, and the AC power flow within the distribution network to ensure both security and efficiency.

2. A strong stability proof of the solution, based on the Nikaido–Isoda reformulation of the lower-level generalized Nash equilibrium problem (GNEP).

Together, these contributions significantly advance demand response strategies, enhancing the functionality and efficiency of local energy markets.

The rest of the paper is organized as follows: Section 2 details the decision-making models for the DSO and prosumers, using a GNEP to model the peer-to-peer market. Moving forward, Section 3 is dedicated to establishing the stability of our proposed solution, profiting from the Nikaido–Isoda reformulation of the lower-level GNEP. Section 4 showcases our numerical findings on instances built on top of standard IEEE distribution grids, illustrating the practicality and effectiveness of our approach. Lastly, Section 5 concludes the paper, summarizing our contributions and suggesting directions for future research.
Recent studies have increasingly explored tariff design as a method for demand-side management, capitalizing on the responsiveness of prosumers to price signals. These approaches are often based on Stackelberg games and bilevel optimization to reflect the complex interplay between the different classes of decision makers involved. Notably, Schittekatte et al. (2018) pioneered the use of a best-response algorithm to achieve a Stackelberg equilibrium between the DSO and the prosumers, a methodology that was expanded by Hoarau and Perez (2019) to encompass electromobility considerations. Pricing techniques based on bilevel optimization have been proposed by Askeland et al. (2020, 2021) and Anjos et al. (2021). Similarly, Tushar et al. (2014) examined a single-leader multi-follower Stackelberg game where the strategies of prosumers were obtained as the solution of a generalized Nash equilibrium problem (GNEP). While Maharjan et al. (2016) applied an iterative approach for a multi-leader multi-follower Stackelberg game. Aussel et al. (2020) propose a trilevel demand-side management model for energy pricing, involving the energy supplier, intermediary agents, and end users. The three decision-making levels are: the supplier optimizing prices, intermediary agents deciding on price and quantity of exchanges, and end users adjusting demand based on the prices they see. Their focus is on the modeling of intermediary agents and on solution selection on the set of generalized Nash equilibria of the second-level. In general, these studies determine the incentive policy as a constant incentive price, whereas our proposed model employs an affine function.

The works mentioned above reveal gaps in incorporating the impact of prosumers' actions on physical variables such as voltage levels, congestion, and reactive power injections. The integration of distributed energy resources (DERs) and the decentralization of power grids have brought to the forefront the critical role of reactive power support (Dall'Anese et al., 2015). DERs are increasingly recognized for their potential to provide reactive power, a capability that is essential for the efficient operation of networks leveraging these resources. Recent literature, including studies by Potter et al. (2023), Bozionek et al. (2022), and Wolgast et al. (2022) underscores the growing emphasis on creating reactive power markets and tariffs designed to address reactive power contributions. We emphasize recent studies that incorporate physical elements, such as flow constraints and reactive power injection, into the tariff design process while accounting for peer-to-peer energy markets, reflecting a shift toward more comprehensive solutions (Pediaditis et al., 2021).

Current literature in tariff design remains in the scope of classical approaches, like the one by Labbé et al. (1998), where the tariff is synthesized as a constant coefficient of a linear expression in the objective function of the prosumers. However, advancements in information systems and local energy management technologies require a critical reassessment of the foundational premises about intelligibility of tariffs (Jin et al., 2017). These technological advancements invite a reevaluation of the assumptions surrounding tariff structures, suggesting the potential for more dynamic and responsive tariff design frameworks that better align with contemporary energy management capabilities.

As a result, we move beyond constant tariffs based on discrete levels, opting for tariff coefficients that are seen by the prosumers as affine functions of their decision variables. Despite this increased complexity, the focus on simplicity is concretized by ensuring that the decision-making model addressed by prosumers is convex and computationally tractable. This approach finds its origins in works by Zheng and Basar (1982) and is based on the control-theoretic view of incentives presented by Ho et al. (1980). There, under a Stackelberg game setting, the leader calculates optimal cost terms on the objective function of the follower by means of announcing its own reaction function in advance—in our case, a tariff function—thereby called reverse Stackelberg game (Groot et al., 2016). Although the finer distinctions within the definitions of the game—theoretic framework are discussed in greater detail by Ho et al. (1981), in the end, the computation of solutions boils down to solving a mathematical program with equilibrium constraints (MPECs).

The MPECs are a class of problems, of which bilevel optimization is a special case, that allows for modeling two layers of decision-making, in which an equilibrium problem (the lower-level problem) is nested into an optimization problem (the upper-level problem) (Luo et al., 1996). They have been applied in tariff design and taxation through modeling single-leader multiple-follower games (Hart et al., 2021). The formulation of these programs involves the integration of complementarity conditions into optimization models, enabling the detailed analysis of single-leader multi-follower interactions in competitive markets (Raghunathan and Biegler, 2003; Flegel and Kanzow, 2005; Ye, 2005).

Different notions of equilibrium might be chosen to represent the lower-level problem of a MPEC. We focus on the generalized Nash equilibrium problem (GNEP) for its ability to represent a peer-to-peer energy market equilibrium, as done by Le Cadre et al. (2020). This solution concept extends that of the classical Nash equilibrium problem, representing the equilibrium of multiple optimization problems coupled through their objective functions and feasibility sets (Facchinei and Kanzow, 2007; Fischer et al., 2014). Jointly-convex GNEPs appear as an important family of problems for which there are well-known and useful reformulations: (i) the KKT reformulation, in the case
## Notation

Boldface letters are used for matrices and vectors, (·)ᵀ for transposition, and |·| denotes the absolute value of a number or the cardinality of a set. Let **A** be a matrix, ||**A**||_F and ||**A**||_2 denote the Frobenius norm and the 2-norm, respectively.

Let $\mathcal{N}$ be the set of players of a game; we use classical game theoretic notation to define $\mathbf{z}_{-i} := (\mathbf{z}_j)_{j \in \mathcal{N} \setminus \{i\}}$, i.e., a vector stacking the entries of $\mathbf{z}$ that correspond to all players except $i$.

## Nomenclature

### Sets and set-valued maps

$\mathcal{N}^0$ Set of prosumers including the slack prosumer ($\mathcal{N}^0 = \mathcal{N} \cup \{0\}$).

$\mathcal{N}$ Set of prosumers.

$\mathcal{K}$ Set of time steps.

$\mathcal{G}_i$ Feasible generator dispatch set of prosumer $i$.

$\Omega_i(\cdot)$ Feasibility set of prosumer $i$. $\Omega_i := \mathbb{R}^{|\mathcal{N}^0|-n_i} \Rightarrow \mathbb{R}^{n_i}$.

### Parameters

**B** Incidence matrix associating prosumers to corresponding buses.

**H** Matrix of linear coefficients associating bus injections to voltage magnitudes.

$\tilde{\mathbf{H}}$ Matrix of linear coefficients associating prosumer injections to voltage magnitudes.

**h** Power flow affine model intercept.

**R** Matrix of resistive linear coefficients of the *LinDisstFlow* model.

**X** Matrix of reactive linear coefficients of the *LinDisstFlow* model.

$E_i$ Connection capacity of prosumer $i$.

$c_{ij}^k$ Bilateral energy trade price from prosumer $i$ to $j$ at time step $k \in \mathcal{K}$, for $(i, j) \in \mathcal{E}$.

$\hat{v}_i, \check{v}_i$ Upper and lower limits of the squared voltage magnitude at bus $i \in \mathcal{B}$, respectively.

$\hat{s}_{G,i}$ Distributed generator capacity of prosumer $i \in \mathcal{N}^0$.

$\hat{p}_{G,i}^k$ Maximum available active power for the distributed generator of prosumer $i$ at time step $k$.

$\check{\theta}_{G,i}$ Minimum power factor of prosumer $i$.

$\eta_{B,i}$ Square root of the battery's round trip efficiency of prosumer $i$.

$\eta_{G,i}$ Loss coefficient of the distributed generator of prosumer $i$.

$\hat{s}_{B,i}$ Capacity of the battery of prosumer $i$.

$\hat{p}_{B,i}$ Nominal power of the battery of prosumer $i$.

$n_i$ Number of decision variables of prosumer $i$.

### Leader variables

$\mathbf{x}$ Vector of leader variables stacking the tariff function coefficients $\mathbf{x} := (\mathbf{\Phi}_i^k, \mathbf{\phi}_i^k)_{i \in \mathcal{N}^0, k \in \mathcal{K}}$.

$\mathbf{\Phi}_i^k$ Matrix in $\mathbb{R}^{2 \times 2}$ representing the affine coefficients of the tariff function assigned to prosumer $i$ at time step $k$.

$\mathbf{\phi}_i^k$ Vector in $\mathbb{R}^2$ representing the constant coefficients of the tariff function assigned to prosumer $i$ at time step $k$.

$\mathbf{\gamma}_i^k$ Affine function of tariff coefficients such that $\mathbf{\gamma}_i^k(\mathbf{u}_i^k) = \mathbf{\Phi}_i^k \mathbf{u}_i^k + \mathbf{\phi}_i^k$.

### Followers variables

$\mathbf{u}_i^k$ Vector of variables of prosumer $i$ at time step $k$.

$p_i^k$ Net active power injection of prosumer $i$ at time step $k$.

$q_i^k$ Net reactive power injection of prosumer $i$ at time step $k$.

$p_{G,i}^k$ Active power generated by prosumer $i$ at time step $k$.

$q_{G,i}^k$ Reactive power generated by prosumer $i$ at time step $k$.

$p_{E,ij}^k$ Energy offered by prosumer $i$ to prosumer $j$ at time step $k$.

## 2. Problem formulation

The proposed Stackelberg model involves two distinct types of agents: (i) prosumers and (ii) the DSO. Prosumers participate in a peer-to-peer market, where their interactions are driven by self-interest, potentially causing adverse effects on the physical integrity of the underlying electrical grid. To address this challenge, we consider a DSO with the capability to enforce economic incentives aimed at enhancing the grid's operational performance.

Let $(\mathcal{N}^0, \mathcal{E})$ be a symmetric simple directed graph modeling the financial network, where $\mathcal{N}^0 := \mathcal{N} \cup \{0\}$, $\mathcal{N}$ is the set of prosumers, and $\mathcal{E}$ is the set of edges representing possible peer-to-peer interactions. We define prosumer 0 as a slack prosumer representing the external grid at the substation. As usual, the distribution network is a directed rooted tree with $\mathcal{B}^0$ being the set of nodes (buses). The substation node (root) is labeled as 0, and $\mathcal{B} := \mathcal{B}^0 \setminus \{0\}$, where $\mathcal{B}$ is the set of buses that may host prosumers. We define a map $i \mapsto A_i$, from $\mathcal{B}$ to $\mathcal{B}^0$ that returns the parent of each node.

The link between the market network and the electrical grid is illustrated in Figure 1. Each prosumer is connected to exactly one bus. We define a mapping $\mathcal{B} : \mathcal{N}^0 \to \mathcal{B}^0$ that assigns the hosting bus to each prosumer, thereby $\mathcal{B}(2) = 1$ in the example shown in Figure 1. The slack prosumer (external grid) is always connected to the substation bus; therefore, $\mathcal{B}(0) = 0$ is always true.

In the following sections, we will outline the models that capture the decision-making processes of the two classes of players: the prosumer and the DSO.
standard constraint qualifications are met, and (ii) the Nikaido–Isoda reformulation; both reformulations are detailed
by Facchinei and Kanzow (2007). The Nikaido–Isoda reformulation links jointly convex GNEPs and traditional
optimization problems, allowing the analysis of GNEPs at the lower–level of MPECs as parametric optimization
problems. Consequently, this allows the application of the notion of stability in parametric optimization, as developed
by Dempe (2002), whose importance in an application context has been highlighted by Caruso et al. (2020).

Our paper pioneers the application of reverse Stackelberg games to demand response in power distribution systems, exploring the dynamics of peer-to-peer markets alongside voltage and reactive power management. We address this through a bilevel programming framework with second-order conic constraints that uniquely overcomes the common issue of non-uniqueness in lower-level solutions. Notably, our research is the first within the domain of electricity market demand response to offer a proof of stability for the solutions. Our main contributions are:

1. The development of an inclusive demand-response model that integrates a peer-to-peer market, battery operations, and the AC power flow within the distribution network to ensure both security and efficiency.

2. A strong stability proof of the solution, based on the Nikaido–Isoda reformulation of the lower–level generalized Nash equilibrium problem (GNEP).

Together, these contributions significantly advance demand response strategies, enhancing the functionality and efficiency of local energy markets.

The rest of the paper is organized as follows: Section 2 details the decision–making models for the DSO and prosumers, using a GNEP to model the peer–to–peer market. Moving forward, Section 3 is dedicated to establishing the stability of our proposed solution, profiting from the Nikaido–Isoda reformulation of the lower–level GNEP. Section 4 showcases our numerical findings on instances built on top of standard IEEE distribution grids, illustrating the practicality and effectiveness of our approach. Lastly, Section 5 concludes the paper, summarizing our contributions and suggesting directions for future research.

**Notation**

Boldface letters are used for matrices and vectors, (·)ᵀ for transposition, and |·| denotes the absolute value of a number or the cardinality of a set. Let **A** be a matrix, ||**A**||_F and ||**A**||_2 denote the Frobenius norm and the 2-norm, respectively.

Let $\mathcal{N}$ be the set of players of a discrete-time dynamic game; we use classical game theoretic notation to define $\mathbf{u}_{-i} := (\mathbf{z}_j)_{j \in \mathcal{N} \setminus \{i\}}$, i.e., a vector stacking the entries of $\mathbf{u}$ that correspond to all players except $i$. Superindices denote the discrete time index, e.g., $\mathbf{u}_i^k$ is the vector of decision variables associated to player $i$ at time $k$.

**Nomenclature**

Sets and set-valued maps

$\mathcal{N}^0$ Set of prosumers including the slack prosumer ($\mathcal{N}^0 = \mathcal{N} \cup \{0\}$).

$\mathcal{N}$ Set of prosumers.

$\mathcal{K}$ Set of time steps.

$\mathcal{G}_i$ Feasible generator dispatch set of prosumer $i$.

$\Omega_i(\cdot)$ Feasibility set of prosumer $i$. $\Omega_i : \mathbb{R}^{|\mathcal{N}|-n_i} \Rightarrow \mathbb{R}^{n_i}$.

Parameters

**B** Incidence matrix associating prosumers to corresponding buses.

**H** Matrix of linear coefficients associating bus injections to voltage magnitudes.

**\tilde{H}** Matrix of linear coefficients associating prosumer injections to voltage magnitudes.

**h** Power flow affine model intercept.

**R** Matrix of resistive linear coefficients of the *LinDistFlow* model.

**X** Matrix of reactive linear coefficients of the *LinDistFlow* model.

$E_i$ Connection capacity of prosumer $i$.

$c_{ij}^k$ Bilateral energy trade price from prosumer $i$ to $j$ at time step $k \in \mathcal{K}$, for $(i, j) \in \mathcal{E}$.
<figure><img src="image_2.png" /><figcaption>Figure 1: The financial network and the supporting electrical grid.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Distributed generator capacity set (depicted in purple) for different operating modes in the active/reactive power plane.</figcaption></figure>

## 2.1. The prosumer model

Each prosumer owns and manages a single photovoltaic generator and a single battery, and may buy or sell energy in a peer-to-peer fashion with its neighbors in ($\mathcal{N}$, $\mathcal{E}$). The prosumer schedules its decisions in order to maximize their utility function. To do that, each prosumer solves a dynamic optimization problem over a finite horizon with time step indices $k \in \mathcal{K} := \{1, ..., K\}$.

To begin, we establish a model for the production capacity of individual prosumers. In this context, we assume that each prosumer possesses a photovoltaic distributed generator, with $p_{G,i}^k$ and $q_{G,i}^k$ representing the active and reactive power output of prosumer $i \in \mathcal{N}$ at time step $k \in \mathcal{K}$. The feasible dispatch set in the active and reactive power space is modeled as follows

$$C_i^k := \{(p_{G,i}, q_{G,i}) \in \mathbb{R}^2 : \qquad (1a)$$

$$\check{p}_{G,i}^k \le p_{G,i} \le \hat{p}_{G,i}^k \qquad (1b)$$

$$p_{G,i}^2 + q_{G,i}^2 \le s_{G,i}^2 \qquad (1c)$$

$$q_{G,i} \le \check{\theta}_i p_{G,i} \qquad (1d)$$

The definition provided in (1) comprehensively encompasses various standard operating modes, as detailed in Dall'Anese et al. (2015). These modes are illustrated in Figure 2, each corresponding to distinct selections of the parameters $\check{p}_{G,i}^k$ and $\theta_i$. Specifically, the modes include the Maximum Power Point (MPP), Reactive Power Control (RPC), Active Power Control (APC), and Optimal Inverter Dispatch (OID).

Apart from the distributed generator, each prosumer owns a battery, modeled as a 1-dimensional dynamic system such that the state of charge $s_{B,i}^k$ satisfies

$$s_{B,i}^k = s_{B,i}^{k-1} + \eta_{B,i} p_{B,i}^{+,k} - \frac{1}{\eta_{B,i}} p_{B,i}^{-,k}, \quad \forall k \in \mathcal{K} \qquad (2)$$

where $p_{B,i}^k = p_{B,i}^{+,k} - p_{B,i}^{-,k}$ is the energy absorbed by the battery of prosumer $i$, over time step $k$, and $p_{B,i}^{+,k} \ge 0$ and $p_{B,i}^{-,k} \ge 0$ are auxiliary variables that help modeling the battery energy losses in the charging and discharging cycle. The positive parameter $\eta_{B,i}$ is the square root of the round cycle battery efficiency. In addition, we impose the following capacity bounds

$$0 \le s_{B,i}^k \le \hat{s}_{B,i}, \quad \forall k \in \mathcal{K} \qquad (3)$$

$$-\hat{p}_{B,i} \le p_{B,i}^k \le \hat{p}_{B,i}, \quad \forall k \in \mathcal{K} \qquad (4)$$

Having defined the models for the distributed generator and battery, we now delve into the elements that shape the preferences of prosumers. We assume there are three components: (i) the tariff/incentive cost terms, (ii) the income obtained via peer-to-peer exchanges, and (iii) the distributed generator losses. Combining the three components results in the following prosumer disutility function

$$J_{F,i}(T(\cdot), \mathbf{u}_i, \mathbf{u}_{-i}) = \sum_{k \in \mathcal{K}} \left( T_i^k(p_i^k, q_i^k) - U_i^k \left( \mathbf{p}_{E,i}^k, (p_{E,ji}^k)_{j \in \mathcal{N}_i^0} \right) + l_i^k(p_{G,i}^k, q_{G,i}^k) \right) \qquad (5)$$

where $\mathbf{u}_i^k := (p_i^k, q_i^k, p_{G,i}^k, q_{G,i}^k, p_{B,i}^{+,k}, p_{B,i}^{-,k}, s_{B,i}^k, \mathbf{p}_{E,i}^k)$ comprises the decision variables of prosumer $i \in \mathcal{N}$ at time step $k \in \mathcal{K}$, with $\mathbf{p}_{E,i}^k := (p_{E,0i}^k, (p_{E,ij}^k)_{j \in \mathcal{N}_i^0})$. The function $T_i^k(\cdot)$ is the tariff/incentive function, $U_i^k(\cdot)$ characterizes the gain obtained in the peer-to-peer market, and $l_i^k(\cdot)$ is the distributed generator losses. Note that minimizing the disutility function (5) is the same as maximizing the utility.

The profit from peer-to-peer energy transactions can be expressed as follows:

$$U_i^k \left( \mathbf{p}_{E,i}^k, (p_{E,ji}^k)_{j \in \mathcal{N}_i^0} \right) := \sum_{j \in \mathcal{N}_i^0} c_{ij}^k p_{E,ij}^k - c_{ji}^k p_{E,ji}^k \qquad (6)$$

where $c_{ij}^k$ is the price of energy sold from prosumer $i$ to prosumer $j$ at time step $k$, $p_{E,ij}^k \ge 0$ represents the amount of energy that prosumer $i$ is willing to offer to prosumer $j$. Hence, the actual energy transfer from $i$ to $j$, is given by

$$\Delta p_{E,ij}^k = p_{E,ij}^k - p_{E,ji}^k, \quad \forall j \in \mathcal{N}_i^0 \qquad (7)$$

which is a joint decision between prosumer $i$ and prosumer $j$. Note that the sign constraint forcing $p_{E,ij}^k$ to be non-negative allows modeling the difference between buying and selling price associated with the transfer $\Delta p_{E,ij}^k$.

We introduce a 1-norm regularization term that accounts for the losses of the photovoltaic generator. It is defined as

$$l_i^k(p_{G,i}^k, q_{G,i}^k) := \eta_{G,i} (p_{G,i}^k + |q_{G,i}^k|), \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (8)$$

This captures the preference of prosumers in favor of not dispatching active or reactive power unless it benefits them, whether through incentives, market conditions, or self-consumption.
$\hat{v}_i$, $\tilde{v}_i$ Upper and lower limits of the squared voltage magnitude at bus $i \in B$, respectively.

$\hat{s}_{G,i}$ Distributed generator capacity of prosumer $i \in \mathcal{N}$.

$\hat{p}_{G,i}^k$ Maximum available active power for the distributed generator of prosumer $i$ at time step $k$.

$\tilde{v}_{G,i}$ Minimum power factor of prosumer $i$.

$\eta_{B,i}$ Square root of the battery's round trip efficiency of prosumer $i$.

$\eta_{G,i}$ Loss coefficient of the distributed generator of prosumer $i$.

$\hat{s}_{B,i}$ Capacity of the battery of prosumer $i$.

$\hat{p}_{B,i}$ Nominal power of the battery of prosumer $i$.

$n_i$ Number of decision variables of prosumer $i$.

## Leader variables

$\mathbf{x}$ Vector of leader variables stacking the tariff function coefficients $\mathbf{x} := (\mathbf{\Phi}_i^k, \mathbf{\phi}_i^k)_{i \in \mathcal{N}', k \in \mathcal{K}}$.

$\mathbf{\Phi}_i^k$ Matrix in $\mathbb{R}^{2 \times 2}$ representing the affine coefficients of the tariff function assigned to prosumer $i$ at time step $k$.

$\mathbf{\phi}_i^k$ Vector in $\mathbb{R}^2$ representing the constant coefficients of the tariff function assigned to prosumer $i$ at time step $k$.

$\gamma_i^k$ Affine function of tariff coefficients such that $\gamma_i^k(\mathbf{u}_i^k) = \mathbf{\Phi}_i^k \mathbf{u}_i^k + \mathbf{\phi}_i^k$.

## Followers variables

$\mathbf{u}_i^k$ Vector of variables of prosumer $i$ at time step $k$.

$p_i^k$ Net active power injection of prosumer $i$ at time step $k$.

$q_i^k$ Net reactive power injection of prosumer $i$ at time step $k$.

$p_{G,i}^k$ Active power generated by prosumer $i$ at time step $k$.

$q_{G,i}^k$ Reactive power generated by prosumer $i$ at time step $k$.

$p_{E,ij}^k$ Energy offered by prosumer $i$ to prosumer $j$ at time step $k$.

$s_{B,i}^k$ State of charge of the battery of prosumer $i$ at time step $k$.

$p_{B,i}^k$ Absorbed power by the battery of prosumer $i$ at time step $k$.

## 2. Problem formulation

The proposed Stackelberg model involves two distinct types of agents: (i) prosumers and (ii) the DSO. Prosumers participate in a peer-to-peer market, where their interactions are driven by self-interest, potentially causing adverse effects on the physical integrity of the underlying electrical grid. To address this challenge, we consider a DSO with the capability to enforce economic incentives aimed at enhancing the grid's operational performance.

The proposed incentive mechanism is designed to be implemented likewise to how current grid usage tariffs are enforced. In practice, once users are connected to the grid, they are automatically subject to the current tariff rules. Similarly, in our model, prosumers (followers) are bound by the incentive mechanism as soon as they are connected to the grid. The incentive functions depend on their grid usage. Prosumers are also assumed to have prior knowledge of these incentive rules, allowing them to adjust their behavior strategically. In terms of regulation, the proposed approach would rely on similar frameworks already in place to support traditional grid usage tariffs at the distribution network level. Also, the information necessary to compute the incentives is already collected by smart meters in practice.

Let $(\mathcal{N}^0, \mathcal{E})$ be a symmetric simple directed graph modeling the financial network, where $\mathcal{N}^0 := \mathcal{N} \cup \{0\}$, $\mathcal{N}$ is the set of prosumers, and $\mathcal{E}$ is the set of edges representing possible peer-to-peer interactions. We define prosumer 0 as a slack prosumer representing the external grid at the substation. As usual, the distribution network is a directed rooted tree with $\mathcal{B}^0$ being the set of nodes (buses). The substation node (root) is labeled as 0, and $\mathcal{B} := \mathcal{B}^0 \setminus \{0\}$,
In the following, we present the optimization problem to be solved by each prosumer

$$ \min_{\mathbf{u}_i} J_{F,i}(T(\cdot), \mathbf{u}_i, \mathbf{u}_{-i}) \qquad (9a) $$

$$ \text{s.t.} \quad p_i^k = \sum_{j \in \mathcal{N}_i^0} (p_{E,ij}^k - p_{E,ji}^k), \quad \forall k \in \mathcal{K} \quad (\mu_i^k) \quad (9b) $$

$$ p_i^k = p_{G,i}^k - p_{B,i}^k - p_{D,i}^k, \quad \forall k \in \mathcal{K} \qquad (v_{P,i}^k) \qquad (9c) $$

$$ q_i^k = q_{G,i}^k - q_{D,i}^k, \quad \forall k \in \mathcal{K} \qquad (v_{Q,i}^k) \qquad (9d) $$

$$ s_{B,i}^k = s_{B,i}^{k-1} + \eta_{B,i} p_{B,i}^{+,k} - \frac{p_{B,i}^{-,k}}{\eta_{B,i}}, \quad \forall k \in \mathcal{K} \quad (v_{B,i}^k) \quad (9e) $$

$$ (p_{G,i}^k, q_{G,i}^k) \in \mathcal{G}_i^k, \quad \forall k \in \mathcal{K} \qquad (\lambda_{G,i}^k) \qquad (9f) $$

$$ 0 \le s_{B,i}^k \le \hat{s}_{B,i}, \quad \forall k \in \mathcal{K} \qquad (\lambda_{SL,i}^k, \lambda_{SU,i}^k) \qquad (9g) $$

$$ 0 \le p_{B,i}^{+,k} \le \hat{p}_{B,i}, \quad \forall k \in \mathcal{K} \qquad (\lambda_{BL,i}^{+,k}, \lambda_{BU,i}^{+,k}) \qquad (9h) $$

$$ 0 \le p_{B,i}^{-,k} \le \hat{p}_{B,i}, \quad \forall k \in \mathcal{K} \qquad (\lambda_{BL,i}^{-,k}, \lambda_{BU,i}^{-,k}) \qquad (9i) $$

$$ 0 \le p_{E,ij}^k, \quad \forall (j, k) \in \mathcal{N}_i^0 \times \mathcal{K} \qquad (\lambda_{EL,ij}^k) \qquad (9j) $$

where (9b) represents the balance between the net active power injection and the peer-to-peer market trade, constraints (9c) and (9d) are the active and reactive power balance with parameters $p_{D,i}^k$ and $q_{D,i}^k$ representing the active and reactive power demand, respectively. The dynamics of the battery are modeled via (9e), with capacity constraints (9g), (9h), and (9i). Finally, the sign constraints imposed on the bilateral offers are enforced through (9j).

The computational complexity of (9) drastically depends on the selected $T(\cdot)$. If $T(\cdot)$ is a convex quadratic function, problem (9) forms a convex quadratic program with second-order cone constraints.

**Remark 1.** $p_{E,0i}^k$ is a decision variable of prosumer i representing the energy they are willing to buy from the DSO (slack prosumer indexed by 0) at time step k. Therefore the energy transfer $\Delta p_{E,i0}^k = p_{E,i0}^k - p_{E,0i}^k$ is a decision of prosumer i. In fact, we could replace $(p_{E,i0}^k - p_{E,0i}^k)$ by $\Delta p_{E,i0}^k$ whenever it appears in problem (9). As a consequence, we have that if T is convex, then (9) is always feasible because prosumer can always satisfy their power consumption needs as $\Delta p_{E,i0}^k$ is assumed unbounded.

## 2.2. The peer-to-peer market model

The objective of this section is to leverage a model for the peer-to-peer market equilibrium. We consider two relevant classes of market designs: with and without a local market operator (LMO). The corresponding solution concepts representing the output of the market in each of the cases mentioned above are: (i) generalized Nash equilibria (GNEs) Facchinei and Kanzow (2007) and (ii) welfare optima. In the following paragraphs, we characterize the two solution concepts. Then, in Appendix A, we establish explicit connections between them in terms of dual variables and parameters of the prosumer model.

### 2.2.1. The set of GNEs

We assume that if there is no LMO in place and prosumers are free to trade with each other, then the set of market equilibria is characterized by the GNEs. In the context of this work, a GNE is a solution where no prosumer can improve their objective function by unilaterally changing their decision, given the decisions of all other prosumers. Prosumers behave selfishly and independently, minimizing their objective function under local and coupling constraints, based on (10). It extends the classical Nash equilibrium by also considering coupling constraints between the players.

The optimization problem (9) is coupled to the decision making of the other prosumers by the objective function (9a) and constraints (9b), since $\mathbf{p}_{E,ji}$, with $j \in \mathcal{N}_i^j$, is a vector of parameters in prosumer i's decision problem, and simultaneously a vector of variables for prosumer j. We explicitly denote the coupling between the optimization problems by writing (9) in compact form as

$$ \min_{\mathbf{u}_{F,i} \in \Omega_{F,j}(\mathbf{u}_{F,-i})} J_{F,i}(T_i(\cdot), \mathbf{u}_{F,i}, \mathbf{u}_{F,-i}) \qquad (10) $$

where $\mathbf{u}_{F,i} := (\mathbf{u}_{F,i}^k)_{k \in \mathcal{K}}$ is the action profile of prosumer i over the whole planning horizon, and $\Omega_{F,i}(\cdot)$ is a point-to-set mapping—also referred to as set-valued map (see (Rockafellar and Wets, 1998, Chapter 5)),—representing the feasible set of prosumer i as a function of the action profiles of the other prosumers.

More precisely, consider a game with a set $\mathcal{N}$ of players, and the tuple $(\Omega_{F,i}, J_{F,i}, T_i)_{i \in \mathcal{N}}$ as the data.

**Definition 1.** The set of GNEs, denoted by $S(T)$, is the set of elements $\bar{\mathbf{u}} \in \Pi_{i \in \mathcal{N}} \mathbb{R}^{n_i}$ such that

$$ \bar{\mathbf{u}}_i \in \arg \min_{\mathbf{u}_i \in \Omega_i(\mathbf{u}_{-i})} J_{F,i}(T_i(\cdot), \mathbf{u}_i, \mathbf{u}_{-i}), \quad \forall i \in \mathcal{N} \qquad (11) $$

The concept of GNEs is also referred to as pseudo-Nash equilibria; see Harker (1991). In our case, the set $S(T)$ can be interpreted as modeling the behavior of prosumers uncoordinately making bilateral contracts of energy, subject to tariff T. The fundamental task of this work is to design T such that all elements in $S(T)$ represent operationally secure and efficient action profiles. A representation of $S(T)$ in terms of a system of equations is studied in Appendix A for convex quadratic T.

### 2.2.2. The set of welfare optima

On the other hand, when there is a market clearing mechanism driving the system to an economically efficient equilibrium, i.e., one that maximizes the social welfare of all prosumers or, equivalently, minimizes the total disutility. The following definition captures the corresponding solution concept

**Definition 2.** The set of economically efficient market outputs, denoted by $W(T)$, is the set of elements $\bar{\mathbf{u}} \in \Pi_{i \in \mathcal{N}} \mathbb{R}^{n_i}$ such that

$$ \bar{\mathbf{u}} \in \arg \min_{\mathbf{u} \in \Omega} \sum_{i \in \mathcal{N}} J_{F,i}(T_i(\cdot), \mathbf{u}_i, \mathbf{u}_{-i}) \qquad (12) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: The financial network and the supporting electrical grid.</figcaption></figure>

where $\mathcal{B}$ is the set of buses that may host prosumers. We define a map $i \mapsto A_i$, from $\mathcal{B}$ to $\mathcal{B}^0$ that returns the parent of each node. The distribution network graph is defined to have a tree structure as Guerrero et al. (2019), though other types of distribution networks can also be considered; see the discussion in Section 2.3.

The link between the market network and the electrical grid is illustrated in Figure 1. Each prosumer is connected to exactly one bus. We define a mapping $\mathcal{B} : \mathcal{N}^0 \to \mathcal{B}^0$ that assigns the hosting bus to each prosumer, thereby $\mathcal{B}(2) = 1$ in the example shown in Figure 1. The slack prosumer (external grid) is always connected to the substation bus; therefore, $\mathcal{B}(0) = 0$ is always true.

In the following sections, we will outline the models that capture the decision-making processes of the two classes of players: the prosumer and the DSO.

## 2.1. The prosumer model

Prosumers are expected to solve a dynamic optimization problem considering a planning horizon, which is crucial for effective battery management. Advanced energy management systems for device operation and automated trading utilize predictions of future energy prices and consumption (Finnah et al., 2022; Jin et al., 2017). In our model, each prosumer owns and manages a single photovoltaic generator and a single battery and may buy or sell energy in a peer-to-peer fashion with its neighbors in $(\mathcal{N}, \mathcal{E})$. To maximize their utility function, each prosumer schedules decisions by solving a dynamic optimization problem over a finite horizon with time step indices $k \in \mathcal{K} := \{1, ..., K\}$. In practice, this involves a rolling horizon optimization strategy, where decisions are planned several steps ahead but implemented incrementally, allowing for real-time adjustments based on new information or market changes.

To begin, we establish a model for the production capacity of individual prosumers. In this context, we assume that each prosumer possesses a photovoltaic distributed generator, with $p_{G,i}^k$ and $q_{G,i}^k$ representing the active and reactive power output of prosumer $i \in \mathcal{N}$ at time step $k \in \mathcal{K}$. The feasible dispatch set in the active and reactive power space is modeled as follows

$$ G_i^k := \{(p_{G,i}^k, q_{G,i}^k) \in \mathbb{R}^2 : \qquad (1a) $$

$$ \check{p}_{G,i} \le p_{G,i}^k \le \hat{p}_{G,i}^k \qquad (1b) $$

$$ (p_{G,i}^k)^2 + (q_{G,i}^k)^2 \le \hat{s}_{G,i}^2 \qquad (1c) $$

$$ q_{G,i}^k \le \theta_i p_{G,i}^k \qquad (1d) $$

The definition provided in (1) comprehensively encompasses various standard operating modes, as detailed in Dall'Anese et al. (2015). These modes are illustrated in Figure 2, each corresponding to distinct selections of the parameters $\check{p}_{G,i}^k$ and $\theta_i$. Specifically, the modes include the Maximum Power Point (MPP), Reactive Power Control (RPC), Active Power Control (APC), and Optimal Inverter Dispatch (OID).

Apart from the distributed generator, each prosumer owns a battery, modeled as a 1-dimensional dynamic system such that the state of charge $s_{B,i}^k$ satisfies

$$ s_{B,i}^k = s_{B,i}^{k-1} + \eta_{B,i} p_{B,i}^{+,k} - \frac{1}{\eta_{B,i}} p_{B,i}^{-,k}, \quad \forall k \in \mathcal{K} \qquad (2) $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Electric branch associated with bus $i \in B$.</figcaption></figure>

We can interpret (12) as representing the decision-making of all prosumers jointly minimizing the total disutility.

**Remark 2.** It is natural to ask if there exist conditions under which individual prosumers, if acting selfishly, might want to deviate from $W(T)$. This is the same as asking if there exists an instance of the problem for which $\bar{u}_F \in W(T) \setminus S(T)$. Similarly, we may pose questions regarding the reciprocal difference $S(T) \setminus W(T)$, i.e., if there exists an output of a market without LMO that is not economically efficient. Following a similar framework as in Le Cadre et al. (2020), in Appendix A, we derive explicit expressions characterizing the relationship between $S$, $W$ in terms of dual variables and parameters of (9).

## 2.3. The effect of prosumers on the electrical grid

This section focuses on modeling the impact of prosumers' behavior on the operation of the distribution system. Specifically, we aim to establish the relationship between the net active and reactive power injected by prosumers and the resulting changes in voltage levels and system losses. We drop the time-step indices unless said otherwise because all relations presented here are static.

The optimal power flow (OPF) model is a non-convex quadratically constrained quadratic problem (QCQP) which is NP-Hard in general. Because of the difficulty of solving the OPF, several approximations have been made. Between them, convex relaxations are of the families: SDP, SOCP, and linear; see Low (2014a). SDP relaxations present the difficulty of extracting solutions if the rank constraint is not met. SOCP relaxations are not exact when voltage upper bound constraints are active Low (2014b). Therefore, we believe linear approximations are better suited for situations with high penetration of renewable energy resources since the system is exposed to overvoltage conditions. Among the linear OPF models, we highlight the linear OPF model presented in Coffrin and Van Hentenryck (2014), the two models for unbalanced networks presented in Bernstein et al. (2018), and the widely referenced LinDistFlow model proposed in Baran and Wu (1989). In the proposed work, we adopt the LinDistFlow model which is detailed in the following.

For each bus $i \in B$, the following equations hold

$$ \tilde{P}_i = \sum_{j \in C_i} \tilde{P}_j - \tilde{P}_i \quad (13a) $$

$$ \tilde{Q}_i = \sum_{j \in C_i} \tilde{Q}_j - \tilde{Q}_i \quad (13b) $$

$$ v_{A_i} = v_i - 2\tilde{r}_i \tilde{P}_i - 2\tilde{x}_i \tilde{Q}_i \quad (13c) $$

Here, $\tilde{P}_i$ and $\tilde{Q}_i$ represent the active and reactive power flows, respectively, on the upstream branch of bus $i \in B$, measured in the downstream direction. The parameters $\tilde{r}_i$ and $\tilde{x}_i$ denote the resistance and reactance of that branch, while $v_i$ corresponds to the square of the voltage magnitude at bus $i \in B$. The set $C_i$ represents the set of successor buses connected to $i$. A summary of the notation is provided in Figure 3. The LinDistFlow is limited to radial (tree graph) networks with balanced voltage and currents. If the reader requires a linear model capturing a more general setting, such as unbalanced and/or non-radial networks, they may refer to the first-order Taylor linearization in Bernstein et al. (2018).

To calculate the system losses, define matrices $\tilde{\mathbf{R}} \in \mathbb{R}^{|B| \times |B|}$ and $\tilde{\mathbf{X}} \in \mathbb{R}^{|B| \times |B|}$ with coefficients

$$ \tilde{\mathbf{R}}_{ij} := \sum_{k \in \mathcal{L}(i) \cap \mathcal{L}(j)} \tilde{r}_k, \quad \tilde{\mathbf{X}}_{ji} := \sum_{k \in \mathcal{L}(i) \cap \mathcal{L}(j)} \tilde{x}_k \quad (14) $$

where $\mathcal{L}(i)$ is the set of lines between bus $i$ and the slack bus 0, including $i$. Both $\tilde{\mathbf{R}}$ and $\tilde{\mathbf{X}}$ are positive definite matrices, and the losses of the system can be expressed as a convex quadratic function

$$ L(\tilde{\mathbf{p}}, \tilde{\mathbf{q}}) = \tilde{\mathbf{p}}^\top \tilde{\mathbf{R}} \tilde{\mathbf{p}}^k + \tilde{\mathbf{q}}^\top \tilde{\mathbf{R}} \tilde{\mathbf{q}} \quad (15) $$

The operational security constraints are voltage magnitude bounds at each bus. Recall that in (13), our proxy for the voltage magnitudes are the voltage magnitudes squared. Therefore, we express the security constraints limits as the following linear constraints

$$ \hat{v}_i^2 \le v_i \le \hat{v}_i^2, \quad \forall i \in B \quad (16) $$

Now, it is necessary to translate the previous expressions into prosumer related-quantities. The correspondence between prosumers and buses can be described by a relation $B \times \mathcal{N} \rightarrow \{0, 1\}$, represented by the incidence matrix $\mathbf{B} \in \{0, 1\}^{|B| \times |\mathcal{N}|}$, such that

$$ B_{ij} := \begin{cases} 1, & \text{if prosumer } j \text{ is connected to bus } i \\ 0, & \text{otherwise} \end{cases} \quad (17) $$

Let us consider again the time step indices and obtain the voltage magnitude and losses in terms of the prosumers' injections. Recall that $\mathbf{p}^k := (p_i^k)_{i \in \mathcal{N}}$ and $\mathbf{q}^k := (q_i^k)_{i \in \mathcal{N}}$, for all $k \in \mathcal{K}$, and aggregate the effect of prosumers by bus by applying the following transformation based on (17)

$$ \tilde{\mathbf{p}}^k = \mathbf{B} \mathbf{p}^k \quad (18a) $$

$$ \tilde{\mathbf{q}}^k = \mathbf{B} \mathbf{q}^k \quad (18b) $$

Such that

$$ \mathbf{R} := \mathbf{B}^\top \tilde{\mathbf{R}} \mathbf{B} \quad (19a) $$

$$ \mathbf{X} := \mathbf{B}^\top \tilde{\mathbf{X}} \mathbf{B} \quad (19b) $$

The linear transformation $\mathbf{B}$ preserves positive definiteness; therefore, again, we obtain a convex quadratic function for the
<figure><img src="image_4.png" /><figcaption>**Figure 2:** Distributed generator capacity set (depicted in purple) for different operating modes in the active/reactive power plane.</figcaption></figure>

where $p_{B,i}^k = p_{B,i}^{+,k} - p_{B,i}^{-,k}$ is the energy absorbed by the battery of prosumer $i$, over time step $k$, and $p_{B,i}^{+,k} \ge 0$ and $p_{B,i}^{-,k} \ge 0$ are auxiliary variables that help modeling the battery energy losses in the charging and discharging cycle. The positive parameter $\eta_{B,i}$ is the square root of the round cycle battery efficiency. In addition, we impose the following capacity bounds

$$0 \le s_{B,i}^k \le \hat{s}_{B,i}, \forall k \in \mathcal{K} \qquad (3)$$

$$-\hat{p}_{B,i} \le p_{B,i}^k \le \hat{p}_{B,i}, \forall k \in \mathcal{K} \qquad (4)$$

Having defined the models for the distributed generator and battery, we now delve into the elements that shape the preferences of prosumers. We assume there are three components: (i) the tariff/incentive cost terms, (ii) the income obtained via peer-to-peer exchanges, and (iii) the distributed generator losses. Combining the three components results in the following prosumer disutility function

$$J_{F,i}(T_i(\cdot), \mathbf{u}_i, \mathbf{u}_{-i}) := \sum_{k \in \mathcal{K}} \left( T_i^k(p_i^k, q_i^k) - U_i^k(\mathbf{p}_{E,i}^k, (p_{E,ji}^k)_{j \in \mathcal{N}_i^0}) + l_i^k(p_{G,i}^k, q_{G,i}^k) \right) \qquad (5)$$

where $\mathbf{u}_i^k := (p_i^k, q_i^k, p_{G,i}^k, q_{G,i}^k, p_{B,i}^{+,k}, p_{B,i}^{-,k}, s_{B,i}^k, \mathbf{p}_{E,i}^k)$ comprises the decision variables of prosumer $i \in \mathcal{N}$ at time step $k \in \mathcal{K}$, with $\mathbf{p}_{E,i}^k := (p_{E,0i}^k, (p_{E,ij}^k)_{j \in \mathcal{N}_i^0})$. The function $T_i^k(\cdot)$ is the tariff/incentive function, $U_i^k(\cdot)$ characterizes the gain obtained in the peer-to-peer market, and $l_i^k(\cdot)$ is the distributed generator losses. Note that minimizing the disutility function (5) is the same as maximizing the utility.

The profit from peer-to-peer energy transactions can be expressed as follows:

$$U_i^k(\mathbf{p}_{E,i}^k, (p_{E,ji}^k)_{j \in \mathcal{N}_i^0}) := \sum_{j \in \mathcal{N}_i^0} c_{ij}^k p_{E,ij}^k - c_{ji}^k p_{E,ji}^k \qquad (6)$$

where $\mathcal{N}_i^0$ represents the set of prosumers that are connected to prosumer $i$, including the slack prosumer; parameter $c_{ij}^k$ is the price of energy sold from prosumer $i$ to prosumer $j$ at time step $k$; and $p_{E,ij}^k \ge 0$ represents the amount of energy that prosumer $i$ is willing to offer to prosumer $j$. Hence, the actual energy transfer from $i$ to $j$, is given by

$$\Delta p_{E,ij}^k = p_{E,ij}^k - p_{E,ji}^k, \forall j \in \mathcal{N}_i^0 \qquad (7)$$

which is a joint decision between prosumer $i$ and prosumer $j$. Transactions only occur between prosumers connected with each other in the financial network. Note that the sign constraint forcing $p_{E,ij}^k$ to be non-negative allows modeling the difference between buying and selling price associated with the transfer $\Delta p_{E,ij}^k$.

We introduce a 1-norm regularization term that accounts for the losses of the photovoltaic generator. It is defined as

$$l_i^k(p_{G,i}^k, q_{G,i}^k) := \eta_{G,i} (p_{G,i}^k + |q_{G,i}^k|), \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (8)$$

This captures the preference of prosumers in favor of not dispatching active or reactive power unless it benefits them, whether through incentives, market conditions, or self-consumption.
system losses in terms of the prosumers' active and reactive power injections

$$L(\mathbf{p}^k, \mathbf{q}^k) = \mathbf{p}^k \mathbf{R} \mathbf{p}^k + \mathbf{q}^k \mathbf{R} \mathbf{q}^k \quad (20)$$

We further compactify notation by defining $\mathbf{H} = (\mathbf{R}^T, \mathbf{X}^T)$ to obtain a linear relation between squared voltage magnitude and the net active and reactive power injections made by prosumers

$$v^k = \mathbf{H} \begin{pmatrix} \mathbf{p}^k \\ \mathbf{q}^k \end{pmatrix} + \mathbf{1} \quad (21)$$

Finally, recall that at the substation (bus 0), the distribution network is connected to an external grid (transmission system). We may express the net systemic power output export as the concave function

$$P(\mathbf{p}^k, \mathbf{q}^k) = \mathbf{1}^T \mathbf{p}^k - L(\mathbf{p}^k, \mathbf{q}^k), \quad (22)$$

interpreted as the difference between the total amount of energy produced within the distribution network and the losses.

## 2.4. The incentive design problem

In this section, we will leverage a reverse Stackelberg model to determine the incentive functions $T_i^k(\cdot)$, for each pair $(i, k) \in \mathcal{N} \times \mathcal{K}$. We refer to it as the incentive design problem, which is to be solved assuming global knowledge of the parameters that characterize the electrical grid and the prosumers. The objective is to drive the behavior of prosumers towards maximizing a predefined operational performance metric considering security constraints. We interpret the incentive terms $T_i^k(p_i^k, q_i^k)$ as a local perturbation on the utility function of each prosumer. We refer to them as “local” because the domain corresponds to the decisions of only prosumer $i$. Ultimately, we expect optimal global operation to emerge by locally perturbing the utility functions of the players.

The objective function of the regulatory entity (the leader) is an operational performance metric. For simplicity of exposition, we define it as the total exported active power

$$J_L(\mathbf{p}, \mathbf{q}) := \sum_{k \in \mathcal{K}} P(\mathbf{p}^k, \mathbf{q}^k) \quad (23)$$

which is a strictly concave quadratic function to be maximized.

It is worth noting that the leader's objective function solely relies on the variables $\mathbf{p}$ and $\mathbf{q}$. Therefore, the applied incentive is specifically targeted toward these variables. For computational tractability reasons, we restrict ourselves to a parametric family of functions on $(\mathbf{p}, \mathbf{q})$, such that

$$T_i^k(p, q) := \gamma_i^k(p, q)^T \begin{pmatrix} p \\ q \end{pmatrix} \quad (24)$$

where $\gamma_i^k: \mathbb{R}^2 \to \mathbb{R}^2$ is a semidefinite positive affine transformation. To be more precise, we state that $\gamma_i^k \in \mathcal{A}$ for all $(i, k) \in \mathcal{N} \times \mathcal{K}$ iff. there exist coefficients $\mathbf{\Phi}_i^k \in \mathbb{R}^{2 \times 2}$ and $\mathbf{\phi}_i^k$, such that $\mathbf{\Phi}_i^k = \mathbf{\Phi}_i^{k\top} \ge 0$ and

$$\gamma_i^k(p, q) = \frac{1}{2} \mathbf{\Phi}_i^k \begin{pmatrix} p \\ q \end{pmatrix} + \mathbf{\phi}_i^k \quad (25)$$

where

$$\mathbf{\Phi}_i^k := \begin{pmatrix} \mathbf{\Phi}_{pp,i}^k & \mathbf{\Phi}_{pq,i}^k \\ \mathbf{\Phi}_{pq,i}^k & \mathbf{\Phi}_{qq,i}^k \end{pmatrix}, \quad \mathbf{\phi}_i^k := \begin{pmatrix} \mathbf{\phi}_{p,i}^k \\ \mathbf{\phi}_{q,i}^k \end{pmatrix}$$

are the incentive function coefficients meant to be determined by the DSO. Hence, a tuple $(\mathbf{\Phi}, \mathbf{\phi})$ uniquely determines the incentive function $\gamma$.

We write the vector of tariff costs as a function of the vector of control actions for each time step $k$ by stacking the entries corresponding to each prosumer

$$T^k(\mathbf{u}^k) = (\gamma^k(\mathbf{u}^k)^\top \mathbf{u}^k \quad (26)$$

Accordingly, we define the matrix $\tilde{\mathbf{\Phi}}^k$ and vector $\tilde{\mathbf{\phi}}^k$ such that

$$\gamma^k(\mathbf{u}^k) = \tilde{\mathbf{\Phi}}^k \mathbf{u}^k + \tilde{\mathbf{\phi}}^k \quad (27)$$

Prosumers are assumed to behave rationally, solving problem (9) to minimize their disutility (maximize their profit). This is the main driver for the proposed approach. To that end, some design choices over the shape of functions $T_i$ should be put in place so that (9) is simple enough for prosumers to solve. First, the incentive cost term, $T_i^k(p_i^k, q_i^k)$, is known in advance to each prosumer, and it is a function of their own local active and reactive power output. Second, $T_i^k$ is a quadratic convex function.

Operation performance metrics of distribution networks ultimately capture the cost of operation. It might consider Ohmic losses on electrical branches<fnref n="1" />, total exported energy, voltage violations, line congestion, etc.

### 2.4.1. The reverse Stackelberg game model

Determining the incentive functions requires characterizing how they influence the behavior of prosumers and, ultimately, the system's operation. There is a need for a different kind of equilibrium concept than the one utilized to model the peer-to-peer market since the roles of the DSO and the prosumers are not symmetric.

In the peer-to-peer market equilibrium model described in Section 2.2, prosumers announced their strategies simultaneously. Here, the DSO has to announce the incentive functions first and then let the prosumers react accordingly. To capture this asymmetry, we rely on a particular case of Stackelberg equilibrium described in the following paragraphs.

As in Ho et al. (1981), we acknowledge a difference between announcing a strategy and performing an action. It is the order of announcing strategies rather than the order of actions that distinguishes the leader from the follower. The reverse Stackelberg equilibrium solution concept considers the following timeline:

1. the leader announces their strategy

2. knowing the leader's strategy, the followers choose their action

<footnote n="1">¹Lines and transformers.</footnote>
In the following, we present the optimization problem that models the decision-making process of each prosumer:

$$ \min_{\mathbf{u}_i} J_{F,i}(T(\cdot), \mathbf{u}_i, \mathbf{u}_{-i}) \qquad (9a) $$

$$ \text{s.t.} \quad p_i^k = \sum_{j \in \mathcal{N}_i^0} (p_{E,ij}^k - p_{E,ji}^k), \qquad \forall k \in \mathcal{K} \qquad (9b) $$

$$ p_i^k = p_{G,i}^k - p_{B,i}^k - p_{D,i}^k, \qquad \forall k \in \mathcal{K} \qquad (9c) $$

$$ q_i^k = q_{G,i}^k - q_{D,i}^k, \qquad \forall k \in \mathcal{K} \qquad (9d) $$

$$ s_{B,i}^k = s_{B,i}^{k-1} + \eta_{B,i} p_{B,i}^{+k} - \frac{p_{B,i}^{-k}}{\eta_{B,i}}, \qquad \forall k \in \mathcal{K} \qquad (9e) $$

$$ (p_{G,i}^k, q_{G,i}^k) \in \tilde{\mathcal{G}}_i^k, \qquad \forall k \in \mathcal{K} \qquad (9f) $$

$$ 0 \le s_{B,i}^k \le \hat{s}_{B,i}, \qquad \forall k \in \mathcal{K} \qquad (9g) $$

$$ 0 \le p_{B,i}^{+k} \le \hat{p}_{B,i}, \qquad \forall k \in \mathcal{K} \qquad (9h) $$

$$ 0 \le p_{B,i}^{-k} \le \hat{p}_{B,i}, \qquad \forall k \in \mathcal{K} \qquad (9i) $$

$$ 0 \le p_{E,ij}^k, \qquad \forall (j,k) \in \mathcal{N}_i^0 \times \mathcal{K} \qquad (9j) $$

where (9b) represents the balance between the net active power injection and the peer-to-peer market trade, constraints (9c) and (9d) are the active and reactive power balance with parameters $p_{D,i}^k$ and $q_{D,i}^k$ representing the active and reactive power demand, respectively. The dynamics of the battery are modeled via (9e). Capacity constraints of the distributed generator is considered in (9f) by means of a polyhedral inner approximation of (1), $\tilde{\mathcal{G}}_i^k$. The capacity of the battery is enforced by (9g), (9h), and (9i). Finally, the sign constraints imposed on the bilateral offers are enforced through (9j).

The computational complexity of (9) drastically depends on the selected $T(\cdot)$. If $T(\cdot)$ is a convex quadratic function, problem (9) forms a convex quadratic program with second-order cone constraints.

**Remark 1.** $p_{E,0i}^k$ is a decision variable of prosumer i representing the energy they are willing to buy from the DSO (slack prosumer indexed by 0) at time step k. Therefore the energy transfer $\Delta p_{E,i0}^k = p_{E,i0}^k - p_{E,0i}^k$ is a decision of prosumer i. In fact, we could replace $(p_{E,i0}^k - p_{E,0i}^k)$ by $\Delta p_{E,i0}^k$ whenever it appears in problem (9). As a consequence, we have that if $T$ is convex, then (9) is always feasible because prosumer can always satisfy their power consumption needs as $\Delta p_{E,i0}^k$ is assumed unbounded.

## 2.2. The peer-to-peer market model

The objective of this section is to leverage a model for the peer-to-peer market equilibrium. We consider two relevant classes of market designs: with and without a local market operator (LMO). The corresponding solution concepts representing the output of the market in each of the cases mentioned above are: (i) generalized Nash equilibria (GNEs) (Facchinei and Kanzow, 2007) and (ii) welfare optima. In the following paragraphs, we characterize the two solution concepts.

### 2.2.1. The set of GNEs

We assume that if there is no LMO in place and prosumers are free to trade with each other, then the set of market equilibria is characterized by the GNEs. In the context of this work, a GNE is a solution where no prosumer can improve their objective function by unilaterally changing their decision, given the decisions of all other prosumers. Prosumers behave selfishly and independently, minimizing their objective function under local and coupling constraints, based on (10). It extends the classical Nash equilibrium by also considering coupling constraints between the players.

The optimization problem (9) is coupled to the decision making of the other prosumers by the objective function (9a) and constraints (9b), since $p_{E,ji}$, with $j \in \mathcal{N}_i$, is a vector of parameters in prosumer i's decision problem, and simultaneously a vector of variables for prosumer j. We explicitly denote the coupling between the optimization problems by writing (9) in compact form as

$$ \min_{\mathbf{u}_{F,i} \in \Omega_{F,i}(\mathbf{u}_{F,-i})} J_{F,i}(T_i(\cdot), \mathbf{u}_{F,i}, \mathbf{u}_{F,-i}) \qquad (10) $$
3. the action of the leader is determined based on their previously announced strategy as a function of the actions of the followers

In the proposed model, the strategy of the leader is the tariff function $\gamma(\cdot)$, the action of the follower is the vector $\mathbf{u}$ whose elements are the active and reactive power injections, and the action of the leader is $\gamma(\mathbf{u})$. The reverse Stackelberg optimization problem is given as

$$ \text{\textbf{\textit{\textbf{max}}} } \gamma(\cdot) \quad J_{\text{L}}(\gamma(\cdot), \mathbf{u}) \qquad (28a) $$

$$ \text{s.t.} \quad \gamma(\cdot) \in \mathcal{A} \qquad (28b) $$

$$ \gamma(\cdot) \in \mathcal{D} \qquad (28c) $$

$$ \hat{\nu} \leq \mathbf{H}\mathbf{u}^k + \mathbf{1} \leq \hat{\nu} \quad \forall k \in \mathcal{K} \qquad (28d) $$

$$ \mathbf{u} \in S(\gamma(\cdot)) \qquad (28e) $$

where (28b) forces $\gamma(\cdot)$ to belong to the family of semidefinite positive affine functions from $\mathbb{R}^2$ to $\mathbb{R}^2$, denoted by $\mathcal{A}$, constraint (28c) imposes the budget condition limiting the action of the leader, constraint (28d) enforces the voltage levels, and (28e) is the lower-level equilibrium problem. The quoted max{$\cdot$} operator is utilized to emphasize that (28) might be ill-posed as the lower level problem might have a multiplicity of solutions (28e), see Dempe and Dutta (2012). Uniqueness of solutions in the space of active and reactive net power injection is studied in Section 3.

Depending on the focus, the lower-level (28e) might be replaced by

$$ \mathbf{u} \in \mathcal{W}(\gamma(\cdot)) \qquad (29) $$

to solve the problem as if the market clearing mechanism were efficient in the economic sense, i.e., if the equilibrium maximizes the welfare or, equivalently, minimizes the total disutility.

### 2.4.2. The incentive budget constraints

For each $(i, k) \in \mathcal{N} \times \mathcal{K}$, we assume that there is a physical bound on $\mathbf{w}_i^k := (p_i^k, q_i^k)$ given by the prosumer-to-grid connection capacity, which takes the form

$$ ||\mathbf{w}_i^k||_2 \leq E_i \qquad (30) $$

where $E_i$ is the connection capacity measured in the same physical units of apparent power, i.e. per unit or volt-ampere (VA) units. Therefore, we restrict $(\mathbf{\Phi}_i^k, \mathbf{\phi}_i^k)$ so that we ensure that for any $\mathbf{w}_i^k$ satisfying (30), the incentive cost terms are bounded. More precisely, we impose the following budget condition

$$ (\mathbf{\Phi}_i^k, \mathbf{\phi}_i^k) \in \mathcal{D}_i, \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (31) $$

such that

$$ D_i^k := \left\{ (\mathbf{\Phi}_i^k, \mathbf{\phi}_i^k) \in \mathbb{S}_2 \times \mathbb{R}^2 \mid \begin{array}{l} \forall \mathbf{w}_i^k \in \mathbb{R}^2 : \|\mathbf{w}_i^k\|_2 \le E_i, \\ -M_i^k \le (\mathbf{w}_i^k)^\top \mathbf{\Phi} \mathbf{w}_i^k + \mathbf{\phi}^\top \mathbf{w}_i^k \le M_i^k \end{array} \right\} \qquad (32) $$

where $\mathbb{S}_2$ is the set of symmetric positive semidefinite matrices of size 2.

**Proposition 1.** For each $(i, k) \in \mathcal{N} \times \mathcal{K}$, if $(\mathbf{\Phi}_i^k, \mathbf{\phi}_i^k)$ satisfies

$$ E_i ||\mathbf{\Phi}_i^k||_F + ||\mathbf{\phi}_i^k||_2 \le \frac{M_i}{E_i} \qquad (33) $$

then $(\mathbf{\Phi}_i^k, \mathbf{\phi}_i^k) \in \mathcal{D}_i$.

**PROOF.** For simplicity of notation, ignore the indices $i$ and $k$. We write (33) as

$$ E ||\mathbf{\Phi}||_F + ||\mathbf{\phi}||_2 \le \frac{M}{E} $$

Recall that the Frobenius norm is an upper bound on the matrix 2-norm. Then,

$$ E^2 ||\mathbf{\Phi}||_2 + E ||\mathbf{\phi}||_2 \le M $$

The connection capacity condition (30) states that $E$ is an upper bound on $||\mathbf{w}||_2$, therefore the following is true

$$ ||\mathbf{w}||_2^2 ||\mathbf{\Phi}||_2 + ||\mathbf{\phi}||_2 ||\mathbf{w}||_2 \le M $$

Then, we apply the the Cauchy-Schwarz inequality to obtain

$$ ||\mathbf{w}||_2 ||\mathbf{\Phi w}||_2 + ||\mathbf{\phi}||_2 ||\mathbf{w}||_2 \le M $$

followed by

$$ ||\mathbf{w}^\top \mathbf{\Phi} \mathbf{w} + \mathbf{\phi}^\top \mathbf{w}||_2 \le M $$

Consequently,

$$ -M \le \mathbf{w}^\top \mathbf{\Phi} \mathbf{w} + \mathbf{\phi}^\top \mathbf{w} \le M $$

□

Problem (28) can be more explicitly written such that $\gamma(\cdot)$ is parameterized in terms of the coefficients $(\mathbf{\Phi}, \mathbf{\phi})$.

Problem (28) can be more explicitly written such that $\gamma(\cdot)$ is parameterized in terms of the coefficients $(\mathbf{\Phi}, \mathbf{\phi})$.

The formulation of Problem (28) can be further refined by explicitly parameterizing $\gamma(\cdot)$ in terms of the coefficients $(\mathbf{\Phi}, \mathbf{\phi})$.

$$ \max_{\mathbf{\Phi}, \mathbf{\phi}} J_{\text{L}}(\mathbf{\Phi}, \mathbf{\phi}, \mathbf{u}) \qquad (34a) $$

$$ \text{s.t.} \quad \mathbf{\Phi}_{pp,i}^k + \mathbf{\Phi}_{qq,i}^k \ge 0 \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (34b) $$

$$ (\mathbf{\Phi}_{pq,i}^k)^2 \le \mathbf{\Phi}_{pp,i}^k \mathbf{\Phi}_{qq,i}^k \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (34c) $$

$$ E_i ||\mathbf{\Phi}_i^k||_F + ||\mathbf{\phi}_i^k||_2 \le \frac{M_i}{E_i} \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (34d) $$

$$ \hat{\nu} \le \mathbf{H}^k \mathbf{u}^k + \mathbf{h}^k \le \hat{\nu} \qquad \forall k \in \mathcal{K} \qquad (34e) $$

$$ \mathbf{u} \in S(\mathbf{\Phi}, \mathbf{\phi}) \qquad (34f) $$

where the objective function (34a) is concave quadratic, constraints (34b) and (34f) ensure that each matrix $\mathbf{\Phi}_i^k \in \mathbb{R}^{2 \times 2}$
<figure><img src="image_5.png" /><figcaption>Figure 3: Electric branch associated with bus $i \in \mathcal{B}$.</figcaption></figure>

where $\mathbf{u}_{F,i} := (\mathbf{u}_{F,i}^k)_{k \in \mathcal{K}}$ is the action profile of prosumer $i$ over the whole planning horizon, and $\Omega_{F,i}(\cdot)$ is a point-to-set mapping—also referred to as set-valued map (Rockafellar and Wets, 1998, Chapter 5),—representing the feasible set of prosumer $i$ as a function of the action profiles of the other prosumers.

More precisely, consider a game with a set $\mathcal{N}$ of players, and the tuple $(\Omega_{F,i}, J_{F,i}, T_i)_{i \in \mathcal{N}}$ as the data.

**Definition 1.** The set of GNEs, denoted by $S(T)$, is the set of elements $\bar{\mathbf{u}} \in \prod_{i \in \mathcal{N}} \mathbb{R}^{n_i}$ such that

$$ \bar{\mathbf{u}}_i \in \arg\min_{\mathbf{u}_i \in \Omega_i(\mathbf{u}_{-i})} J_{F,i}(T_i(\cdot), \mathbf{u}_i, \mathbf{u}_{-i}), \quad \forall i \in \mathcal{N} \qquad (11) $$

The concept of GNEs is also referred to as pseudo-Nash equilibria (Harker, 1991). In our case, the set $S(T)$ can be interpreted as modeling the behavior of prosumers uncoordinately making bilateral contracts of energy, subject to tariff $T$. We assume the market output can be modeled as the GNE that prosumers reach if they have perfect knowledge of each other's strategies. This equilibrium model reflects the market outcome based on competitive behavior with shared information, rather than the actual execution of individual decisions. The fundamental task of this work is to design $T$ such that all elements in $S(T)$ represent operationally secure and efficient action profiles.

### 2.2.2. The set of welfare optima

On the other hand, when there is a market clearing mechanism driving the system to an economically efficient equilibrium, i.e., one that maximizes the social welfare of all prosumers or, equivalently, minimizes the total disutility. The following definition captures the corresponding solution concept

**Definition 2.** The set of economically efficient market outputs, denoted by $\mathcal{W}(T)$, is the set of elements $\bar{\mathbf{u}} \in \prod_{i \in \mathcal{N}} \mathbb{R}^{n_i}$ such that

$$ \bar{\mathbf{u}} \in \arg\min_{\mathbf{u} \in \Omega} \sum_{i \in \mathcal{N}} J_{F,i}(T_i(\cdot), \mathbf{u}_i, \mathbf{u}_{-i}) \qquad (12) $$

We can interpret (12) as representing the decision-making of all prosumers jointly minimizing the total disutility.

**Remark 2.** It is natural to ask if there exist conditions under which individual prosumers, if acting selfishly, might want to deviate from $\mathcal{W}(T)$. This is the same as asking if there exists an instance of the problem for which $\bar{\mathbf{u}}_F \in \mathcal{W}(T) \setminus S(T)$. Similarly, we may pose questions regarding the reciprocal difference $S(T) \setminus \mathcal{W}(T)$, i.e., if there exists an output of a market without LMO that is not economically efficient. For further detail on the relation between the sets $S(T)$ and $\mathcal{W}(T)$ please refer to the Appendix A.

## 2.3. The effect of prosumers on the electrical grid

This section focuses on modeling the impact of prosumers' behavior on the operation of the distribution system. Specifically, we aim to establish the relationship between the net active and reactive power injected by prosumers and the resulting changes in voltage levels and system losses. In this section, we drop the time-step indices unless said otherwise because all relations presented here are static.

The optimal power flow (OPF) model is a non-convex quadratically constrained quadratic problem (QCQP) which is NP-Hard in general. Because of the difficulty of solving the OPF, several approximations have been made. Between them, convex relaxations are of the families: SDP, SOCP, and linear (Low, 2014a). SDP relaxations present the difficulty of extracting solutions if the rank constraint is not met. SOCP relaxations are not exact when voltage upper bound constraints are active (Low, 2014b). Therefore, we believe linear approximations are better suited for situations with high penetration of renewable energy resources since the system is exposed to overvoltage conditions. Among
is positive semidefinite, equation (34d) represents the budget conditions, (34e) are the voltage limits, and (34f) models the lower-level. Note that the complexities of problem (28) arise primarily from (34f), as it would otherwise be a simpler second-order conic problem.

The problem can be further compactly written as

$$ \gamma(\cdot) \in \arg \min_{\gamma(\cdot) \in \Gamma(\mathbf{u})} J_L(\gamma(\cdot), \mathbf{u}) \quad (35a) $$

$$ \mathbf{u}_i \in \arg \min_{\tilde{\mathbf{u}}_i \in \Omega_{F,i}(\mathbf{u}_{-i})} J_F(\gamma_i(\cdot), \tilde{\mathbf{u}}_i, \mathbf{u}_{-i}), \forall i \in \mathcal{N} \quad (35b) $$

### 3. Properties

In this section, our concern is to study how the lower-level equilibrium might respond to perturbations in the parameters defining the objective function of prosumers. Specifically, our concern is on studying stability in the sense of Dempe (2002) for the lower level as a parametric optimization problem. First, we show that the GNEP can be reformulated as an optimization problem, classically known as the Nikaido–Isoda (NI) reformulation, see Facchinei and Kanzow (2007). Then we show that the set of solutions to the GNEP, $S(T)$, is strongly stable towards perturbations in the parameters that define it by proving that the NI reformulation is strongly stable. As corollaries, we leverage guarantees of uniqueness of solutions in the net active and reactive power output space.

**Definition 3.** (Jointly convex GNEP) We say that (11) is a jointly convex GNEP if, for every prosumer $i \in \mathcal{N}$ and every $\mathbf{u}_{-i}$, the objective function $J_F(T(\cdot), \cdot, \mathbf{u}_{-i})$ is convex, the set $\Omega_i(\mathbf{u}_{-i})$ is closed and convex, and for some closed convex $\Omega \subset \mathbb{R}^n$ and all $i \in \mathcal{N}$ we have

$$ \Omega_i(\mathbf{u}_{-i}) = \{ \mathbf{u}_i \in \mathbb{R}^{n_i} \mid (\mathbf{u}_i, \mathbf{u}_{-i}) \in \Omega \} $$

As remarked in Facchinei and Kanzow (2007), if the feasible set of each player of a GNEP is described by inequalities and if all the shared constraints, i.e., the ones that involve variables from more than one player are the same for every player involved in the constraint, then the GNEP in question is a jointly convex GNEP. By this remark, it is clear that the GNEP in (11) is jointly convex. By the discussion above we derive the following result

**Lemma 1.** The game defined in (11) is a jointly convex GNEP.

First, we establish a classical result that bridges jointly convex GNEPs and optimization problems

**Lemma 2.** If $T_i$ is convex for each $i \in \mathcal{N}$, a vector $\bar{\mathbf{u}}$ is a solution to (11) if and only if $\bar{\mathbf{u}}$ is a global minimum of the optimization problem

$$ \min_{\mathbf{u}} \hat{V}(T, \mathbf{u}) \quad (36a) $$

$$ \text{s.t.} \quad \mathbf{u} \in \Omega, \quad (36b) $$

with zero objective function value, where

$$ \hat{V}(T, \mathbf{u}) := \sup_{\mathbf{y} \in \Omega} \psi(T, \mathbf{u}, \mathbf{y}) \quad (37) $$

and

$$ \psi(T, \mathbf{u}, \mathbf{y}) := \sum_{i \in \mathcal{N}} (J_{F,i}(T_i, \mathbf{u}_i, \mathbf{u}_{-i}) - J_{F,i}(T_i, \mathbf{y}_i, \mathbf{u}_{-i})) \quad (38) $$

**Lemma 2** corresponds to the classical Nikaido–Isoda (NI) reformulation in the case of jointly convex GNEPs, hence the proof is omitted, see Facchinei and Kanzow (2007). The NI reformulation is an important analytical tool that allows to transform the equilibrium in (11) into an optimization problem.

**Lemma 3.** If $T_i(\cdot)$ is convex for each $i \in \mathcal{N}$, problem (36) is a convex optimization problem with polyhedral feasible set

**PROOF.** Note that the objective function of each follower can be expressed as the sum

$$ J_F(T_i(\cdot), \mathbf{u}_i, \mathbf{u}_{-i}) = h_i(\mathbf{u}_i, T_i(\cdot)) + g_i(\mathbf{u}_{-i}), \forall i \in \mathcal{N} \quad (39) $$

with

$$ h_i(\mathbf{u}_i, T_i(\cdot)) := \sum_{k \in \mathcal{K}} \left( T_i^k(p_i^k, q_i^k) + l_i(p_{G,i}^k, q_{G,i}^k) \right) $$

and

$$ g_i(\mathbf{u}_{-i}) := \sum_{k \in \mathcal{K}} U_i^k (\mathbf{p}_{E,ij}, (p_{E,ij})_{j \in \mathcal{N}_i^0}) $$

Then the NI–function becomes

$$ \psi(T(\cdot), \mathbf{u}, \mathbf{y}) = \sum_{i \in \mathcal{N}} (h_i(\mathbf{u}_i, T_i(\cdot)) - h_i(\mathbf{y}_i, T_i(\cdot))) \quad (40) $$

For fixed $T(\cdot)$, function $\psi(T(\cdot), \mathbf{u}, \mathbf{y})$ is convex in $\mathbf{u}$, for all $\mathbf{y} \in \Omega$, therefore $\hat{V}(T(\cdot), \mathbf{u})$ is convex in $\mathbf{u}$ as the supremum operation in (37) preserves convexity. $\square$

**Theorem 1.** If $T(\cdot)$ is convex, the set of GNEs $S(T)$ of (11) is convex. Moreover, if $T(\cdot)$ is strictly convex, the set of GNEs of (11), $S(T)$, is uniquely defined in $(\mathbf{p}, \mathbf{q})$.

**PROOF.** By Lemma 3, problem (36) is a convex optimization problem. Therefore, the set of optimal solutions is a convex set. By Lemma 2, the solution set of (36) corresponds to $S(T)$ if the objective value of (36) is zero at optimality. Therefore, to prove Theorem 1, it suffices to prove that the optimal value of (36) is zero for any given $T(\cdot)$. For that, we express the NI-function as in (40), and considering $\Omega$ is closed, the optimal value of (36) can be written as

$$ \hat{V}^* = \min_{\mathbf{u} \in \Omega} \left\{ \max_{\mathbf{y} \in \Omega} \sum_{i \in \mathcal{N}} (h_i(\mathbf{u}_i, T_i(\cdot)) - h_i(\mathbf{y}_i, T_i(\cdot))) \right\} \quad (41) $$

which is zero. Therefore $S(T)$ is convex.

If $T(\cdot)$ is strictly convex, (40) is strictly convex with respect to $(\mathbf{p}, \mathbf{q})$, therefore the NI reformulation (36) is unique in $(\mathbf{p}, \mathbf{q})$, thus, by Lemma 2, so does the set $S(T)$. $\square$

The DSO exclusively considers the net active and reactive power output, denoted as $(\mathbf{p}, \mathbf{q})$, as the operational performance metrics hinge on these variables. Hence, our focus is
the linear OPF models, we highlight the linear OPF model presented by Coffrin and Van Hentenryck (2014), the two models for unbalanced networks presented by Bernstein et al. (2018), and the widely referenced *LinDistFlow* model proposed by Baran and Wu (1989).

In this work, we adopt the *LinDistFlow* model, detailed below. We choose this model primarily for its simplicity of exposition, as it is easy to present and construct. Despite being a linear approximation of the Optimal Power Flow (OPF), it provides sufficient accuracy for our purposes, although it is only valid for radial networks. For non-radial networks, the proposed demand response model can be adapted to incorporate the linear models presented by Coffrin and Van Hentenryck (2014) and Bernstein et al. (2018).

For each bus $i \in B$, the following equations model the relationship between active power, reactive power, and voltage

$$p_i = \sum_{j \in C_i} \tilde{P}_j - \tilde{P}_i \qquad (13a)$$

$$q_i = \sum_{j \in C_i} \tilde{Q}_j - \tilde{Q}_i \qquad (13b)$$

$$v_{A_i} = v_i + 2\tilde{r}_i \tilde{P}_i + 2\tilde{x}_i \tilde{Q}_i \qquad (13c)$$

Here, $\tilde{P}_i$ and $\tilde{Q}_i$ represent the active and reactive power flows, respectively, on the upstream branch of bus $i \in B$, measured in the downstream direction. The parameters $\tilde{r}_i$ and $\tilde{x}_i$ denote the resistance and reactance of that branch, while $v_i$ corresponds to the square of the voltage magnitude at bus $i \in B$. The set $C_i$ represents the set of successor buses connected to $i$. A summary of the notation is provided in Figure 3. The *LinDistFlow* is limited to radial (tree graph) networks with balanced voltage and currents. If the reader requires a linear model capturing a more general setting, such as unbalanced and/or non-radial networks, they may refer to the first-order Taylor linearization proposed by Bernstein et al. (2018).

To calculate the system losses, define matrices $\tilde{\mathbf{R}} \in \mathbb{R}^{|B| \times |B|}$ and $\tilde{\mathbf{X}} \in \mathbb{R}^{|B| \times |B|}$ with coefficients

$$\tilde{\mathbf{R}}_{ij} := \sum_{k \in \mathcal{L}(i) \cap \mathcal{L}(j)} \tilde{r}_k, \quad \tilde{\mathbf{X}}_{ji} := \sum_{k \in \mathcal{L}(i) \cap \mathcal{L}(j)} \tilde{x}_k \qquad (14)$$

where $\mathcal{L}(i)$ is the set of lines between bus $i$ and the slack bus 0, including $i$. Both $\tilde{\mathbf{R}}$ and $\tilde{\mathbf{X}}$ are positive definite matrices, and the losses of the system can be expressed as a convex quadratic function

$$L(\tilde{\mathbf{p}}, \tilde{\mathbf{q}}) = \tilde{\mathbf{p}}^{\top} \tilde{\mathbf{R}} \tilde{\mathbf{p}}^{k} + \tilde{\mathbf{q}}^{\top} \tilde{\mathbf{R}} \tilde{\mathbf{q}} \qquad (15)$$

The operational security constraints are voltage magnitude bounds at each bus. Recall that in (13), our proxy for the voltage magnitudes are the voltage magnitudes squared. Therefore, we express the security constraints limits as the following linear constraints

$$\hat{v}_i^2 \le v_i \le \hat{v}_i^2, \quad \forall i \in B \qquad (16)$$

Now, it is necessary to translate the previous expressions into prosumer related-quantities. The correspondence between prosumers and buses can be described by a relation $\mathbf{B} \times \mathcal{N} \to \{0, 1\}$, represented by the incidence matrix $\mathbf{B} \in \{0, 1\}^{|\mathcal{B}| \times |\mathcal{N}|}$, such that

$$B_{ij} := \begin{cases} 1, & \text{if prosumer } j \text{ is connected to bus } i \\ 0, & \text{otherwise} \end{cases} \qquad (17)$$

Let us consider again the time step indices and obtain the voltage magnitude and losses in terms of the prosumers' injections. Recall that $\mathbf{p}^k := (p_i^k)_{i \in \mathcal{N}}$ and $\mathbf{q}^k := (q_i^k)_{i \in \mathcal{N}}$, for all $k \in \mathcal{K}$, and aggregate the effect of prosumers by bus by applying the following transformation based on (17)

$$\tilde{\mathbf{p}}^k = \mathbf{B} \mathbf{p}^k \qquad (18a)$$

$$\tilde{\mathbf{q}}^k = \mathbf{B} \mathbf{q}^k \qquad (18b)$$
on ensuring the robustness of the solution, specifically within the (p, q)-space, against variations in lower-level model parameters. Thus, providing a robust stability guarantee in this projected space suffices.

Expanding on this concept, it worth noting that the objective function (34a) and the coupling constraints (34e) in the tariff design problem (34) solely depend on (p, q). Consequently, establishing a uniqueness guarantee for $S(T(\cdot))$ in (p, q) circumvents challenges associated with non-uniqueness of solutions in mathematical programs with equilibrium constraints. This guarantee is easily attained using the tools we have leveraged thus far.

To ensure robustness of the solution against parameter variations in the lower-level problem, we extend our approach beyond merely ensuring uniqueness to focus the notion of stability. Our approach involves drawing upon the concept of strong stability for local solutions, as defined in Dempe (2002), which we paraphrase as follows

**Definition 4 (Dempe (2002)).** A local solution $(\bar{x}, \bar{y})$ of a parametric optimization problem

$$ \min_{\mathbf{x}} f(\mathbf{x}, \mathbf{y}) \qquad (42a) $$

$$ \text{s.t.} \quad g(\mathbf{x}, \mathbf{y}) \le 0 \qquad (42b) $$

is said to be strongly stable if there exist neighborhoods $V_\delta(\bar{y})$ with $\delta > 0$ of $\bar{y}$ and $V_\epsilon(\bar{x})$ with $\epsilon > 0$ of $\bar{x}$, and a uniquely determined continuous vector-valued function $x: V_\delta(\bar{y}) \to V_\epsilon(\bar{x})$ such that $x(y)$ is the unique local optimal solution of problem (42) in $V_\epsilon(\bar{x})$ for all $y \in V_\delta(\bar{y})$.

First, we introduce a projection operator. Let $\bar{u} \in \Omega$ represent a solution at the lower level, with $\bar{u}_i^k := (\bar{p}_i^k, \bar{q}_i^k, \bar{p}_{G,i}^k, \bar{q}_{G,i}^k, \bar{p}_{B,i}^{+k}, \bar{p}_{B,i}^{-k}, \bar{s}_{B,i}^k, \bar{p}_{E,i}^k)$. We define $\Pi_{pq}(\bar{u}_i^k) := (\bar{p}_i^k, \bar{q}_i^k)$ and extend the notation to aggregate variables, such that $\Pi_{pq}(\bar{u}) := (\bar{p}, \bar{q})$.

**Proposition 2.** If $(\bar{\Phi}, \bar{\phi}, \bar{u})$ is an optimal solution to (34) and $\bar{w} = \Pi_{pq}(\bar{u})$, then the projected solution $(\bar{\Phi}, \bar{\phi}, \bar{w})$ is strongly stable.

**PROOF.** From (41), $V(T(\cdot), \mathbf{u})$ can be written as

$$ V(T(\cdot), \mathbf{u}) = \sum_{i \in \mathcal{N}'} h_i(T(\cdot), \mathbf{u}) + C \qquad (43) $$

where C is constant with respect to $\mathbf{u}$, given by

$$ C := \max_{y \in \Omega} \left\{ - \sum_{i \in \mathcal{N}'} h_i(T(\cdot), y) \right\} $$

thus, $V(T(\cdot), \mathbf{u})$ is differentiable in $\mathbf{u}$. In fact, the partial_hesian matrix of the lagrangian function of the NI reformulation (36) with respect to $\mathbf{w}$ is a block-diagonal matrix

$$ \frac{\partial^2 L}{\partial \mathbf{w}^2} = \mathrm{diag} \left( \left( \bar{\Phi}_i^k \right)_{i \in \mathcal{N}, k \in \mathcal{K}} \right) $$

Therefore $\frac{\partial^2 L}{\partial \mathbf{w}^2}$ is positive definite because each block diagonal element is positive definite as $\epsilon > 0$; thus $(\bar{\Phi}, \bar{\phi}, \bar{w})$ satisfies the strong sufficient optimality condition of second order Dempe (2002). In addition, (36) satisfies the Mangasarian-Fromowitz constraint qualification (MFCQ). Then, Proposition 2 follows directly from Theorem 4.4 in Dempe (2002). $\square$

Proposition 2 stands as a robustness guarantee, ensuring that the tariff solution remains reliable in the face of parameter variations of the lower-level. Note that stability defined in this sense implies uniqueness of the market output with respect to $(\mathbf{p}, \mathbf{q})$.

Furthermore, the stability property not only applies for the GNEP-based market setup. It remains true for both when the market is modeled as a GNEP and when it is modeled as an optimization problem that aims to minimize total disutility, as presented in Definition 2. To see this, note that if $T(\cdot)$ is strictly convex, then (12) is a strictly convex problem with respect to $(\mathbf{p}, \mathbf{q})$, and the same line of analysis developed in this section can be applied. This implies that the stability condition we have established is equally applicable to both managed and unmanaged market designs, solidifying the versatility of the proposed model.

We highlight that the NI reformulation is a valuable tool for characterizing the stability of jointly convex GNEPs. Note that by characterizing the stability of the solution of the NI reformulation, we bypass the notions of stability in variational systems such as Rockafellar (2023) which are commonly used in more general variational systems when there is no optimization problem reformulation at hand.

# 4. Numerical results

An instance of the problem comprises a distribution network, the prosumers with their devices, and the financial grid. For each instance, the computational experiments follow a scheme in which the tariff is calculated by the DSO and passed to the prosumers in a rolling horizon fashion<fnref n="2" />. In this case, tariff functions for all prosumers are calculated considering a planning horizon of 4 hours and applied only for 1 hour. We run the experiments utilizing a 10-minute time step setting; therefore, each instance of (34) to be solved has $|\mathcal{K}| = 6 \times 4 = 24$ time steps.

## 4.1. A toy problem

This section is dedicated to leveraging a toy problem, with the objective of gaining deeper insights into the proposed tariff design and offering a tangible proof of concept. First, we show how the DSO is able to induce safe operation by prompting the voltage profiles to remain within nominal operation bounds. Then, we establish the main differences between the proposed approach and the conventional constant tariff, which involves only bilinear terms between leader and follower variables.

The experimental setup for this toy problem involves a set $\mathcal{N} = \{1, 2, 3\}$ of prosumers connected to a 4 bus electrical

<footnote n="2">²Also called receding horizon.</footnote>
Such that

$$ \mathbf{R} := \mathbf{B}^{\top} \tilde{\mathbf{R}} \mathbf{B} \qquad (19a) $$

$$ \mathbf{X} := \mathbf{B}^{\top} \tilde{\mathbf{X}} \mathbf{B} \qquad (19b) $$

The linear transformation **B** preserves positive definiteness; therefore, again, we obtain a convex quadratic function for the system losses in terms of the prosumers' active and reactive power injections

$$ L(\mathbf{p}^k, \mathbf{q}^k) = \mathbf{p}^k \mathbf{R} \mathbf{p}^k + \mathbf{q}^k \mathbf{R} \mathbf{q}^k \qquad (20) $$

We further compactify notation by defining $\mathbf{H} = (\mathbf{R}^{\top}, \mathbf{X}^{\top})$ to obtain a linear relation between squared voltage magnitude and the net active and reactive power injections made by prosumers

$$ \nu^k = \mathbf{H} \begin{pmatrix} \mathbf{p}^k \\ \mathbf{q}^k \end{pmatrix} + \mathbf{1} \qquad (21) $$

Finally, recall that at the substation (bus 0), the distribution network is connected to an external grid (transmission system). We may express the net systemic power output export as the concave function

$$ P(\mathbf{p}^k, \mathbf{q}^k) = \mathbf{1}^{\top} \mathbf{p}^k - L(\mathbf{p}^k, \mathbf{q}^k), \qquad (22) $$

interpreted as the difference between the total amount of energy produced within the distribution network and the losses.

## 2.4. The incentive design problem

In this section, we will leverage a reverse Stackelberg model to determine the incentive functions $T_i^k(\cdot)$, for each pair $(i, k) \in \mathcal{N} \times \mathcal{K}$. We refer to it as the incentive design problem, which is to be solved by the DSO assuming global knowledge of the parameters that characterize the electrical grid and the prosumers. Likewise, the DSO announces its incentive policy such that each prosumer $i$ solves (9) with full knowledge of $T_i^k(\cdot)$ for all $k \in \mathcal{K}$. This is similar to how taxation rules are announced in advance to ensure transparency and predictability in real-life scenarios and follows the spirit of the taxation design framework proposed by Labbé et al. (1998).

The objective of the incentive design problem is to drive the distribution system towards maximizing a predefined operational performance metric by influencing prosumer behavior through their responses to the incentive functions. We interpret the incentive terms $T_i^k(p_i^k, q_i^k)$ as a local perturbation on the utility function of each prosumer, termed “local” because they only affect the decisions of the specific prosumer $i$. Ultimately, we aim for the optimal operation of the entire distribution system to emerge from locally perturbing the utility functions of individual players.

Operation performance metrics of distribution networks ultimately capture the cost of operation. It might consider Ohmic losses on electrical branches<fnref n="1" />, total exported energy, voltage violations, line congestion, etc. For simplicity of exposition, we define it as the total exported active power

$$ J_{\mathrm{L}}(\mathbf{p}, \mathbf{q}) := \sum_{k \in \mathcal{K}} P(\mathbf{p}^k, \mathbf{q}^k) \qquad (23) $$

which is a strictly concave quadratic function to be maximized.

It is worth noting that the leader's objective function solely relies on the variables **p** and **q**. Therefore, the applied incentive is specifically targeted toward these variables. We restrict ourselves to a parametric family of functions on $(\mathbf{p}, \mathbf{q})$, such that

$$ T_i^k(p_i^k, q_i^k) = \gamma_i^k(p_i^k, q_i^k)^{\top} \begin{pmatrix} p_i^k \\ q_i^k \end{pmatrix} \qquad (24) $$

where $\gamma_i^k: \mathbb{R}^2 \to \mathbb{R}^2$ is a semidefinite positive affine transformation. To be more precise, we define the space of positive definite affine functions $\mathcal{A}$ such that $\gamma_i^k$ belongs to $\mathcal{A}$ for all $(i, k) \in \mathcal{N} \times \mathcal{K}$ if and only if there exist coefficients $\Phi_i^k \in \mathbb{R}^{2 \times 2}$ and $\phi_i^k$, such that $\Phi_i^k = \Phi_i^{k\top} \ge 0$ and

$$ \gamma_i^k(p_i^k, q_i^k) = \frac{1}{2} \Phi_i^k \begin{pmatrix} p_i^k \\ q_i^k \end{pmatrix} + \phi_i^k \qquad (25) $$

<footnote n="1">¹Lines and transformers.</footnote>
grid where the set of buses is $B^0 = \{0, 1, 2, 3\}$, with parent buses $A_1 = 0, A_2 = 1$, and $A_3 = 0$. The set of prosumers is $\mathcal{N} = \{1, 2, 3\}$ such that $\mathbf{B} = \text{diag}(1, 1, 1)$ as depicted in Figure 1. The distribution system corresponds to the IEEE 4-bus test feeder IEEE (1991) which has been adapted to this setting by balancing the line segments using Carson's equation and Kron reduction. Distributed loads have been converted into balanced three-phase loads. The resistance of each electrical branch has been increased so that the resistance to reactance ratio is exactly 2.3—a normal ratio in distribution systems.—In this case, each prosumer model represents a group of actual prosumers. This is because the IEEE 4-bus system is a simplified electrical equivalent of a whole distribution system. Therefore, to analyze the effect of the prosumers in the operation of the grid, it is necessary to aggregate the effect of a large group of them.

The time series data for the problem is depicted in Figure 4. On top is the maximum available solar power, and at the bottom is the active power load. The signals have been obtained from one-minute time series data from a set of individual homes from the U.S. Austin region Pecan Street (2021). The time-series have been upsampled, averaging the values over periods of 10 minutes to accommodate them to the current time setting. A constant power factor of 0.85 is considered; therefore, the reactive power load is proportional to the active power load by a constant factor of $\tan(\cos^{-1}(0.85)) \approx 0.61$.

Naturally, the high solar irradiance period occurs around midday. During this period, there are high chances of violating the upper voltage limit because the net active power injections increase and $\mathbf{R}$ in (19a) is a positive definite matrix that maps the increment in power injection into an increment of voltage magnitudes. Without the tariff incentive, the voltage magnitude of electrical Node 2 would have been out of bounds. Figure 5 depicts the voltage profiles considering the tariff. Note that the voltage profile gets flattened to remain within bounds during the high irradiance interval. This is achieved by inducing the active and reactive net power profiles shown in figures 6 and 8. A similar process occurs during the high-demand periods in the evening, where voltages tend to lower, and they get flattened as a result of the tariff to avoid a lower limit voltage violation.

In Figure 6, we observe that in time periods near midday, the induced net power injection depicted in orange has less noise than the blue line representing the maximum available net active power. Recall that the tariff gets updated and passed to the prosumers hourly. The red vertical lines show the moments in which the tariff is recomputed and passed to the prosumers. Notably, within these intervals, the net active power injection maintains a near-constant level despite the underlying potential for variation. This observation suggests a denoising effect associated with the tariff updates, which likely contributes to enhancing the predictability of the system's operation—a beneficial characteristic for operating the distribution network.

## 4.2. Experiments on larger instances

In evaluating our proposed method, we look at various metrics and extend our tests to larger scenarios based on the IEEE 4-bus and 34-bus systems. We adopt the same rolling horizon approach as we did in the earlier section solving 168

<figure><img src="image_5.png" /><figcaption>Figure 4: Time series data of maximum available power (top), and active power load (bottom).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Voltage magnitude per bus.</figcaption></figure>

instances of the tariff design problem, one each hour over a week. We considered the following performance metrics:

* **Hourly Average Number of Voltage Violations per Bus (NVV):** This is the count of average voltage magnitude deviations from the acceptable range per bus, calculated over 10-minute intervals.

* **Hourly Average Net Energy Export per Prosumer (NEE) in kW:** We measure the total energy exported to the grid, divided by the evaluation period and the number of prosumers.

* **Hourly Average Distribution Grid Losses (DGL) in kW:** This represents the total resistive losses in the grid, averaged over the evaluation period.

* **Hourly Average Disutility per Prosumer (DPP) in $:** This is the sum of each prosumer's disutility, averaged over the time of the evaluation.

* **Hourly Average Tariff Cost per Prosumer (TC) in $:** We accumulate the tariff costs from the prosumers'
where

$$ \mathbf{\Phi}_i^k := \begin{pmatrix} \mathbf{\Phi}_{pp,i}^k & \mathbf{\Phi}_{pq,i}^k \\ \mathbf{\Phi}_{pq,i}^k & \mathbf{\Phi}_{qq,i}^k \end{pmatrix}, \quad \mathbf{\Phi}_i^k := \begin{pmatrix} \phi_{p,i}^k \\ \phi_{q,i}^k \end{pmatrix} $$

are the incentive function coefficients meant to be determined by the DSO. Hence, a tuple $(\mathbf{\Phi}, \mathbf{\phi})$ uniquely determines the incentive function $\gamma$.

In compact notation, we write the vector of tariff costs as a function of the vector of control actions for each time step $k$ by stacking the entries corresponding to each prosumer

$$ T^k(\mathbf{u}^k) = (\gamma^k(\mathbf{u}^k))^\top \mathbf{u}^k \qquad (26) $$

Accordingly, we define the matrix $\tilde{\mathbf{\Phi}}^k$ and vector $\tilde{\mathbf{\phi}}^k$ such that

$$ \gamma^k(\mathbf{u}^k) = \tilde{\mathbf{\Phi}}^k \mathbf{u}^k + \tilde{\mathbf{\phi}}^k \qquad (27) $$

We define $\gamma(\cdot)$ as an affine function with a positive semidefinite coefficient matrix for two main reasons. First, it ensures that problem (9) remains a convex quadratic problem, supporting the assumption that prosumers behave rationallymeaning they solve problem (9) to global optimality, thereby maximizing their profit. Second, the quadratic form introduced in the prosumers' cost function is crucial for establishing the stability of problem (9) in the $(\mathbf{p}, \mathbf{q})$ space, as a parametric optimization problem, which is addressed in Section 3. In addition, while other families of functions might also be suitable, we have chosen affine functions due to their simplicity and well-studied properties in the context of reverse Stackelberg games (Groot et al., 2016).

Note that the definition in (24) allows prosumers to experience null tariff costs if they choose to rely solely on self-consumption. This means that a prosumer $i$ might decide not to interact with the grid by setting their net active and reactive power to zero, i.e., $p_i^k = 0$ and $q_i^k = 0$, for some period $k$. As a result, prosumers would neither pay any tariff nor receive any monetary incentive for that period.

### 2.4.1. The reverse Stackelberg game model

Determining the incentive functions requires characterizing how they influence the behavior of prosumers and, ultimately, the system's operation. There is a need for a different kind of equilibrium concept than the one utilized to model the peer-to-peer market since the roles of the DSO and the prosumers are not symmetric. In the peer-to-peer market equilibrium model described in Section 2.2, prosumers announced their strategies simultaneously. Here, the DSO has to announce the incentive functions first and then let the prosumers react accordingly. To capture this asymmetry, we rely on a particular case of Stackelberg equilibrium described in the following paragraphs.

As Ho et al. (1981), we acknowledge a difference between announcing a strategy and performing an action. It is the order of announcing strategies rather than the order of actions that distinguishes the leader from the follower. The reverse Stackelberg equilibrium solution concept considers the following timeline:

1. the leader announces their strategy $\gamma(\cdot)$

2. knowing the leader's strategy, the followers choose their action $\mathbf{u}$

3. the action of the leader $\gamma(\mathbf{u})$ is determined based on their previously announced strategy as a function of the actions of the followers

Step 3 is predetermined given the actions of the followers and does not involve further choice, rather it can be interpreted as the evaluation of $\gamma(\mathbf{u})$ which is already considered in the decision problem of prosumers.

In the proposed model, the strategy of the leader is the tariff function $\gamma(\cdot)$, the action of the follower is the vector $\mathbf{u}$ whose elements are the active and reactive power injections, and the action of the leader is $\gamma(\mathbf{u})$. The reverse Stackelberg optimization problem is given as

$$ \max_{\gamma(\cdot)} \quad J_L(\gamma(\mathbf{u}), \mathbf{u}) \qquad (28a) $$

$$ \text{s.t.} \quad \gamma(\cdot) \in \mathcal{A} \qquad (28b) $$

$$ \gamma(\cdot) \in \mathcal{D} \qquad (28c) $$
<figure><img src="image_7.png" /><figcaption>**Figure 6:** Net active power injection at maximum power point (MPP) vs. the induced signal.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7:** Reactive power generation of prosumer 2.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8:** Reactive power generation of prosumer 2.</figcaption></figure>

objective function over each period and then average it out.

Upon examining the NVV and NEE columns in Table 1, a comparison between the proposed approach and the constant tariff policy reveals a notable observation. While both strategies ensure the system operates safely (NVV = 0), the proposed tariff design stands out by achieving this objective alongside maintaining higher levels of net energy export (NEE). However, it is important to note that the proposed approach does incur higher levels of ohmic losses. This can be attributed to the

<figure><img src="image_10.png" /><figcaption>**Figure 9:** Cumulative probability density function objective function of the DSO, IEEE-4bus test feeder.</figcaption></figure>

increased amounts of energy being exported, which naturally leads to greater loss in the transmission process.

With the focus on operation performance, we did not consider the tariff budget in the objective function of the leader. Note that in Table 1, the TC metric changes signs from the proposed tariff to the constant tariff. The results are likely to be different if considering the budget upperbound $M_i$ for $i \in \mathcal{N}$ in the objective function of the leader, which can be done without changing the structure of the problem.

While these metrics offer insights into different aspects of performance, our primary focus is on minimizing the DSO's objective function, as it is indicative of the overall effectiveness of our model. The graph in Figure 9 illustrates the cumulative probability distribution of the DSO's objective function values across the three strategies we have put to the test.

## 5. Conclusions

In our research, we have successfully demonstrated the efficacy of a reverse Stackelberg game approach to demand response in power distribution networks. The proposed framework positions the DSO as the leader with the ability to impose a tariff function, to which the prosumers, interpreted as followers, react. This framework also encompasses a peer-to-peer market among the prosumers, modeled as a GNEP, and takes into account both active and reactive power generation and their effects on the network's voltage levels. This comprehensive approach addresses a commonly neglected aspect in current research, namely, the impact of market outcomes on the operation of the distribution network.

The numerical experiments indicate a marked improvement in maximizing the net power output and maintaining the voltage levels within safe margins, compared to the conventional constant time-of-use tariff, particularly when perturbations in the followers' objective function parameters are introduced. To understand the reason behind the effectiveness in the presence of perturbations in the parameters of the lower-level, we provide theoretical results ensuring that the obtained solutions are always strongly stable. Moreover, our method stands out in its ability to denoise net active and reactive power injection of prosumers, which is also attributed to the stability property mentioned above. These findings underscore the potential of our bilevel programming solution to not only streamline the
$$ \tilde{v} \leq \mathbf{H}\mathbf{u}^k + 1 \leq \hat{v} \qquad \forall k \in \mathcal{K} \qquad (28d) $$

$$ \mathbf{u} \in S(\gamma(\cdot)) \qquad (28e) $$

where (28b) forces $\gamma(\cdot)$ to belong to the family of semidefinite positive affine functions from $\mathbb{R}^2$ to $\mathbb{R}^2$, denoted by $\mathcal{A}$, constraint (28c) imposes the budget condition limiting the action of the leader, constraint (28d) enforces the voltage levels, and (28e) is the lower-level equilibrium problem. The quoted max{$\cdot$} operator is utilized to emphasize that (28) might be ill-posed as the lower level problem might have a multiplicity of solutions (28e) (Dempe and Dutta, 2012). Uniqueness of solutions in the space of active and reactive net power injection is studied in Section 3.

Following the style of Groot et al. (2016), (28) could be more compactly written as

$$ \bar{\gamma}(\cdot) \in \arg \max_{\gamma(\cdot) \in \Gamma(\bar{\mathbf{u}})} J_L(\gamma(\bar{\mathbf{u}}), \bar{\mathbf{u}}) \qquad (29a) $$

$$ \bar{\mathbf{u}}_i \in \underset{\mathbf{u}_i \in \Omega_{F,i}(\bar{\mathbf{u}}_{-i})}{\arg \min} J_F(\bar{\gamma}_i(\mathbf{u}_i), \mathbf{u}_i, \bar{\mathbf{u}}_{-i}), \forall i \in \mathcal{N} \qquad (29b) $$

In the literature of reverse Stackelberg games, obtaining a solution is often divided into two stages (Groot et al., 2016). In the first stage, a target equilibrium, is defined, typically corresponding to the team problem solution, i.e., the equilibrium of the game if all players (both the leader and the followers) act as if they are on the same team, aiming to maximize a shared objective function, usually the leader's. The team problem solution is found by solving a centralized optimization problem where the decisions of both the leader and the followers are optimized together. In the second stage, the leader designs an incentive strategy (often an affine or other type of function) to induce the followers to act in a way that achieves or approximates this team problem solution in a decentralized manner.

However, a limitation of this approach is that the team problem solution may not always be inducible, meaning that there may not exist an incentive function strategy $\gamma(\cdot)$ that can drive the followers to achieve the team problem solution $\bar{\mathbf{u}}$. In fact, questions about existence of solutions developed in works by Worku et al. (2023) and Groot et al. (2016) are formulated around this issue.

To address this limitation, we propose a solution approach in which both target equilibrium and the incentive functions are found simultaneously. In our model, the target equilibrium is defined as the best equilibrium that can be induced by an affine function. If the team problem solution is inducible by an affine strategy, our target equilibrium naturally coincides with that solution. However, if the team problem is not inducible, our target equilibrium represents the best achievable outcome. Therefore, the solution approach proposed in this work enhances the applicability of reverse Stackelberg games to more practical settings with constrained solution spaces. In this case, questions about existence of solutions.

### 2.4.2. The incentive budget constraints

For each $(i, k) \in \mathcal{N} \times \mathcal{K}$, we assume that there is a physical bound on $\mathbf{w}_i^k := (p_i^k, q_i^k)$ given by the prosumer-to-grid connection capacity, which takes the form

$$ \|\mathbf{w}_i^k\|_2 \leq E_i \qquad (30) $$

where $E_i$ is the connection capacity measured in the same physical units of apparent power, i.e. per unit or volt-ampere (VA) units. Therefore, we restrict $(\mathbf{\Phi}_i^k, \mathbf{\phi}_i^k)$ so that we ensure that for any $\mathbf{w}_i^k$ satisfying (30), the incentive cost terms are bounded. More precisely, we impose the following budget condition

$$ (\mathbf{\Phi}_i^k, \mathbf{\phi}_i^k) \in \mathcal{D}_i^k, \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (31) $$

such that

$$ \mathcal{D}_i^k := \left\{ (\mathbf{\Phi}_i^k, \mathbf{\phi}_i^k) \in \mathbb{S}_2 \times \mathbb{R}^2 \mid \begin{array}{l} \forall \mathbf{w}_i^k \in \mathbb{R}^2 : \|\mathbf{w}_i^k\|_2 \leq E_i, \\ -M_i^k \leq (\mathbf{w}_i^k)^\top \mathbf{\Phi} \mathbf{w}_i^k + \mathbf{\phi}^\top \mathbf{w}_i^k \leq M_i^k \end{array} \right\} \qquad (32) $$

where $M_i^k$ are positive constants and $\mathbb{S}_2$ is the set of symmetric positive semidefinite matrices of size 2.

**Proposition 1.** For each $(i, k) \in \mathcal{N} \times \mathcal{K}$, if $(\mathbf{\Phi}_i^k, \mathbf{\phi}_i^k)$ satisfies

$$ E_i \|\mathbf{\Phi}_i^k\|_F + \|\mathbf{\phi}_i^k\|_2 \leq \frac{M_i^k}{E_i}, \qquad (33) $$
<table><caption>**Table 1**
Mean value of selected performance metrics.</caption><thead><tr><th rowspan="2">Strategies</th><th colspan="4">IEEE 4-bus test feeder</th><th colspan="4">IEEE 34-bus test feeder</th></tr><tr><th>NVV</th><th>NEE</th><th>DGL</th><th>DPP</th><th>TC</th><th>NVV</th><th>NEE</th><th>DGL</th><th>DPP</th><th>TC</th></tr></thead><tbody><tr><td>No tariff</td><td>0.087</td><td>0.351</td><td>0.028</td><td>-0.300</td><td>0.</td><td>0.101</td><td>0.118</td><td>0.089</td><td>-0.106</td><td>0.</td></tr><tr><td>Constant tariff</td><td>0.</td><td>0.218</td><td>0.022</td><td>-0.119</td><td>0.088</td><td>0.</td><td>0.073</td><td>0.065</td><td>-0.030</td><td>0.078</td></tr><tr><td>Proposed</td><td>0.</td><td>0.339</td><td>0.027</td><td>-0.292</td><td>-0.022</td><td>0.</td><td>0.108</td><td>0.077</td><td>-0.126</td><td>-0.056</td></tr></tbody></table>

objective function of the DSO but also to bring consistency and predictability to prosumer behavior, thereby reinforcing the overall efficiency and reliability of local energy markets.

We highlight the application of the Nikaido–Isoda reformulation to the lower–level GNEP to characterize the set of GNES. Looking forward, we believe this reformulation holds significant potential for practical implementation, particularly in solving Stackelberg games with a GNEP at the lower level through a branch–and–cut strategy.

One limitation of the current approach is the assumption of fixed and known peer–to–peer prices, which should be the endogenous result of the bidding process among agents in a local market platform. Future work should relax this assumption, enabling a more dynamic and realistic modeling of price mechanisms in these markets.

Another critical area for future investigation is the incorporation of inherent uncertainties, such as solar irradiance and demand fluctuations, into our model. Addressing this right–hand side uncertainty, particularly in the lower level of the bilevel problem, is paramount. Drawing inspiration from robust optimization methods successfully applied in voltage regulation problems within active distribution networks, we propose exploring these methods within bilevel robust optimization. This approach can potentially enhance the resilience and efficacy of our model in the face of uncertainty.

## References

Anjos, M. F., Brotcorne, L., and Gomez-Herrera, J. A. (2021). Optimal setting of time-and-level-of-use prices for an electricity supplier. *Energy*, 225:120517.

Askeland, M., Backe, S., Bjarghov, S., and Korpäs, M. (2021). Helping end-users help each other: Coordinating development and operation of distributed resources through local power markets and grid tariffs. *Energy Economics*, 94:105065.

Askeland, M., Burandt, T., and Gabriel, S. A. (2020). A stochastic mpec approach for grid tariff design with demand-side flexibility. *Energy Systems*.

Aussel, D., Brotcorne, L., Lepaul, S., and von Niederhäusern, L. (2020). A trilevel model for best response in energy demand-side management. *European Journal of Operational Research*, 281(2):299–315.

Baran, M. and Wu, F. (1989). Network reconfiguration in distribution systems for loss reduction and load balancing. *IEEE Transactions on Power Delivery*, 4(2):1401–1407.

Ben-Tal, A. and Nemirovski, A. (2001). On polyhedral approximations of the second-order cone. *Mathematics of Operations Research*, 26(2):193–205.

Bernstein, A., Wang, C., Dall'Anese, E., Le Boudec, J.-Y., and Zhao, C. (2018). Load flow in multiphase distribution networks: Existence, uniqueness, non– singularity and linear models. *IEEE Transactions on Power Systems*, 33(6):5832–5843.

Bozionek, J., Wolgast, T., and Nieße, A. (2022). Design and evaluation of a multi-level reactive power market. *Energy Informatics*, 5(1):6.

Caruso, F., Lignola, M. B., and Morgan, J. (2020). *Regularization and Approximation Methods in Stackelberg Games and Bilevel Optimization*, pages 77–138. Springer International Publishing, Cham.

Coffrin, C. and Van Hentenryck, P. (2014). A linear-programming approximation of ac power flows. *INFORMS Journal on Computing*, 26(4):718–734.

Dall'Anese, E., Dhople, S. V., Johnson, B. B., and Giannakis, G. B. (2015). Optimal Dispatch of Residential Photovoltaic Inverters Under Forecasting Uncertainties. *IEEE Journal of Photovoltaics*, 5(1):350–359.

Dempe, S. (2002). *Parametric Optimization*, pages 61–118. Springer US, Boston, MA.

Dempe, S. and Dutta, J. (2012). Is bilevel programming a special case of a mathematical program with complementarity constraints? *Mathematical Programming*, 131(1):37–48.

Facchinei, F. and Kanzow, C. (2007). Generalized nash equilibrium problems. *4OR*, 5(3):173–210.

Fischer, A., Herrich, M., and Schönefeld, K. (2014). Generalized nash equilibrium problems - recent advances and challenges. *Pesquisa Operacional*, 34:521–558.

Flegel, M. and Kanzow, C. (2005). Abadie-type constraint qualification for mathematical programs with equilibrium constraints. *Journal of Optimization Theory and Applications*, 124:595–614.

Groot, N., De Schutter, B., and Hellendoorn, H. (2016). Optimal Affine Leader Functions in Reverse Stackelberg Games: Existence Conditions and Characterization. *Journal of Optimization Theory and Applications*, 168(1):348–374.

Harker, P. T. (1991). Generalized nash games and quasi-variational inequalities. *European Journal of Operational Research*, 54(1):81–94.

Hart, W., Laird, C., Watson, J., Woodruff, D. L., Hackebeil, G., Nicholson, B., and Siirola, J. D. (2021). *Mathematical Programs with Equilibrium Constraints*. Pyomo — Optimization Modeling in Python.

Ho, Y.-C., Luh, P., and Muralidharan, R. (1981). Information structure, stackelberg games, and incentive controllability. *IEEE Transactions on Automatic Control*, 26(2):454–460.

Ho, Y.-c., Luh, P. B., and Olsder, G. J. (1980). A control-theoretic view on incentives. In *1980 19th IEEE Conference on Decision and Control including the Symposium on Adaptive Processes*, pages 1160–1170.

Hoarau, Q. and Perez, Y. (2019). Network tariff design with prosumers and electromobility: Who wins, who loses? *Energy Economics*, 83:26–39.

IEEE (1991). IEEE tests feeders. https://site.ieee.org/pes-testfeeders/resources/.

Jin, X., Baker, K., Christensen, D., and Isley, S. (2017). Foresee: A user–centric home energy management system for energy efficiency and demand response. *Applied Energy*, 205:1583–1595.

Labbé, M., Marcotte, P., and Savard, G. (1998). A bilevel model of taxation and its application to optimal highway pricing. *Management Science*, 44(12 PART 1):1608–1622.

Le Cadre, H., Jacquot, P., Wan, C., and Alasseur, C. (2020). Peer-to-peer electricity market analysis: From variational to Generalized Nash Equilibrium. *European Journal of Operational Research*, 282(2):753–771.

Low, S. H. (2014a). Convex relaxation of optimal power flow—part i: Formulations and equivalence. *IEEE Transactions on Control of Network Systems*, 1(1):15–27.

Low, S. H. (2014b). Convex relaxation of optimal power flow—part ii: Exactness. *IEEE Transactions on Control of Network Systems*, 1(2):177–189.
then $(\mathbf{\Phi}_i^k, \mathbf{\phi}_i^k) \in \mathcal{D}_i^k$.

PROOF. For simplicity of notation, ignore the indices $i$ and $k$. We write (33) as

$$E||\mathbf{\Phi}||_F + ||\mathbf{\phi}||_2 \le \frac{M}{E}$$

Recall that the Frobenius norm is an upper bound on the matrix 2-norm. Then,

$$E^2||\mathbf{\Phi}||_2 + E||\mathbf{\phi}||_2 \le M$$

The connection capacity condition (30) states that $E$ is an upper bound on $||\mathbf{w}||_2$, therefore the following is true

$$||\mathbf{w}||_2^2||\mathbf{\Phi}||_2 + ||\mathbf{\phi}||_2||\mathbf{w}||_2 \le M$$

Then, we apply the the Cauchy-Schwarz inequality to obtain

$$||\mathbf{w}||_2||\mathbf{\Phi}\mathbf{w}||_2 + ||\mathbf{\phi}||_2||\mathbf{w}||_2 \le M$$

followed by

$$||\mathbf{w}^\top\mathbf{\Phi}\mathbf{w} + \mathbf{\phi}^\top\mathbf{w}||_2 \le M$$

Consequently,

$$-M \le \mathbf{w}^\top\mathbf{\Phi}\mathbf{w} + \mathbf{\phi}^\top\mathbf{w} \le M$$

□

The formulation of Problem (28) can be further clarified by explicitly parameterizing $\gamma(\cdot)$ in terms of the coefficients $(\mathbf{\Phi}, \mathbf{\phi})$

$$\max_{\mathbf{\Phi},\mathbf{\phi}} J_L(\mathbf{\Phi}, \mathbf{\phi}, \mathbf{u}) \qquad (34a)$$

$$\text{s.t.} \quad \Phi_{pp,i}^k + \Phi_{qq,i}^k \ge \varepsilon \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (34b)$$

$$(\Phi_{pq,i}^k)^2 + \varepsilon \le \Phi_{pp,i}^k \Phi_{qq,i}^k \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (34c)$$

$$E_i ||\mathbf{\Phi}_i^k||_F + ||\mathbf{\phi}_i^k||_2 \le \frac{M_i^k}{E_i} \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (34d)$$

$$\hat{\mathbf{v}} \le \mathbf{H}^k \mathbf{u}^k + \mathbf{h}^k \le \hat{\mathbf{v}} \qquad \forall k \in \mathcal{K} \qquad (34e)$$

$$\mathbf{u}_i \in \arg \min_{\tilde{\mathbf{u}}_i \in \Omega_{F,i}(\mathbf{u}_{-i})} J_F(\mathbf{\Phi}_i, \mathbf{\phi}_i, \tilde{\mathbf{u}}_i, \mathbf{u}_{-i}), \qquad \forall i \in \mathcal{N} \qquad (34f)$$

where the objective function (34a) is concave quadratic, constraints (34b) and (34c) ensure that each matrix $\mathbf{\Phi}_i^k \in \mathbb{R}^{2 \times 2}$ is positive semidefinite, equation (34d) represents the budget conditions, (34e) are the voltage limits, and (34f) models the lower-level.

The difficulty of problem (34) comes from its bilevel structure, due to the inclusion of constraint (34f). To address this, problem (34) can be reformulated into a single-level optimization problem by replacing equation (34f) with its KKT reformulation; see Appendix A for further details. Since the KKT conditions define a closed set, if the feasible set of (34) is not empty, either an optimal solution exists or it is an unbounded optimization problem. As the objective function of the leader depends solely on $(\mathbf{p}, \mathbf{q})$, the objective function is bounded because $(\mathbf{p}, \mathbf{q})$ is bounded by construction. Therefore, if (34) is feasible, then a solution exists.

### 3. Properties

In this section, our concern is to study how the lower-level equilibrium might respond to perturbations in the parameters defining the objective function of prosumers. Specifically, our concern is on studying stability in the sense of Dempe (2002) for the lower level as a parametric optimization problem. First, we show that the GNEP can be reformulated as an optimization problem, classically known as the Nikaido–Isoda (NI) reformulation (Facchinei and Kanzow, 2007). Then we show that the set of solutions to the GNEP, $S(T)$, is strongly stable towards perturbations in the parameters that define it by proving that the NI reformulation is strongly stable. As corollaries, we leverage guarantees of uniqueness of solutions in the net active and reactive power output space.
Luo, Z.-Q., Pang, J.-S., and Ralph, D. (1996). *Mathematical Programs with Equilibrium Constraints*. Cambridge University Press.

Maharjan, S., Zhu, Q., Zhang, Y., Gjessing, S., and Başar, T. (2016). Demand response management in the smart grid in a large population regime. *IEEE Transactions on Smart Grid*, 7(1):189–199.

Pecan Street (2021). Pecan Street Dataport. https://www.pecanstreet.org/dataport/.

Pediaditis, P., Papadaskalopoulos, D., Papavasiliou, A., and Hatziargyriou, N. (2021). Bilevel optimization model for the design of distribution use-of-system tariffs. *IEEE Access*, 9:132928–132939.

Potter, A., Haider, R., Ferro, G., Robba, M., and Annaswamy, A. M. (2023). A reactive power market for the future grid. *Advances in Applied Energy*, 9:100114.

Raghunathan, A. and Biegler, L. (2003). Mathematical programs with equilibrium constraints (mpecs) in process engineering. *Comput. Chem. Eng.*, 27:1381–1392.

Rockafellar, R. T. (1970). *Convex analysis*. Princeton Mathematical Series. Princeton University Press, Princeton, N. J.

Rockafellar, R. T. (2023). Generalized nash equilibrium from a variational analysis perspective. *Set-Valued and Variational Analysis: Theory and Applications*.

Rockafellar, R. T. and Wets, R. J. B. (1998). *Variational Analysis*. Springer Berlin Heidelberg.

Schittekatte, T., Momber, I., and Meeus, L. (2018). Future-proof tariff design: Recovering sunk grid costs in a world where consumers are pushing back. *Energy Economics*, 70:484–498.

Tushar, W., Zhang, J. A., Smith, D. B., Poor, H. V., and Thiébaux, S. (2014). Prioritizing consumers in smart grid: A game theoretic approach. *IEEE Transactions on Smart Grid*, 5(3):1429–1438.

Wolgast, T., Ferenz, S., and Nieße, A. (2022). Reactive power markets: A review. *IEEE Access*, 10:28397–28410.

Ye, J. (2005). Necessary and sufficient optimality conditions for mathematical programs with equilibrium constraints. *Journal of Mathematical Analysis and Applications*, 307:350–369.

Zheng, Y.-P. and Basar, T. (1982). Existence and derivation of optimal affine incentive schemes for stackelberg games with partial information: a geometric approach. *International Journal of Control*, 35(6):997–1011.

## A. The KKT reformulation of the lower-level

In this section, we construct the KKT reformulation of the lower-level GNEP (11). It is well known that if the optimization problem of each player satisfies standard constraint qualifications, solving the GNEP is equivalent to solving the concatenation of the KKT systems of all players Facchinei and Kanzow (2007). We start by constructing an inner polyhedral approximation the feasible dispatch set of each DG (1), allowing us to represent (9) as a polyhedron. In the case $\Omega(\mathbf{u}_{-i})$ is a polyhedron, no matter the choice of $\mathbf{u}_{-i}$, and the objective function of each prosumer is convex quadratic as long as $\Phi_i^k \ge 0$, the problem of each prosumer satisfies the so-called refined Slater constraint qualification (RSCQ), presented as Theorem 28.2 in Rockafellar (1970). As a consequence, the stacked KKT conditions of all followers become necessary and sufficient for a solution $\bar{\mathbf{u}}$ to be a solution to (11), i.e., $\bar{\mathbf{u}} \in S(T)$.

First, we construct a polyhedral inner approximation of (1c) following Ben-Tal and Nemirovski (2001). Therefore, replacing (1c) by the following set of inequalities for all $j = 0, ..., n_{SOC}$

$$-\alpha_{G,ij}p_G - \beta_{G,ij} \le q_G \quad (44a)$$

$$\alpha_{G,ij}p_G + \beta_{G,ij} \ge q_G \quad (44b)$$

Where

$$\beta_{G,ij} = q_{G,i,j}^0 - \alpha_{G,ij}p_{G,ij}^0 \quad (45)$$

$$\alpha_{G,ij} = \frac{q_{G,i,j+1}^0 - q_{G,i,j}^0}{p_{G,i,j+1}^0 - p_{G,i,j}^0} \quad (46)$$

$$(p_{G,ij}^0, q_{G,ij}^0) = \bar{s}_{G,i} \left( \sin\left(\frac{j\pi}{4n}\right), \cos\left(\frac{j\pi}{4n}\right) \right) \quad (47)$$

The KKT conditions of the GNEs are the system of equations (9b) to (9e), (48), and (49).

Similarly, we can leverage the KKT conditions for the welfare equilibria; thereby, obtaining necessary and sufficient conditions for (12). The KKT conditions for the welfare equilibria are the system of equations: (9b) to (9e), (48a) to (48g), (48i), (49), and

$$-\mu_i^k + \mu_j^k - \lambda_{EL,ij}^k = 0 \quad \forall (i, j, k) \in \mathcal{N} \times \mathcal{N}_i \times \mathcal{K} \quad (50a)$$

$$\mu_i^k + c_{E,i0}^k + \lambda_{EL,i0}^k = 0 \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (50b)$$

In other words, the KKT system of the welfare optimization problem is the same as of the GNEs, except that (48h) have to be replaced by (50).

**Proposition 3.** The following is a necessary condition for a solution to (1) to be a solution to (2).

$$\lambda_{EL,ij}^k + \lambda_{EL,ij}^k = 0 \quad \forall (i, j, k) \in \mathcal{N} \times \mathcal{N}_i \times \mathcal{K} \quad (51)$$

*PROOF.* The financial network is a symmetric directed graph, therefore if $(i, j)$ belongs to the edges, so do $(j, i)$. Then, from (50a) we obtain

$$-\mu_j^k + \mu_i^k - \lambda_{EL,ji}^k = 0 \quad \forall (i, j, k) \in \mathcal{N} \times \mathcal{N}_i \times \mathcal{K} \quad (52)$$

Summing (50a) and (52), we obtain (51). □

**Proposition 4.** The following are sufficient conditions for a GNE to be a Welfare equilibrium

$$\mu_i^k = \mu_j^k \quad (53a)$$

$$\lambda_{EL,ij}^k = \lambda_{EL,ji}^k = 0 \quad (53b)$$

*PROOF.* Reorder the terms in (53a) to obtain

$$\mu_i^k - \mu_j^k = 0$$

Add a zero-valued $\lambda_{EL,ij}^k$ on the left hand side to obtain

$$\mu_i^k - \mu_j^k + \lambda_{EL,ij}^k = 0$$

Thereby, we have arrived to (50a), hence the GNE is a welfare equilibrium. □
**Definition 3.** (Jointly convex GNEP) We say that (11) is a jointly convex GNEP if, for every prosumer $i \in \mathcal{N}$ and every $\mathbf{u}_{-i}$, the objective function $J_F(T(\cdot), \cdot, \mathbf{u}_{-i})$ is convex, the set $\Omega_i(\mathbf{u}_{-i})$ is closed and convex, and for some closed convex $\Omega \subset \mathbb{R}^n$ and all $i \in \mathcal{N}$ we have

$$ \Omega_i(\mathbf{u}_{-i}) = \{ \mathbf{u}_i \in \mathbb{R}^{n_i} \mid (\mathbf{u}_i, \mathbf{u}_{-i}) \in \Omega \} $$

As remarked by Facchinei and Kanzow (2007), if the feasible set of each player of a GNEP is described by inequalities and if all the shared constraints, i.e., the ones that involve variables from more than one player are the same for every player involved in the constraint, then the GNEP in question is a jointly convex GNEP. By this remark, it is clear that the GNEP in (11) is jointly convex. By the discussion above we derive the following result

**Lemma 1.** The game defined in (11) is a jointly convex GNEP.

First, we establish a classical result that bridges jointly convex GNEPs and optimization problems

**Lemma 2.** If $T_i$ is convex for each $i \in \mathcal{N}$, a vector $\bar{\mathbf{u}}$ is a solution to (11) if and only if $\bar{\mathbf{u}}$ is a global minimum of the optimization problem

$$ \min_{\mathbf{u}} \hat{V}(T, \mathbf{u}) \qquad (35a) $$

$$ \text{s.t.} \quad \mathbf{u} \in \Omega, \qquad (35b) $$

with zero objective function value, where

$$ \hat{V}(T, \mathbf{u}) := \sup_{\mathbf{y} \in \Omega} \psi(T, \mathbf{u}, \mathbf{y}) \qquad (36) $$

and

$$ \psi(T, \mathbf{u}, \mathbf{y}) := \sum_{i \in \mathcal{N}} (J_{F,i}(T_i, \mathbf{u}_i, \mathbf{u}_{-i}) - J_{F,i}(T_i, \mathbf{y}_i, \mathbf{u}_{-i})) \qquad (37) $$

Lemma 2 corresponds to the classical Nikaido–Isoda (NI) reformulation in the case of jointly convex GNEPs, hence the proof is omitted (Facchinei and Kanzow, 2007). The NI reformulation is an important analytical tool that allows to transform the equilibrium in (11) into an optimization problem.

**Lemma 3.** If $T_i(\cdot)$ is convex for each $i \in \mathcal{N}$, problem (35) is a convex optimization problem [ ].

**PROOF.** Note that the objective function of each follower can be expressed as the sum

$$ J_F(T_i(\cdot), \mathbf{u}_i, \mathbf{u}_{-i}) = h_i(\mathbf{u}_i, T_i(\cdot)) + g_i(\mathbf{u}_{-i}), \quad \forall i \in \mathcal{N} \qquad (38) $$

with

$$ h_i(\mathbf{u}_i, T_i(\cdot)) := \sum_{k \in \mathcal{K}} \left( T_i^k(p_i^k, q_i^k) + l_i^k(p_{G,i}^k, q_{G,i}^k) + c_{0i}^k p_{E,0i}^k - \sum_{j \in \mathcal{N}_i^0} c_{ij}^k p_{E,ij} \right) $$

and

$$ g_i(\mathbf{u}_{-i}) := \sum_{k \in \mathcal{K}} \sum_{j \in \mathcal{N}} c_{ji}^k p_{E,ji}^k $$

Then the NI–function becomes

$$ \psi(T(\cdot), \mathbf{u}, \mathbf{y}) = \sum_{i \in \mathcal{N}} (h_i(\mathbf{u}_i, T_i(\cdot)) - h_i(\mathbf{y}_i, T_i(\cdot))) \qquad (39) $$

For fixed $T(\cdot)$, function $\psi(T(\cdot), \mathbf{u}, \mathbf{y})$ is convex in $\mathbf{u}$, for all $\mathbf{y} \in \Omega$, therefore $\hat{V}(T(\cdot), \mathbf{u})$ is convex in $\mathbf{u}$ as the supremum operation in (36) preserves convexity. $\square$
Stationarity conditions:

$$ \Phi_{pp,i}^k p_i^k + \Phi_{pq,i}^k q_i^k + \phi_{p,i}^k + \mu_i^k + \nu_{P,i}^k = 0 \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (48a) $$

$$ \Phi_{qq,i}^k q_i^k + \Phi_{pq,i}^k p_i^k + \phi_{q,i}^k + \nu_{Q,i}^k = 0 \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (48b) $$

$$ \nu_{P,i}^k + \alpha_{G,i} \sum_{j \in W} (\lambda_{GL,i,j}^k + \lambda_{GU,i,j}^k) = 0 \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (48c) $$

$$ \nu_{P,i}^k - \alpha_{G,i} \sum_{j \in W} (\lambda_{GU,i,j}^k - \lambda_{GL,i,j}^k) = 0 \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (48d) $$

$$ \nu_{P,i}^k + \eta_{B,i} \nu_{B,i}^k + \lambda_{PBU,i}^{+k} - \lambda_{PBL,i}^{+k} = 0 \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (48e) $$

$$ \nu_{P,i}^k - \frac{1}{\eta_{B,i}} \nu_{B,i}^k + \lambda_{PBU,i}^{-k} - \lambda_{PBL,i}^{-k} = 0 \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (48f) $$

$$ \nu_{B,i}^k - \nu_{B,i}^{k+1} + \lambda_{SBU,i}^k - \lambda_{SBL,i}^k = 0 \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K}' \quad (48g) $$

$$ c_{E,i,j}^k + \mu_i^k + \lambda_{EL,i,j}^k = 0 \quad \forall (i, j, k) \in \mathcal{N} \times \mathcal{N}_i^0 \times \mathcal{K} \quad (48h) $$

$$ c_{E,0i}^k + \mu_i^k - \lambda_{EL,0i}^k = 0 \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (48i) $$

Complementarity slackness conditions:

$$ \lambda_{GL,i,j}^k \ge 0 \perp -\alpha_{G,i,j} p_G - \beta_{G,i,j} \le q_G \quad \forall (i, j, k) \in \mathcal{N} \times W_i \times \mathcal{K} \quad (49a) $$

$$ \lambda_{GU,i,j}^k \ge 0 \perp \alpha_{G,i,j} p_G + \beta_{G,i,j} \ge q_G \quad \forall (i, j, k) \in \mathcal{N} \times W_i \times \mathcal{K} \quad (49b) $$

$$ \lambda_{SL,i}^k \ge 0 \perp s_{B,i}^k \ge 0 \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (49c) $$

$$ \lambda_{SU,i}^k \ge 0 \perp s_{B,i}^k \le \hat{s}_{B,i} \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (49d) $$

$$ \lambda_{PBL,i}^{+k} \ge 0 \perp p_{B,i}^{+k} \ge 0 \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (49e) $$

$$ \lambda_{PBL,i}^{-k} \ge 0 \perp p_{B,i}^{-k} \ge 0 \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (49f) $$

$$ \lambda_{PBU,i}^{+k} \ge 0 \perp p_{B,i}^{+k} \le \hat{p}_{B,i} \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (49g) $$

$$ \lambda_{PBU,i}^{-k} \ge 0 \perp p_{B,i}^{-k} \le \hat{p}_{B,i} \quad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (49h) $$

$$ \lambda_{EL,i,j}^k \ge 0 \perp p_{E,i,j}^k \ge 0 \quad \forall (i, j, k) \in i \in \mathcal{N}_i^0 \times \mathcal{N}_i^0 \times \mathcal{K} \quad (49i) $$

**Proposition 5.** The following are sufficient conditions for a GNE to be a Welfare equilibrium

$$ c_{E,i,j}^k = c_{E,ji}^k \quad (54a) $$

$$ \lambda_{EL,i,j}^k = \lambda_{EL,ji}^k = 0 \quad (54b) $$

PROOF. From (48h) and (54a) we have

$$ \mu_i^k + \lambda_{EL,i,j}^k = \mu_j^k + \lambda_{EL,ji}^k $$

Therefore

$$ \mu_i^k = \mu_j^k $$

and by Proposition 4, the GNE is also a welfare equilibrium.

Propositions 3 to 5 characterize, to some extent, the two classes of market designs considered in this work.
**Theorem 1.** If $T(\cdot)$ is convex, the set of GNEs $S(T)$ of (11) is convex. Moreover, if $T(\cdot)$ is strictly convex, the set of GNEs of (11), $S(T)$, is uniquely defined in $(\mathbf{p}, \mathbf{q})$.

**PROOF.** By Lemma 3, problem (35) is a convex optimization problem. Therefore, the set of optimal solutions is a convex set. By Lemma 2, the solution set of (35) corresponds to $S(T)$ if the objective value of (35) is zero at optimality. Therefore, to prove Theorem 1, it suffices to prove that the optimal value of (35) is zero for any given $T(\cdot)$. For that, we express the NI-function as in (39), and considering $\Omega$ is closed, the optimal value of (35) can be written as

$$ \hat{V}^* = \min_{\mathbf{u} \in \Omega} \left\{ \max_{\mathbf{y} \in \Omega} \sum_{i \in \mathcal{N}'} (h_i(\mathbf{u}_i, T_i(\cdot)) - h_i(\mathbf{y}_i, T_i(\cdot))) \right\} \quad (40) $$

which is zero. Therefore $S(T)$ is convex.

If $T(\cdot)$ is strictly convex, (39) is strictly convex with respect to $(\mathbf{p}, \mathbf{q})$, therefore the NI reformulation (35) is unique in $(\mathbf{p}, \mathbf{q})$, thus, by Lemma 2, so does the set $S(T)$. $\square$

The DSO exclusively considers the net active and reactive power output, denoted as $(\mathbf{p}, \mathbf{q})$, as the operational performance metrics hinge on these variables. Hence, our focus is on ensuring the robustness of the solution, specifically within the $(\mathbf{p}, \mathbf{q})$-space, against variations in lower-level model parameters. Thus, providing a robust stability guarantee in this projected space suffices.

To ensure robustness of the solution against parameter variations in the lower-level problem, we extend our approach beyond merely ensuring uniqueness to focus the notion of stability. Our approach involves drawing upon the concept of strong stability for local solutions, as defined by Dempe (2002), which we paraphrase as follows

**Definition 4 (Dempe (2002)).** A local solution $(\bar{\mathbf{x}}, \bar{\mathbf{y}})$ of a parametric optimization problem

$$ \min_{\mathbf{x}} f(\mathbf{x}, \mathbf{y}) \quad (41a) $$

$$ \text{s.t.} \quad g(\mathbf{x}, \mathbf{y}) \le 0 \quad (41b) $$

is said to be strongly stable if there exist neighborhoods $\mathcal{U}_{\delta}(\bar{\mathbf{y}})$ with $\delta > 0$ of $\bar{\mathbf{y}}$ and $\mathcal{V}_{\epsilon}(\bar{\mathbf{x}})$ with $\epsilon > 0$ of $\bar{\mathbf{x}}$, and a uniquely determined continuous vector-valued function $x: \mathcal{U}_{\delta}(\bar{\mathbf{y}}) \to \mathcal{V}_{\epsilon}(\bar{\mathbf{x}})$ such that $x(y)$ is the unique local optimal solution of problem (41) in $\mathcal{V}_{\epsilon}(\bar{\mathbf{x}})$ for all $y \in \mathcal{U}_{\delta}(\bar{\mathbf{y}})$.

First, we introduce a projection operator. Let $\bar{\mathbf{u}} \in \Omega$ represent a solution at the lower level, with $\bar{\mathbf{u}}_i^k := (\bar{p}_i^k, \bar{q}_i^k, \bar{p}_{G,i}^k, \bar{q}_{G,i}^k, \bar{p}_{B,i}^{+k}, \bar{p}_{B,i}^{-k}, \bar{s}_{B,i}^k, \bar{p}_{E,i}^k)$. We define $\Pi_{pq}(\bar{\mathbf{u}}_i^k) := (\bar{p}_i^k, \bar{q}_i^k)$ and extend the notation to aggregate variables, such that $\Pi_{pq}(\bar{\mathbf{u}}) := (\bar{\mathbf{p}}, \bar{\mathbf{q}})$.

**Proposition 2.** If $(\bar{\mathbf{\Phi}}, \bar{\mathbf{\phi}}, \bar{\mathbf{u}})$ is an optimal solution to (34) and $\bar{\mathbf{w}} = \Pi_{pq}(\bar{\mathbf{u}})$, then the projected solution $(\bar{\mathbf{\Phi}}, \bar{\mathbf{\phi}}, \bar{\mathbf{w}})$ is strongly stable.

**PROOF.** First, we combine (36) and (39) to write

$$ \hat{V}(T(\cdot), \mathbf{u}) = \sup_{\mathbf{y} \in \Omega} \sum_{i \in \mathcal{N}'} h_i(T(\cdot), \mathbf{u}) - h_i(T(\cdot), \mathbf{y}) \quad (42) $$

Then, note that the $\sup\{\cdot\}$ taken only over $y$, therefore we can write

$$ \hat{V}(T(\cdot), \mathbf{u}) = \sum_{i \in \mathcal{N}} h_i(T(\cdot), \mathbf{u}) + C \quad (43) $$

where $C$ is constant with respect to $\mathbf{u}$, given by

$$ C := \max_{\mathbf{y} \in \Omega} \left\{ - \sum_{i \in \mathcal{N}} h_i(T(\cdot), \mathbf{y}) \right\} $$
thus, $\hat{V}(T(\cdot), u)$ is differentiable in $u$. In fact, by applying (43), the partial hessian matrix of the lagrangian function of the NI reformulation (35) with respect to $w$ is a block-diagonal matrix

$$ \frac{\partial^2 L}{\partial w^2} = \mathrm{diag} \left( \left( \overline{\Phi}_i \right)^k \right)_{i \in \mathcal{N}, k \in \mathcal{K}} $$

Therefore $\frac{\partial^2 L}{\partial w^2}$ is positive definite because each block diagonal element is positive definite as $\epsilon > 0$; thus $(\overline{\Phi}, \overline{\phi}, \overline{w})$ satisfies the strong sufficient optimality condition of second order (Dempe, 2002). In addition, (35) satisfies the Mangasarian-Fromowitz constraint qualification (MFCQ). Then, Proposition 2 follows directly from Theorem 4.4 in Dempe (2002). □

**Proposition 2** stands as a robustness guarantee, ensuring that the tariff solution remains reliable in the face of parameter variations of the lower-level. Note that stability defined in this sense implies uniqueness of the market output with respect to $(p, q)$.

Furthermore, the stability property not only applies for the GNEP-based market setup. It remains true for both when the market is modeled as a GNEP and when it is modeled as an optimization problem that aims to minimize total disutility, as presented in Definition 2. To see this, note that if $T(\cdot)$ is strictly convex, then (12) is a strictly convex problem with respect to $(p, q)$, and the same line of analysis developed in this section can be applied. This implies that the stability condition we have established is equally applicable to both managed and unmanaged market designs, solidifying the versatility of the proposed model.

We highlight that the NI reformulation is a valuable tool for characterizing the stability of jointly convex GNEPs. Note that by characterizing the stability of the solution of the NI reformulation, we bypass the notions of stability in variational systems such as Rockafellar (2023) which are commonly used in more general variational systems when there is no optimization problem reformulation at hand.

Finally, it is important to highlight that the objective function (34a) and the coupling constraints (34e) in the tariff design problem (34) depend exclusively on $w = (p, q)$. Consequently, by Theorem 1, problem (34) is well defined, leading to the following proposition:

**Proposition 3.** If the solution set of (29) is non-empty, then every solution yields the same objective function value for the leader.

**PROOF.** By (23), the objective function of the leader is expressed in terms of $w$. Furthermore, by Theorem 1, and representing $\gamma$ using parameters $(\Phi, \phi)$, problem (29) can be equivalently formulated as the following MPEC

$$ \max_{(\Phi, \phi, w)} \{ J_L(w) : (\Phi, \phi, w) \text{ satisfy (34b) to (34e), } w = w^*(\Phi, \phi) \}, \quad (44) $$

where $w^*(\Phi, \phi)$ corresponds to the projection on $w$ of the set of GNEs, which defines a singleton written as

$$ w_i^*(\Phi_i, \phi_i) = \Pi_{pq} \left( \arg \min_{u_i \in \Omega_{F,i}} J_F(\Phi_i, \phi_i, u_i, \bar{u}_{-i}) \right), \quad \forall i \in \mathcal{N}. \quad (45) $$

Since (44) is an optimization problem, it is clear that any optimal solution yields the same objective value. □

## 4. Numerical results

An instance of the problem comprises a distribution network, the prosumers with their devices, and the financial grid. For each instance, the computational experiments follow a scheme in which the tariff is calculated by the DSO and passed to the prosumers in a rolling horizon fashion<fnref n="2" />. With our approach, the leader and followers plan their operations several steps ahead but make decisions incrementally. Tariff functions for all prosumers are calculated considering a planning horizon of 4 hours and applied only for 1 hour. We run the experiments utilizing a 10-minute time step setting; therefore, each instance of (34) to be solved has $|\mathcal{K}| = 6 \times 4 = 24$ time steps. Optimal strategies for both leader and

<footnote n="2">²Also called receding horizon.</footnote>
followers are calculated under the rolling horizon setting, this allows for prosumers to consider future variations in solar irradiance and manage their batteries accordingly.

To compute the optimal affine strategies ($\Phi, \phi$), problem (34) is solved by replacing the lower-level (34f) by its KKT reformulation depicted in Appendix A. This results in a quadratically constrained quadratic program with complementarity constraints. Globally optimal solutions are obtained for the toy problem case in Section (4.1) by applying the spatial branching capabilities of Gurobi version 10.0.1 (Gurobi Optimization, LLC, 2023); whereas in Section 4.2, only locally optimal solutions are obtained utilizing IpOpt (Wächter and Biegler, 2006). The experiments were conducted on a computer equipped with an Intel(R) Xeon(R) CPU E7-4809 v3 @ 2.0 GHz, with 32 physical CPUs and 64 GB of RAM.

To evaluate the proposed approach against a more conventional one, we also consider an alternative constant tariff design. This involves solving (34) with $\Phi_i^k = 0$ and removing the positive definiteness constraints (34b) and (34c). As a result, applying the KKT reformulation simplifies the problem to a quadratic program with complementarity constraints. Since the constant tariff lacks the stability property of the proposed approach, this allows us to compare a stable approach with one that does not hold that property.

## 4.1. A toy problem

This section is dedicated to leveraging a toy problem, with the objective of gaining deeper insights into the proposed tariff design and offering a tangible proof of concept. First, we show how the DSO is able to induce safe operation by prompting the voltage profiles to remain within nominal operation bounds. Then, we establish the main differences between the proposed approach and the conventional constant tariff, which involves only bilinear terms between leader and follower variables.

The experimental setup for this toy problem involves a set $\mathcal{N} = \{1, 2, 3\}$ of prosumers connected to a 4 bus electrical grid where the set of buses is $\mathcal{B}^0 = \{0, 1, 2, 3\}$, with parent buses $A_1 = 0, A_2 = 1$, and $A_3 = 0$. The set of prosumers is $\mathcal{N} = \{1, 2, 3\}$ such that $\mathbf{B} = \text{diag}(1, 1, 1)$ as depicted in Figure 1. The distribution system corresponds to the IEEE 4-bus test feeder (IEEE, 1991) which has been adapted to this setting by balancing the line segments using Carsons equation and Kron reduction. Distributed loads have been converted into balanced three-phase loads. The resistance of each electrical branch has been increased so that the resistance to reactance ratio is exactly 2.3—a normal ratio in distribution systems.—In this case, each prosumer model represents a group of actual prosumers. This is because the IEEE 4-bus system is a simplified electrical equivalent of a whole distribution system. Therefore, to analyze the effect of the prosumers in the operation of the grid, it is necessary to aggregate the effect of a large group of them.

The time series data for the problem is depicted in Figure 4. On top is the maximum available solar power, and at the bottom is the active power load. The signals have been obtained from one-minute time series data from a set of individual homes from the U.S. Austin region (Pecan Street, 2021). The time-series have been upsampled, averaging the values over periods of 10 minutes to accommodate them to the current time setting. A constant power factor of 0.85 is considered; therefore, the reactive power load is proportional to the active power load by a constant factor of $\tan(\cos^{-1}(0.85)) \approx 0.61$.

Naturally, the high solar irradiance period occurs around midday. During this period, there are high chances of violating the upper voltage limit because the net active power injections increase and $\mathbf{R}$ in (19a) is a positive definite matrix that maps the increment in power injection into an increment of voltage magnitudes. Without the tariff incentive, the voltage magnitude of electrical Node 2 would have been out of bounds. Figure 5 depicts the voltage profiles considering the tariff. Note that the voltage profile gets flattened to remain within bounds during the high irradiance interval. This is achieved by inducing the active and reactive net power profiles shown in figures 6 and 7a. A similar process occurs during the high-demand periods in the evening, where voltages tend to lower, and they get flattened as a result of the tariff to avoid a lower limit voltage violation.

In Figure 6, we observe that in time periods near midday, the induced net power injection depicted in orange has less noise than the blue line representing the maximum available net active power. Recall that the tariff gets updated and passed to the prosumers hourly. The red vertical lines show the moments in which the tariff is recomputed and passed to the prosumers. Notably, within these intervals, the net active power injection maintains a near-constant level despite the underlying potential for variation. This observation suggests a denoising effect associated with the tariff updates, which likely contributes to enhancing the predictability of the system's operation—a beneficial characteristic for operating the distribution network.
<figure><img src="image_6.png" /><figcaption>**Figure 4:** Time series data of maximum available power (top), and active power load (bottom).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5:** Voltage magnitude per bus.</figcaption></figure>

## 4.2. Experiments on larger instances

In evaluating our proposed method, we look at various metrics and extend our tests to larger scenarios based on the IEEE 4-bus and 34-bus systems. We adopt the same rolling horizon approach as we did in the earlier section solving 168 instances of the tariff design problem, one each hour over a week. We considered the following performance metrics:

* **Hourly Average Number of Voltage Violations per Bus (NVV):** This is the count of average voltage magnitude deviations from the acceptable range per bus, calculated over 10-minute intervals. Ideally, NVV should be as low as possible, with $NVV = 0$ indicating that the given policy effectively regulates voltage within the acceptable range.

* **Hourly Average Net Energy Export per Prosumer (NEE) in kW:** We measure the total energy exported to the grid, divided by the evaluation period and the number of prosumers. A higher NEE is preferable, as it indicates that the DSO is effectively minimizing the amount of energy imported from the external grid.

* **Hourly Average Distribution Grid Losses (DGL) in kW:** This represents the total resistive losses in the grid, averaged over the evaluation period. Lower DGL values are preferable as they indicate reduced energy losses and associated transmission costs.
<figure><img src="image_8.png" /><figcaption>**Figure 6:** Net active power injection at maximum power point (MPP) vs. the induced signal.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7:** Reactive power generation of a single prosumer under proposed and constant tariffs.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8:** Cumulative probability density function objective function of the DSO, IEEE-4bus test feeder.</figcaption></figure>

* **Hourly Average Disutility per Prosumer (DPP) in $:** This is the sum of each prosumer's disutility, averaged over the time of the evaluation. Lower DPP values are desirable as they reflect lower costs for prosumers, indicating higher satisfaction.

* **Hourly Average Tariff Cost per Prosumer (TC) in $:** We accumulate the tariff costs from the prosumers' objective function over each period and then average it out. Lower TC values are beneficial for prosumers as they reduce their costs, though they may be less favorable for the decision maker. Nevertheless, note that the objective function of the leader does not consider incentive costs.
<table><caption>**Table 1**
Mean value of selected performance metrics.</caption><thead><tr><th rowspan="2">Strategies</th><th colspan="5">IEEE 4-bus test feeder</th><th colspan="5">IEEE 34-bus test feeder</th></tr><tr><th>NVV</th><th>NEE</th><th>DGL</th><th>DPP</th><th>TC</th><th>NVV</th><th>NEE</th><th>DGL</th><th>DPP</th><th>TC</th></tr></thead><tbody><tr><td>No tariff</td><td>0.087</td><td>0.351</td><td>0.028</td><td>-0.300</td><td>0.</td><td>0.101</td><td>0.118</td><td>0.089</td><td>-0.106</td><td>0.</td></tr><tr><td>Constant tariff</td><td>0.</td><td>0.218</td><td>0.022</td><td>-0.119</td><td>0.088</td><td>0.</td><td>0.073</td><td>0.065</td><td>-0.030</td><td>0.078</td></tr><tr><td>Proposed</td><td>0.</td><td>0.339</td><td>0.027</td><td>-0.292</td><td>-0.022</td><td>0.</td><td>0.108</td><td>0.077</td><td>-0.126</td><td>-0.056</td></tr></tbody></table>

Upon examining the NVV and NEE columns in Table 1, a comparison between the proposed approach and the constant tariff policy reveals a notable observation. While both strategies ensure the system operates safely (NVV = 0), the proposed tariff design stands out by achieving this objective alongside showcasing levels of net energy export that surpass the constant tariff approach in about 50%. However, it is important to note that the proposed approach does incur higher levels of ohmic losses. This is due to the higher energy exports, which naturally increase transmission losses, as the tariff design problem itself favors exporting over minimizing losses.

With the focus on operation performance, we did not consider the tariff budget in the objective function of the leader. Note that in Table 1, the TC metric changes signs from the proposed tariff to the constant tariff. The results are likely to be different if considering the budget upperbound $M_i$ for $i \in \mathcal{N}$ in the objective function of the leader, which can be done without changing the structure of the problem. Also, we computed the average resolution time across all IEEE 34-bus instances. The results show that the proposed approach, while improving system performance, is approximately 10 times slower than the constant tariff method due to its increased computational complexity.

While these metrics offer insights into different aspects of performance, our primary focus is on minimizing the DSO's objective function, as it is indicative of the overall effectiveness of our model. The graph in Figure 8 illustrates the cumulative probability distribution of the DSO's objective function values across the three strategies we have put to the test.

The profiles of reactive power power injection depicted in Figure 7 serve as an explanation for the difference in performance observed between the conventional constant tariffs and the proposed ones. Note that the proposed tariff showcases more degrees of freedom in inducing desired values for reactive power injection, as shown in Figure 7a, whereas the points that are inducible by the constant tariffs belong to a discrete set. In the case of constant tariffs, values that do not belong to either $\{\hat{s}_{B,i}^k, 0, \hat{s}_{B,i}^k\}$ are unstable equilibria. It is also important to note that Figure 7 was generated using a box inner approximation of $C_i^k$, so the lower and upper bounds, shown in red, are constant.

## 5. Conclusions

In our research, we have successfully demonstrated the efficacy of a reverse Stackelberg game approach to demand response in power distribution networks. The proposed framework positions the DSO as the leader with the ability to impose a tariff function, to which the prosumers, interpreted as followers, react. This framework also encompasses a peer-to-peer market among the prosumers, modeled as a GNEP, and takes into account both active and reactive power generation and their effects on the network's voltage levels. This comprehensive approach addresses a commonly neglected aspect in current research, namely, the impact of market outcomes on the operation of the distribution network.

The numerical experiments indicate a marked improvement in maximizing the net power output and maintaining the voltage levels within safe margins, compared to the conventional constant time-of-use tariff, particularly when perturbations in the followers' objective function parameters are introduced. To understand the reason behind the effectiveness in the presence of perturbations in the parameters of the lower-level, we provide theoretical results ensuring that the obtained solutions are always strongly stable. Moreover, our method stands out in its ability to denoise net active and reactive power injection of prosumers, which is also attributed to the stability property mentioned above. These findings underscore the potential of our bilevel programming solution to not only streamline the objective function of the DSO but also to bring consistency and predictability to prosumer behavior, thereby reinforcing the overall efficiency and reliability of local energy markets.
We highlight the application of the Nikaido–Isoda reformulation to the lower–level GNEP to characterize the set of GNEs. Looking forward, we believe this reformulation holds significant potential for practical implementation, particularly in solving Stackelberg games with a GNEP at the lower level through a branch–and–cut strategy.

One limitation of the current approach is the assumption of fixed and known peer-to-peer prices, which should be the endogenous result of the bidding process among agents in a local market platform. Future work should relax this assumption, enabling a more dynamic and realistic modeling of price mechanisms in these markets.

Another critical area for future investigation is the incorporation of inherent uncertainties, such as solar irradiance and demand fluctuations, into our model. Addressing this right-hand side uncertainty, particularly in the lower level of the bilevel problem, is paramount. Drawing inspiration from robust optimization methods successfully applied in voltage regulation problems within active distribution networks, we propose exploring these methods within bilevel robust optimization. This approach can potentially enhance the resilience and efficacy of our model in the face of uncertainty.

## References

Anjos, M. F., Brotcorne, L., and Gomez-Herrera, J. A. (2021). Optimal setting of time-and-level-of-use prices for an electricity supplier. *Energy*, 225:120517.

Askeland, M., Backe, S., Bjarghov, S., and Korpås, M. (2021). Helping end-users help each other: Coordinating development and operation of distributed resources through local power markets and grid tariffs. *Energy Economics*, 94:105065.

Askeland, M., Burandt, T., and Gabriel, S. A. (2020). A stochastic mpec approach for grid tariff design with demand-side flexibility. *Energy Systems*.

Aussel, D., Brotcorne, L., Lepaul, S., and von Niederhäusern, L. (2020). A trilevel model for best response in energy demand-side management. *European Journal of Operational Research*, 281(2):299–315.

Baran, M. and Wu, F. (1989). Network reconfiguration in distribution systems for loss reduction and load balancing. *IEEE Transactions on Power Delivery*, 4(2):1401–1407.

Bernstein, A., Wang, C., Dall'Anese, E., Le Boudec, J.-Y., and Zhao, C. (2018). Load flow in multiphase distribution networks: Existence, uniqueness, non-singularity and linear models. *IEEE Transactions on Power Systems*, 33(6):5832–5843.

Bozionek, J., Wolgast, T., and Nieße, A. (2022). Design and evaluation of a multi-level reactive power market. *Energy Informatics*, 5(1):6.

Caruso, F., Lignola, M. B., and Morgan, J. (2020). *Regularization and Approximation Methods in Stackelberg Games and Bilevel Optimization*, pages 77–138. Springer International Publishing, Cham.

Coffrin, C. and Van Hentenryck, P. (2014). A linear-programming approximation of ac power flows. *INFORMS Journal on Computing*, 26(4):718–734.

Dall'Anese, E., Dhople, S. V., Johnson, B. B., and Giannakis, G. B. (2015). Optimal Dispatch of Residential Photovoltaic Inverters Under Forecasting Uncertainties. *IEEE Journal of Photovoltaics*, 5(1):350–359.

Dempe, S. (2002). *Parametric Optimization*, pages 61–118. Springer US, Boston, MA.

Dempe, S. and Dutta, J. (2012). Is bilevel programming a special case of a mathematical program with complementarity constraints? *Mathematical Programming*, 131(1):37–48.

Facchinei, F. and Kanzow, C. (2007). Generalized nash equilibrium problems. *4OR*, 5(3):173–210.

Finnah, B., Gönsch, J., and Ziel, F. (2022). Integrated day-ahead and intraday self-schedule bidding for energy storage systems using approximate dynamic programming. *European Journal of Operational Research*, 301(2):726–746.

Fischer, A., Herrich, M., and Schönefeld, K. (2014). Generalized nash equilibrium problems - recent advances and challenges. *Pesquisa Operacional*, 34:521–558.

Flegel, M. and Kanzow, C. (2005). Abadie-type constraint qualification for mathematical programs with equilibrium constraints. *Journal of Optimization Theory and Applications*, 124:595–614.

Groot, N., De Schutter, B., and Hellendoorn, H. (2016). Optimal Affine Leader Functions in Reverse Stackelberg Games: Existence Conditions and Characterization. *Journal of Optimization Theory and Applications*, 168(1):348–374.

Guerrero, J., Chapman, A. C., and Verbi, G. (2019). Decentralized p2p energy trading under network constraints in a low-voltage network. *IEEE Transactions on Smart Grid*, 10(5):5163–5173.

Gurobi Optimization, LLC (2023). Gurobi Optimizer Reference Manual.

Harker, P. T. (1991). Generalized nash games and quasi-variational inequalities. *European Journal of Operational Research*, 54(1):81–94.

Hart, W., Laird, C., Watson, J., Woodruff, D. L., Hackebeil, G., Nicholson, B., and Siirola, J. D. (2021). *Mathematical Programs with Equilibrium Constraints*. Pyomo Optimization Modeling in Python.

Ho, Y.-C., Luh, P., and Muralidharan, R. (1981). Information structure, stackelberg games, and incentive controllability. *IEEE Transactions on Automatic Control*, 26(2):454–460.

Ho, Y.-c., Luh, P. B., and Olsson, G. J. (1980). A control-theoretic view on incentives. In *1980 19th IEEE Conference on Decision and Control including the Symposium on Adaptive Processes*, pages 1160–1170.

Hoarau, Q. and Perez, Y. (2019). Network tariff design with prosumers and electromobility: Who wins, who loses? *Energy Economics*, 83:26–39.
*IEEE* (1991). IEEE tests feeders. https://site.ieee.org/pes-testfeeders/resources/.

Jin, X., Baker, K., Christensen, D., and Isley, S. (2017). Foresee: A user-centric home energy management system for energy efficiency and demand response. *Applied Energy*, 205:1583–1595.

Labbé, M., Marcotte, P., and Savard, G. (1998). A bilevel model of taxation and its application to optimal highway pricing. *Management Science*, 44(12 PART 1):1608–1622.
Le Cadre, H., Jacquot, P., Wan, C., and Alasseur, C. (2020). Peer-to-peer electricity market analysis: From variational to Generalized Nash Equilibrium. *European Journal of Operational Research*, 282(2):753-771.

Low, S. H. (2014a). Convex relaxation of optimal power flowpart i: Formulations and equivalence. *IEEE Transactions on Control of Network Systems*, 1(1):15-27.

Low, S. H. (2014b). Convex relaxation of optimal power flowpart ii: Exactness. *IEEE Transactions on Control of Network Systems*, 1(2):177-189.

Luo, Z.-Q., Pang, J.-S., and Ralph, D. (1996). *Mathematical Programs with Equilibrium Constraints*. Cambridge University Press.

Maharjan, S., Zhu, Q., Zhang, Y., Gjessing, S., and Baar, T. (2016). Demand response management in the smart grid in a large population regime. *IEEE Transactions on Smart Grid*, 7(1):189-199.

Pecan Street (2021). Pecan Street Dataport. https://www.pecanstreet.org/dataport/.

Pediaditis, P., Papadaskalopoulos, D., Papavasiliou, A., and Hatziargyriou, N. (2021). Bilevel optimization model for the design of distribution use-of-system tariffs. *IEEE Access*, 9:132928-132939.

Potter, A., Haider, R., Ferro, G., Robba, M., and Annaswamy, A. M. (2023). A reactive power market for the future grid. *Advances in Applied Energy*, 9:100114.

Raghunathan, A. and Biegler, L. (2003). Mathematical programs with equilibrium constraints (mpecs) in process engineering. *Comput. Chem. Eng.*, 27:1381-1392.

Rockafellar, R. T. (1970). *Convex analysis*. Princeton Mathematical Series. Princeton University Press, Princeton, N. J.

Rockafellar, R. T. (2023). Generalized nash equilibrium from a variational analysis perspective. *Set-Valued and Variational Analysis: Theory and Applications*.

Rockafellar, R. T. and Wets, R. J. B. (1998). *Variational Analysis*. Springer Berlin Heidelberg.

Schittekatte, T., Momber, I., and Meeus, L. (2018). Future-proof tariff design: Recovering sunk grid costs in a world where consumers are pushing back. *Energy Economics*, 70:484-498.

Tushar, W., Zhang, J. A., Smith, D. B., Poor, H. V., and Thiébaux, S. (2014). Prioritizing consumers in smart grid: A game theoretic approach. *IEEE Transactions on Smart Grid*, 5(3):1429-1438.

Wächter, A. and Biegler, L. T. (2006). On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. *Mathematical Programming*, 106(1):25-57.

Wolgast, T., Ferenz, S., and NieSSe, A. (2022). Reactive power markets: A review. *IEEE Access*, 10:28397-28410.

Worku, S. B., Tsegaw, B. B., and Kassa, S. M. (2023). Existence and computations of best affine strategies for multilevel reverse stackelberg games. *Mathematical Methods of Operations Research*, 97(3):339-366.

Ye, J. (2005). Necessary and sufficient optimality conditions for mathematical programs with equilibrium constraints. *Journal of Mathematical Analysis and Applications*, 307:350-369.

Zheng, Y.-P. and Basar, T. (1982). Existence and derivation of optimal affine incentive schemes for stackelberg games with partial information: a geometric approach. *International Journal of Control*, 35(6):997-1011.

## A. The KKT reformulation of the lower-level

In this section, we construct the KKT reformulation of the lower-level GNEP (11). It is well known that if the optimization problem of each player satisfies standard constraint qualifications, solving the GNEP is equivalent to solving the concatenation of the KKT systems of all players (Facchinei and Kanzow, 2007). We start by constructing an inner polyhedral approximation the feasible dispatch set of each DG (1), allowing us to represent (9) as a polyhedron. In the case $\Omega(\mathbf{u}_{-i})$ is a polyhedron, no matter the choice of $\mathbf{u}_{-i}$, and the objective function of each prosumer is convex quadratic as long as $\Phi_i^k \ge 0$, the problem of each prosumer satisfies the so-called refined Slater constraint qualification (RSCQ), presented as Theorem 28.2 in Rockafellar (1970). As a consequence, the stacked KKT conditions of all followers become necessary and sufficient for a solution $\bar{\mathbf{u}}$ to be a solution to (11), i.e., $\bar{\mathbf{u}} \in S(T)$. Moreover, because the KKT conditions define a closed set, an optimal solution to (34) is guaranteed to exist as long as the feasible set is non-empty.

First, we construct a polyhedral inner approximation of (1c) following. Therefore, replacing (1c) by the following set of inequalities:

$$ -\alpha_{G,ij} p_{G,i}^k - \beta_{G,ij} \le q_{G,i}^k, \quad j = 0, \dots, n_{\text{SOC}} \qquad (46a) $$

$$ \alpha_{G,ij} p_{G,i}^k + \beta_{G,ij} \ge q_{G,i}^k, \quad j = 0, \dots, n_{\text{SOC}} \qquad (46b) $$

Where

$$ \beta_{G,ij} = q_{G,i,j}^0 - \alpha_{G,ij} p_{G,ij}^0 \qquad (47) $$

$$ \alpha_{G,ij} = \frac{q_{G,i,j+1}^0 - q_{G,ij}^0}{p_{G,i,j+1}^0 - p_{G,ij}^0} \qquad (48) $$
$$ (p_{G,ij}^0, q_{G,ij}^0) = \bar{s}_{G,i} \left( \sin\left(\frac{j\pi}{4n_{\text{SOC}}}\right), \cos\left(\frac{j\pi}{4n_{\text{SOC}}}\right) \right) \quad (49) $$

By doing so, problem (9) can be explicitly written as a quadratic optimization problem with polyhedral feasibility set as

$$ \min_{\mathbf{u}_i} J_{F,i} (T_i(\cdot), \mathbf{u}_i, \mathbf{u}_{-i}) \quad (50a) $$

$$ \text{s.t.} \quad p_i^k = \sum_{j \in \mathcal{N}_i^0} (p_{E,ij}^k - p_{E,ji}^k), \quad \forall k \in \mathcal{K} \qquad (\mu_i^k) \quad (50b) $$

$$ p_i^k = p_{G,i}^k - p_{B,i}^k - p_{D,i}^k, \quad \forall k \in \mathcal{K} \qquad (v_{P,i}^k) \quad (50c) $$

$$ q_i^k = q_{G,i}^k - q_{D,i}^k, \quad \forall k \in \mathcal{K} \qquad (v_{Q,i}^k) \quad (50d) $$

$$ s_{B,i}^k = s_{B,i}^{k-1} + \eta_{B,i} p_{B,i}^{+k} - \frac{p_{B,i}^{-,k}}{\eta_{B,i}}, \quad \forall k \in \mathcal{K} \qquad (v_{B,i}^k) \quad (50e) $$

$$ -\alpha_{G,ij} p_{G,i}^k - \beta_{G,ij} \le q_{G,i}^k, \quad \forall j = 1, \dots, n_{\text{SOC}}, \forall k \in \mathcal{K} \qquad (\lambda_{\text{GL},ij}^k) \quad (50f) $$

$$ \alpha_{G,ij} p_{G,i}^k + \beta_{G,ij} \ge q_{G,i}^k, \quad \forall j = 1, \dots, n_{\text{SOC}}, \forall k \in \mathcal{K} \qquad (\lambda_{\text{GU},ij}^k) \quad (50g) $$

$$ \tilde{p}_{G,i}^k \le p_{G,i}^k \le \hat{p}_{G,i}^k, \quad \forall k \in \mathcal{K} \qquad (\lambda_{\text{PGL},i}^k, \lambda_{\text{PGU},i}^k) \quad (50h) $$

$$ 0 \le s_{B,i}^k \le \hat{s}_{B,i}, \quad \forall k \in \mathcal{K} \qquad (\lambda_{\text{SL},i}^k, \lambda_{\text{SU},i}^k) \quad (50i) $$

$$ 0 \le p_{B,i}^{+,k} \le \hat{p}_{B,i}, \quad \forall k \in \mathcal{K} \qquad (\lambda_{\text{BL},i}^{+,k}, \lambda_{\text{BU},i}^{+,k}) \quad (50j) $$

$$ 0 \le p_{B,i}^{-,k} \le \hat{p}_{B,i}, \quad \forall k \in \mathcal{K} \qquad (\lambda_{\text{BL},i}^{-,k}, \lambda_{\text{BU},i}^{-,k}) \quad (50k) $$

$$ 0 \le p_{E,ij}^k, \quad \forall (j, k) \in \mathcal{N}_i^0 \times \mathcal{K} \qquad (\lambda_{\text{EL},ij}^k) \quad (50l) $$

The KKT conditions for the lower-level GNEP are comprised of the system of equations (50b) to (50e), for all $i \in \mathcal{N}$, along with (51), and (52), as presented below.

Stationarity conditions:

$$ \Phi_{pp,i}^k p_i^k + \Phi_{pq,i}^k q_i^k + \phi_{p,i}^k + \mu_i^k + v_{P,i}^k = 0 \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (51a) $$

$$ \Phi_{qq,i}^k q_i^k + \Phi_{pq,i}^k p_i^k + \phi_{q,i}^k + v_{Q,i}^k = 0 \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (51b) $$

$$ v_{P,i}^k + \sum_{j \in \mathcal{W}} \alpha_{G,ij} (\lambda_{\text{GL},ij}^k + \lambda_{\text{GU},ij}^k) + \lambda_{\text{PGL},i}^k - \lambda_{\text{PGU},i}^k = 0 \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (51c) $$

$$ v_{Q,i}^k - \sum_{j \in \mathcal{W}} (\lambda_{\text{GU},ij}^k - \lambda_{\text{GL},ij}^k) = 0 \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (51d) $$

$$ v_{P,i}^k + \eta_{B,i} v_{B,i}^k + \lambda_{\text{PBU},i}^{+,k} - \lambda_{\text{PBL},i}^{+,k} = 0 \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (51e) $$

$$ v_{P,i}^k - \frac{1}{\eta_{B,i}} v_{B,i}^k + \lambda_{\text{PBU},i}^{-,k} - \lambda_{\text{PBL},i}^{-,k} = 0 \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (51f) $$

$$ v_{B,i}^k - v_{B,i}^{k+1} + \lambda_{\text{SBU},i}^k - \lambda_{\text{SBL},i}^k = 0 \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K}' \quad (51g) $$

$$ c_{E,ij}^k + \mu_i^k + \lambda_{\text{EL},ij}^k = 0 \qquad \forall (i, j, k) \in \mathcal{N} \times \mathcal{N}_i^0 \times \mathcal{K} \quad (51h) $$

$$ c_{E,0i}^k + \mu_i^k - \lambda_{\text{EL},0i}^k = 0 \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \quad (51i) $$

Complementarity slackness conditions:

$$ \lambda_{\text{GL},ij}^k \ge 0 \perp -\alpha_{G,ij} p_{G,i}^k - \beta_{G,ij} \le q_{G,i}^k \qquad \forall (i, j, k) \in \mathcal{N} \times \mathcal{W}_i \times \mathcal{K} \quad (52a) $$
$$ \lambda_{GU,ij}^{k} \geq 0 \perp \alpha_{G,ij} p_{G,i}^{k} + \beta_{G,ij} \geq q_{G,i}^{k} \qquad \forall (i, j, k) \in \mathcal{N} \times \mathcal{W}_{i} \times \mathcal{K} \qquad (52b) $$

$$ \lambda_{SL,i}^{k} \geq 0 \perp s_{B,i}^{k} \geq 0 \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (52c) $$

$$ \lambda_{SU,i}^{k} \geq 0 \perp s_{B,i}^{k} \leq \hat{s}_{B,i} \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (52d) $$

$$ \lambda_{PBL,i}^{+k} \geq 0 \perp p_{B,i}^{+k} \geq 0 \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (52e) $$

$$ \lambda_{PBL,i}^{-k} \geq 0 \perp p_{B,i}^{-k} \geq 0 \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (52f) $$

$$ \lambda_{PBU,i}^{+k} \geq 0 \perp p_{B,i}^{+k} \leq \hat{p}_{B,i} \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (52g) $$

$$ \lambda_{PBU,i}^{-k} \geq 0 \perp p_{B,i}^{-k} \leq \hat{p}_{B,i} \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (52h) $$

$$ \lambda_{EL,ij}^{k} \geq 0 \perp p_{E,ij}^{k} \geq 0 \qquad \forall (i, j, k) \in i \in \mathcal{N}_{i}^{0} \times \mathcal{N}_{i}^{0} \times \mathcal{K} \qquad (52i) $$

Similarly, we can leverage the KKT conditions for the welfare equilibria; thereby, obtaining necessary and sufficient conditions for (12). The KKT conditions for the welfare equilibria are the system of equations: (9b) to (9e), (51a) to (51g), (51i), (52), and

$$ - \mu_{i}^{k} + \mu_{j}^{k} - \lambda_{EL,ij}^{k} = 0 \qquad \forall (i, j, k) \in \mathcal{N} \times \mathcal{N}_{i} \times \mathcal{K} \qquad (53a) $$

$$ \mu_{i}^{k} + c_{E,i0}^{k} + \lambda_{EL,i0}^{k} = 0 \qquad \forall (i, k) \in \mathcal{N} \times \mathcal{K} \qquad (53b) $$

In other words, the KKT system of the welfare optimization problem is the same as of the GNEs, except that (51h) have to be replaced by (53).

**Proposition 4.** The following is a necessary condition for a solution to (1) to be a solution to (2).

$$ \lambda_{EL,ij}^{k} + \lambda_{EL,ij}^{k} = 0 \qquad \forall (i, j, k) \in \mathcal{N} \times \mathcal{N}_{i} \times \mathcal{K} \qquad (54) $$

**PROOF.** The financial network is a symmetric directed graph, therefore if (i, j) belongs to the edges, so do (j, i). Then, from (53a) we obtain

$$ - \mu_{j}^{k} + \mu_{i}^{k} - \lambda_{EL,ji}^{k} = 0 \qquad \forall (i, j, k) \in \mathcal{N} \times \mathcal{N}_{i} \times \mathcal{K} \qquad (55) $$

Summing (53a) and (55), we obtain (54). □

**Proposition 5.** The following are sufficient conditions for a GNE to be a Welfare equilibrium

$$ \mu_{i}^{k} = \mu_{j}^{k} \qquad (56a) $$

$$ \lambda_{EL,ij}^{k} = \lambda_{EL,ji}^{k} = 0 \qquad (56b) $$

**PROOF.** Reorder the terms in (56a) to obtain

$$ \mu_{i}^{k} - \mu_{j}^{k} = 0 $$

Add a zero-valued $\lambda_{EL,ij}^k$ on the left hand side to obtain

$$ \mu_{i}^{k} - \mu_{j}^{k} + \lambda_{EL,ij}^{k} = 0 $$

Thereby, we have arrived to (53a), hence the GNE is a welfare equilibrium. □

**Proposition 6.** The following are sufficient conditions for a GNE to be a Welfare equilibrium

$$ c_{E,ij}^{k} = c_{E,ji}^{k} \qquad (57a) $$

$$ \lambda_{EL,ij}^{k} = \lambda_{EL,ji}^{k} = 0 \qquad (57b) $$

**PROOF.** From (51h) and (57a) we have

$$ \mu_{i}^{k} + \lambda_{EL,ij}^{k} = \mu_{j}^{k} + \lambda_{EL,ji}^{k} $$

Therefore

$$ \mu_{i}^{k} = \mu_{j}^{k} $$

and by Proposition 5, the GNE is also a welfare equilibrium. □
In summary, we have shown the relationship between the KKT conditions for welfare equilibria and those for GNE. The necessary and sufficient conditions outlined in the propositions highlight how a GNE can also satisfy the conditions for welfare maximization under specific circumstances. By leveraging the symmetry properties of the financial network and the interactions between prosumers, we establish that a GNE can be a welfare equilibrium when certain dual variable conditions hold.