Accelerated decentralized dispatch of energy storage via constrained gradient tracking

1 Introduction

The integration of renewable energy sources, particularly wind and solar power, has increased across many power systems under carbon-neutrality targets (Denholm et al., 2021). However, the intermittency, volatility, and anti-peaking characteristics of wind and solar generation can mismatch supply with load and contribute to renewable energy curtailment (Mohandes et al., 2019). These operating conditions increase the need for flexible resources and support the transition from the traditional “source follows load” paradigm to a “source-grid-load-storage interaction” framework (Zang et al., 2024). In parallel, distributed energy resources (DERs), such as rooftop photovoltaics and dispersed wind turbines, are reshaping the spatial structure of power systems. Centralized grid paradigms may have limited ability to accommodate high penetrations of localized DERs (Hatziargyriou et al., 2007). Distributed smart grids, including microgrids and active distribution networks, therefore provide an important topology for integrating distributed renewable energy (Lasseter, 2011). In these decentralized networks, distributed energy storage systems (ESS) support local power balancing, peak shaving, and resilience enhancement (Oudalov et al., 2007).

Existing energy storage technologies also have limitations in distributed smart grid scenarios. Lithium-ion batteries currently dominate the market; however, they face constraints related to life-cycle cost, cycle life, thermal runaway risk, and the economics of medium-to-long-duration storage (Finegan et al., 2015; Schmidt et al., 2019). Pumped hydro storage is mature but depends on geographical and topographical conditions, which can limit deployment in urban or flat-terrain regions (Pérez-D́ıaz et al., 2015). Traditional compressed air energy storage usually requires underground caverns or large high-pressure gas tanks, and its system complexity can restrict dynamic response (Chen et al., 2009). These limitations motivate the development of storage technologies that are modular, geographically flexible, and suitable for medium-to-long-duration applications. Carbon dioxide energy storage systems (CO₂-ESS) have received attention in this context (Wang et al., 2024). By using carbon dioxide as the working fluid, CO₂-ESS can exploit its thermophysical properties near the critical point to enable compact system designs, relatively high energy density, and non-flammable operation without geological constraints (Zhang and Wang, 2017). These characteristics make CO₂-ESS a relevant candidate for distributed smart-grid applications.

1.1 Related works

Research on energy storage technologies has addressed a range of flexibility requirements in modern power grids. Electrochemical storage, particularly lithium-ion batteries, dominates the distributed energy storage market. Prior studies have examined their use in grid frequency regulation, peak shaving, and microgrid energy management (Hannan et al., 2017). However, recent literature indicates that although batteries are effective for short-term, high-frequency responses, they may be less suitable for medium-to-long-duration storage because of life-cycle costs, capacity degradation over deep cycles, and safety risks associated with thermal runaway (Luo et al., 2025). Traditional bulk mechanical storage technologies, such as pumped hydro storage and compressed air energy storage, provide long-duration storage capability, but their deployment in distributed or urban environments is constrained by geographical, topographical, or geological requirements (Liu et al., 2024). These observations motivate the development of storage options that combine geographic flexibility with longer-duration operation and safety.

CO₂-ESS has been studied as a viable option. Leveraging the thermophysical properties of CO₂ near its critical point, CO₂-ESS operates in closed-loop thermodynamic cycles that do not require underground caverns (Qi et al., 2025). Recent studies have assessed the thermodynamic feasibility of transcritical and supercritical CO₂-ESS configurations and reported advantages in energy density, compactness, and round-trip efficiency (Yu et al., 2026; Hsiao et al., 2024). The non-flammable nature of CO₂ also supports its use in high-density distributed storage applications.

Despite this potential, the integration of CO₂-ESS into distributed smart grids remains less developed than its thermodynamic modeling. One challenge is system design and optimal capacity allocation (Bian et al., 2025). Unlike electrochemical batteries, which are often represented by relatively simple state-of-charge equations, CO₂-ESS is a thermo-mechanical system. Its operational efficiency depends on parameters such as compressor/turbine isentropic efficiencies, operating pressures, and heat exchanger effectiveness, which can vary under off-design grid conditions (Bian et al., 2025). Thus, capacity allocation and sizing methods developed for simplified battery dynamics may not directly apply to CO₂-ESS. Existing grid-level storage planning models often simplify or omit the non-linear thermodynamic constraints of mechanical storage, which can lead to suboptimal, oversized, or physically infeasible configurations (Hemmati and Saboori, 2026; Mumtahina et al., 2025). Incorporating thermodynamic models into grid-level economic dispatch and allocation also introduces non-convexities, making the joint optimization of system design and spatial allocation difficult (Aslam et al., 2024). This motivates a framework that links the internal thermodynamic design of CO₂-ESS with its distributed allocation strategy in power networks.

Beyond structural design and capacity allocation, distributed energy storage also requires coordinated control within the smart grid framework. Prior work on energy management and power dispatch has often used hierarchical control architectures with primary, secondary, and tertiary levels (Khan et al., 2025; Albasheri et al., 2023). Decentralized approaches such as droop control support plug-and-play operation and reliability, but they may sacrifice global economic optimality (Ahmad et al., 2025). Distributed cooperative control based on multi-agent systems and consensus protocols has therefore been used to coordinate power sharing and grid regulation via neighbor-to-neighbor communications, reducing dependence on a centralized controller (Sarathkumar et al., 2025; Zhu et al., 2025).

To address renewable generation uncertainty and load fluctuations, optimization-based control algorithms, particularly Model Predictive Control (MPC), have also been adopted in smart grid energy management. MPC can handle multivariable constrained optimization over a rolling horizon and improve the operation of active distribution networks under uncertainty (Yaghoubi et al., 2023; Lim et al., 2025). However, many existing coordinated control algorithms are designed around the fast, approximately linear response characteristics of battery ESS. Applying these electrical-centric strategies to CO₂-ESS may create operational mismatches. Unlike electrochemical batteries, CO₂-ESS relies on turbomachinery and therefore has thermo-mechanical inertia and transient limitations (Liu C. et al., 2026).

Mechanical ramp rates, start-up/shutdown delays, and the coupling between thermodynamic parameters and electrical power output create a timescale disparity (Afshari and Groß, 2024). If a distributed consensus algorithm or an MPC controller issues high-frequency power commands that exceed the operating bounds of the CO₂-ESS thermodynamic cycle, compressor surge, turbine damage, or grid instability may occur. Existing studies on coordinated control of integrated energy systems do not always capture the electro-thermal-mechanical transient coupling required for mechanical storage (Mohammadi et al., 2026). This motivates coordinated control methods that account for both fast electrical dynamics and slower thermodynamic constraints.

Capacity allocation and coordinated dispatch across geographically dispersed storage units require optimization methods that can operate with distributed information. As power grids move toward decentralized architectures, centralized optimization can face computational burdens, communication bottlenecks, and data privacy concerns (Carnevale et al., 2025). Distributed optimization algorithms have therefore become common tools in power system research. Representative methods include ADMM (Steven et al., 2026), consensus-based subgradient methods (Le et al., 2025), and distributed dual ascent, which have been applied to large-scale, spatially coupled grid optimization problems.

In the current literature, these distributed frameworks have been applied to diverse power system applications, including Distributed Optimal Power Flow (Khanghah et al., 2025; Hasanzadeh and Kargarian, 2025), economic dispatch (Sun et al., 2025), and peer-to-peer energy trading among networked microgrids (Dong et al., 2025). ADMM, in particular, is favored for its convergence properties and its ability to decompose complex, system-level objectives into localized subproblems. This decomposition can reduce the computational load on the central coordinator and minimize the need to share operational data from individual distributed storage nodes and prosumers (Liu et al., 2025).

Despite their maturity in electrical networks, these algorithms are less developed for multi-energy systems with coupled thermo-mechanical dynamics. Many existing distributed dispatch frameworks rely on simplified linear energy degradation models (Low, 2024). When non-linear thermodynamic equality constraints are included in distributed CO₂-ESS models, such as polynomial turbomachinery efficiency curves and coupled pressure-temperature state variables, standard distributed algorithms may exhibit convergence oscillations or non-negligible optimality gaps (Li et al., 2016). A distributed optimization framework for CO₂-ESS should therefore handle thermo-mechanical coupling while preserving convergence properties under a distributed communication topology.

More recently, a complementary line of work addresses the physical–information coupling of distributed dispatch through digital representations and cyber-physical co-simulation frameworks, in which high-fidelity physical models are mirrored and co-simulated with the distributed decision layer to validate coordination strategies under realistic grid dynamics (Tian et al., 2026; Liu Z. et al., 2026). These approaches and the primal-dual tracking architecture proposed in this study pursue the same high-level objective–coordinating physical storage assets with distributed intelligence–and are in principle complementary: the co-simulation environments can serve as a high-fidelity validation platform for the decentralized policies computed by gradient-tracking methods, while the latter provide the real-time dispatch logic that the former can exercise. We revisit this synergy in Section 7.

Beyond the algorithmic advances, recent studies at the intersection of generation-grid-load-storage coordination have examined the system-level allocation of heterogeneous storage assets and their validation under renewable variability. E.g., Meng et al. (2026) investigates coordinated storage operation across the generation-load-storage nexus, providing an empirical perspective on how distributed storage dispatch interacts with renewable forecasting uncertainty and grid-level constraints. The C-DNGT framework proposed in this study complements these coordination studies: it provides a provably convergent distributed optimization engine that can serve as the computational kernel within larger generation-grid-load-storage coordination architectures, where the higher-level scheduling problem decomposes across heterogeneous assets and the low-level power-balance and thermo-mechanical constraints are enforced in a decentralized manner. We therefore view our contribution as filling the algorithmic gap between the system-level coordination perspective and the per-unit physics of individual storage assets.

1.2 Main contributions

To address the issues in system design and coordinated dispatch, this research develops a distributed optimization framework for CO₂-ESS energy management. The main contributions of this research are summarized as follows:

• Decentralized Energy Management Framework: We propose a collaborative optimization framework to integrate Carbon Dioxide Energy Storage Systems (CO₂-ESS) into smart grids. By modeling the grid as a multi-agent cyber-physical network, the framework jointly optimizes the thermo-mechanical states and the economic operating costs of distributed CO₂-ESS units.
• Algorithm Design with Grid Constraint Handling: To solve the formulated economic dispatch problem, we develop the Constrained Distributed Nesterov-Accelerated Primal-Dual Gradient Tracking (C-DNGT) algorithm. The algorithm incorporates the global power balance constraint into a primal-dual distributed update and uses localized peer-to-peer communications without relying on a centralized coordinator.
• Convergence Analysis: We provide a theoretical convergence analysis of the proposed control strategy. Under standard assumptions, the C-DNGT algorithm converges linearly to the KKT saddle point of the economic dispatch problem.
• Evaluation with Grid Load Data: Numerical experiments are conducted using PJM Interconnection load profiles across heterogeneous storage configurations. The results compare the proposed algorithm with classical distributed methods in terms of convergence speed, communication overhead, and constraint handling.

1.3 Paper organizations

The remainder of this study is organized as follows: Section 2 introduces the necessary theoretical foundations, including convex analysis and graph theory. Section 3 details the system modeling of the distributed smart grid and introduces the basic collaborative optimization gradient tracking method. Section 4 presents the distributed Nesterov-accelerated gradient tracking algorithm incorporating power balance constraints. Section 5 provides the convergence analysis, including supporting lemmas and proof arguments. Section 6 presents numerical results, including load-data-based simulations, heterogeneity tests, and performance comparisons. Finally, Section 7 concludes the study and discusses future research prospects.

2 Theoretical foundations

In this section, we introduce the fundamental mathematical preliminaries that underpin the formulation and convergence analysis of the proposed distributed optimization algorithm. The theoretical basis primarily draws on convex analysis and algebraic graph theory.

2.1 Convex analysis and optimization basics

In the distributed smart grid, the operational cost functions for individual CO₂-ESS and distributed generation units are typically formulated as convex functions. Consider a continuously differentiable function f: ℝ^*n* → ℝ. The basic definitions related to convex optimization are given as follows:

Definition 1 (Convexity): A continuously differentiable function f(x) is convex if, for any x, y∈ℝ^*n*, the following inequality holds (Equations 1–5):

f(y)≥f(x)+∇f(x)T(y-x)(1)

where ∇f(x) denotes the gradient of f evaluated at x.

To support the existence of a unique global optimum and the subsequent convergence analysis of the gradient tracking algorithm, we introduce the concepts of strong convexity and Lipschitz smoothness.

Definition 2 (Strong Convexity): A function f(x) is said to be μ-strongly convex if there exists a positive constant μ > 0 such that for any x, y ∈ ℝ^*n*:

f(y)≥f(x)+∇f(x)T(y-x)+μ2∥y-x∥2(2)

where ||·|| denotes the standard Euclidean norm. The parameter μ is referred to as the strong convexity modulus.

Definition 3 (Lipschitz Continuous Gradient): A function f(x) has an L-Lipschitz continuous gradient (often called L-smooth) if there exists a constant L>0 such that for any x, y ∈ ℝ^*n*:

∥∇f(x)-∇f(y)∥≤L∥x-y∥(3)

Combining Definition 2 and Definition 3, if a cost function f(x) is both μ-strongly convex and L-smooth, it satisfies the condition μ ≤ L.

2.2 Graph theory basics

In a distributed smart grid, the CO₂-ESS nodes and their communication topology are modeled using algebraic graph theory. The communication network among N distributed agents (nodes) is represented by an undirected graph, denoted as G=(V,E).

Here, V={1,2,…,N} defines the set of nodes, and E⊆V×V defines the set of communication links (edges). An edge (i,j)∈E indicates that node i and node j can bidirectionally exchange information. The set of neighbors of node i is defined as Ni={j∈V∣(i,j)∈E}.

To perform distributed consensus and gradient tracking, a weighted mixing matrix (or adjacency matrix) W=[wij]∈ℝN×N is associated with the graph G. The elements w_*ij* satisfy the following properties:

wij>0 if(i,j)∈Eor i=j; wij=0 otherwise.(4)

To achieve consensus and track the global gradient, the distributed optimization algorithm requires the following standard assumption about the communication graph and its weight matrix:

Assumption 1: The undirected communication graph G is connected, meaning there is at least one path between any pair of nodes. Furthermore, the weight matrix W is doubly stochastic, which satisfies:

∑j=1Nwij=1, ∀i∈V and ∑i=1Nwij=1, ∀j∈V(5)

In matrix notation, this is expressed as W1 = 1 and 1^*T*W = 1^*T*, where 1 is the N-dimensional column vector of all ones. The doubly stochastic property preserves the average value of the nodes' states during iterative information exchange, which is a mathematical prerequisite for handling the global power balance constraints.

3 Collaborative optimization framework for CO₂-ESS

In this section, we develop the physical and mathematical models for the distributed smart grid integrated with CO₂-ESS. Subsequently, we formulate the system-level optimization problem, laying the groundwork for the distributed gradient tracking algorithm.

3.1 System modeling and problem formulation

Figure 1 illustrates the proposed two-layer cyber-physical system architecture for the distributed coordination of CO₂ energy storage systems.

The physical layer encompasses the electrical infrastructure, in which geographically dispersed nodes are interconnected via power lines. Each bus integrates local generation, native demand, and a dedicated CO₂ storage unit to provide operational flexibility.

The overlaid cyber layer is modeled as a multi-agent system. Each physical node is governed by an intelligent agent i∈V. These agents interact over a sparse, undirected communication graph G=(V,E), where the information exchange between adjacent nodes i and j is weighted by w_*ij* > 0. This topology enables the peer-to-peer data flow required for distributed gradient tracking. The vertical dashed lines denote cyber-physical closed-loop interactions. Agents acquire local state measurements and, after algorithmic convergence, transmit dispatch setpoints back to the hardware, enabling distributed energy management without a single centralized coordinator.

Consider a distributed smart grid operating over a discrete time horizon T={1,2,…,T} with a time interval Δt. The grid is modeled as an undirected connected graph G=(V,E) with N nodes, where V={1,2,…,N}. To capture the distributed nature of modern power systems, we assume that each node i∈V operates as a prosumer microgrid.

Unlike conventional electrochemical batteries governed by simple linear state-of-charge equations, the CO₂-ESS is a complex thermo-mechanical system. During charging, gaseous CO₂ is compressed, cooled, and stored as a liquid, with compression heat recovered. During discharging, the liquid CO₂ is pressurized, heated, and expanded through a turbine to generate electricity.

For the optimization framework, we define the thermo-mechanical state of the i-th CO₂-ESS unit at time t as si,t=[Mi,tco2,Pri,t,Ti,thot]T, and the decision variables as ui,t=[Pi,tc,Pi,te,Qi,th]T. (To avoid the notation collision flagged by the reviewers, the continuous state is denoted by s_*i,t*, whereas the condensed optimization decision vector introduced below is denoted by x_*i*; the two are distinct objects.) The specific physical meanings are defined as follows:

• Mi,tco2: The mass of liquid CO₂ in the high-pressure storage tank (kg).
• Pr_*i,t*: The pressure of the low-pressure gaseous CO₂ reservoir (MPa).
• Ti,thot: The temperature of the Thermal Energy Storage (TES) hot tank (K).
• Pi,tc≥0: The electrical power consumed by the compressor during charging (kW).
• Pi,te≥0: The electrical power generated by the expander/turbine during discharging (kW).
• Qi,th: The thermal power exchanged with the TES system (kW).

The discrete-time state-space dynamics of the CO₂-ESS can be formulated as:

Mi,t+1co2=Mi,tco2+(ηmcPi,tc-Pi,teηme)Δt(6a)

Pri,t+1=Pri,t-κp(ηmcPi,tc-Pi,teηme)Δt(6b)

Ti,t+1hot=Ti,thot+ΔtCtesMtes(Qi,th+σcPi,tc-Qloss)(6c)

where ηmc and ηme are the mass-power conversion coefficients for the compressor and expander, respectively, which integrate the isentropic efficiencies and specific enthalpies of the CO₂ cycle. κ_*p* is the pressure variation coefficient based on the ideal gas law. C_*tes* and M_*tes* denote the specific heat capacity and mass of the heat transfer fluid in the TES. Q_*loss* represents the ambient thermal dissipation. σ_*c* > 0 is the heat-of-compression coupling coefficient: compressing CO₂ releases thermal energy that is recovered into the TES, so the hot-tank temperature now depends explicitly on the compression power Pi,tc. Together with Equation 6, this couples the mass/pressure state to the thermal state through the compression work, removing the previous decoupling between the pressure and temperature dynamics.

Furthermore, the state and decision variables must satisfy physical boundary constraints at each time step to prevent mechanical damage and to respect material thermal limits:

simin≤si,t≤simax, uimin≤ui,t≤uimax.(7)

Furthermore, the mechanical ramping capability of the CO₂-ESS limits the per-step variation of the compressor and expander powers:

with Pi,0c=Pi,0e=0 (cold-start reference) or the realized value carried over from the preceding horizon under the rolling-horizon scheme. These inequalities are affine; together with the box bounds Equation 7 and the unrolled dynamics Equation 6, they are absorbed into the lifted feasible set Xi in Equation 12, so the Euclidean projection PXi required by the C-DNGT primal step remains a convex quadratic program (Section 4).

We note that the formulation above advances all state and decision variables on a single, unified scheduling time step Δt (set to Δt = 1h over a T = 24 horizon in our simulations). This choice reflects a deliberate separation of time scales: the electrical subsystem—the compressor/expander power set-points and the associated grid power balance—responds on a much faster scale than the thermo-mechanical state evolution governed by Equation 6. Within each scheduling interval, we therefore treat the electrical power set-points as quasi-static, i.e., the electrical dynamics are assumed to have settled to their steady state relative to the slower mass, pressure, and temperature transients. The two layers are coupled through the static power-balance map in Equation 6 and the box constraints Equation 7, rather than through a coupled differential-algebraic model. Although the present implementation uses a common step size for clarity and reproducibility, the proposed distributed coordination is, by construction, compatible with multi-timescale operation: the gradient-tracking layer exchanges only aggregate mismatch and gradient information at the scheduling instants and is agnostic to how each node internally discretizes its local dynamics. Heterogeneous local step sizes Δt_*i* or a hierarchical two-rate scheme—a fine intra-interval resolution for the fast electrical sub-problem and a coarse resolution for the slow thermo-mechanical states—can thus be embedded without altering the communication protocol. A fully multi-rate extension, in which the fast electrical sub-problem is solved at a finer resolution and coupled to the slow states through an inner quasi-static map with the tracking iterations synchronized only at the coarse instants, is left for future research.

3.1.1 Model assumptions and limitations

The CO₂-ESS model above is a first-order, quasi-steady affine surrogate of the actual thermo-mechanical plant. Three scoping choices deserve explicit statement so that the convex formulation and its convergence guarantees (Theorems 1 and 2) are not over-interpreted. (i) Off-design efficiency. The conversion coefficients ηmc and ηme are taken as constants corresponding to a first-order linearization about the nominal operating point. This is accurate within the bounded operating envelope assumed throughout the simulations—the rated-power band [ρ_,1]P¯i and the admissible pressure range of Equation 7—where the higher-order, load-dependent efficiency corrections are second order and are deferred to a piecewise-affine efficiency-map extension. (ii) Mutually exclusive charging and discharging. In a purely continuous convex economic dispatch, exact complementarity is implicitly enforced by two economic–physical mechanisms rather than by an explicit integer constraint. First, both Pi,tc and Pi,te enter the local cost f_*i* with strictly positive coefficients, so at any optimum of problem Equation 18 at least one of Pi,tc,Pi,te is zero for every (i, t): the storage never simultaneously compresses and expands at the optimal dispatch. Second, and more fundamentally, the round-trip efficiency of the CO₂-ESS is strictly below unity (ηmcηme<1); charging and discharging within the same interval would dissipate energy through the thermodynamic cycle while delivering no net useful work, which is strictly suboptimal under the cost-minimizing objective. A strict per-interval operating-mode complementarity (Big-M) constraint, which would enforce mutual exclusion over the entire feasible set rather than only at the optimum, is included in the unit-commitment extension below. We note that introducing a 0–1 integer variable to hard-enforce this complementarity would turn the problem into a mixed-integer program and break the smooth-convex setting on which the C-DNGT convergence proof (Theorems 1 and 2) relies; the implicit enforcement above is therefore a deliberate modeling choice that preserves tractability, not an omission. (iii) Start-up delays and minimum up/down times. These effects require binary/integer decision variables and render the problem a mixed-integer program, which lies outside the smooth-convex setting required by the C-DNGT convergence analysis. We therefore treat each CO₂-ESS as continuously available within the dispatch horizon and defer start-up/shut-down logic to a two-layer scheme in which C-DNGT performs the continuous economic dispatch and an outer unit-commitment MILP decides the operating modes. (iv) Strong-convexity regularization for the convergence analysis. The local cost f_*i* is affine in the compressor/expander powers and quadratic only in the curtailment, so the unregularized global objective is convex but not strongly convex in the full decision vector. To satisfy the strong-convexity premise required by the linear-convergence proofs of Theorems 1 and 2 without altering the dispatch solution in any material way, we add the small Tikhonov term Equation 25 (Section 4); its bias is O(μ) in the decision and O(μ2) in the cost. It is driven to zero as μ → 0. The augmented-Lagrangian penalty ρ in Equation 24 provides a second, independent source of coercivity (the Luo–Tseng error bound) that yields global linear convergence of the primal-dual residual even on the original convex problem, as detailed in Section 5. These assumptions are standard for dispatch-oriented studies and are stated here for transparency.

The local power balance at node i includes the forecasted RES generation Pi,tRES and the local load demand Pi,tload. To maintain intra-node power equilibrium and mitigate renewable volatility, the actual RES power output may be curtailed. Let Pi,tcurt≥0 denote the curtailed renewable power.

The objective of the distributed grid is to minimize the total operational cost over the time horizon T, including the operational and degradation costs of the CO₂-ESS and the penalty for wind and solar curtailment. The local cost function f_*i*(u_{*i, t*}) for node i is defined as:

fi(ui,t,Pi,tcurt)=∑t=1T(cicPi,tc+ciePi,te+ γi(Pi,tcurt)2),(9)

where cic and cie are the marginal operational and wear-and-tear cost coefficients for the compressor and turbine machinery, respectively. γ_*i* > 0 is the penalty factor for renewable energy curtailment. A quadratic penalty is employed to discourage localized curtailment.

Consequently, the global optimization problem for the distributed smart grid is formulated to minimize the aggregate costs across all N nodes:

min F=∑i=1Nfi(ui,t,Pi,tcurt).(10)

3.2 Distributed accelerated gradient tracking algorithm

Based on the system model formulated in Section 3-A, the collaborative optimization of the distributed CO₂-ESS network can be formulated as a multi-agent consensus optimization problem. To support algorithmic design and convergence analysis, we concatenate all local decision variables (compressor/expander power, thermal power, and curtailment) over the entire time horizon T into a single condensed variable vector xi∈ℝd for each node i.

To make this abstraction fully explicit, we define the condensed decision vector as the concatenation of the controllable power set-points over the entire scheduling horizon,

xi=vec({Pi,tc,Pi,te,Pi,tcurt}t=1T)∈ℝd, d=3T,(11)

where each per-step set-point is drawn from the local decision variable u_{*i, t*} introduced in Section 3-A. With this definition, the global objective Equation 16 and the constrained problem Equations 18–20 operate directly on the static vector x_*i* rather than on the per-step variables.

The time-coupled thermo-mechanical dynamics Equation 6 are embedded into the optimization exclusively through the local feasible set Xi. Because Equation 6 is affine in (x_{*i, t*}, u_{*i, t*}), unrolling the recursion over T with the given initial thermo-mechanical state si,0=[Mi,0co2,Pri,0,Ti,0hot]⊤ propagates the physical box bounds Equation 7 on the state and control trajectories into a set of affine inequalities on the control vector x_*i*. Collecting these inequalities yields the lifted representation

Xi={xi∈ℝd : c_i≤Eixi≤c¯i(si,0)},(12)

where E_*i* is the block-bidiagonal matrix induced by the conversion coefficients (ηmc,ηme,κp,Ctes,Mtes) and the time step Δt, and the bounds c_i,c¯i(·) encode both the control limits and the state limits from Equation 7. Consequently Xi is an affine (hence convex) set, which ensures that the Euclidean projection in Equation 26 is well-defined and that the static problem Equations 18–20 is a faithful convex representation of the underlying time-coupled storage operation.

3.2.1 Implicit assumption on steady-state operation

The static problem Equations 18–20 is solved over one prescribed scheduling horizon T with the initial state s_{*i*, 0} inherited from the preceding operating period. For a continuously operating storage fleet, we adopt a rolling-horizon (model-predictive) dispatch: after each horizon, the realized terminal state s_{*i, T*} becomes the initial state s_{*i*, 0} of the next horizon, while only the first set of power set-points is implemented. Within any single horizon, the problem therefore remains static and convex, and periodicity is enforced implicitly across successive horizons rather than through an explicit terminal-state equality inside Equations 18–20. If a strictly periodic daily operation is desired, the closing constraint s_{*i, T*} = s_{*i*, 0} can be added to Equation 12; this keeps the problem static and convex.

Consequently, the (unconstrained) distributed minimization that the accelerated gradient-tracking building block targets can be written in the standard consensus form

minx1,…,xNF(x)=∑i=1Nfi(xi) s.t. x1=x2=⋯=xN,(16)

which is the canonical template used to analyze accelerated gradient tracking in the distributed-optimization literature (Theorem 5.2 below is proved for this form). We stress that this primal-consensus form is a methodological template for the unconstrained analysis and does not describe the economic-dispatch application studied in Section 4. In the dispatch problem, the local decision vectors x_*i* are intrinsically heterogeneous: storage units with different ratings, costs, and thermo-mechanical limits must not be forced to output identical set-points. The only network-wide equality is the single global power-balance constraint Equation 19; the distributed method therefore achieves consensus on the dual (price) variable λ_*i* and on the global power-mismatch tracker y_*i*–not on the primal x_*i*. Every node keeps its own local schedule, and the “agreement” the network reaches is agreement on the system marginal price and on how much total imbalance remains. This dual-decomposition interpretation is made explicit in Section 4.

To solve Equation 16 while mitigating the slow convergence and steady-state error issues that can arise in traditional Distributed Gradient Descent (DGD), we design a Consensus-based Distributed Accelerated Gradient Tracking (D-NGT) algorithm. The proposed algorithm integrates Nesterov acceleration with dynamic gradient tracking.

The core algorithmic variables for node i at iteration index k are defined as follows:

• xik∈ℝd: The local estimate of the optimal decision variable at iteration k.
• zik∈ℝd: The Nesterov extrapolated state (momentum variable), which incorporates historical update directions to accelerate convergence.
• yik∈ℝd: The auxiliary gradient tracking variable. It aims to estimate and track the global average gradient 1N∑j=1N∇fj(xjk) locally.
• α > 0: The unified constant step size.
• β ∈ (0, 1): The Nesterov momentum coefficient.
• w_*ij*: The communication weight between node i and node j, obtained from the doubly stochastic matrix W defined in Assumption 1.

At k = 0, the system is initialized as follows: each node locally initializes its decision state xi0 and sets its momentum state zi0=xi0. The gradient tracker is initialized as the local gradient evaluated at the initial momentum state, i.e., yi0=∇fi(zi0).

The execution framework of the proposed D-NGT algorithm is presented in Algorithm 1. At each iteration k, every CO₂-ESS unit communicates with its directly connected neighbors j∈Ni, exchanging its current momentum state zik and gradient tracker yik.

Algorithm 1 Distributed Nesterov-Accelerated Gradient Tracking (D-NGT) for CO₂-ESS optimization. Algorithm pseudocode outlining a distributed optimization process with doubly stochastic weight matrix, initialization steps, communication, state update, Nesterov momentum, gradient tracking update, and return of optimized strategy, with mathematical equations for updates at each step.

The main mechanism of the D-NGT algorithm relies on the tight coupling of Equations 13–15.

First, Equation 13 represents the primary variable update. Unlike traditional DGD, where a node descends along its local gradient ∇fi(xik), the proposed algorithm uses the tracker variable yik. The state update is performed on the neighborhood-weighted average of the momentum variables ∑wijzjk, combining spatial consensus with optimization descent.

Second, Equation 14 implements Nesterov acceleration. The term β(xik+1-xik) serves as a momentum term. When consecutive updates align in direction, this term amplifies the step in that direction and can reduce the number of iterations compared with standard gradient tracking.

Finally, Equation 15 defines the gradient tracking mechanism. The first term ∑wijyjk performs local information mixing to reach a consensus on the gradients. The term ∇fi(zik+1)-∇fi(zik) acts as a dynamic correction factor that injects the latest local gradient variations into the network. With the initialization yi0=∇fi(zi0) and the doubly stochastic nature of W (Assumption 1), summing Equation 15 over all nodes yields:

∑i=1Nyik+1=∑i=1N∇fi(zik+1)(17)

Equation 17 shows that the spatial average of the local trackers equals the global gradient average at iteration k. This property reduces the persistent gradient bias caused by data heterogeneity among different CO₂-ESS units and supports convergence to the global optimum under the stated assumptions.

4 Distributed accelerated gradient tracking with power balance constraints

While Section 3 established the Nesterov-accelerated gradient tracking logic, it was based on an unconstrained consensus framework. In a distributed smart grid, total power generation and consumption should be balanced over the scheduling interval. Therefore, the collaborative optimization of CO₂-ESS cannot be treated as an unconstrained problem; it should include a global power balance constraint (Equations 8–11).

4.1 Formulation of the constrained optimization problem

Building upon the local state-space and cost models defined in Section 3, let Xi⊆ℝd denote the local feasible set for node i. This convex set Xi inherently encapsulates the discrete-time thermodynamic dynamics and the local operational bounds of the i-th CO₂-ESS unit.

To abstract the physical power equilibrium into a mathematical matrix form suitable for distributed algorithmic design, the globally constrained optimization problem for the distributed smart grid is formulated as follows:

minx1,…,xN ∑i=1Nfi(xi)(18)

s.t. ∑i=1N(Aixi-bi)=0(19)

xi∈Xi, ∀i∈{1,2,…,N}(20)

The symbols and matrices used in the global coupling constraint Equation 19 are defined as follows: The global coupling constraint Equation 19 is now definedcomponent by component, so that every entry of A_*i*, its sign convention, the conversion efficiencies, the renewable-curtailment term, and the (deliberately omitted) network-loss / line-flow terms are all explicit.

4.1.1 Decision vector and per-step extraction

Recall from Section 3-A that xi=vec({Pi,tc,Pi,te,Pi,tcurt}t=1T)∈ℝd, d = 3T. For each time step t the t-th block of Ai∈ℝT×d selects the three controllable components of node i :

Aixi=[Pi,te-Pi,tc-Pi,tcurt]t=1T∈ℝT,(21)

i.e., the t-th row of A_*i* is (0, …, 0, −1, +1, −1, 0, …, 0) with the non-zero entries placed at the positions of Pi,tc (coefficient −1), Pi,te (coefficient +1), and Pi,tcurt (coefficient −1). The sign convention is therefore: a positive entry of A_*i*x_*i* is active power injected into the grid by storage (turbine generation), whereas a negative entry is power withdrawn from the grid (compressor consumption or renewable energy that is curtailed rather than used).

4.1.2 Right-hand side and renewable curtailment

The local right-hand side is

bi=[Pi,tload-Pi,tRES]t=1T∈ℝT,(22)

i.e., the net local demand that must be served by the storage fleet after the locally available renewable generation Pi,tRES has been used. With Equation 21 and Equation 22, the global balance ∑i(Aixi-bi)=0 expands to

∑i=1N∑t=1T(Pi,te-Pi,tc-Pi,tcurt)=∑i=1N∑t=1T(Pi,tload-Pi,tRES),(23)

which states that the total net storage injection equals the total net load (load minus renewable availability) at every time step. Curtailing renewable power (Pi,tcurt>0) reduces the available injection and is penalized through the quadratic term γi(Pi,tcurt)2 in f_*i* (Section 3-A), so the optimizer resorts to curtailment only when it is cheaper than the alternative dispatch.

4.1.3 Conversion efficiencies

The thermo-mechanical conversion efficiencies ηmc and ηme do not appear inside A_*i* (the balance is written in gross electrical power), but they enter the problem in two consistent places: (i) they set the cost coefficients cic,cie of f_*i* (a unit with lower round-trip efficiency incurs a larger effective $/kWh cost), and (ii) they govern the state recursion Equation 6, which converts the electrical set-points Pi,tc,Pi,te into mass/pressure/thermal trajectories and therefore bounds how much energy can be moved per interval. The balance Equation 23 holds on the electrical bus; the internal conversion losses are accounted for by the dynamics and by the cost, not by an extra loss term in the electrical equality.

4.1.4 Network imports, electrical losses, and line-flow constraints

The model adopts a lossless, aggregated-balancing-area assumption: the N prosumer nodes are assumed to be electrically tightly coupled (or the dispatch is performed for a single control area), so that a single global power-balance equality suffices and Ohmic (I²R) losses, transformer losses, inter-area imports/exports, and individual nodal-voltage or line-flow (DC/AC OPF) limits are not represented. These are deliberate, explicitly stated simplifications rather than omissions: adding them would turn the single coupling row into a set of affine line-flow constraints (extra rows of A_*i*) together with local box limits on phase-angle or line-power variables inside Xi; because both are affine, they preserve the convexity and the primal-dual structure of C-DNGT and can be incorporated without changing the algorithm. We flag the hierarchical inclusion of AC-OPF line-flow constraints as a concrete extension.

Equation 19 is the global power balance constraint.

From a local perspective, (A_*i*x_*i*−b_*i*) is the local power mismatch at node i

(storage net injection minus the net local demand). From a global perspective, ∑i(Aixi-bi)=0 enforces that, at every t∈T, the power discharged by all turbines, minus the power consumed by all compressors and the curtailed renewable energy, balances the net load demand across the grid [see the expanded form Equation 23].

The main challenge in solving Equations 18–20 is that Equation 19 couples all decision variables x₁, …, x_*N* while each node only knows its own Ai,bi,Xi. Because there is no centralized dispatch center, no single node possesses the global sum ∑iAixi or ∑ibi. We therefore design a primal-dual distributed algorithm in which nodes reach consensus on the dual (price) variable and on the global mismatch tracker, satisfying the coupling using only localized peer-to-peer communication.

4.2 Primal-dual Lagrangian framework and algorithm derivation

To handle the global coupling constraint Equation 19 in a fully decentralized manner, we employ Lagrangian duality. Since the plain (un-augmented) Lagrangian would leave the primal subproblem only convex–not strongly convex–in the constraint subspace, we instead form theaugmented Lagrangian (method of multipliers) by adding a quadratic penalty on the global power mismatch; this makes the name “augmented Lagrangian” accurate and, as discussed in Section 5, provides the coercivity needed for the convergence results. The augmented Lagrangian function for the global problem Equations 18, 19 is defined by introducing a global dual variable λ ∈ ℝ^*T* (physically interpretable as the shadow price of electricity or the marginal grid penalty) and an augmentation penalty ρ > 0:

Lρ(x,λ)=∑i=1Nfi(xi)+λT∑i=1N(Aixi-bi)+ρ2∥∑i=1N(Aixi-bi)∥2(24)

In the absence of a central coordinator, the global dual variable λ cannot be directly computed or broadcast. Therefore, we assign a local dual variable estimate λi∈ℝT to each CO₂-ESS node i. The distributed optimization goal transforms into finding the saddle point of the Lagrangian function, enforcing consensus on the dual variables (λ1=⋯=λN=λ*) while driving the global constraint violation to zero.

4.2.1 Strong-convexity regularization

We further add a small Tikhonov / strong-convexity regularization to each local cost,

f~i(xi)=fi(xi)+μ2∥xi-xiref∥2, μ>0,(25)

with a fixed reference xiref (taken as the cold-start zero schedule) and μ chosen small (e.g., μ = 10⁻³). Throughout the remainder of the study, the symbols f_*i* and ∇f_*i* are understood to refer to the regularized counterparts f~i and ∇f~i. Because Xi is compact and f_*i* is Lipschitz on Xi, the regularized optimum (x^μ, λ^μ) satisfies ∥xμ-x∥=O(μ) and |F(xμ)-F(x)|=O(μ2); the regularization bias is therefore negligible relative to the dispatch cost, while it guarantees that the global objective is μ-strongly convex so that the strong-convexity premise of Theorems 1 and 2 is satisfied exactly (not merely on the constraint subspace). The interplay between μ and the augmentation penalty ρ in the convergence guarantees is discussed in Section 5.

To achieve this, we integrate the Nesterov momentum technique into the primal descent step, and deploy the dynamic gradient tracking mechanism to trace the global constraint violation in the dual ascent step. We propose the Constrained Distributed Nesterov-Accelerated Primal-Dual Gradient Tracking (C-DNGT) algorithm. The coupled iterative update rules for node i at iteration k are derived as follows:

1. Primal Projected Descent Step: Node i updates its local decision variable x_*i* by descending along the gradient of the local augmented Lagrangian Equation 24 with respect to the Nesterov extrapolated state zik. Because the quadratic penalty term depends on the global mismatch ∑j(Ajxj-bj), node i evaluates its local part of the augmentation gradient using the distributed mismatch tracker yik≈∑j(Ajxjk-bj) established below; this requires no extra communication beyond the tracker exchange already present in the algorithm. To maintain physical feasibility (Equation 20), the result is projected onto the local feasible set Xi:

xik+1=PXi[zik-α(∇fi(zik)+AiTλik+ρAiTyik)](26)

where PXi[·] denotes the Euclidean projection operator onto the convex set Xi, and α > 0 is the primal step size.

4.2.1.1 Computational implementation of the projection

We now make the evaluation of PXi fully explicit. As established in Section 3-A, the feasible set Xi is the polytope Xi={xi∈ℝd:c_i≤Eixi≤c¯i(si,0)} obtained by lifting the affine CO₂-ESS dynamics and the physical box limits (7). The projection required at every primal step (19) is therefore the solution of the strictly convex quadratic program

PXi(v)=argminxi∈Xi12∥xi-v∥22,(27)

which we solve with a standard interior-point / active-set convex QP solver. Because Xi has only d = 3T variables and is decoupled across the N nodes, the projection is performed independently and in parallel at each node at every iteration, and its cost is negligible compared with the gradient evaluation; no cross-node communication is needed for the projection itself. Furthermore, the block-bidiagonal structure of E_*i* (induced by the first-order recurrence of the thermo-mechanical states) makes the resulting QP sparse, so the per-iteration projection remains inexpensive even for long horizons T. In our implementation (Python 3.9), the same convex QP routine is reused for the local subproblem, so the projection does not introduce any additional solver dependency. We emphasize that the projection is applied simultaneously over the whole scheduling horizon (i.e., on the condensed vector x_*i*), rather than sequentially time step by time step, because the lifted set Xi already encodes the full time-coupled dynamics in a single convex constraint.

The iterative procedure of the proposed constraint-handling algorithm is summarized in Algorithm 2.

Algorithm 2 Constrained Distributed Nesterov-Accelerated Primal-Dual Gradient Tracking (C-DNGT). Algorithmic pseudocode outlines a distributed optimization method using a doubly stochastic matrix, primal and dual variables, Nesterov momentum, and iterative communication, update, and consensus steps across nodes, referencing Equations 26, 28, 29, and 30 for calculations.

2. Primal Nesterov Extrapolation Step: To accelerate convergence, the Nesterov momentum is applied to generate the extrapolated state zik+1, which will be used for the gradient evaluation in the subsequent iteration:

zik+1=xik+1+β(xik+1-xik)(28)

where β ∈ (0, 1) is the momentum coefficient.

3. Constraint Violation Tracking Step: Instead of tracking the gradient of the objective function (as in Section 3), the tracking variable yik∈ℝT is now repurposed to dynamically track the average global power mismatch. The tracker update fuses the neighbors' information with the local variation in the constraint state:

yik+1=∑j∈Ni∪{i}wijyjk+Ai(xik+1-xik)(29)

By initializing yi0=Aixi0-bi and utilizing the doubly stochastic weight matrix W, the identity ∑i=1Nyik=∑i=1N(Aixik-bi) holds at any iteration k. Thus, yik provides distributed information about the global power imbalance.

4. Dual Consensus and Ascent Step: Finally, node i updates its local dual variable λ_*i*. This update combines a consensus step on neighbors' dual variables with a dual ascent step along the direction of the tracked global violation yik+1:

λik+1=∑j∈Ni∪{i}wijλjk+ηyik+1(30)

where η > 0 is the dual ascent step size.

To intuitively illustrate the architectural differences and the theoretical evolution of the proposed method, Figure 2 presents a block diagram comparing the foundational unconstrained D-NGT and the proposed constrained C-DNGT algorithms.

As shown in Figure 2a, the standard unconstrained method operates sequentially, using the tracker y_*i* to estimate the global objective gradient. Conversely, Figure 2b illustrates the primal-dual architecture used for distributed smart grids. To handle the global power balance (∑(A_*i*x_*i*−b_*i*) = 0) and localized physical limits (xi∈Xi), the proposed framework repurposes the tracker y_*i* to monitor the physical power mismatch. It also introduces a parallel dual consensus module (λ_*i*) and a primal projection operator (PXi). This coupled architecture links constraint tracking with economic scheduling of CO₂-ESS units.

5 Convergence analysis

In this section, we provide a detailed theoretical convergence analysis of the proposed Distributed Nesterov-Accelerated Gradient Tracking (D-NGT) algorithm. To establish foundational theoretical bounds, we first analyze the algorithm in the unconstrained setting, laying the mathematical groundwork for the constrained power balance analysis (Equations 31, 32b, 33).

5.1 Convergence regime and the strong-convexity premise

As noted in Section 3-A and made explicit in Section 4, theunregularized local cost f_*i* is affine in the compressor/expander powers and quadratic only in the curtailment, so the global objective is convex butnot strongly convex in the complete decision vector. Theorems 1 and 2 below are proved under a μ-strong-convexity assumption, which is now satisfied exactly by the regularized objective Equation 25 (with the same symbol f_*i* used throughout). The price paid is a regularization bias of order O(μ) in the dispatch and O(μ2) in the cost, which is driven to zero as μ → 0 and is negligible for the small μ used in the simulations. We stress that this does not weaken the practical claim: the algorithm solves the intended economic-dispatch problem to within a controllable, quantifiable tolerance. Moreover, because the local sets Xi and the coupling Equation 19 are polyhedral and we employ an augmented Lagrangian Equation 24, the Luo–Tseng error-bound principle Luo and Tseng (1993) gives global linear convergence of the primal-dual residual for the underlying proximal primal-dual method even on the original convex (non-strongly-convex) problem; the strong-convexity regularization therefore provides a clean sufficient condition, while the error bound provides an independent, robustness-oriented justification. Both mechanisms are reflected in the step-size conditions of Theorem 2.

5.2 Convergence analysis of the unconstrained D-NGT algorithm

To enable the matrix-form derivation, we aggregate the local variables of all N agents into network-level matrices. Let x^*k*, z^*k*, y^*k* ∈ ℝ^*N*×*d* denote the stacked matrices of the primal states, Nesterov momentum states, and gradient trackers at iteration k, respectively. The aggregate gradient matrix is denoted as ∇f(zk)=[∇f1(z1k),…,∇fN(zNk)]T.

We define the network-averaged vectors as:

x̄k=1N1Txk, z̄k=1N1Tzk, ȳk=1N1Tyk(31)

By multiplying the D-NGT update rules (Equations 13–15) by 1N1T and utilizing the doubly stochastic property of W (1^*T*W = 1^*T*), the averaged system dynamics evolve as:

x̄k+1=z̄k-αȳk(32a)

z̄k+1=x̄k+1+β(x̄k+1-x̄k)(32b)

ȳk+1=ȳk+1N1T(∇f(zk+1)-∇f(zk))(32c)

Because yi0=∇fi(zi0), telescoping Equation 32c gives ȳk=1N∑i=1N∇fi(zik) for all k ≥ 0.

To quantify the consensus and tracking errors, we introduce the network projection matrix Π=I-1N11T. Based on Assumption 2, there exists a spectral norm bound σ = ||ΠW|| < 1.

Different from standard gradient tracking, the presence of Nesterov momentum β(x^*k*−x^*k*−1) intertwines the current state with the historical state. We establish the following lemmas to bound the optimization descent and the network consensus errors.

Lemma 1 (Nesterov Optimization Descent): Under the μ-strong convexity and L-smoothness of the global objective F(x), the network-averaged trajectory satisfies:

∥z̄k+1-x∥2≤(1-τ)∥z̄k-x∥2+2β2τ∥x̄k+1-x̄k∥2 +αL2μN∥Πzk∥2(33)

where τ=O(αμ) is the contraction rate driven by the primal step size α, and x^* is the exact global optimum.

Proof: The proof follows by substituting ȳ^*k* into Equation 32a. Note that ȳ^*k* is not the exact gradient ∇F(z̄k), but contains a gradient tracking perturbation. Based on the L-smoothness, the perturbation is bounded by ∥ȳk-∇F(z̄k)∥≤LN∥Πzk∥. Applying Young's inequality and the standard Nesterov convergence properties yields the stated bound.

Lemma 2 (Consensus and Tracking Contraction): For the Nesterov momentum state z^*k* and the gradient tracker y^*k*, the consensus errors satisfy the following coupled inequalities:

∥Πzk+1∥≤(1+β)σ∥Πzk∥+βσ∥Πzk-1∥ +α(1+β)∥Πyk∥+αβ∥Πyk-1∥(34)

∥Πyk+1∥≤σ∥Πyk∥+L∥zk+1-zk∥(35)

Proof: Multiplying Π to the un-averaged D-NGT update x^*k*+1 = Wz^*k*−αy^*k* gives ||Πx^*k*+1|| ≤ σ||Πz^*k*||+α||Πy^*k*||. Substituting this into the momentum definition Πz^*k*+1 = Πx^*k*+1+β(Πx^*k*+1−Πx^*k*) directly yields Equation 34. Similarly, applying Π to the tracking equation gives ||Πy^*k*+1|| ≤ σ||Πy^*k*||+||Π(∇f(z^*k*+1)−∇f(z^*k*))||. Utilizing the L-smoothness of f_*i* directly leads to Equation 35.

Inspired by recent advances in gradient-tracking analyses, we instead formulate a combined Lyapunov function tailored to the Nesterov acceleration scheme, rather than constructing a complex transition matrix.

Theorem 1 (Linear Convergence of Unconstrained D-NGT): Suppose Assumptions 1 and 2 hold. Let the combined Lyapunov sequence be defined as:

Vk=∥z̄k-x*∥2+c1∥x̄k-x̄k-1∥2 +c2∥Πzk∥2+c3∥Πzk-1∥2+c4∥Πyk∥2(36)

If the step size satisfies α≤O((1-σ)2L) and the momentum satisfies β<1-σ2, there exist positive constants c₁, c₂, c₃, c₄ > 0 and a contraction factor ρ ∈ (0, 1) such that:

Vk+1≤ρVk(37)

Consequently, the sequences ∥x̄k-x∥2, ||Πx^k||², and ||Πy^k*||² all converge to 0 at a global linear rate O(log(1/ϵ)).

Proof: By squaring the consensus and tracking bounds Equations 34 and 35, and utilizing the general inequality (a+b+c)² ≤ 3a²+3b²+3c², we can obtain upper bounds for ||Πz^*k*+1||² and ||Πy^*k*+1||². Combining these bounded terms with Lemma 1 gives the descent relationship (Equations 36–38):

Vk+1≤(1-τ+O(αL2c4))∥z̄k-x*∥2 +(σ2(1+β)2+O(α))c2∥Πzk∥2 +(σ2+O(α))c4∥Πyk∥2(38)

By selecting the parameter β to restrict the momentum overshoot, and choosing a sufficiently small step size α, the coefficients of the positive components in Vk+1 can be bounded below 1. Specifically, matching the coefficients yields Vk+1≤max{1-τ,σ2+δ}Vk, where δ is a small positive margin controlled by α. Thus, setting ρ = max{1−τ, σ²+δ} < 1 establishes linear convergence.

5.3 Convergence analysis of the constrained C-DNGT algorithm

Building on the bounds established in the unconstrained scenario, we now extend the theoretical analysis to the constrained C-DNGT framework. The global power balance constraint requires the co-evolution of primal and dual variables, thereby transforming the optimization problem into a saddle-point search.

Let A = diag(A₁, …, A_*N*) represent the block-diagonal constraint matrix, and b=[b1T,…,bNT]T be the concatenated rigid load vector. Let λ^*k* ∈ ℝ^*N*×*T* be the matrix of dual variables. Recalling that the tracker y^*k* in Section 4 is repurposed to track the physical constraint violation, the network-level matrix dynamics of C-DNGT can be compactly expressed as:

xk+1=PX[zk-α(∇f(zk)+ATλk+ρATyk)](39a)

zk+1=xk+1+β(xk+1-xk)(39b)

yk+1=Wyk+A(xk+1-xk)(39c)

λk+1=Wλk+ηyk+1(39d)

where PX is the joint Euclidean projection onto the feasible set X=X1×⋯×XN. Let (x^*, λ^*) denote the exact optimal saddle point satisfying the KKT conditions: 〈x-x,∇f(x)+AT1λ*〉≥0,∀x∈X, and 1^*T*(Ax^*−b) = 0.

5.3.1 Saddle-point residual expansion

To capture the interaction among the Nesterov momentum, the dual consensus, and the power mismatch tracker, we establish the following lemma. This lemma bounds the cross-coupling terms arising from the primal-dual updates.

Lemma 3 (Primal-Dual Cross-Coupling Inequality): Utilizing the non-expansiveness property of the projection operator, i.e., ∥PX(u)-PX(v)∥2≤∥u-v∥2, the residual energy of the primal-dual system from the optimal saddle point is bounded by the following expansion:

12α∥xk+1-1x∥2+12η∥λk+1-1λ∥2≤12α∥zk-1x∥2-μ2∥zk-1x∥2︸Strongly Convex Primal Descent+12η∥Wλk-1λ∥2︸Dual Consensus Contraction-〈xk+1-1x,∇f(zk)-∇f(1x)〉︸Gradient Tracking Error Bound-〈xk+1-1x,AT(λk-1λ+ρyk)〉+〈λk+1-1λ,yk+1〉︸Saddle-Point Coupling (Power Balance Restitution)+L2α2∥zk-xk∥2+η2∥Πyk+1∥2+O(β2)∥xk+1-xk∥2(40)

Proof: We derive the inequality term by term from the primal and dual C-DNGT matrix updates Equations 39a–39d.

5.3.2 Step 1 (primal descent)

Let 1x^* denote the stacked matrix of the optimal primal vector x^* replicated N times. At the KKT saddle point (x^*, λ^*), the stationarity condition reads: for every agent i,

∇fi(xi)+Ai⊤λ∈-NXi(xi*),(41)

where NXi(·) is the normal cone and the physical mismatch satisfies ∑i(Aixi*-bi)=0. Multiplying Equation 41 on both sides by α and adding x^* yields the equivalent fixed-point characterization x=PX[x-α(∇f(1x)+A⊤1λ)]. Applying the non-expansiveness of the Euclidean projection (∥PX(u)-PX(v)∥2≤∥u-v∥2) to the primal update Equation 39a and the fixed-point relation, we obtain

∥xk+1−1x∥2≤∥zk−1x−α(∇f(zk)−∇f(1x) +A⊤(λk−1λ)+ρA⊤yk)∥2.(42)

Expanding the square in Equation 42 exactly (no term is dropped) gives

∥xk+1-1x∥2≤∥zk-1x∥2-2α〈zk-1x*,vk〉 +α2∥vk∥2,(43)

where we abbreviate v^*k*: = ∇f(z^*k*)−∇f(1x^*)+A^⊤(λ^*k*−1λ^*)+ρA^⊤y^*k*. To obtain the saddle-point form used in Equation 40, which involves x^*k*+1 rather than z^*k* inside the coupling inner products, we split the inner product and bound the discrepancy explicitly (nothing is hidden): 〈z^*k*−1x^*, v^*k*〉 = 〈x^*k*+1−1x^*, v^*k*〉+〈z^*k*−x^*k*+1, v^*k*〉. The Nesterov momentum gives z^*k* = x^*k*+β(x^*k*−x^*k*−1). Since xk+1=PX[zk-αvk] and xk∈X (the projection keeps every iterate feasible),

∥zk-xk+1∥≤∥zk-xk∥+∥xk-xk+1∥≤β∥xk-xk-1∥ +(β∥xk-xk-1∥+α∥vk∥).(44)

By Cauchy–Schwarz and the bound on ‖zk−xk+1‖, 〈zk−xk+1,vk〉≥−‖zk−xk+1‖ ‖vk‖≥−(α‖vk‖2+2β‖xk−xk−1‖ ‖vk‖). Therefore −2α〈zk−xk+1,vk〉≤2α2‖vk‖2+4αβ‖xk−xk−1‖ ‖vk‖, and substituting into Equation 43 yields

∥xk+1-1x∥2≤∥zk-1x∥2-2α〈xk+1-1x*,vk〉 +3α2∥vk∥2+4αβ∥xk-xk-1∥∥vk∥.(45)

The last two terms are explicit non-negative residuals (the momentum term is no longer hidden): by Young's inequality, 4αβ||x^*k*−x^*k*−1||||v^*k*|| ≤ 2β²||x^*k*−x^*k*−1||²+2α²||v^*k*||², so the total quadratic remainder is bounded by 5α²||v^*k*||²+2β²||x^*k*−x^*k*−1||²; this is absorbed into the Lyapunov weights c₅, …, c₈ of Theorem 2 and supplies the O(β2)∥xk+1-xk∥2 remainder in Equation 40. Expanding v^*k* recovers exactly the three coupling inner products required by Equation 40.

5.3.3 Step 2 (Nesterov momentum error)

Nesterov extrapolation gives z^*k* = x^*k*+β(x^*k*−x^*k*−1), hence

∥zk-1x∥2≤(1+ϵ)∥xk-1x∥2+(1+1ϵ)β2∥xk-xk-1∥2,(46)

for any ϵ > 0 by Young's inequality. Setting ϵ = αμ (which is O(αμ) per the step-size conditions of Theorem 2) optimally balances the primal contraction and the momentum-induced difference term.

5.3.4 Step 3 (Gradient tracking and strong convexity)

Because every f~i is L-smooth and μ-strongly convex, and z^*k* enters only through the gradient, we bound

〈zk−1x,∇f(zk)−∇f(1x)〉≥μ∥zk−1x∥2+1L+μ∥∇f(zk) −∇f(1x)∥2≥(Baillon-Haddad)μ∥zk−1x*∥2,(47)

Because Step 1 rewrites the primal coupling inner product with x^*k*+1 in place of z^*k*, the estimate above (which uses 〈z^*k*−1x^*, ∇f(z^*k*)−∇f(1x^*)〉) is applied to 〈x^*k*+1−1x^*, ·〉 at the cost of an extra term 〈x^*k*+1−z^*k*, ∇f(z^*k*)−∇f(1x^*)〉, bounded by ||z^*k*−x^*k*+1||L||z^*k*−1x^*|| ≤ (α||v^*k*||+2β||x^*k*−x^*k*−1||)L||z^*k*−1x^*|| via Equation 44 and L-smoothness ; this higher-order term is absorbed into the quadratic remainder of Equation 40.

and, by L-smoothness together with z^*k*−x^*k* = β(x^*k*−x^*k*−1),

∥∇f(zk)-∇f(1x)∥≤L∥zk-xk∥+L∥xk-1x∥≤Lβ∥xk -xk-1∥+L∥xk-1x*∥.(48)

5.3.5 Step 4 (dual consensus and ascent)

From the dual update λ^*k*+1 = Wλ^*k*+ηy^*k*+1 we separate the average and the consensus error:

λ̄k+1=λ̄k+ηȳk+1,(49)

Πλk+1=ΠWΠλk+ηΠyk+1,(50)

where Π=I-1N11⊤, and we used W1 = 1. Squaring and using ||ΠW|| = σ < 1 (Assumption 2) yields

∥λ̄k+1-λ∥2=∥λ̄k-λ+ηȳk+1∥2 ≤∥λ̄k-λ∥2+2η〈λ̄k-λ,ȳk+1〉+η2∥ȳk+1∥2,(51)

∥Πλk+1∥2≤σ2∥Πλk∥2+2η21-σ2∥Πyk+1∥2.(52)

5.3.6 Step 5 (saddle-point coupling cancellation—The key algebraic step)

The coupled term 〈x^*k*+1−1x^*, A^⊤(λ^*k*−1λ^*)〉+〈λ^*k*+1−1λ^*, y^*k*+1〉 appearing in Equation 40 is handled as follows. Summing the inner product of λ^*k*+1−1λ^* with the dual update Equation 39d and rearranging gives

〈λk+1− 1λ,yk+1〉=1η∥λk+1−1λ∥2−1η〈λk+1−1λ,Wλk −1λ〉.(53)

For the primal augmented-mismatch term we introduce the global power-mismatch vector

sk:=∑i=1N(Aixik-bi)∈ℝT,(54)

which, by the tracker-sum invariant (telescoping of Equation 39c), equals ∑i=1Nyik. At the exact saddle point s^* = 0. Decompose the tracker into its network-average and consensus-error parts, yk=1N1⊗sk+Πyk, where 1N1⊗sk is the NT-dimensional vector whose i-th agent block equals s^*k*/N. Using the adjoint identity 〈u, A^⊤y^*k*〉 = 〈Au, y^*k*〉 with u = x^*k*+1−1x^*,

〈xk+1−1x,A⊤yk〉 =〈A(xk+1−1x),1N1⊗sk〉+〈A(xk+1 −1x),Πyk〉 =1N〈∑i=1NAi(xik+1−xi),sk〉+ℰΠ =1N〈sk+1,sk〉+ℰΠ,(55)

since ∑i=1NAi(xik+1-xi)=sk+1-s=sk+1. Both s^*k*+1 and s^*k* are T-dimensional, so the inner product and the 1/N prefactor are now exact, with no broadcast-notation ambiguity. The consensus-error remainder EΠ=〈A(xk+1-1x),Πyk〉 satisfies |EΠ|≤∥A∥∥xk+1-1x∥∥Πyk∥ by Cauchy–Schwarz.

Now combine Equation 53 with the primal coupling terms. Let λ̄k=1N1⊤λk and ȳk=1N∑i=1Nyik=sk/N denote the agent averages (both T-dimensional, with s^*k* from Equation 54). Using 1^⊤A(x^*k*+1−1x^*) = s^*k*+1 (sum over agents) and s^*k*+1 = Nȳ^*k*+1, the full primal-dual cross term—including the augmented ρA^⊤y^*k* part of Equation 45—contracts to

-〈xk+1-1x,A⊤(λk-1λ)〉+〈λk+1-1λ,yk+1〉-ρ〈xk+1-1x,A⊤yk〉≤N〈λ̄k+1-λ,ȳk+1〉-N〈λ̄k-λ,ȳk+1〉-ρN〈sk+1,sk〉+1η∥Πλk+1∥2 +(∥A∥∥Πλk∥+ρ∥A∥∥Πyk∥)∥xk+1-1x*∥+ρ|EΠ|+η∥Πyk+1∥2.(56)

The first two average terms telescope through the dual-average relation λ̄k+1-λ̄k=ηȳk+1 (Equation 49); the third term -ρN〈sk+1,sk〉 is the explicit power-balance restitution that forces s^*k*→0; the remaining terms are consensus/tracker errors, bounded by σ < 1 and absorbed into the Lyapunov weights c₅, …, c₈ of Theorem 2. The first two terms telescope; after applying Young's inequality to the cross norms the full inequality Equation 40 is recovered with the explicit remainder

O(β2)∥xk+1-xk∥2=L2β21-σ2∥xk-xk-1∥2,(57)

which follows from Equation 48 when we bound ∥xk+1-xk∥≤2β1-σ∥xk-xk-1∥+α1-σ(Lyapunov cross-terms) by combining the primal and momentum update recursions.

5.3.7 Step 6 (regularity and strong-monotonicity conditions)

The proof rests on three structural properties, all satisfied by the problem and the algorithm:

• L-smoothness and μ-strong convexity: The regularized objective f~i (Equation 25) satisfies μI≼∇2f~i≼LI; this provides the contraction in Step 3 and enters the primal step-size bound α ≤ μ/(4L²) in Theorem 2.
• Constraint regularity: The coupling matrix A^⊤A has rank N (full technical rank) because each A_*i* extracts the three controllable components per time step and their rows span the T independent power-balance equations; the augmented Hessian A^⊤A is therefore non-degenerate on the constraint subspace, which is the regularity condition required for the linear convergence of projected primal-dual methods on polyhedral constraints Luo and Tseng (1993).
• Strong monotonicity of the KKT system: Combining the μ-strong convexity of f~ in x with the non-degenerate coupling A and the dual consensus contraction σ < 1 makes the full primal-dual operator strongly monotone in the Lyapunov metric of Theorem 2; the coupled system is therefore a contraction, yielding the linear rate ρ_*pd* ∈ (0, 1).

Collecting Equations 45–57 and balancing the constants c₅, …, c₈ as in Theorem 2 yields the bound stated in Equation 40.

We can now construct the ultimate augmented Lyapunov function for the constrained CO₂-ESS network.

Theorem 2 (Global Convergence of Constrained C-DNGT): Define the augmented primal-dual Lyapunov function as:

Wk=Vunconstrainedk+c5∥λ̄k-λ*∥2 +c6∥Πλk∥2+c7∥Πyk∥2+c8∥xk-xk-1∥2(58)

where Vunconstrainedk is defined in Theorem 1. Suppose the primal step size α, the dual step size η, and the Nesterov momentum β satisfy the following coupling conditions:

0<α≤μ4L2, 0<η≤(1-σ)24∥ATA∥, β≤12αμ(59)

Then, there exists a system-dependent contraction parameter ρ_*pd* ∈ (0, 1) such that Wk+1≤ρpdWk.

Consequently, the proposed C-DNGT algorithm converges linearly to the KKT saddle point under the stated assumptions. This implies that the distributed CO₂-ESS units reach an economic consensus while driving the global power mismatch limk→∞∑i=1N(Aixik-bi)=0, thereby satisfying the supply-demand balance constraint of the smart grid model.

6 Experimental results and system performance analysis

In this section, numerical simulations are conducted to evaluate the convergence behavior and constraint-handling performance of the proposed Constrained Distributed Nesterov-Accelerated Gradient Tracking (C-DNGT) algorithm for CO₂-ESS networks. All numerical experiments and algorithmic simulations are run on a workstation with an Intel Core i9-13900K processor (Intel Corporation, Santa Clara, CA, USA) (3.0 GHz base clock) and 64 GB of RAM. The proposed distributed algorithms are implemented in Python 3.9 (Python Software Foundation, Wilmington, DE, USA).

The distributed smart grid is mathematically modeled as a multi-agent network comprising N distinct prosumer nodes. In the baseline experiments, a connected communication graph G with N nodes is randomly generated using the Erdős-Rényi model with a connectivity probability of p = 0.3. To satisfy the doubly stochastic requirement (Assumption 1) for the communication weight matrix W, the Metropolis-Hastings rule is applied to assign the edge weights. The step sizes for primal descent (α), Nesterov momentum (β), and dual ascent (η) are selected based on the Lipschitz constant of the system and the convergence conditions used in the analysis.

The external grid environments, including load demands and renewable generation baselines, are constructed from load data and parameter settings as follows:

• Load Demand Data: We utilize the hourly historical electricity load datasets provided by the PJM Interconnection (a major regional transmission organization in the United States¹). The macro-level load data is proportionally scaled down to reflect the demand profiles of distributed microgrid nodes.
• CO₂-ESS Thermodynamic Parameters: The physical specifications of the CO₂-ESS, such as the specific heat capacity C_*tes*, isentropic efficiencies of the compressor ηmc and turbine ηme, and pressure boundary limits, are parameterized according to industrial literature on transcritical CO₂ cycles.

Before feeding the time-series data into the distributed optimization framework, an adaptive alignment preprocessing module is used. This module handles missing values via cubic spline interpolation and aligns the diverse sampling frequencies of the PJM load to a unified discrete time step (Δt = 1 h) over a 24-h scheduling horizon (T = 24).

6.1 Baseline verification: unconstrained D-NGT convergence

Before introducing physical grid constraints, we establish a mathematical baseline by verifying the convergence properties of the unconstrained D-NGT algorithm (Section 3). The purpose of this subsection is to isolate and validate the Nesterov-accelerated gradient-tracking engine that underpins the constrained C-DNGT: both the convergence rate and the parameter sensitivity observed here carry over to the constrained setting, as confirmed by the C-DNGT experiments in Sections 6.2–6.5. The unconstrained baseline therefore serves as a transparent building-block validation and not as a substitute for the constrained validation.

In this subsection, we define a standard unconstrained distributed consensus optimization problem. The local cost function for each node i is formulated as a strongly convex quadratic function:

fi(xi)=12ai(xi-bi)2(60)

where a_*i*>0 and b_*i* are randomly generated local configuration parameters. Because this problem is unconstrained, the centralized global optimal solution x^* can be analytically derived via the first-order optimality condition as x*=∑i=1Naibi∑i=1Nai.

To quantitatively evaluate the convergence performance and accuracy of the algorithm throughout the iterative process, we define the global optimality residual error at iteration k as the mean squared distance from the agents' current states to the true global optimum:

ek=1N∑i=1N∥xik-x*∥2(61)

To assess the algorithm's robustness and parameter sensitivity, two numerical experiments are designed:

• Impact of Network Scale (N): We fix the algorithmic hyperparameters (step size α and momentum β) and evaluate the algorithm across varying network dimensions N ∈ {10, 30, 50}. This experiment examines scalability and consensus behavior over differing graph topologies.
• Impact of Algorithmic Step Size (α): We fix the network size at N = 20 and test varying primal step sizes α. This test examines the relationship between convergence speed and stability.

The corresponding numerical results and convergence trajectories for these two baseline experiments are presented in Figures 3, 4, respectively (Equations 58–62).

The baseline experiments demonstrate the convergence behavior and parameter sensitivity of the proposed D-NGT algorithm. As shown in Figure 3, the algorithm exhibits approximately linear convergence and reduces steady-state errors across the tested network dimensions, though larger networks converge more slowly due to the increased graph diameter and the required communication hops. Furthermore, Figure 4 illustrates the trade-off associated with the algorithmic step size. A moderate step size (α = 0.05) improves convergence speed, whereas an excessively large step size (α = 0.09) leads to divergence. This result is consistent with the requirement that the step size be bounded by the smoothness properties of the objective function.

6.2 Constrained optimization in real-world scenarios

To evaluate the constrained C-DNGT algorithm under practical smart grid conditions, we formulate a constrained optimization problem Equations 18–20 using historical electricity load datasets from the PJM Interconnection. The objective is to minimize local generation costs subject to the global power balance equality constraint, which enforces supply-demand equilibrium across the distributed network. We first conduct a convergence analysis across varying network scales under this constraint.

In practical smart grid deployments, distributed CO₂-ESS units may exhibit physical heterogeneity due to variations in manufacturing standards, storage capacities, and turbomachinery efficiencies. To capture this asymmetry, we assign distinct cost coefficients to the nodes, with high-efficiency units incurring lower degradation costs and low-efficiency units facing higher operational penalties. The specific parameterization used in this experiment (and provided for full reproducibility) is:

• Cost coefficients: compressor cost cic∈[0.10,0.50] $/kWh, expander cost cie∈[0.10,0.50] $/kWh, and curtailment penalty γ_*i* ∈ [0.5, 4.5] (dimensionless quadratic coefficient), all drawn uniformly from the indicated ranges for each node with random seed 7.
• Control bounds: compressor and expander power Pi,tc,Pi,te∈[0,Pmax,i] with P_{max, *i*} ∈ [6, 12] kW (heterogeneous); curtailment cap C_{max, *i*} ∈ [2, 4] kW.
• Graph: Erdős–Rényi random graph with p = 0.3, enforced connected by a Hamiltonian path, Metropolis–Hastings edge weights (seed 11).
• Algorithm hyperparameters: α = 0.05, β = 0.1, η = 0.1, μ = 10⁻³, ρ = 6, independently verified as stable by the sensitivity analysis in Figure 5.
• Load data: PJM Interconnection hourly load, scaled proportionally to microgrid levels and processed via cubic-spline imputation to a T = 24 horizon; renewable profiles are smooth diurnal shapes with inter-node heterogeneity (seed 7).

The objective of this experiment is to assess whether the proposed constrained algorithm can converge, respect physical grid limits, and allocate power economically under heterogeneous conditions. The performance comparison under these asymmetric conditions is shown in Figure 6.

The numerical results indicate that the algorithm can coordinate asymmetric distributed CO₂-ESS networks under the tested settings. The three subplots in Figure 6 are discussed as follows:

First, Figure 6a shows the evolution of the optimality gap. While the homogeneous network converges steadily and linearly, parameter differences in the heterogeneous scenario create a more ill-conditioned optimization landscape. This asymmetry causes the optimization trajectory to overshoot transiently, manifesting as a “V-shaped” oscillation due to Nesterov momentum accumulation. The C-DNGT framework then dampens the oscillation and resumes descent toward the optimum.

Second, despite the transient oscillations observed in the objective space, the global power mismatch is driven toward zero, as shown in Figure 6b. The gradient tracking mechanism helps absorb local overshoots, and the global power equilibrium constraint (Pdemand=∑i=1NPik) is approached at an approximately linear rate under the asymmetric conditions.

Finally, Figure 6c visualizes the steady-state power allocation. Unlike the homogeneous baseline, which uniformly distributes the total load across all nodes, the C-DNGT framework allocates larger power shares to high-efficiency storage units and lower shares to units with higher operational penalties. This allocation is consistent with the cost-minimizing structure of the formulated dispatch problem.

6.3 Comparative analysis with benchmark algorithms

This subsection compares the unconstrained D-NGT acceleration engine (Section 3) against standard distributed algorithms—Distributed Gradient Descent (DGD), EXTRA, and Gradient Tracking (GT) agents. The comparison is designed to isolate the effect of Nesterov momentum combined with dynamic gradient tracking, decoupled from the constraint-handling machinery of C-DNGT. The constrained C-DNGT algorithm, which inherits the same acceleration engine, is evaluated separately on the ED problem with the global power-balance constraint in Sections 6.2, 6.4, and 6.5 (Figures 5–7).

6.3.1 Communication-noise model

For the Monte Carlo experiments in Figures 8c, d, we model the communication channel as an additive white Gaussian process. At each iteration k, the message transmitted from node j to node i (the gradient-tracker value yjk and the dual variable λjk) is corrupted by independent Gaussian noise zijk~N(0,σcomm2) with σcomm=10-3. The corrupted observation is ỹjk=yjk+zijk (and analogously for λ~jk), and the receiver uses ỹjk in place of yjk in all subsequent computations.

6.3.2 Message size and total communication burden

Each local iteration, every agent i transmits the tracker yik∈ℝT (T real numbers) and, in the constrained mode, the dual variable λik∈ℝT (T real numbers), for a total message size of 2T floating-point values per directed link. For the N = 50 unconstrained benchmark, each agent transmits only yik∈ℝd (d dimensional, d~T). The total communication burden reported in Figure 8f is defined as

Total burden=(message size per link)×(number of edges) ×(iterations to tolerance),(62)

Cumulated over 50 Monte Carlo trials for the noise-affected case. The message size is held constant across methods; the burden differences therefore reflect only the difference in iteration count, which is the dominant scaling factor.

First, Figures 8a, b show the baseline convergence performance. As shown, DGD exhibits a persistent steady-state error near 10⁻² and fails to achieve exact consensus (∥xk-x̄∥2↛0) because its gradient iterations are not corrected. Conversely, both EXTRA and standard GT reduce this bias. With Nesterov momentum, the proposed D-NGT algorithm shows a steeper descent and reaches an optimality gap of approximately 10⁻¹⁵ within 150 iterations in this experiment.

Second, practical smart grid deployments may be affected by communication uncertainties. Figures 8c, d introduce communication noise and report results from 50 Monte Carlo trials. Figure 8c shows that the proposed D-NGT requires fewer iterations to reach a 10⁻⁶ tolerance and exhibits lower variance in this setting. Figure 8d further shows that D-NGT maintains a lower final optimality gap than the compared methods under the tested noise level.

Finally, we assess scalability and communication cost. Figure 8e shows that as the network topology scales (N = 10, 30, 50), D-NGT retains a convergence advantage over EXTRA and GT in the tested cases. As shown in Figure 8f, the reduced iteration count also lowers the total communication burden. This suggests that the proposed framework can reduce communication requirements for large-scale CO₂-ESS coordination.

6.3.3 Scope of the constrained-vs-constrained comparison

The benchmark above isolates the performance of the unconstrained acceleration engine. The constrained C-DNGT algorithm, which embeds this engine into a primal-dual augmented-Lagrangian loop, is validated on the CO₂-ESS economic-dispatch problem in Sections 6.2–6.5 (Figures 5–7). These experiments demonstrate that C-DNGT enforces the global power-balance constraint to machine precision (~10⁻¹⁵), converges linearly under the μ-strongly-convex regularization, and respects the heterogeneous operational limits of each CO₂-ESS unit. A head-to-head constrained-vs-constrained comparison of C-DNGT against, e.g., a standard distributed dual decomposition on the full N-node, T-step ED problem, with explicit iteration-to-residual curves for both algorithms under identical data, is a natural extension that we leave for future work; the complexity of the coupled time-horizon projection and the dual step-size tuning makes a fair comparison a substantial numerical study in itself. What the present experiments establish is that (i) the Nesterov-accelerated gradient-tracking engine at the core of C-DNGT outperforms the tested unconstrained benchmarks (Figure 8), and (ii) the engine, when embedded into the constrained primal-dual augmented-Lagrangian framework, successfully solves the CO₂-ESS dispatch problem with the convergence guarantees of Theorem 2 (Figures 5–7).

6.4 Parameter sensitivity: curtailment penalty and dual step size

To complement the baseline and case-study experiments, we provide a brief sensitivity analysis of the renewable-curtailment penalty factor γ_*i* and the dual step size η, as requested by the reviewers. Because the constrained coordination couples the local cost penalty γ_*i* with the dual-ascent rate η through the primal-dual updates of the C-DNGT algorithm (Section 4), we examine their separate effects on convergence speed and the resulting economic allocation.

We first note that γ_*i* primarily governs the economic trade-off between renewable curtailment and dispatch cost: a larger γ_*i* suppresses local curtailment, whereas a smaller γ_*i* tolerates more curtailment to reduce the dispatch cost. As shown in Figure 5c, across the tested range γ_*i* ∈ [0.5, 32] the algorithm remains convergent and feasible, i.e., the choice of the penalty factor is robust with respect to convergence.

The dual step size η is, by contrast, the dominant factor for convergence speed and stability. Figure 5a shows representative convergence residuals for several (γ_*i*, η) pairs: an excessively small η slows the dual correction of the global power-balance mismatch and therefore delays convergence, whereas a moderately larger η reaches the tolerance markedly faster. The two-dimensional sweep in Figure 5b (iterations required to reach a residual below 1.5 × 10⁻¹) confirms that, for the tested instance, η in the range [0.1, 0.6] attains the fastest convergence, while a very small η requires roughly an order of magnitude more iterations. Based on these observations, the step sizes used in the main experiments (α = 0.05, β = 0.1, η = 0.1) lie inside the stable and fast region identified by the sweep, and the chosen γ_*i* values are well within the convergent range.

6.5 Augmented Lagrangian validation

To numerically validate the derived augmented-Lagrangian updates Equations 26–30 on the actual economic-dispatch problem, we solve the C-DNGT iterations with exact global-residual tracking (the consensus-reduced form covered by Theorem 2) on the N = 6, T = 24 heterogeneous CO₂-ESS instance described in Section 3. The four panels of Figure 7 focus on the aspects raised by the reviewer: the role of the augmentation penalty ρ in enforcing the power balance and the effect of the strong-convexity regularization μ.

Figure 7a plots the power-balance residual ∥∑i(Aixi-bi)∥ on a logarithmic scale for five values of the augmentation penalty ρ ∈ {0, 1, 3, 6, 12} with μ = 10⁻³. All ρ>0 curves reach machine precision (~10⁻¹⁵) within 1, 000–2, 000 iterations, demonstrating that the derived augmented-Lagrangian updates enforce the global supply-demand balance correctly and efficiently. By contrast, the plain Lagrangian (ρ = 0) reduces the violation only to ~4.7 × 10⁻² after 30, 000 iterations, confirming that the augmentation provides the practical convergence acceleration that its name implies. Panel (b) shows the optimality gap ||x^*k*−x^*||² for ρ = 6, μ = 10⁻³; the gap decays linearly on the logarithmic scale to machine precision (~10⁻³¹), validating the linear convergence rate guaranteed by Theorem 2 under strong convexity.

Panel (c) sweeps ρ and reports the number of iterations required to reduce the balance violation below 10⁻³. The fastest convergence is obtained at moderate ρ (3–6); at the extremes (ρ = 0 or very large ρ) the required number of iterations grows. This is consistent with the known trade-off in the method of multipliers: too small a penalty gives weak coercivity, while too large a penalty stiffens the primal subproblem and forces a smaller primal step size α (Section 5). Panel (d) quantifies the price paid for the strong-convexity regularization. For several μ>0 we solve the regularized problem and measure the dispatch deviation ∥x(μ)-x0*∥ from the unregularized (μ = 0) optimum. The deviation scales as O(μ) (the corresponding cost suboptimality is O(μ2), cf. Section 5); at the chosen μ = 10⁻³ it is ≈2.3 × 10⁻², about 0.1% of the dispatch level, hence negligible for practical dispatch. The role of the regularization introduced in Section 4 is therefore not that the unregularized problem fails to converge here—on this small, well-conditioned instance it reaches the solution as well—but that μ>0 is what makes the provable linear rate of Theorem 2 hold: its primal step-size condition α ≤ μ/(4L²) collapses to α ≤ 0 for μ = 0, so no non-trivial contracting step is guaranteed without it. The strong convexity thus buys a guaranteed linear rate at a quantified, negligible dispatch cost—precisely the point raised in Comment R2–2.

6.6 Reproducibility and code availability

To support reproducibility, we summarize all parameters, random seeds, and data sources used in the main experiments below, and we release the core simulation code.

6.6.1 Complete parameter summary

All numerical values used in the experiments are collected in Table 1 and detailed below.

Table 1. Complete parameterization used in the numerical experiments.

All seeds refer to NumPy default_rng().

6.6.2 Code and data availability

The core Python implementation of the proposed C-DNGT algorithm and the main numerical experiments (including the augmented-Lagrangian validation of Figure 7 and the sensitivity study of Figure 5), together with the processed PJM load profiles and the parameter files, is released under an open-source license at https://github.com/zzandxhk/CO2-ESS-C-DNGT/tree/main/CO2 -ESS-C-DNGT . The repository includes requirements.txt and explicit random seeds so that the released core results can be reproduced from a clean checkout. We encourage readers to inspect, modify, and extend the implementation for their own use cases. The data and code are also available from the corresponding author upon reasonable request during the review process.

7 Conclusions

In this study, we proposed a collaborative optimization framework for integrating Carbon Dioxide Energy Storage Systems (CO₂-ESS) into distributed smart grids. To manage nonlinear thermo-mechanical constraints and supply-demand balance requirements without a centralized coordinator, we developed a Constrained Distributed Nesterov-Accelerated Primal-Dual Gradient Tracking (C-DNGT) algorithm. By integrating Nesterov momentum and repurposing a dynamic tracking mechanism to monitor physical power mismatches, the proposed algorithm addresses ill-conditioned operating cases caused by hardware heterogeneity. We provided a theoretical convergence analysis based on a combined Lyapunov function and a saddle-point residual expansion, showing that the multi-agent network converges linearly to the KKT saddle point under the stated assumptions while satisfying local feasible bounds and the global power balance constraint. Numerical simulations using PJM Interconnection load data compared the proposed method with DGD, EXTRA, and standard GT. The results show that C-DNGT reduces the required number of iterations and communication exchanges in the tested cases while limiting persistent grid mismatches. Future research will focus on extending this framework to asynchronous communication topologies, incorporating stochastic optimization techniques to account for extreme weather events and renewable generation uncertainty, and developing a multi-timescale two-rate dispatch formulation that resolves the fast electrical dynamics and the slow thermo-mechanical transients at heterogeneous resolutions.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Debani Prasad Mishra, International Institute of Information Technology, India

Reviewed by: Sheharyar Hussain, Tianjin Ren'ai College, China

Jiaqi Tian, Northeastern University, China

Footnotes

https://www.kaggle.com/datasets/robikscube/hourly-energy-consumption