1
GradMAP: Gradient-Based Multi-Agent Proximal Learning for Grid-Edge Flexibility
arXiv:2604.24549v1 [cs.LG] 27 Apr 2026
Yihong Zhou*, Hongtai Zeng, and Thomas Morstyn Department of Engineering Science, University of Oxford *Correspondence: [email protected]
Abstract—Coordinating large populations of grid-edge devices requires learning methods that remain fully decentralised in deployment while still respecting three-phase AC distributionnetwork physics. This paper proposes gradient-based multiagent proximal learning (GradMAP) to address this challenge. GradMAP trains independent neural-network policies for each agent without any parameter sharing, and each agent uses only its own local observation for online decision-making without communication. During offline training, GradMAP embeds a differentiable three-phase AC power-flow model in a primal–dual learning loop and uses implicit differentiation to propagate exact network-constraint violations to update the policy parameters. To speed up training, GradMAP reuses expensive environment gradients through a proximal surrogate within a trust region defined in the more direct policy-output (action) space, instead of the probability distribution space used in other works, such as PPO. In case studies with 1,000 agents managing batteries, heat pumps, and controllable generators on the IEEE 123-bus feeder, GradMAP learns decentralised policies that minimise three-phase AC load-flow constraint violations within 15 minutes of training on a single workstation-class NVIDIA RTX PRO 5000 Blackwell 48GB GPU. This is a 3–5× training speed-up over gradientbased self-supervised learning benchmarks and substantially better training efficiency than multi-agent reinforcement-learning benchmarks. In out-of-sample tests, GradMAP also delivers among the lowest operating cost and constraint violations. Index Terms—Multi-agent learning, implicit differentiation, primal-dual learning, three-phase unbalanced AC power flow, grid-edge flexibility.
I. I NTRODUCTION Decarbonisation and the electrification of heating and transportation are driving rapid growth in the number of grid-edge resources, including distributed renewables, battery storage, electric vehicles, and heat pumps. These devices can provide valuable flexibility through controllable power consumption or injection. Collectively, they can offer flexibility on a scale comparable to multiple large conventional power plants, with their global value estimated at USD $270 billion by 2040 [1]. However, traditional centralised dispatch approaches become increasingly impractical as the number of devices grows because of both computational limits (large-scale optimisation under uncertainty) and communication constraints (latency and bandwidth across the grid edge). To address these challenges, existing work has explored two broad classes of methods: traditional model-based coordination methods and, more recently, learning-based methods. Among traditional approaches, distributed optimisation solves the system-wide coordination problem through iterative message passing [2], [3]. Another traditional approach is hierarchical aggregation–disaggregation [4], [5], in which
a system operator can perform smaller-scale dispatch over aggregated feasible power regions submitted by aggregators, who then map aggregate instructions back to individual gridedge devices. However, distributed optimisation often suffers from significant communication overheads, slow convergence, and numerical instability [6]. Meanwhile, aggregation-based methods rely on accurate characterisations of the cost and feasible region of the aggregated flexibility resources, which is challenging due to uncertainty, nonlinear device dynamics, and nonlinear alternating-current (AC) distribution-network constraints [4]. Another promising direction is learning-based decisionmaking, including imitation learning (IL), reinforcement learning (RL), and self-supervised learning. In IL, an agent learns a policy by mimicking a set of expert decision trajectories, which may be obtained from historical operational data or generated by solving offline optimisation problems. This approach has been shown to be effective for power-system dispatch [7], [8]. However, it requires a significant amount of expert data, which is unlikely to be available for new systems and systems undergoing rapid change. Although expert training data can be generated through offline simulations, this is costly when the underlying decision problem is difficult to solve, such as for large-scale AC optimal power flow (OPF). In RL, agents learn through interaction with an environment that can largely be treated as a black box, as long as observations and rewards are returned in response to actions. This broad applicability makes RL appealing for many decisionmaking problems [9]–[11]. Multi-agent (MA) RL (MARL) is a particularly relevant extension for coordinating large numbers of grid-edge devices, because a single-agent RL would require an impractical communication capability that monitors millions of devices and disseminates online dispatch signals to them. Centralised single agent RL would also raise issues related to the autonomy and and privacy of individual device owners. The MA setting addresses these issues by assigning each device (or group) its own agent, which can make online decisions without communication while learning coordination through offline training. Recent studies have explored MARL for grid-edge device coordination [11]–[13]. However, these works follow the standard black-box RL setting, treating the power-system environment as a black box and neglecting the fact that power-system operators have been using wellestablished grid physics and mathematical models for decades [14]–[16]. Ignoring available network modelling methods is unlikely to be information-efficient. Some RL studies have also recognised the benefits of envi-
2
ronment models, especially when they are differentiable. Refs. [17]–[19] showed that exploiting reward gradients w.r.t. actions can improve the efficiency of RL learning and lead to better constraint satisfaction. Here a challenge is that evaluating such environment gradients can be computationally expensive. For example, AC power-flow gradients w.r.t. power injections typically need to be obtained through implicit differentiation. To alleviate this cost, [20] and [21] exploited gradient reuse, where the environment gradient is used to construct improved actions for multi-step supervised policy learning. [22] pointed out their methods may lead to misaligned learning objective that reduces efficiency, and proposed to reuse gradients with a probability trust-region clipping and a KL penalty that leads to improved performance, which is similar to the well-known proximal policy optimisation (PPO) approach in standard black-box RL [23]. However, it remains questionable whether KL- and clipping-based trust-region regularisation defined in the probability distribution space is still the best choice for the differentiable case, because a direct trust region in policyoutput space is more closely aligned with the coordinates in which the gradients are defined. Moreover, these works [20]– [22] were only focused on the single-agent setting. Self-supervised learning trains models on unlabeled data by generating its own supervisory signals. In power-system decision-making, existing works [24]–[27] start from standard power-system optimisation that finds optimal power dispatch, and then extend it to a more nonlinear non-convex problem that finds the optimal parameterised decision policies generating the dispatch, thereby removing the MDP assumption in RL. In this context, (environment) gradients are the commonly exploited self-supervisory signals. However, existing methods require gradient evaluation at every policy update, leading to the same challenge of computationally expensive gradient evaluation previously noted for differentiable RL [17]–[19]. Ref. [28] tried to train a surrogate environment model to accelerate gradient evaluation, but such a data-driven powergrid surrogate may introduce extra bias. Also, existing selfsupervised learning literature [24]–[28] failed to consider the multi-agent setting. As the above discussion shows, several key gaps remain in learning-based decision-making. For IL [7], [8], the main limitation is the difficulty of obtaining expert decision-making data for training. For RL, a major gap is that well-established power-system models and their differentiability are often unexploited, leading to less efficient learning [10]–[13]. Although some RL studies have explored gradient information and the gradient reuse [20]–[22], the trust-region construction is not well aligned with the policy-output space (where the gradients are computed) and existing works have focused only on the single-agent setting. Self-supervised learning exploits powersystem environment gradients and avoids the MDP assumption in RL [24]–[28], but existing studies failed to explore the more efficient gradient reuse and are also limited to the single-agent case. Finally, none of the existing gradient-based methods considered exact AC power-flow constraints, which are essential in power system operation. In this paper, we propose a gradient-based multi-agent proximal learning (GradMAP) framework for large-scale grid-
edge flexibility coordination that addresses all the identified gaps, with the following three main contributions: 1) We develop a gradient-informed, MDP-free multi-agent (MA) learning framework that scales to 1,000 agents with independent neural-network parameters and purely local observations, while remaining trainable on a single workstation-class GPU within only 15 minutes. 2) We introduce a proximal update that mitigates the computational overhead of gradient evaluation by enabling repeated reuse of the same gradients (from the environment to agent policy outputs). We apply direct trust-region regularisation in policy-output space rather than through distribution-space KL and policyprobability clipping used in PPO-style works [22], [23]. We numerically demonstrate that our method performs better in both operation cost, AC network constraints, and training time than the KL-based variant [22], exactgradient-based self-supervised learning approaches related to [17]–[19], [24]–[27], and standard black-box MARL baselines including IPPO and MAPPO [29]. 3) We embed an exact differentiable three-phase unbalanced AC power-flow solver in the training loop and compute voltage gradient w.r.t. nodal power injections via implicit differentiation, allowing AC networkfeasibility to directly shape the agent policy learning. The paper is organised as follows. Section II presents the problem formulation. Section III describes our proposed “GradMA” framework, which extends single-agent selfsupervised learning to the MA setting that re-evaluates gradients for every update. Section IV describes the full GradMAP framework, which enhances GradMA with proximal updates that reuse gradients for stronger efficiency. Section V reports the experimental results and Section VI concludes the paper. II. P ROBLEM F ORMULATION We consider a population of homes connected to a threephase unbalanced distribution feeder. Each home is represented by one agent, and each agent controls one flexible device, such as a battery, a heat pump, or a controllable generator, together with exogenous uncontrollable rooftop PV and local demand. This one-agent-per-home, one-deviceper-home assumption is representative of many residential flexibility settings and is adopted for ease of presentation. The proposed framework could handle homes with multiple controllable assets either by augmenting the local state-action space or by introducing multiple coordinated agents within the same home. The proposed framework could also be straightforwardly extended beyond residential flexibility to a wider set of grid-edge coordination problems (e.g. EV charging networks, commercial battery storage systems). The objective is to learn decentralised policies that minimise the operation cost of all homes while satisfying both local device constraints and feeder-level physical limits. The potential user of this learning framework may be an aggregator or the distribution system operator (DSO). We assume that this user has access to the network model, and this is possible since this information sharing supports network constraint
3
handling. Note that our focus is on the communication-free decision-making framework rather than fair profit-sharing or market design. Although we use the term “flexibility”, demand baselining methodologies are not the focus here, which mainly affect market settlement. With this in mind, we first introduce the device models and then the three-phase network model that couples them. A. Grid-Edge Device Models We consider three representative classes of grid-edge devices that capture important cross-time couplings encountered in grid-edge flexibility: batteries with nonlinear efficiency and energy coupling, heat pumps with thermal coupling, and controllable distributed generators with ramp-rate coupling. Given our home-level abstraction, we next describe the dynamics, constraints, and cost function of each controllable device type. 1) Battery: For a battery managed by agent i, its battery dynamics, constraints, and operating cost are summarised as: bat bat Êi,t+1 = Ei,t + ηi (P̂i,t ) P̂i,t ∆t, (1a) Ei,t+1 = clip(Êi,t+1 , 0, Eimax ), (1b) bat P bat = (Ei,t+1 − Ei,t )/(ηi (P̂i,t ) ∆t), (1c) i,t bat bat bat (1d) −Pi,max ≤ P̂i,t ≤ Pi,max , max 0 ≤ Êi,t+1 ≤ Ei , (1e) bat Ei,T ≥ Ei,target , (1f) deg energy bat (1g) Ci,t = Ci,t + Ci,t , C deg = cdeg |P bat |∆t, (1h) i,t i,t i energy imp net + exp net + Ci,t = ρt [pi,t ] − ρt [−pi,t ] , (1i) pnet = P load − P pv + P bat , (1j) i,t i,t i,t i,t bat where Ei,t denotes the stored energy, and P̂i,t denotes the power command issued by the agent (control variable). t is the time step within the control horizon T := {1, . . . , T }. For interpretability, we define an auxiliary variable Êi,t+1 to represent the implied next energy state if the battery follows exactly bat the power command P̂i,t for the whole time step ∆t without any physical limits. However, in the real world, the battery stops charging or discharging once the stored energy reaches its bounds, which is captured by the clipping operation in (1b). bat Therefore, the actual AC-side power Pi,t , averaged over ∆t, should be back-calculated from the true energy change and the bat efficiency ηi (P̂i,t ), as shown in (1c). The power command bat P̂i,t is constrained by (1d), and the implied next energy state Êi,t+1 is constrained by (1e) (although mathematically redundant, its violation provides valuable learning signals as will be discussed in Section III-A). We also introduce an endstate energy requirement in (1f), which ensures the battery retains sufficient stored energy for the next episode. Instead of making the common assumption of constant charging and discharging efficiency used in much of the battery optimisation literature (which provides computational tractability but can introduce a significant “sim-to-real” gap), we use a set of more accurate nonlinear efficiency curves obtained by exponential fitting to measurement data [30], which captures the reduced efficiency typically observed at low power levels. Our case studies will randomly assign a curve to each battery system.
Eqs. (1g)–(1j) separate throughput-based battery degradation from the electricity cost (with import/export tariffs ρimp t and ρexp t ) at the household meter. Note that in full detail, degradation cost is also nonlinear and consists of both calendar aging and cycle aging [31]. For simplicity we do not include these features, but our inclusion of nonlinear efficiency shows the framework could straightforwardly be extended to also incorporate nonlinear degradation cost. 2)Heat Pump: The model for a heat pump is: ∆t ϑout hp t − ϑi,t , + COPi P̂i,t ϑ̂i,t+1 = ϑi,t + Ci Ri set set ϑi,t+1 = clip ϑ̂i,t+1 , ϑi − δi , ϑiout + δi , ϑt −ϑi,t Ci hp ∆t (ϑi,t+1 − ϑi,t ) − Ri , P = i,t COPi hp hp 0 ≤ P̂i,t ≤ Pi,max , set ϑset − δ i ≤ ϑ̂i,t+1 ≤ ϑi + δi , i ϑi,T ≥ ϑhp i,target , hp energy use Ci,t = Ci,t + Ci,t , C use = cdeg |P hp |∆t, i,t i i,t energy imp net + exp net + C = ρ [p t i,t ] − ρt [−pi,t ] , i,t pnet = P load − P pv + P hp , i,t i,t i,t i,t
(2a) (2b) (2c) (2d) (2e) (2f) (2g) (2h) (2i) (2j)
hp where ϑi,t denotes the indoor temperature, and P̂i,t denotes the electrical power command issued by the agent (control variable). Similar to the battery model, here we have defined an auxiliary variable ϑ̂i,t+1 to represent the implied next indoor temperature if the heat pump operates fully according to hp the power command P̂i,t . Here we assume that the indoor temperature will also stop at the comfort boundaries through a feedback controller, which is captured by the clipping hp operation in (2b). Therefore, the actual AC-side power Pi,t (averaged over ∆t) should be back-calculated from the true temperature change and the RC thermal model [32], as shown hp in (2c). The power command P̂i,t is constrained by (2d), and the implied next indoor temperature ϑ̂i,t+1 is constrained by (2e). We also introduce an end-state temperature requirement by (2f), which ensures thermal flexibility for the next episode. Eqs. (2g)–(2j) models the cost for the agent that includes the simplified wear-and-tear cost (2h) and the electricity cost (2i).
3) Controllable Generator: The model for a controllable generator is summarised as: gen gen gen (3a) ∆P̂i,t = P̂i,t − Pi,t−1 , gen gen (3b) ∆Pi,t = clip ∆P̂i,t , ∆P i , ∆P i , P gen = P gen + ∆P gen , (3c) i,t−1 i,t ( geni,t gen gen Pi,min ≤ P̂i,t ≤ Pi,max , (3d) gen ∆P i ≤ ∆P̂i,t ≤ ∆P i , (3e) gen energy fuel (3f) Ci,t = Ci,t + Ci,t , C fuel = (a P gen + b (P gen )2 )∆t, (3g) i i,t i i,t i,t energy imp net + exp net + Ci,t = ρt [pi,t ] − ρt [−pi,t ] , (3h) pnet = P load − P pv − P gen . (3i) i,t i,t i,t i,t Following the same logic as the battery and heat pump models, gen P̂i,t denotes the generation setpoint command issued by the
4
gen is the auxiliary implied power agent (control variable), ∆P̂i,t gen change, and Pi,t denotes the actual generated power after respecting the ramping constraints, achieved through the ramp clipping in (3b) based on the lower and upper ramping bounds gen ∆P i , ∆P i . The actual generated power Pi,t is obtained by gen is constrained by (3d), and (3c). The setpoint command P̂i,t the implied power change is constrained by (3e). Eqs. (3f)– (3i) model the total agent operation cost as the sum of the quadratic fuel cost and the electricity cost. Due to the clipping, the per-time-step state constraints are satisfied naturally, so the critical constraints that each agent must learn to satisfy are the end-state requirements in (1f) and (2f), as well as the network constraints introduced next.
B. Three-Phase Network Model The household net loads introduced above (device power plus uncertain PV and local demand) are connected to different phases and different nodes of a three-phase unbalanced distribution network. Here we assume one household only connects to a single phase. Following standard formulations for unbalanced distribution feeders, we consider a network with one slack bus and NB three-phase PQ buses, and define: Y 3N • st ∈ C B : Net complex power injections from Wyeconnected components at PQ buses. 3N D • st ∈ C B : Net complex power injections from Deltaconnected components at PQ buses. 3 0 • vt ∈ C : Fixed complex voltage at the slack bus. 3NB • vt ∈ C : Complex nodal voltages at PQ buses. 3N • it ∈ C B : Net phase current injections at PQ buses. 3N D • it ∈ C B : Phase-to-phase currents for PQ buses. The electrical characteristics is represented by the bus B +1)×3(NB +1) admittance matrix Y ∈ C3(N , partitioned as Y00 Y0L Y := , (4) YL0 YLL where the sub-matrices denote self and mutual admittances between the slack bus (0) and non-slack buses (L). The AC power flow equations governing steady-state operation are Y (5a) diag(H ⊤ iD t )vt + st = diag(vt )it , D st = diag(Hvt )iD , (5b) t it = YL0 vt0 + YLL vt , (5c) vmin ≤ |vt | ≤ vmax (5d) where (·) indicates complex conjugation, and H is the transformation matrix relating phase-to-neutral and phase-to-phase quantities (see [33] for details). Eq. (5d) specifies the network voltage magnitude limit, where | · | performs the elementwise norm operation. Real-world distribution networks are typically voltage-constrained and the line current limits can be disregarded [3], especially for distribution networks with long power lines [34]. Incorporating the line current constraints is also straightforward. C. Multi-Agent Learning Formulation Let Nbat , Nhp , and Ngen denote the battery, heat-pump, and generator agent sets, respectively. Let oi,t be the local dev observation of agent i, and let P̂i,t = πi (oi,t ) denote its control action (power command of each of the devices such gen as P̂i,t ). We can formulate the multi-agent learning problem
through an optimisation lens, which helps us get rid of the Markov restriction # " Tin standard RL: X gen X X X hp bat Ci,t E Ci,t + min Ci,t + θ:={θi }N i=1
t=1
i∈Nbat
i∈Nhp
i∈Ngen
(6a) s. t. ∀i, t (6b) battery agents satisfy (1), ∀i ∈ Nbat (6c) heat-pump agents satisfy (2), ∀i ∈ Nhp (6d) generator agents satisfy (3), ∀i ∈ Ngen (6e) with network constraints (5a)–(5d), ∀t. (6f) Problem (6) is very similar to a standard dispatch problem, except that the decision variables are not the power dispatch dev actions P̂i,t , but the parameterised policies {πi (·, θi )}N i=1 (such as neural networks (NNs) with parameters θi ) that generate them from local observations without online communication. The direct cost summation across all agents in the objective (6a) is appropriate in the present cooperative setting because the goal is to minimise social operating cost, or equivalently maximise social welfare under exogenous prices, and we focus on the fundamental MA dispatch problem rather than fair allocation. The key challenge is not only the network coupling in (5a)–(5d), but also the requirement that good dispatch actions must be generated from local information alone. This is the key difference between our MA framework and standard self-supervised learning approaches [24]–[27]. Note that for better learning performance, our agent policies will not directly output device-level power. Instead, let ai,t denote the policy output. It is clipped to [−1, 1] and decoded into device-level power by ai,t + 1 hp hp bat bat P̂i,t = ai,t Pi,max , P̂i,t = Pi,max , 2 ai,t + 1 gen gen gen gen P̂i,t = Pi,min + Pi,max − Pi,min , (7) 2 which naturally enforce the device power limits. dev P̂i,t = πi (oi,t , θi ),
III. G RADIENT-BASED M ULTI -AGENT L EARNING As discussed in Section II-C, unlike single-agent approaches [24]–[27], our parameterised policies {πi (·, θi )}N i=1 are restricted to local observations for each agent in the online phase. Because of this restriction, we adopt the centralised training and decentralised execution (CTDE) framework [29], which trains the agent policies with global information but executing them locally, to extend the exact gradient-based self-supervised learning to the MA case (GradMA). Fig. 1 illustrates the architecture of GradMA. A. Exact Gradient-Based Multi-Agent Primal-Dual Learning Problem (6) is essentially a constrained nonlinear program. A common approach for handling this problem is primal-dual learning [24], where the hard constraints are first incorporated into the objective to form a Lagrangian loss function. Then gradients (or higher-order information) are used to perform loss minimisation, and Lagrangian multipliers are gradually increased to impose heavier penalties for constraint violations. The standard Lagrangian would attach multipliers to every constraint. For the problem scale considered here, this leads to a large dual space and unstable learning dynamics. We
5
therefore aggregate violations into six physically meaningful violation channels in: voltage magnitude VVolt , battery per-step energy VBstp , battery terminal energy VBend , heat-pump per-step temperature VHstp , heat-pump terminal temperature VHend , and generator ramping VGrmp : T X max([|vt | − vmax ]+ + [vmin − |vt |]+ ) VVolt := , 0.5 (vmax − vmin ) T /N t=1 VBstp :=
T X X ([Êi,t − E max ]+ + [−Êi,t ]+ ) i
bat Pi,max T
VHend :=
i
δi T
i∈Nhp
VGrmp :=
δi
costs gradients w.r.t. actions
I V. Proximal L ear ning (Gr adM AP)
I I I . Gr adM A
,
0.5 (∆P i − ∆P i ) T
Here, max(·) returns the maximum element of the voltage vector. The denominator in each channel is a normalisation factor that brings the corresponding violation term, when summed over all time steps in a rollout, to a scale comparable to those of the other terms. The network-constraint violation term is further scaled up by the number of agents so as to preserve its importance in the Lagrangian as the number of agents increases. Although the battery and heat-pump models in Section II-A already enforce per-step state-limit satisfaction through clipping, the violation channels based on the power commands in VBstp and VHstp are empirically found important for effective learning in both gradient-based methods and RL. The rollout Lagrangian h XT is then calculated as follows: tot L(θ, λ) = E Ct + λVolt VVolt + λBstp VBstp + t=1 i λBend VBend + λHstp VHstp + λHend VHend + λGrmp VGrmp (8) with X X hp X gen M bat ( Ci,t + Ci,t + Ci,t ), Cttot := 1 PN i=1 |Pi,max | i∈Nbat i∈Nhp i∈Ngen N (9) where the denominator also scales the total cost to make it unitless, so that it is comparable to the constraint-violation terms. M = 200 is a tunable constant that adjusts the importance of the total-cost term in the Lagrangian. With the Lagrangian defined, the primal-dual learning will consist an outer dual loop with Kdual steps and an inner primal loop with Kprimal steps. In the inner loop, the policy parameters are updated by gradient descent (or other optimiser such as Adam) on the Lagrangian (8). After the inner primal updates, the dual variables are each updated by projected ascent on the batch-mean (˜·) violation channels at a dual learning rate α: λr ← [λr + α Ṽr ]+ , r ∈ {Volt, Bstp, Bend, Hstp, Hend, Grmp}, (10) The algorithm stops when the primal-dual gap is sufficiently small or after a maximum number of iterations. Each of the initial dual values λr can be set to zero.
global network constraint violation
II.A and II.B. Differentiable Environment
gradients w.r.t. NN params
Fully Decentr alised Execution actions 1
actions 2
actions 3
... actions N
...
update all agent neural networks (NNs)
N agents local obs 1
.
per-agent local device constraint violation
chain rule re-parametrisation trick
gradients w.r.t. NN params (change in each reuse)
,
gen gen + T X X [∆P̂i,t − ∆P i ]+ + [∆P i − ∆P̂i,t ] t=1 i∈Ngen
auto-differentiation per-agent
chain rule and re-param trick
i
+ X [ϑhp i,target − ϑi,T ]
I I I .A. L agr angian and Pr imal Dual L ear ning
Trust region on policy outputs
T X X [ϑ̂i,t − (ϑset + δi )]+ + [(ϑset − δi ) − ϑ̂i,t ]+ t=1 i∈Nhp
I I I .B. I mplicit Differentiation
reuse for Kprox times
,
t=1 i∈Nbat bat X [Ei,target − Ei,T ]+ VBend := , bat Ei,target i∈Nbat
VHstp :=
Centr alised Tr aining
local obs 2
local obs 3
local ... obs N
Fig. 1. Architecture diagram of the proposed GradMA and GradMAP framework.
B. Implicit Differentiation for AC Power Flow As illustrated in Fig. 1, a core part in the primal-dual training process is the evaluation of gradients of the Lagrangian (8). Most of the gradient terms can be obtained through automatic differentiation (AD) available in modern deep learning frameworks. The AC network constraint part is different: the voltage profile is obtained from an iterative threephase power-flow solve, so naively applying AD would unroll all solver iterations and backpropagate through the entire loop, which creates computational and memory issues. Instead, we use implicit differentiation on the converged fixed point. −1 Using the Z-Bus formulation [35], let Z = YLL and w = −ZYL0 vt0 . Then the power flow can be written as vt = Φ(vt , st ) := Ziinj (vt , st ) + w, (11) where st stacks the Wye and Delta load injections and iinj (·) is the corresponding current injection map implied by (5a)–(5b). At the converged solution vt⋆ = Φ(vt⋆ , st ), differentiating both sides w.r.t. the load injections gives ∂v ⋆ ∂Φ ⋆ ∂Φ ⋆ t I− (vt , st ) = (v , st ), (12) ∂vt ∂st ∂st t ⋆ ∂v The term ∂stt is the Jacobian matrix of the converged voltage w.r.t. the load injections. Explicitly computing this full Jacobian matrix is computationally expensive and unnecessary. To update the policies, the actual quantity we target to evaluate is the gradient of the voltage-violation scalar term w.r.t. the load ∂VVolt Volt injection vector at time step t, namely ∂V ∂st . Let gv := ∂vt⋆ , which can be obtained directly by AD. Then, by the chain rule and (12), we can re-write our target gradient as ∂v ⋆ ⊤ ∂VVolt t = gv ∂st ∂st ∂Φ ⊤ −⊤ ∂Φ ⋆ = (vt⋆ , st ) I− (vt , st ) gv . (13) ∂st ∂vt −⊤ ∂Φ If we define a vector γt := I − ∂v (vt⋆ , st ) gv , then we t can calculate γt by solving ⊤ ∂Φ ⋆ I− (vt , st ) γt = gv , (14) ∂vt
6
We solve this linear system using the Krylov iterative method BiCGSTAB, which avoids explicitly forming the coefficient matrix in Eq. (14) and can thus reduce memory usage compared with a direct solve (e.g., LU- or QR-based). Then our target gradient is calculated as ∂Φ ⊤ ∂VVolt = (vt⋆ , st ) γt . (15) ∂st ∂st Note that many Krylov solvers cannot directly handle the complex-valued system in (14). Thus we convert the complex variables into stacked real-imaginary vectors, e.g., vt := [Re(vt )⊤ , Im(vt )⊤ ]⊤ . C. Re-parameterisation Trick for Stochastic Policy Learning Standard self-supervised gradient methods are usually presented with deterministic policies [24]–[27]. In practice, however, deterministic policies can make the optimisation landscape brittle and can get trapped in poor local minima. To handle this issue, we introduce stochastic policies that enable exploration and smoothing, and optimise them through the reparameterisation trick [21]. More specifically, for each agent, the policy outputs the parameters of a Gaussian distribution over a normalised action, πθi (· | oi,t ) = N µθi (oi,t ), σθ2i (oi,t ) , (16) Given oi,t , the sampled action is ai,t = µθi (oi,t ) + σθi (oi,t ) ϵi,t , ϵi,t ∼ N (0, 1). (17) ai,t is then clipped to [−1, 1] and decoded into device-level power according to (7). If ∂L/∂ai,t denotes the action gradient returned by the differentiable environment, then the gradients in policy-output coordinates are ∂L ∂L ∂L ∂L = , = ϵi,t , (18) ∂µi,t ∂ai,t ∂σi,t ∂ai,t which can then be used for policy optimisation. IV. G RADIENT-BASED P ROXIMAL L EARNING In exact-gradient-based GradMA (and also existing selfsupervised learning approaches), after each differentiable environment rollout, fresh gradient is computed and updates the policy only once. However, differentiating the environment particularly with the implicit differentiation for the three-phase AC power flow can be dominant computational bottlenecks. This section thus introduces GradMAP, which first collects rollouts, computes the corresponding Lagrangian gradients w.r.t. agents’ policy outputs (rather than policy parameters), and then reuses these cached gradients for multiple inner updates within a trust region in the policy-output space, as illustrated in Fig. 1. Suppose for each primal step a rollout batch has been old generated by an old policy with outputs (µold i,t , σi,t ), and let µ σ (gi,t , gi,t ) be the corresponding environment gradients w.r.t. old the policy outputs (µold i,t , σi,t ) via (18). GradMAP tries to optimise thehfollowing problem in each primal step: i µ σ min E gi,t µi,t (θi ) + gi,t σi,t (θi ) − τ H(πθi (· | oi,t )) θ:={θi }
s.t. T (θ) ≤ ϵtr , (19) where τ H(πθ (· | oi,t )) is an entropy bonus to encourage exploration, which is commonly used in stochastic policy learning [29]. It is apparent that, without the trust region and the entropy bonus, differentiating (19) w.r.t. policy outputs (not
µ σ θ) gives exactly the gradients to be reused (gi,t , gi,t ). The trust region is defined as: 2 old 2 1/2 , (20) T (θ) := E ∥µi,t (θi )−µold i,t ∥2 +∥σi,t (θi )−σi,t ∥2 where E represents the average over all agents, all time steps, and all rollouts in the batch (this is feasible because we have assumed centralised offline training as part of CTDE). The trust region is defined directly in the policy-output coordinates (µ, σ) because the reused first-order information is also defined in these coordinates after re-parameterisation. This is both theoretically and intuitively preferable. Theoretically, the GradMAP surrogate is a first-order approximation in the same output coordinates, so the proximal term directly controls the validity region of that approximation. Intuitively, it limits how far the action mean and standard deviation are allowed to move, rather than indirectly constraining a probability ratio or KL quantity in the probability distribution space that is not the coordinate system in which the cached gradients are defined. In implementation, we do not solve the constrained problem directly. Instead, the trust-region constraint is converted into a quadratic soft h penalty, yielding the GradMAP surrogate loss: µ σ ℓGMAP = E gi,t µi,t (θi ) + gi,t σi,t (θi ) − τ H(πθi (· | oi,t )) i β 2 old 2 ∥ +∥σ (θ )−σ ∥ ) + (∥µi,t (θi )−µold i,t i i,t 2 i,t 2 , (21) 2 where β is an adaptively scaled trust-region penalty coefficient. We can then perform Kprox updates to minimise ℓGMAP . Note that in the Kprox updates, although the same environment gradients w.r.t. the agent policy outputs are re-used, the gradients w.r.t. the agent NN parameters are changing due to both the changing trust region penalty and the update in NN parameters. After each primal step we average T (θ) and then update β by dynamic multiplicative scaling (following the KL-penalty version ofPPO [23]): 1.1 β, T (θ) > ϵtr , β ∈ [50, 104 ]. (22) β ← β/1.1, T (θ) < ϵtr /2, β, otherwise, Algorithm 1 summarises the triple-loop structure of GradMAP: a outer dual loop for constraint penalty updates, a middle primal loop under fixed dual variables, and a fast innermost surrogate-optimisation loop. The innermost updates are computationally cheap because of the gradient reuse.
V. C ASE S TUDIES A. Experimental Setup To demonstrate the performance of GradMAP, we present a case study based on the IEEE 123-bus mid-voltage feeder [36] with 1,000 agents. For scaling comparisons we additionally present case studies with the IEEE 13-bus medium-voltage feeder with 10 agents and the reduced EULV 124-bus lowvoltage feeder with 100 agents [36]. Fig. 2 presents the IEEE 123-bus network together with the connection points of the 1,000 agents used in the main study. For training and testing, three-phase unbalanced AC power flow is solved via the Z-Bus method [35]. Note that when computing the implicit networkconstraint gradient by solving (14) with a Krylov iterative solver, a small number of NaN network gradients may arise at the start of training if the Krylov solve fails to converge. We can set these few NaNs to zero without impacting training,
7
Fig. 2. IEEE123 feeder topology used in the main case study together with the connection locations of the 1,000 agents. Gray nodes denote feeder buses, coloured branches indicate occupied phase connections, and marker types distinguish the 334 battery, 333 heat-pump, and 333 generator agents distributed across the network.
and this issue typically becomes less frequent as training progresses and the policy becomes better behaved. Preconditioning or the more memory-intensive direct solve (such as LU- or QR-based) can also be used to improve robustness. Real-world time-series datasets are used in this paper. The non-controllable household load data are obtained from the residential measurement dataset in [37], from which we extract 2,000 half-hourly household demand time series spanning 15 months after preprocessing. The hourly solar-generation and air-temperature time series are sourced from [38]. The import and export price data are obtained from the UK Octopus Energy Agile tariff [39], a residential tariff scheme that varies on an half-hourly and daily basis. Fig. 3 illustrates these data over seven consecutive days. Linear interpolation is used to increase the temporal resolution. In our case studies, each time step is 15 min, and each episode spans one day with a 24-h horizon, giving T = 96. We use 12 consecutive months of data for training and hyperparameter tuning, and the subsequent month, January (31
Price (£/kWh) Temp (°C) Solar (kW) Load (kW)
0.50
Algorithm 1 GradMAP training procedure 1: Initialise all agents’ policy parameters {θi }, dual variables λ, and trust coefficient β 2: for m = 1, . . . , Kdual do 3: for j = 1, . . . , Kprimal do 4: Interact with the environment and collect a batch of rollouts 5: Store rollout statistics (µold , σ old ) and differentiate the rollout Lagrangian (8) w.r.t. actions 6: Convert action gradients into cached (g µ , g σ ) using the re-parameterisation rule (18) 7: for e = 1, . . . , Kprox do 8: Update θi of all agents by minimising (21) using a deep learning optimiser (e.g. Adam) 9: end for 10: Update β according to (22) 11: end for 12: Update each dual multiplier by (10) 13: end for
Import price Export price
0.25 5
Home 0
Home 249
Home 499
Home 749
Home 999
0 10 0 10 0 Day 18 Day 19 Day 20 Day 21 Day 22 Day 23 Day 24
Fig. 3. Five-day visualisation of the datasets used in the IEEE123 1,000agent case study.
days), for out-of-sample testing. During testing, the days are treated as consecutive, such that the final state of one day is used as the initial state of the next day. All decentralised policies use the same 8-dimensional local observation at the current time step, defined as imp out loc ⊤ oi,t = [t̄, x̄i,t , d¯net , p̄exp (23) t , θ̄t , ūi,t , v̄i,t ] , i,t , p̄t where the entries are the normalised time t̄ := t/T , device state x̄i,t (energy level for batteries, temperature for heat pumps, and power at the previous time step for generators), net uncontrollable demand d¯net i,t (uncontrollable load minus solar and ρ̄exp generation), import and export prices ρ̄imp t , outdoor t out temperature θ̄t , terminal-state urgency ūi,t , and local voltage loc v̄i,t . All features are normalised to the interval [0, 1] using either the device parameters or the training data statistics. The terminal-state urgency feature ūi,t was found to be helpful for satisfying the terminal-state constraints. It is calculated as: hp bat 1{i∈Nbat } (Ei,target − Ei,t ) 1{i∈Nhp } (ϑi,target − ϑi,t ) ūi,t = + . (1 − t̄)(Eimax /2) (1 − t̄)δi (24) The reported performance metrics are the total operating cost (6a), the voltage-violation magnitude, the battery and heat-pump end-state violation magnitudes, the training wallclock time, and the number of primal steps to reach specific performance. We evaluate the proposed GradMAP against several benchmarks. To highlight the advantage of our proposed gradient reuse with a policy-output-space trust region, we include GradMA-RPO that adapts the single-agent PPOstyle probability-distribution-space trust region from [22]. We also compare against GradMA (Section III), our proposed MA extension of existing self-supervised learning methods [24]–[27]. Finally, we include standard black-box MARL algorithms, namely IPPO and MAPPO [29]. Following the MAPPO implementation from [29], all agents share a single global critic network. During offline training, this shared critic receives the full set of observations from all agents, and returns an agent-specific value estimate. Note that [29] also proposed using the global state together with each agent’s local observation rather than concatenating all observations, in order to maintain computational scalability. However, our empirical tests show that this performs worse. For this reason, our MAPPO can only be tested in the 10-agent and 100-agent cases, but is not computationally feasible in the 1,000-agent case. Another computationally efficient variant of MAPPO outputs the same value estimate for all agents [29], but we
8
Fig. 4. All 1,000 agents’ normalised power on the held-out test day 18. Positive values represent demand while negative values represent generation. Grey dashed lines represent the 0.95–1.05 p.u. voltage limits. In the left panel (without network constraint signals), the resulting phase voltage ranges (in p.u.) are 0.9680–1.0810, 1.0000–1.0489, and 0.9878–1.0862. In the right panel, the corresponding ranges are 0.9702–1.0581, 0.9950–1.0550, and 0.9865–1.0595.
Power State Demand / (kW) (kWh/°C)solar (kW)
Solar
Bat 1
Bat 2
5 0
Power State Demand / (kW) (kWh/°C)solar (kW)
10 0 10 0 10
10
HP 1
HP 2
Day 18
Day 18
Solar
Demand
Bat 1
Bat 2
5 0
Day 18
NA
NA
Day 18
Day 18
0
Day 18
Day 18
Without network constraint signals 10
5 0
Gen 2
10
HP 1
HP 2
Gen 1
Gen 2
NA
NA
Day 18
Day 18
5 0
0 21 20 19 Day 18
Gen 1
5 0
0 21 20 19
5 0
10 0 10 0 10
With network constraint signals
Demand
10 Day 18
Day 18
0
Fig. 5. Zoomed-in illustration for Fig. 4, where we randomly selected 2 agents for each device type. Bat refers to batteries, HP refers to heat pumps, and Gen refers to generators. Grey dotted lines refer to the min/max power/energy limits, while grey solid lines represent the temperature set point ϑset i .
found this performs even worse than IPPO. We also include a naive baseline to demonstrate that we achieve meaningful learning outcomes. In this baseline, the batteries remain idle, the heat pumps always maintain the indoor temperature at the reference level, and the generators always operate at 100% power. Table I summarises the main hyperparameter settings. All methods use primal–dual learning, batches of 500 episodes per primal step, and the Adam optimiser. GradMA uses a larger learning rate because it evaluates the exact gradient at each update. The policy network for each agent is a small NN with a single hidden layer of 16 neurons (so two layers in total) and tanh activation, which performs better empirically than deeper or wider fully-connected networks. Although each network is small, it should be noted that there are 1,000 independent agents in total with no parameter sharing. The number of training iterations for each method is chosen based on the
TABLE I M AIN TRAINING HYPERPARAMETERS . Hyperparameter GradMAP GradMA-RPO GradMA IPPO MAPPO Actor NN 2 layers with 16 neurons same same same same Critic NN N/A N/A N/A 2×64 2×64 Kdual 20 20 20 40 40 Kprimal 10 10 60 30 30 Kprox 80 80 N/A 80 80 NN learning rate 5×10−4 5×10−4 2×10−3 5×10−4 5×10−4 Dual lr α 150 150 150 150 150 Init. std 0.14 0.14 0.14 1.0 1.0 Entropy coeff. 0.01 0.01 0.01 0.005 0.005 adaptive prox ratio gate [0.2, 2.0] Trust / clipping T (θ) ≤ 0.03 + KL 0.2 none PPO clip 0.2 PPO clip 0.2
training convergence of the cost and constraint violations. Compared with GradMAP and GradMA-RPO, GradMA requires more primal steps to converge because it does not reuse gradient information. IPPO and MAPPO reuse environment rollout data multiple times with clipping, but their sample efficiency is lower; we therefore increase both their number of dual steps (Kdual ) and primal steps (Kprimal ). In addition, because IPPO and MAPPO are sample-based methods, their initial policy standard deviation is set to 1 (0 in log scale), which is a standard choice. The other methods have gradient information available and can therefore use a smaller initial standard deviation of approximately 0.14 (-2 in log scale). To improve model performance for consecutive testing, during training we randomly scale the initial states of all gridedge devices. The script is written in JAX 0.4.38 on Python 3.10.19, and all experiments are run on a single NVIDIA RTX PRO 5000 Blackwell 48 GB GPU (non-data-centre GPU). B. Single-Day Dispatch To illustrate the effectiveness of the proposed GradMAP in handling AC network constraints, Fig. 4 shows the normalised power trajectories for all the 1,000 agents in a single representative day in both the unconstrained case (where the voltage-violation signal is set to zero during training) and the constrained case. As shown, without the voltage-violation signal, the agents learn solely from the price signal, and reduce demand or increase discharge/generation during the afternoon price peak (see Fig. 3 Day 18). In particular, the generators produce power profiles that closely follow the price signal because of the quadratic cost structure. In contrast, when
2000 1500
100 50 0
0
1000
Wall time (s)
2000
0
500
Primal step
1000
Fig. 6. Training convergence plots for the 1,000-agent case. The primal step corresponds to Kprimal in Algorithm 1. Each primal step represents one batched environment-rollout interaction. The horizontal dashed line in the top-left panel marks the lowest training cost achieved by GradMA. Here, the darker colour represents rolling-window average values, while the lighter colour represents per-step values.
the constraint signal is included during training, GradMAP effectively reduces the maximum phase-voltage violation from 0.0362 p.u. to 0.0095 p.u. by reducing generator output and/or increasing battery charging during constrained periods (a clear upward shift in the right panel of Fig. 4). Fig. 5 also provides a zoomed-in visualisation of selected agents. The overall trend is similar. It is noteworthy that one heat pump (HP1) instead reduces its power consumption, possibly due to the three-phase network coupling effect. Satisfying the network constraints also involves a cost trade-off, with the total cost increasing from $4,452 to $4,724 on this day. C. Training Convergence Fig. 6 compares the proposed GradMAP with other learning-based benchmarks in terms of training convergence on the same 1,000-agent case. GradMAP not only converges faster, within only 15 minutes, but also converges to a better solution than all other benchmarks in terms of both operating cost and constraint satisfaction. IPPO is unable to learn effectively in this setting, likely because its critic networks only take local observation and cannot capture the complex network coupling. Although GradMA-RPO also shows fast convergence in constraint satisfaction, its operating cost starts to increase after around 100 primal steps, indicating an inability to effectively balance constraint satisfaction and cost minimization. This supports our argument that our proposed policyoutput-space trust region in GradMAP is more principled and therefore leads to better learning performance than the PPOstyle trust region defined in the probability distribution space [22]. The exact-gradient GradMA is the only method that achieves operating-cost and constraint performance comparable to that of GradMAP, but it is much slower in both wallclock time and primal steps. In particular, focusing on the operating cost in the first row of plots, GradMAP is approximately 3–5× faster in the wall-clock time than GradMA in converging to the same cost level.
GradMA-RPO
GradMA
IPPO
Naive
20
25
30
Battery end
HP end
IPPO
MAPPO
100000 0
Daily max voltage violation (p.u.)
20
0
5
10
0.03 0.02 0.01 0.00
15
Evaluation day
Phase A Phase B Phase C
3 2 1 0
Fig. 7. Out-of-sample benchmark comparison for the 1,000-agent case over 31 consecutive test days. Error bars denote 95-percentile bootstrap confidence intervals.
104
Wall time (s)
Bend + Hend (%)
Volt (%)
40
GradMAP
Daily end-state violation (%)
IPPO
103
GradMAP
GradMA-RPO
GradMA
1.20
Cost ratio
GradMA
1.15 1.10 1.05
102
Daily max voltage violation (p.u.)
Cost (£)
GradMA-RPO
0.03 0.02 0.01 0.00
1.00
Daily end-state violation (%)
GradMAP
2500
Accumulated cost (£)
9
3 2 1
0 IEEE13 EULV-red. IEEE123 IEEE13 EULV-red. IEEE123 (10) (100) (1000) (10) (100) (1000) Fig. 8. Scaling comparison across three different number of agents on different distribution networks.
D. Out-of-Sample Benchmark Comparison Fig. 7 compares the out-of-sample performance over 31 consecutive held-out days. GradMAP always achieves the lowest cumulative cost, remaining about 12% cheaper than the naive baseline, while also maintaining among the lowest voltage and end-state violations. E. Scaling Across Available System Sizes Fig. 8 compares GradMAP and other benchmarks across different problem scales: the 10-, 100-, and 1,000-agent cases. As can be seen, GradMAP requires the lowest learning time across all problem scales due to effective gradient reuse, and can be 4.5× faster than GradMA and 10× faster than IPPO and MAPPO. GradMAP is always the lowest-operating-cost method and consistently maintains nearly the lowest constraint violations. VI. C ONCLUSION This paper proposes GradMAP, a gradient-based multiagent proximal learning framework for coordinating largescale grid-edge flexibility in a fully decentralised manner. By
10
exploiting the differentiability of power-system environments and reusing exact gradients multiple times within a policyoutput-space trust region (rather than a PPO-style distributionspace trust region), GradMAP achieves both computational efficiency and effective learning. Case studies with up to 1,000 agents on the three-phase unbalanced IEEE 123-bus feeder show that GradMAP achieves lower operating cost and better AC power-flow constraint satisfaction than all benchmark methods, including other gradient-based learning approaches, such as our proposed multi-agent extensions of proximal gradient learning and self-supervised learning, as well as multi-agent reinforcement-learning methods. Compared with the best-performing benchmark, GradMAP reduces training time by a factor of 3–5. Using only a single workstation-class NVIDIA RTX PRO 5000 Blackwell 48 GB GPU, GradMAP achieves the best performance on the 1,000-agent case with only 15-minutes of training. Several directions remain for future work. One is to go beyond first-order gradients and incorporate higher-order derivative reuse into GradMAP, as commonly exploited in optimisation-based methods. Another is to extend GradMAP to discrete-action settings. One direct approach is to use the Straight-Through Estimator, which is biased but performs well empirically [40]. The proposed GradMAP method could also be promising for other decision-making domains where there are good differentiable environment simulators but evaluating gradients or higher-order derivatives is expensive, such as robotics and autonomous driving. R EFERENCES [1] The International Energy Agency (IEA), “Digitalisation and Energy,” www.iea.org/reports/digitalisation-and-energy, The International Energy Agency (IEA), Technical Report, 2017, Accessed: 2025-05-22. [2] Z. Chen, C. Guo, S. Dong, Y. Ding, and H. Mao, “Distributed robust dynamic economic dispatch of integrated transmission and distribution systems,” IEEE Transactions on Industry Applications, 2021. [3] J. Zhai et al., “Distributionally robust joint chance-constrained dispatch for integrated transmission-distribution systems via distributed optimization,” IEEE Transactions on Smart Grid, 2022. [4] Y. Zhou, C. Essayeh, and T. Morstyn, “Aggregated feasible active power region for distributed energy resources with a distributionally robust joint probabilistic guarantee,” IEEE Transactions on Power Systems, 2025. [5] H. Früh et al., “Coordinated vertical provision of flexibility from distribution systems,” IEEE Transactions on Power Systems, 2023. [6] Z. Tan, Z. Yan, H. Zhong, and Q. Xia, “Non-iterative solution for coordinated optimal dispatch via equivalent projection—part i: Theory,” IEEE Transactions on Power Systems, vol. 39, no. 1, pp. 890–898, 2024. [7] H. Zhou et al., “LSTM-based energy management for electric vehicle charging in commercial-building prosumers,” Journal of Modern Power Systems and Clean Energy, vol. 9, no. 5, pp. 1205–1216, 2021. [8] S. Gao, C. Xiang, M. Yu, K. T. Tan, and T. H. Lee, “Online optimal power scheduling of a microgrid via imitation learning,” IEEE Transactions on Smart Grid, vol. 13, no. 2, pp. 861–876, 2021. [9] R. S. Sutton and A. G. Barto, Reinforcement Learning: An Introduction, 2nd ed. Cambridge, MA: MIT Press, 2015. [10] J. Cao et al., “Deep reinforcement learning-based energy storage arbitrage with accurate lithium-ion battery degradation model,” IEEE Transactions on Smart Grid, 2020. [11] D. Qiu, Y. Wang, W. Hua, and G. Strbac, “Reinforcement learning for electric vehicle applications in power systems: A critical review,” Renewable and Sustainable Energy Reviews, vol. 173, p. 113052, 2023. [12] F. Charbonnier, T. Morstyn, and M. D. McCulloch, “Scalable multiagent reinforcement learning for distributed control of residential energy flexibility,” Applied energy, vol. 314, p. 118825, 2022. [13] F. Charbonnier et al., “Centralised rehearsal of decentralised cooperation: Multi-agent reinforcement learning for the scalable coordination of residential energy flexibility,” Applied Energy, 2025.
[14]
“Network topology optimisation (nia2 neso087),” National Energy System Operator, Tech. Rep., Dec. 2025. [Online]. Available: https://www.neso.energy/document/373681/download [15] A. Kaya, A. J. Conejo, and S. Rebennack, “Fifty years of power systems optimization,” European Journal of Operational Research, 2025. [16] H.-D. Chiang, Z.-Y. Wang, and L. Zeng, “Dynamic relationship between kkt saddle solutions and optimal solutions in ac opf problems,” IEEE Transactions on Power Systems, vol. 39, no. 1, pp. 1637–1646, 2023. [17] I. Clavera, V. Fu, and P. Abbeel, “Model-augmented actor-critic: Backpropagating through paths,” arXiv preprint arXiv:2005.08068, 2020. [18] J. Liu, H. Guo, and Q. Chen, “Deep differentiable reinforcement learning for generic energy storage day-ahead bidding with transformer neural networks,” in 2024 PESGM. IEEE, 2024, pp. 1–5. [19] Y. Ye et al., “Physics-guided safe policy learning with enhanced perception for real-time dynamic security constrained optimal power flow,” Journal of Modern Power Systems and Clean Energy, 2025. [20] M. A. Z. Mora, M. Peychev, S. Ha, M. Vechev, and S. Coros, “Pods: Policy optimization via differentiable simulation,” in International Conference on Machine Learning. PMLR, 2021, pp. 7805–7817. [21] S. Son, L. Zheng, R. Sullivan, Y.-L. Qiao, and M. Lin, “Gradient informed proximal policy optimization,” Advances in Neural Information Processing Systems, vol. 36, pp. 8788–8814, 2023. [22] H. Zhong, X. Wang, Z. Li, and L. Huang, “Reparameterization proximal policy optimization,” arXiv preprint arXiv:2508.06214, 2025. [23] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” arXiv preprint arXiv:1707.06347, 2017. [24] S. Park and P. Van Hentenryck, “Self-supervised learning for largescale preventive security constrained dc optimal power flow,” IEEE Transactions on Power Systems, vol. 40, no. 3, pp. 2205–2216, 2024. [25] G. Chen and J. Qin, “Neural risk limiting dispatch in power networks: Formulation and generalization guarantees,” IEEE Transactions on Power Systems, 2025. [26] K. Chen, B. Knueven, and W. Jones, “A hard-constrained nn learning framework for rapidly restoring ac-opf from dc-opf,” arXiv preprint arXiv:2602.06255, 2026. [27] A. Anrrango, A. Quisaguano, G. E. Constante-Flores, and C. Li, “Selfsupervised learning of parametric approximation for security-constrained dc-opf,” arXiv preprint arXiv:2601.13486, 2026. [28] G. Chen, J. Qin, and H. Zhang, “Model-free self-supervised learning for dispatching distributed energy resources,” IEEE Transactions on Smart Grid, vol. 16, no. 2, pp. 1287–1300, 2024. [29] C. Yu et al., “The surprising effectiveness of ppo in cooperative multiagent games,” Advances in neural information processing systems, vol. 35, pp. 24 611–24 624, 2022. [30] A. Thingvad, C. Ziras, J. Hu, and M. Marinelli, “Assessing the energy content of system frequency and electric vehicle charging efficiency for ancillary service provision,” in UPEC. IEEE, 2017, pp. 1–6. [31] M. Naumann, F. B. Spingler, and A. Jossen, “Analysis and modeling of cycle aging of a commercial lifepo4/graphite cell,” Journal of Power Sources, vol. 451, p. 227666, 2020. [32] H. Hao, B. M. Sanandaji, K. Poolla, and T. L. Vincent, “Aggregate flexibility of thermostatically controlled loads,” IEEE Transactions on Power Systems, vol. 30, no. 1, pp. 189–198, 2014. [33] A. Bernstein and E. Dall’Anese, “Linear power-flow models in multiphase distribution networks,” in ISGT-Europe. IEEE, 2017. [34] C. Heinrich, P. Fortenbacher, A. Fuchs, and G. Andersson, “Pvintegration strategies for low voltage networks,” in ENERGYCON, 2016. [35] M. Bazrafshan and N. Gatsis, “Comprehensive modeling of three-phase distribution systems via the bus admittance matrix,” IEEE Transactions on Power Systems, vol. 33, no. 2, pp. 2015–2029, 2018. [36] K. P. Schneider et al., “Analytic considerations and design basis for the ieee distribution test feeders,” IEEE Trans on power systems, 2017. [37] R. Wardle et al., “Dataset (TC1a): basic profiling of domestic smart meter customers,” Accessed: Sep, vol. 27, 2020. [38] S. Pfenninger and I. Staffell, “Renewables.ninja.” [Online]. Available: https://www.renewables.ninja/ [39] Octopus Energy, “Agile historical data portal,” https://agile. octopushome.net/historical-data, 2026, accessed: 2026-03-26. [40] Y. Bengio, N. Léonard, and A. Courville, “Estimating or propagating gradients through stochastic neurons for conditional computation,” arXiv preprint arXiv:1308.3432, 2013.