Distributed robust optimization of airborne wind energy-integrated energy systems with virtual battery-based demand response

Abstract

Weak-grid rural settlements in remote areas often face fragile grid infrastructure and high fuel transportation costs, which severely limit both the economic performance and operational reliability of local energy systems. While the rapid growth of renewables helps reduce reliance on fossil fuels and strengthens resilience through multi-energy complementarity, the natural variability of wind and solar power becomes much more problematic in weak-grid settings, driving up energy supply risk costs. To tackle these issues, this paper presents a two-stage scheduling framework for weak-grid integrated energy systems (IES). The framework brings demand-side airborne wind energy (AWE) into the system to capture more stable high-altitude wind resources, with power output regulated by adjusting altitude. A double-loop modeling approach allows both the AWE operating height and its power output to take part directly in optimal scheduling. At the same time, a virtual-battery-based demand response mechanism is used to represent load flexibility. Uncertainties are addressed with a Wasserstein distance-driven distributionally robust optimization (DRO) model that includes Conditional Value-at-Risk (CVaR) to quantify tail risks and balance costs against extreme scenarios. Validation on a township-level energy system in Hami, Xinjiang, China, based on five comparative experiments and sensitivity analysis, shows that the proposed framework cuts total operating cost by 3.67% and reduces tail risk loss (CVaR, α = 0.90) by 7.52% compared with the baseline, with the virtual battery-based demand response contributing an extra 10.3% cost saving. Overall, the framework offers a practical and scalable approach for integrating high shares of renewables in weak-grid regions, markedly improving renewable accommodation and the system’s ability to withstand disturbances.

1 Introduction

Against the overlapping background of global low-carbon targets and rapid energy dig-italization, regional energy systems are undergoing a fundamental transition from traditional single-energy infrastructures toward integrated multi-energy and multi-network architec-tures (Hao et al., 2022; Wen et al., 2025). According to IRENA data, global installed renewable power capacity increased from 2,179 GW at the end of 2017 to 4,448 GW at the end of 2024. In 2024 alone, the newly added capacity reached 585 GW, accounting for 92.5% of the global expansion in newly in-stalled power generation capacity (International Renewable Energy Agency (IRENA), 2025). However, the effectiveness of this transition exhibits substantial regional heterogeneity. In China, the central and eastern regions have achieved large-scale deployment of IES through coordinated operation across power generation, transmission, consumption, and energy storage. In contrast, western regions characterized by weak grid infrastructure, such as Xinjiang, continue to face prominent challenges, including limited grid connectivity, high transportation costs for primary energy resources, and insuf-ficient local absorption capacity for renewable electricity (Sun et al., 2021; Wu et al., 2020). Consequently, establishing a resilient and cost-effective regional integrated energy supply system that efficiently ex-ploits local renewable resources under weak grid conditions remains a critical challenge in both engineering practice and academic research.

Current research on this issue primarily follows two technical approaches: First, eco-nomic studies based on market mechanisms and business models, which focus on designing rational market mechanisms and business models to guide investment and consumption be-havior through economic signals (Ge et al., 2021; Schwidtal et al., 2023). Second is the physical research approach based on system modeling and operational optimization. This approach focuses inward on the system, centering on establishing precise mathematical models to characterize the physical proper-ties and coupling relationships among generation, grid, load, and storage. It employs ad-vanced algorithms such as robust optimization and distributed optimization to address dis-patch challenges under the dual uncertainties of high renewable energy penetration and fluctuating load demand (Sun et al., 2022; Ma et al., 2021). These two approaches complement each other: the physical optimization model provides cost and constraint benchmarks for market mechanism design, while market signals can be incorporated as cost factors into optimization objectives, par-ticipating in the planning and operation of physical systems. Together, they form a technical and economic closed-loop system, offering solutions to the energy transition challenges faced by weak grid regions.

In terms of renewable energy supply, wind power has become the preferred choice due to its scalability, low-carbon characteristics, and mature technology (McKenna et al., 2025). However, near-ground wind power is susceptible to the combined effects of wind speed fluctuations and power uncertainty in weak grid scenarios, making it difficult to mitigate the systemic risks posed by grid weaknesses (Chawda et al., 2024). In contrast, AWE utilizes teth-ered aircraft to capture stable high-altitude wind resources at elevations ranging from hun-dreds to thousands of meters, opening new technical pathways for high-penetration renewa-ble energy integration in areas with weak grid infrastructure. Khursh et al. (2025) con-ducted a comprehensive review of AWE research from 2014 to 2024, summarizing progress across dimensions including system architecture, control strategies, energy conversion, and economic evaluations. They noted that at altitudes of 300–1500 m, wind energy density is higher and wind speed profiles are typically smoother than near-surface layers, providing a physical basis for enhancing annual utilization hours and reducing output fluctuations. Zhang et al. (2025) proposes an evaluation method for high-altitude wind resources and pro-vides a quantitative framework tailored to engineering site selection, offering a reliable ba-sis for localized development in areas with weak grid infrastructure. Schmehl et al. (Fagiano et al., 2022) summarized the major achievements and challenges of autonomous AWE systems in attitude control, trajectory planning, and fault tolerance. They noted that as the technology matures, AWEs have begun evolving from experimental platforms into grid-connected, dispatchable power sources. In September 2025, China’s megawatt-scale commercial airborne wind energy system completed its test flight in Xinjiang. Its material and cost per kilowatt-hour demonstrate significant advantages over land-based wind power (Cheng, 2025). However, existing lit-erature has limited research on integrating this technology as a distributed power source into regional integrated energy systems to achieve coordinated optimization and dispatch. Therefore, this study proposes a dual-loop AWE modeling approach. The inner loop outputs the maximum available power upper envelope, while the outer loop incorporates height as a decision variable into the integrated energy optimization process.

To address the uncertain fluctuations of renewable energy and measure associated risks, scholars primarily employ two traditional approaches. One is stochastic optimization, where decision-makers aim to minimize expected costs or risk-adjusted expected costs by solving two-stage or multi-stage stochastic programming models (Herding et al., 2024). For example, Zhang et al. (2020) developed a two-stage stochastic dispatch model that accounts for the spatiotemporal char-acteristics of wind power. The day-ahead stage determines unit start-up/shutdown and base power output plans, while the real-time stage performs output adjustments under various wind scenarios. This ultimately achieves the economic operation of the multi-energy system in an expected sense. Such approaches can intuitively reveal the probability distribution characteristics of renewable energy output while balancing the expected cost and risk cost of system operation. Another approach is robust optimization. It does not rely on precise assumptions about probability distributions. Instead, by constructing confidence intervals or polyhedral uncertainty sets for uncertain variables such as wind and solar power generation and load, it transforms the problem into a min–max or min–max–min form to find the opti-mal solution under worst-case conditions (Yang et al., 2024). For example, Zong and Yuan (2023) employed a hi-erarchical modeling approach for regional electricity and gas supply networks coupled with distributed energy resources. At the lower level, they developed a two-stage robust optimi-zation model incorporating photovoltaic and load uncertainties. The upper and lower levels were solved collaboratively using ADMM. Results demonstrated that robust optimization significantly reduces the dependence of operational strategies on single forecasting scenar-ios under high renewable energy penetration, thereby enhancing the system’s disturbance resilience. Research indicates (Lin et al., 2025; Du et al., 2024; Li X. et al., 2024) that the results of stochastic optimization are highly dependent on the specified empirical probability distribution. When sample sizes are limited or statistical characteristics drift over time, optimal solutions become extremely sensitive to empirical distributions, leading to distorted evaluation outcomes. Robust opti-mization often requires reserving substantial redundant capacity for extremely low-probability extreme scenarios, resulting in significantly increased investment and oper-ational costs and exhibiting overly conservative characteristics.

Yet, while stochastic and robust optimization methods can address uncertainty, they fail to jointly consider probabilistic distributions and extreme tail risks—particularly in data-limited weak-grid scenarios, where solutions tend to be overly conservative or highly sensitive. If the empirical distribution is directly treated as the true distribution under such conditions, a significant deviation may exist between the expected optimal dispatch and the actual operational risk. Therefore, DRO seeks to re-define uncertainty on this foundation. DRO introduces the concept of a fuzzy set-a collec-tion encompassing all possible probability distributions centered around an empirical dis-tribution or benchmark distribution, with a radius defined by a certain probability distance. By leveraging available historical data and statistical information to anchor fuzzy sets, it avoids excessive conservatism. Simultaneously, it acknowledges data limitations and model errors, ensuring decision reliability by optimizing for worst-case scenarios within the fuzzy sets. Against this backdrop, the choice of distance metric becomes central to the research, with the Wasserstein distance emerging as the preferred option due to its exceptional statis-tical properties. Grounded in the optimal transport theory of probability distributions, it leverages the geometric information of data samples to construct mathematical models that are easy to handle and can even be reformulated as convex optimization problems. Li (Li et al., 2023) et al. proposed a distribution-robust optimization scheduling model for electrothermal IES based on an improved Wasserstein distance. By incorporating distribution information of extreme scenarios into traditional fuzzy sets, they demonstrated that DRO achieves a smoother trade-off between expected cost and worst-case cost under conditions of source-load uncertainty and multi-functional zoning. Furthermore, to effectively quantify tail risks under extreme scenarios, Xuan et al. (2021) introduced the CVaR metric into IES planning. By embedding it into the objective function, they avoided the risk myopia that arises from optimizing solely for expected costs, thereby enhancing the economic robustness of the planning scheme. In summary, this study employs a two-stage framework to construct distribution fuzzy sets centered on empirical distributions using the L1-norm-based Wasserstein-1 distance. This approach characterizes the distributional un-certainty of renewable power output under data-driven conditions. Furthermore, by incor-porating CVaR as a risk metric into the model, the objective function includes the expected value of system operating costs exceeding the VaR threshold at confidence level α. This ap-proach penalizes and constrains low-probability, high-loss scenarios. Consequently, the proposed method effectively suppresses the impact of tail risks on the security and eco-nomic efficiency of power supply in weak grids while ensuring scheduling feasibility, thereby enhancing the resilience of system operations.

Correspondingly, virtual energy storage has become a major research hotspot in the field of demand-side response in recent years. In parallel with the energy structure transformation on the supply side, and to better mitigate the strong fluctuations caused by high-proportion renewable energy, the load side must also evolve from the traditional rigid consumption model toward an active, adjustable, and bidirectional flexible regime (Cabot and Villavicencio, 2024). Demand response mecha-nisms provide critical flexibility to the system by guiding users to adjust their electricity consumption through price signals or incentives. Conventional demand response typically classifies loads into three types: shiftable loads, interruptible loads, and inelastic loads (Vahedipour-Dahraie et al., 2018). To further unlock flexibility, current research is gradually shifting from simple price-based demand response toward treating adjustable loads as variable energy units of equal im-portance to physical energy storage. Wang et al. (Hassankashi et al., 2026) unified the modeling of temperature comfort ranges for adjustable thermal loads with transferable time windows for electrical loads in the operational optimization of regional integrated energy systems, achieving peak shaving and valley filling through coordinated electricity-heat demand response. Liu et al. (Liu et al., 2022) treated interruptible and transferable loads as flexible resources within a two-level dispatch framework for regional integrated energy systems. By coordinating different stakeholders and flexible load management through upper-level game theory and low-er-level optimization, they demonstrated the value of demand-side flexibility in enhancing system economic efficiency and security margins. However, most of these models remain constrained to the peak and trough values of load power curves. They either approximate response magnitude using price elasticity coefficients or characterize load peak-shaving behavior through time windows and switching variables. Consequently, they struggle to continuously and precisely describe the system’s power imbalance in energy terms, nor can they provide equivalent representations directly compatible with physical energy storage across multiple types of flexible loads. To address this issue, some scholars have begun in-corporating the concept of virtual energy storage. They model electric heating and air condi-tioning loads as virtual energy storage units with charge/discharge power and capacity con-straints, introducing a unified energy state variable to participate in demand response and scheduling optimization within integrated energy systems. For instance, Bu et al. (2024) de-veloped a Virtual Energy Storage (VES) model within a community integrated energy system by leveraging building thermal mass. This model collaborates with physical storage systems to achieve demand response control across multiple time scales through rolling optimization. Fu et al. (2024) summarized the modeling approach and key metrics for uniformly representing sources, grids, and loads as virtual energy storage systems, noting that this energy stor-age-based characterization provides a powerful tool for uniformly quantifying demand-side flexibility. Building upon this foundation, this paper further proposes a virtual battery-based demand response model: Equivalent to the shifting behavior of adjustable loads as the charging and discharging processes of a virtual battery, the state of charge (SoC) character-izes its cumulative peak-shifting capacity and sustainable response space. Through power, capacity, and mutual exclusion constraints, the virtual battery is coordinated with physical electrical and thermal energy storage within a unified optimization framework. This ap-proach maximizes the buffering capacity of demand-side flexibility against fluctuating floating wind power in weak grid scenarios.

Recent advances in the optimal scheduling of integrated energy systems from 2024 to 2026 are summarized in Table 1. Gao et al. (2024) proposed a low-carbon scheduling model for IESs considering offshore wind power-to-hydrogen production and dynamic hydrogen blending, addressing the challenges of offshore wind power accommodation, low-carbon utilization of hydrogen energy, and the nonconvexity-induced computational difficulty of hydrogen-blending models. Zhang et al. (2024) developed a two-stage robust operation model for electricity–gas–heat multi-energy microgrids, aiming to coordinate system economy and robustness under heterogeneous uncertainties associated with renewable energy, loads, and multi-energy networks. Fan et al. (2024) proposed a distributionally robust adaptive model predictive control scheduling method, which improves rolling optimization and real-time scheduling adaptability of IESs under renewable energy forecasting errors and distributional uncertainty. Li F. et al. (2024) established a data-driven DRO model considering refined hydrogen utilization and stepped carbon trading, addressing the insufficient exploitation of the low-carbon value of hydrogen energy and the conservatism of scheduling decisions under renewable energy uncertainty. Gao et al. (2025) proposed a three-stage stochastic–robust optimization model for renewable-energy-based IESs, balancing the lack of extreme-scenario protection in stochastic optimization and the excessive conservatism of robust optimization. Wang et al. (2026) developed a multistage robust scheduling model for electricity–hydrogen systems considering nonanticipativity and robust feasibility, addressing multistage decision consistency and operational feasibility under renewable energy uncertainty. Chen et al. (2025) proposed a risk-averse stochastic scheduling model for electricity–heat IESs considering microgrid interfaces, in which CVaR is introduced to quantify the financial risks caused by wind and solar uncertainty, and the model is solved using SOCP. In summary, existing research on AWE primarily focuses on its physical characteristics and low-level control strategies, with few studies integrating it into IES optimal scheduling frameworks, thereby failing to fully exploit its active dispatch potential. Conventional demand response models are largely confined to peak-valley constraints and struggle to provide continuous and precise quantification of power imbalance, while mature modeling methods tailored to weak-grid scenarios remain underdeveloped. Moreover, traditional robust optimization approaches either suffer from excessive conservatism or fail to account for the probabilistic characteristics of uncertainties, rendering them unable to achieve precise trade-offs between scheduling costs and operational risks. To address these limitations, this paper proposes a distributionally robust optimal scheduling method for weak-grid integrated energy systems that incorporates airborne wind energy and virtual battery-based demand response. The main innovations are as follows. A comprehensive physical model of the IES including AWE is constructed, along with an upper envelope model of available power output that accounts for altitude adjustment and energy consumption, treating the platform operating height as an optimizable decision variable to fully unlock its active dispatch potential. On the demand side, a virtual battery-based demand response model is introduced together with a fallback mechanism for extreme wind power shortage scenarios, enabling demand response to contribute regulation capacity while ensuring system feasibility and achieving precise power balance characterization suitable for weak grids. We further develop a two-stage distributionally robust optimization framework based on Wasserstein distance, incorporating CVaR into the objective function for tail risk quantification and solving the model efficiently via the column-and-constraint generation (C&CG) algorithm, thereby realizing a precise trade-off between scheduling costs and extreme operational risks.

2 Mathematical modeling

2.1 Structure of the weak grid integrated energy system

This study investigates a weak-grid regional integrated energy system located in Hami, Xinjiang, China. The region has abundant wind, solar, and geothermal resources. However, due to limited transmission capability and sparse load distribution, relying only on the supply from the upstream power grid is neither economical nor reliable. We construct an IES energy flow model as shown in Figure 1. The main energy inputs are electricity and natural gas, and energy conversion is realized through various energy devices, providing electricity and meeting heating demand for the town on an hourly basis. Airborne wind energy provides a more stable and altitude-adjustable renewable output for the weak-grid integrated energy system and reduces the dependence on fossil fuels; virtual-energy-storage-based demand response realizes peak shaving and valley filling and fast reserve by characterizing and controlling adjustable loads, thereby enhancing the economic performance and resilience of the system under multiple uncertainties.

2.1.1 Gas turbine and waste heat boiler

The energy-conversion and operating constraints of the gas turbine are formulated in Equations 1–3.

where denotes the lower heating value of natural gas, denotes the electrical generation efficiency of the gas turbine, and represents the natural-gas input of the GT in period .

Unit Start-up/Shut-down and Ramping Constraints:where and denote the minimum stable output and the maximum output limit of the gas turbine, respectively; denotes the on/off status of the gas turbine at time , with a value of 1 indicating that the unit is online; and and are the start-up and shut-down variables, respectively.where and denote the ramp-up and ramp-down limits of the gas turbine under normal operation, respectively.

2.1.2 Waste heat recovery

The waste-heat recovery relationships are formulated in Equation 4.

where is the fraction of GT exhaust heat absorbed by the organic Rankine cycle (ORC), and is the ORC electrical efficiency. denotes the heat dissipation loss coefficient; is the fraction of GT exhaust heat absorbed by the waste heat boiler (WHB), is the WHB thermal efficiency, and is the WHB heat output power.

2.1.3 Geothermal well heating supply

The geothermal well, photovoltaic, and electric-boiler models are given in Equations 5–7, respectively.

where denotes the geothermal heat output at time t. is the overall efficiency of the geothermal heat exchanger. is the density of the geothermal fluid. is the specific heat capacity of the geothermal fluid. is the fluid volumetric flow rate at time t. and denote the production temperature and reinjection temperature.

PV:

Coupling between the electric boiler power and the heat supplywhere is the heat output power of the EB in period t. is the input electric power of the EB, and is the electricity-to-heat efficiency of the EB.

2.1.4 Power balance

The electrical power balance constraint is given in Equation 8.

where: , , and denote the electricity generation of the gas turbine, the airborne wind energy system, and the PV plant, respectively; and denote the ORC power generation and the EB power input at period , respectively; and denote the discharging and charging power of the BESS, respectively; and denote the electricity purchased from and sold to the upstream power grid, respectively; and denotes the net electricity demand of the town.

2.1.5 Heat balance

The thermal power balance constraint is given in Equation 9.

where, denotes the geothermal heat output; denotes the heat recovered by the WHB from the GT exhaust gas; denotes the heat output of the EB; and denote the heat discharged from and absorbed by the TESS, respectively; and denote the heat purchased from and sold to the external district heating network, respectively; and denotes the thermal load demand at time .

2.2 Airborne wind energy system

The introduction of the Airborne Wind Energy (AWE) platform is one of the innovations of this study. From a physical perspective, it is a tethered rotorcraft or wing that captures high-altitude wind energy through complex aeroelastic and electromechanical interactions. A complete model includes multi-degree-of-freedom platform motion, tower-top deformation, generator dynamics, and control inputs such as generator torque and blade pitch angle. However, if such a high-dimensional and nonlinear model is directly embedded into the intraday scheduling problem, the computational burden will be very large. In addition, this study focuses on intraday scheduling of the regional integrated energy system at the hourly scale rather than flight control design of the AWE system. Therefore, we establish a double-loop airborne wind energy model. The overall double-loop modeling framework of the AWE system is illustrated in Figure 2. The inner-loop control layer regulates generator torque and pitch angle through appropriate control laws, so that the AWE system operates within a pre-designed power envelope under given altitude and wind speed profiles. The outer-loop scheduling layer no longer solves high-dimensional state-space equations, but only optimizes a few representative decision variables, such as the altitude factor and the net output power. The specific model is given as follows:

The inner-loop control layer uses a continuous-time formulation to describe the platform motion.

Considering the main motions of the airborne platform in the wind-direction plane and the vertical plane, the degrees of freedom can be analogous to the surge, heave, and roll motions of a floating platform. Meanwhile, the equivalent bending deformation at the tower top and the generator rotational speed are introduced. To describe the coupling of these slow variables, the state vector is defined aswhere and denote the surge displacement and the heave displacement of the platform in the horizontal plane along the wind direction, denotes the roll angle of the platform about the longitudinal axis, and denotes the equivalent bending displacement at the tower top. Under the small-disturbance assumption, the structural dynamics can be written as a mass–damping–stiffness system:where and denote the equivalent mass matrix and stiffness matrix, denotes the structural damping matrix, which includes material internal damping and equivalent viscous damping due to wind-field effects and fluid viscosity. denotes the projections of aerodynamic forces such as lift and drag on each degree of freedom, and denotes the components of tether tension in the corresponding degrees of freedom.

To facilitate coupling with the control inputs, the above second-order equations can be written in a first-order state-space form:where is the control input, is the external disturbance (e.g., cable disturbance and turbulence), and contains information such as structural mass, stiffness, and viscous damping coefficients.

In this study, the physical effects of viscous damping are mainly reflected in two aspects. First, part of the platform vibration energy is dissipated through the linear damping term , which limits the amplitudes of surge, heave, and roll. Second, it represents an equivalent reduction of aerodynamic power in real scenarios, thereby reducing the effective wind energy that can be transmitted to the generator side.

Airborne wind energy essentially follows the aerodynamic energy conversion law of wind turbines. Assume that the platform operates under an instantaneous relative wind speed , which is obtained by superposition of the altitude wind speed and the platform motion. The aerodynamic torque can be expressed as:where denotes the air density, R is the rotor radius, denotes the torque coefficient, which depends on the tip-speed ratio and the pitch angle , and is one of the control inputs of the pitch system.

On the generator side, angular momentum conservation is satisfied:where denotes the generator rotational inertia, denotes the generator electromagnetic torque controlled by the converter, and is the equivalent viscous damping coefficient.

The electrical output power of the airborne wind energy system is:where denotes the equivalent efficiency of the power generation chain. In the model of this study, by adjusting and the system operates within a safe power region under given altitude and wind speed conditions.

The outer-loop scheduling layer does not directly solve the high-dimensional state-space equations, but represents them as a power upper envelope and records altitude regulation and energy consumption.

The outer-loop formulation, including altitude mapping, wind-speed extrapolation, power-envelope construction, linear interpolation, altitude-rate constraints, and regulation energy consumption, is given in Equations 16–25.

First, the altitude factor is defined, and the corresponding physical altitude is linearly mapped as follows:where denotes the normalized operating altitude, denotes the actual physical altitude, and denote the maximum and minimum operating altitudes of the airborne platform, respectively.where denotes the wind speed at the reference height . denotes the effective wind speed at this height, and denotes the power exponent related to surface roughness.

In a continuous-time interval , the inner loop receives the altitude decision given by the outer loop and the wind condition information under the scenario , and then obtains the power upper envelope curve based on Equations 10–15:where is the set of all curves that satisfyEqs. (10)–(15).

The outer-loop power output is constrained by the upper envelope curve.

Since comes from the inner-loop optimization definition, it is difficult to directly call it in the outer-loop MILP. Therefore, we evaluate the envelope definition offline and then adopt piecewise linear interpolation in the outer loop.

Since is defined by the inner-loop optimization, it is difficult to use it directly in the outer-loop MILP. Therefore, we evaluate the envelope offline and then adopt piecewise linear interpolation in the outer loop.

Select K altitude nodes and define the interpolation coefficients .

SOS2 constraints are added so that for each ,only two adjacent coefficients are nonzero.

The altitude is obtained by interpolation as:

For each , , , the envelope value equals

The upper bound constraint of available power in the outer loop is:

Considering that platform altitude changes require additional energy to overcome gravity and aerodynamic drag, and to keep a linear structure in the optimization framework, the variation of the altitude factor is introduced, and auxiliary climbing and descending variables, as well as the maximum climbing and descending rates, are defined.where denotes the maximum rate of altitude change per unit time.

The energy consumption of platform altitude regulation is linearly calculated as:where and denote the equivalent energy consumption coefficients per unit increase/decrease of the altitude factor.

2.3 Energy storage modeling

The Battery Energy Storage System (BESS) regulates the power balance of the grid through charging and discharging. It is a key unit to buffer the output fluctuations of renewable energy, such as airborne wind energy, and to improve system flexibility. This study assumes that the BESS adopts fast-response energy storage technology such as lithium batteries. The specific model is as follows: The BESS and TESS models, including state transitions, operating limits, charging/discharging exclusivity, and degradation costs, are formulated in Equations 26–33.where denotes the energy of the BESS in period t, and are the charging power and discharging power, and are the charging and discharging efficiencies, and is the self-discharge coefficient.

To ensure that the battery operates in a safe range, the total stored energy shall satisfy the upper and lower bounds.

The per-unit state of charge is introduced to represent the ratio of the related energy to the maximum capacity, and the charging and discharging power are limited by the converter capacity:

The Thermal Energy Storage System (TESS) shifts thermal energy in time to achieve off-peak operation of heating demand and efficient utilization of waste heat. In this study scenario, the TESS can be regarded as a centralized heat storage device equipped with water tanks or phase change materials, with a certain energy capacity and upper limits of charging and discharging heat power. The specific model is as follows:

The total stored thermal energy shall satisfy the upper and lower bounds:

The per-unit state of charge of thermal energy storage is introduced, and the upper and lower limits of charging and discharging thermal power are given:

To avoid simultaneous charging and discharging of electricity and heat in the same time period, binary variables are used to model the mutual exclusion relationship:

To accurately represent the lifetime degradation caused by aging of the energy storage system in the economic dispatch model, a linear degradation cost function based on energy throughput is adopted:where is the converted degradation cost coefficient per unit charging/discharging electrical energy, and is the comprehensive cost corresponding to per unit charging/discharging thermal energy.

2.4 Modeling of the power grid and district heating network

The electricity and heat purchase/sale constraints are formulated in Equations 34, 35, respectively.

To avoid large-scale electricity purchase and electricity sale at the same time in a physical sense, a binary variable is introduced to describe the mutual exclusion between purchasing and selling electricity. The electricity purchase and sale model with the power grid is as follows:where: when = 1, only electricity purchase is allowed at the node; when = 0, only electricity sale is allowed. In weak-grid scenarios, reflects the limits of the interconnection line capacity and the power supply capability of the upstream grid.

The local heating system is connected to the district heating network through the heat exchange station. and denote the thermal power purchased from and sold to the external heating network, respectively, A binary variable is introduced to describe the mutual exclusion between purchasing and selling heat.

The electricity purchase and sale cost and the heat purchase and sale cost are as follows: The corresponding electricity and heat transaction costs are calculated using Equation 36.where and are the electricity purchase price and electricity selling price at time t, respectively, and denote the heat purchase price and the heat selling price at time , in turn.

2.5 Uncertainty and scenario generation

Under the weak-grid scenario with high renewable energy penetration considered in this study, multiple sources of uncertainty exist in the operation of the integrated energy system, mainly including three aspects. First, the uncertainty in high-altitude wind speed and the available power of airborne wind energy is considered. The airborne wind power output is affected by the high-altitude wind speed profile, turbulence intensity, and platform control strategy, and shows significant randomness and temporal correlation, which is the most important random factor in this study. Second, the uncertainty in electric load and thermal load is considered. Due to the diversity of resident behaviors and fluctuations of weather conditions, electric and thermal loads have random disturbances and day-type differences at the hourly scale. Third, the uncertainty in the electricity price is considered. In weak-grid regions, the peak-period electricity purchase price may be significantly higher than that in off-peak periods.

In this framework, the above uncertain quantities are uniformly represented as scenario-based random vectors. The uncertainty vectors, scenario-generation models, and empirical distribution are defined in Equations Equations 37–41. Considering a rolling dispatch period with a 24 h horizon, the random parameter vector of scenario s is defined as:

Considering that high-altitude wind speed has clear autocorrelation and intraday variation, a first-order autoregressive process is used at the reference height to describe the temporal correlation of wind speed, and a deterministic intraday mean wind speed curve is superposed:where is the wind speed of scenario s at the reference height, and is the intraday mean wind speed obtained by fitting historical data. is the wind speed lag correlation coefficient, is the standard deviation of process noise, and is an independent standard normal random variable. By adjusting and , the smoothness and fluctuation intensity among scenarios are controlled. The initial wind speed is sampled from historical data, and the procedure is repeated times to obtain an original scenario pool of wind speed at the reference height.

The uncertainty of electric load and thermal load mainly comes from differences in user behavior and fluctuations in weather conditions. Based on the baseline load curves of representative days, random disturbances are added to construct load scenarios.where and are the intraday baseline curves of electric and thermal loads obtained by clustering historical data, and are a zero-mean disturbance coefficient, and are used to control the uncertainty level.

In weak-grid regions, the electricity purchase price from the upstream grid usually shows clear peak–valley characteristics and may fluctuate with load level and market conditions. This study constructs electricity price scenarios by superposing random disturbances on a given intraday baseline electricity price:where is a zero-mean random disturbance coefficient of electricity price, and controls the fluctuation amplitude of the price spread.

With the above scenario generation mechanism, given a sample size N, N independent joint scenario trajectories can be obtained to form the sample pool for subsequent distributionally robust optimization, and the empirical distribution is defined as:where is the Dirac measure concentrated on sample . Under the empirical distribution, the probability of each scenario is 1/N.

In the traditional stochastic optimization framework, is directly used as the probability distribution of uncertain parameters and participates in the calculation of the expected cost with fixed weights. However, under weak-grid scenarios with a limited number of samples and where the true distribution is difficult to accurately obtain, relying only on the empirical distribution is prone to overfitting. Therefore, to address the above issue, this study constructs a distributional ambiguity set under the Wasserstein distance and optimizes the worst-case distribution in the upper level of the two-stage optimization model.

3 Methodology

3.1 Decision hierarchy and two-stage structure

The operation decisions of the integrated energy system are divided into two levels according to time scale and information availability: day-ahead decisions and scenario-response decisions.

First of all,without fully knowing future high-altitude wind speed, loads, and electricity price, the system operator needs to determine a set of decision variables in advance, including the baseline output trajectory of the gas turbine, the altitude factor of airborne wind energy and its safe operating boundary, the intraday baseline charging/discharging plan of BESS/TESS and the SOC constraints, the baseline adjustment trajectory of demand response, and the electricity/heat purchase and sale strategies with the upstream power grid and the district heating network. These decisions are denoted by a vector , and the feasible region is described by the linear constraint set .

After the uncertain quantity is revealed, the system can make scenario-dependent corrections to unit outputs, energy storage regulation, and demand response according to the available power of airborne wind energy, loads, and electricity price under a specific scenario, so as to eliminate energy imbalance and reduce the total cost and risk as much as possible. For scenario , the second-stage decision vector is denoted by , and its feasible region is jointly determined by the first-stage decisions and the scenario parameters, denoted by .

The two-stage cost decomposition, recourse function, Wasserstein ambiguity set, and distributionally robust objective are formulated in Equations 42–47.

The total cost of the integrated energy system can be decomposed into the first-stage cost and the second-stage scenario cost:where includes gas turbine start-up/shut-down, some fixed operating costs, and capacity-related costs; reflects the actual operating cost under scenario .

Given a scenario and the first-stage decisions , the second-stage problem can be written as:where denotes, under scenario , the minimum operating cost required to satisfy all balance constraints and equipment operating constraints given the first-stage decisions.

Wasserstein ambiguity set and worst-case expected cost.

In the discrete-scenario framework, centered at , by treating as shifts of scenario weights , the Wasserstein ambiguity set is approximated as a probability set with -norm constraints:where is the Wasserstein radius, which determines the expansion degree of the ambiguity set relative to the empirical distribution. Within this ambiguity set, scenario weights are re-adjusted under a distance constraint to approach a potential worst-case distribution and achieve risk control.

Based on the above ambiguity set, the distributionally robust two-stage optimization objective of this study can be expressed as:

In the discrete-scenario setting, the problem of taking the supremum of expected loss over the ambiguity set in the original distributionally robust optimization model can be equivalently transformed into finding a worst-case probability weight vector in a convex polytope defined by a probability simplex and additional constraints. This vector maximizes the linear combination of scenario losses, and its mathematical form is a solvable linear programming problem. is the weight set induced by the Wasserstein radius , and the parameter has a one-to-one correspondence with the original Wasserstein radius . Based on the probability simplex, this set restricts the total deviation of the weight vector relative to the uniform distribution through an -norm constraint. Its physical meaning is that, within a total deviation not exceeding , the probability mass is allowed to be reallocated from low-cost “moderate” scenarios to high-cost “adverse” scenarios, thereby testing the robustness performance of the dispatch strategy under worst-case probability weighting.

Thus, the two-stage problem can be written as:

The first stage aims to select a solution that is relatively safe for all possible scenario-weight distributions, while the inner problem selects the worst-case scenario weights within the Wasserstein ambiguity set.

3.2 CVaR embedding of the risk measure

Only making the expected cost distributionally robust may still fail to directly control the tail risk of extremely expensive scenarios. Therefore, this study further introduces Conditional Value-at-Risk (CVaR) as a tail-risk measure to explicitly penalize the cost exposure under extreme scenarios.

Han et al. (2023) developed a two-stage distributionally robust conditional value-at-risk (DRCVaR) model based on a moment ambiguity set for scheduling zero-carbon island integrated energy systems. By simultaneously capturing intra-day uncertainty losses and dynamically revising day-ahead schedules, their approach effectively balances economic performance with extreme risk control, providing a valuable reference for CVaR-based risk quantification in island scenarios (Han et al., 2023).

The scenario cost, CVaR formulation, and its tractable robust reformulation are given in Equations 48–52.

Scenario total cost is:

Under a given distribution ,the CVaR of at the confidence level is defined as follows. For the discrete-scenario setting under with scenario weights , the CVaR can be written as:where: is VaR an approximation of the VaR level, for scenario exceeding the VaR in scenario VaR the nonnegative excess cost beyond the VaR in each scenario.

Within this framework, we do not measure CVaR, but instead Wasserstein the worst-case upper bound over the Wasserstein ambiguity set for CVaR an upper bound. This not only means avoiding anchoring the decision to a potentially misleading tail of the empirical distribution, but also constructing a robust risk barrier that strictly covers all distributions in a statistical neighborhood, so as to actively constrain extreme losses under the most adverse cases.where: is the risk-aversion coefficient, which controls the weight of tail risk in the overall objective. The larger is, the more the model tends to sacrifice economic performance to reduce tail risk; the larger is, the more severe the weight shift that the adversary can impose, and the more sensitive the model is to out-of-sample extreme cases.

To solve the above min-max problem in a mixed-integer linear programming framework, we introduce the scalar variable θ as an upper bound of the value of this combined objective under the worst-case distribution. The original problem is equivalently transformed into minimizing θ,and satisfying that for all distributions , the corresponding weighted objective value does not exceed θa single optimization problem.

For all weight vectors we have:

Rather than being a fully equivalent reformulation of the continuous Wasserstein-DRO problem, this method should be regarded as a tractable approximation based on a finite set of representative scenarios.

3.3 Solution algorithm based on column-and-constraint generation

Column-and-Constraint Generation (C&CG) is a decomposition algorithm for two-stage robust and distributionally robust optimization. Its core idea is to keep only a limited number of representative scenarios and their weights in the master problem, while the subproblem, given the master solution, searches for a new worst-case weight distribution over the full scenario pool. If the resulting overall risk exceeds the current upper bound, the corresponding weight vector is added to the master problem. This process repeats to iteratively update the scenario set and constraints until convergence. Within this framework, the C & CG solution procedure is illustrated in Figure 3.

First, we input all parameters of the IES and generate wind speed sequences based on the Monte Carlo method. After wind-shear conversion, the available upper envelope of airborne wind energy and the available PV output are obtained, forming the original scenario pool.

The master problem, worst-case subproblem, and convergence criteria of the C & CG algorithm are given in Equations 53–56.

Then, the iteration counter k = 0, in the -th iteration, given the current scenario set and its associated Wasserstein ambiguity-set approximation constraints, we solve the master problem:where is the set constructed by the explicitly enumerated weight vectors in the current iteration. Solving the master problem yields the first-stage decisions .

Then, the subproblem is solved to search for the worst-case weights. After obtaining , it is necessary to check whether, under the full scenario pool and the Wasserstein ambiguity set, there still exists a more adverse weight distribution that makes the distributionally robust cost exceed the current upper bound . Therefore, the following subproblem is constructed:

This problem on is a linear programming problem and can be transformed into an equivalent dual problem; the optimal value is denoted by , and the optimal weight vector is ,If:

Then, under all allowable weight distributions, the current solution is sufficiently robust, and the algorithm can be considered convergent.

If:

Then, this means that there still exists a more adverse distribution not considered by the master problem; therefore, several scenario indices with large weights are added to the current scenario set, and the corresponding distributionally robust inequalities are added to the master problem in linear form, updating the constraint set.

3.4 Virtual-battery demand response

In a regional integrated energy system with a weak grid and limited transmission capability, relying only on conventional units and energy storage devices to absorb renewable fluctuations, such as airborne wind energy, often requires large capacities, leading to high investment and operating costs. On the other hand, residential and public-building loads have a certain degree of time shiftability and compressibility. If part of the load can be shifted upward and downward while ensuring basic comfort, system flexibility can be enhanced without significantly increasing hardware investment.

However, directly describing demand response as “how much is reduced” and “how much is compensated” makes the time-coupling relationship unclear, and it is difficult to ensure that“the reduced load will definitely be compensated in the future,” and it is also difficult to uniformly characterize the flexibility of different types of loads at the optimization-model level. Therefore, this study adopts the Virtual Battery modeling idea and treats adjustable loads as a virtual energy storage device with capacity, power, efficiency, and terminal constraints. In this framework, load reduction is regarded as charging of the virtual battery, corresponding to reduced electricity consumption in the current period, while upward shifting power is regarded as discharging of the virtual battery, corresponding to increased electricity consumption in future periods. The energy state of the virtual battery represents the accumulated shifted energy, and it should return to the initial state at the terminal time as much as possible to keep the total daily energy use nearly unchanged. However, in extreme scenarios, a certain proportion of reducible load is allowed, constrained by a high penalty term, to ensure that real load shedding is triggered only in extreme cases. The detailed model is as follows:

The virtual-battery demand-response model, including load decomposition, state dynamics, operating constraints, response costs, and the revised power balance, is given in Equations 57–67.

The total electric load is decomposed into a rigid load part and an adjustable load part:where is the rigid load and does not participate in demand response; is the adjustable part. According to the load curves given in Chapter2, the adjustable proportion is represented by the flexibility coefficient as follows:

Under the virtual-battery framework, two continuous variables are introduced: denotes the reduced adjustable load in this period, denotes the upward-shifted power in this period. The equivalent electric load participating in demand response is:where denotes the interruptible load power that is reduced and not compensated under extreme conditions, and denote the upper limits of downward-shifted power and upward-shifted power, respectively, which are determined by together with the total load, is the proportion of real reducible load allowed.

Binary constraints are introduced to prohibit simultaneous upward shifting and downward shifting power:

To ensure energy conservation of demand response over time, the energy state of the virtual battery is introduced, which represents the cumulative energy change from the initial time to the current time. Following the modeling of battery energy storage, the state update equation of the virtual battery is:where: is the equivalent efficiency of downward shifting power, is the equivalent efficiency of upward shifting power, which is used to characterize the comfort loss caused by user behavior adjustment.

The energy state of the virtual battery is subject to capacity constraints:

For users, downward shifting cannot be delayed indefinitely, and the reduced load needs to be compensated within a limited time window. Let be the maximum allowable shifting duration, and the coupling between downward shifting and upward shifting can be described by cumulative constraints:

The downward shifting action performed in period t must be fully compensated within the next time steps.

In addition, to avoid frequent switching between upward and downward shifting in a short time, a smoothing constraint is imposed on the rate of change of and , as follows:where and are the allowable upper limits of changes in upward and downward shifting power between adjacent periods, thereby avoiding, to some extent, abrupt regulation on the user side.

Under extreme scenarios of airborne wind power shortage or high electricity price, it may still occur that all units and energy storage operate at full load, but the total load cannot be fully satisfied. Therefore, this study allows a certain proportion of interruptible load and uses a high penalty coefficient to account for its cost.where: is the total interrupted load power, is the penalty coefficient per unit interrupted load.

The demand response cost is:where: and are the cost coefficients of upward shifting and downward shifting, respectively.

By replacing the original load with the equivalent load after demand response correction, the new power balance is:

Thus, virtual-battery-based demand response acts as a flexible load storage on the power side and participates in load balancing and cost optimization together with the energy storage system. Since its decision variables and constraints are all in linear form, it can be naturally embedded into the distributionally robust optimization framework constructed above and serve as a part of the second-stage scenario response , and under different scenarios, by adjusting , ,and to buffer the risks caused by uncertainties such as airborne wind energy.

4 Case study

This section validates the proposed integrated energy system scheduling framework with airborne wind energy, virtual-battery demand response, and Wasserstein-metric distributionally robust optimization. A township-level case study is conducted in Hami, Xinjiang, China, which is considered a representative weak-grid scenario characterized by sparse loads, limited network support, and high marginal costs for fuel-based generation and energy storage deployment. The basic data, including wind speed and electricity prices, were provided by a collaborating partner under a data-sharing agreement and were subsequently processed by our team. The scheduling horizon is 24 h with an hourly resolution, and the simulation parameters are listed in Table 2.

To ensure reproducibility, all simulations are conducted using the same parameter set, scenario generation pipeline, and solver configuration across cases. The results are reported through a unified set of tables and figures that jointly demonstrate: the realism and diversity of renewable/load scenarios; the convergence behavior of the nested C&CG procedure; the energy-balance feasibility in representative scenarios; and the operational value of AWE height adaptation, energy storage, and VB-DR under distributional uncertainty. In particular, the operating altitude of the AWE is optimized within 500–1000 m, and the energy consumption of altitude regulation is included in the objective function to avoid unrealistic frequent climbing or descending.

In this study, Monte Carlo sampling is adopted on the load side to generate intraday trajectories of wind speed, PV output, and electric/thermal loads. Then, a clustering method is used to extract representative scenarios and assign probabilities, so as to characterize uncertainty. Figure 4 shows the mean curves of the corresponding typical scenarios and the upper-envelope power curve of airborne wind energy for the typical day. It can be observed that, in addition to the intermittency of renewable energy, there is also a tail-coupling phenomenon between low-wind periods and high-load periods.

To separate the contributions of different modeling methods, five comparative cases are set up:

Table 3 provides the cost decomposition and key operating indicators of each case for quantitative comparison. C1 is the first-stage operating cost, and C2 is the worst adjustment cost in the second stage. After introducing CVaR, a risk premium appears. Case 1 adopts the basic dispatch scheme. Case 2 modifies the wind power modeling in Case 1 to the inner–outer loop AWE proposed in this paper, where wind power can provide more stable electricity through altitude regulation. Both C1 and C2 decrease, indicating that the adjustable altitude decision of AWE can improve effective output within the power upper envelope, replace part of the high-cost energy supply with wind power in the day-ahead stage, and reduce cost correction in the second stage. Case 3 adds the demand response mechanism proposed in this paper on the basis of Case 2, and the total cost decreases by about 10.3%. This is because the mechanism can effectively smooth the net load curve and improve the time matching between energy supply and demand. Case 4 introduces Wasserstein-DRO (ε = 0.05) on the basis of Case 3. At this time, C2 is the worst-case expected recourse cost in the uncertainty set, which drives the first stage to adopt a more conservative resource allocation plan. Case 5 further introduces CVaR on the basis of Case 4 and adds a penalty term in the objective function. For extreme scenarios, tail-risk loss is suppressed, thereby improving the model’s resilience under weak-grid scenarios.

Comparing Case 4 and Case 5, the CVaR decreases from CNY 14,206.7 to CNY 13,138.2, indicating that the average additional loss of the system under the worst 5%–10% high-cost scenarios is reduced by approximately CNY 1,068 per occurrence. In practical operation, even under extreme conditions such as a sudden drop in wind speed or a sharp increase in load demand, the operator does not need to purchase electricity at temporarily high prices, because the financial impact is significantly mitigated. For daily operation, the required risk reserve can also be reduced accordingly, suggesting that tail risk has been effectively controlled.

Table 4 compares the proposed method with two external benchmark methods. The typical-scenario scheduling method optimizes the system under only a single representative scenario, while the conventional two-stage optimization method mitigates the impact of uncertainty by introducing recourse decisions. In contrast, the proposed Wasserstein-DRO-CVaR method reduces the out-of-sample average cost by 2.44%, decreases the CVaR by 6.72%, and lowers the expected load curtailment from 35.8 kWh to 18.6 kWh. These results demonstrate that the proposed method not only reduces the operating cost, but more importantly, suppresses tail-risk losses and improves the operational resilience of the system under renewable energy uncertainty.

Figure 5 shows the convergence characteristics of the proposed model under the column-and-constraint generation algorithm. The horizontal axis denotes the iteration number, the left vertical axis is the objective value of the master problem, and the right vertical axis is the worst violation value. The red dashed line is the convergence threshold. From the initial iterations, the number of typical scenarios is small, and the master problem cannot sufficiently describe the worst case, so the optimization result is optimistic. As the iteration number increases, key scenarios with large weights are added, the feasible region of the master problem is tightened, and the objective value shows a trend of rapid increase first and then a slow rise. During the iterations, the master problem obtains more information about the worst distribution, and the violation value gradually decreases. When the violation value approaches the threshold, the subsequent improvement is very limited, indicating that the solution tends to be stable. The curve decreases smoothly without oscillation or abrupt jumps, which fully verifies the stability and robustness of the iterative optimization process.

Figure 6 shows the hourly electric power balance under the typical scenario. The positive stacked bars denote the supply-side outputs, and the negative stacked bars denote the electricity consumption. During the night and early morning (t = 0–5), the load is relatively low; the system is mainly supplied by airborne wind energy and the gas turbine. The demand response mechanism appropriately increases the effective load and actively shifts transferable loads to valley periods. Part of the surplus electricity is stored by the energy storage system to improve resilience against uncertainty and high loads. During the morning (t = 6–10), the base load increases significantly and residential electricity use increases; the energy storage system starts discharging to offset the pressure of load growth. Under weak-grid conditions, the grid structure is weak, and electricity purchase/sale behavior is constrained, so electricity purchase occurs only in a short period to provide a small amount of supplementary electricity during the steepest load growth. The demand response mechanism shifts out a small amount of load at (t = 9), which simulates the discharge behavior of the virtual battery. PV starts generating at (t = 8), and the output increases significantly, effectively alleviating the system pressure during the morning load increase. The airborne wind unit adjusts its altitude appropriately at (t = 8–9) to maintain stable output. During midday (t = 11–13), PV output reaches its peak and solar and wind resources are abundant. Surplus electricity mainly has three destinations: part is stored by the energy storage system, part is consumed by the simulated discharge behavior of the demand response mechanism, and part is converted into heat by the electric boiler; a small amount of electricity is sold. During the afternoon (t = 13–15), the system adjusts the energy storage plan to store a small amount of surplus electricity to cope with the evening peak-load pressure. During the evening and night (t = 16–24), PV output gradually decreases to zero; the airborne wind unit performs altitude adjustments and continues to maintain stable output. The energy storage system and the simulated discharge behavior of demand response cooperate to reduce the load peak and relieve the system pressure during the most stressed periods.

Figure 7 shows the hourly thermal power balance under the typical scenario. The positive stacked bars denote the supply-side outputs, and the negative stacked bars denote heat consumption. Hami, Xinjiang has abundant geothermal resources and provides heating throughout the day. The marginal cost of geothermal heat is lower than that of electricity-to-heat conversion by the electric boiler and external heat purchase in most periods, so geothermal output accounts for a large share of the supply side and is relatively stable. The waste heat boiler absorbs the hot exhaust gas of the gas turbine for heating and cooperates with geothermal heat to form a stable output block. At night (t = 0–5), the thermal load does not fluctuate strongly, and a small amount of surplus heat is stored in the thermal storage system. In the morning (t = 6–11), the thermal load gradually increases; the heat stored at night is released in time, and the system purchases part of the heat to alleviate the pressure of supply-side power ramping. During midday (t = 11–13), renewable resources are abundant, and a small amount of electricity is converted into heat by the electric boiler to participate in the heat supply. In the afternoon (t = 14–17), the load level gradually decreases; geothermal supply remains relatively stable; surplus heat is partly sold and partly absorbed by the thermal storage system. During evening and night (t = 18–23), the thermal storage system releases a small amount of heat; geothermal heat and the waste heat boiler basically meet the heat demand, so there is no heat purchase in this period.

Figure 8 show the charging/discharging power and SoC trajectories of the electric energy storage system and the thermal energy storage system under the typical scenario, respectively. For the electric storage system, during (t = 0–5), the supply is sufficient, and the energy-use level is low; while satisfying power balance, the system stores surplus electricity. During (t = 6–9), the load level increases, but PV and wind resources are abundant and supplement the supply gap of conventional units, so only a small amount of discharge occurs. During (t = 11–13), the storage system absorbs surplus renewable electricity. During (t = 16–21), PV output decreases to zero, and the outputs of airborne wind energy and the gas turbine are lower than the demand; the storage system undertakes the main peak-shaving role to ensure stable operation. Compared with electric storage, the thermal storage system has smaller power fluctuation and a smoother operating curve. This is because the main heat supply is provided by abundant geothermal resources, and thermal load and heating processes have stronger thermal inertia. Compared with electrochemical battery storage, the adjustable power ramp rate of thermal storage is more gradual. Overall, the energy storage system basically follows the operating mechanism of absorbing energy in load valley periods and releasing energy in load peak periods. The two storage systems cooperate to reduce the dependence on external energy, relieve peak pressure, and improve operation robustness under uncertainty.

Figure 9 shows the operating power and SoC trajectory of the virtual-battery-based demand response mechanism under the typical scenario. Load Shifting-in indicates that transferable electricity use is arranged in advance or compensated when the system pressure is low; Load Shifting-out indicates that transferable flexible electricity is temporarily reduced or delayed when the pressure is high, thereby reducing the effective load. The introduction of virtual SoC gives an energy-consistent physical interpretation for demand response behaviors: load shifting-in corresponds to battery charging, and load shifting-out corresponds to battery discharging. From the overall curve, the variation of operating power is consistent with the peak-shaving and valley-filling logic of the energy storage system. Load shifting-in is mainly distributed at the beginning of the day and near windows with renewable surplus, driving the virtual SoC to increase, indicating that flexible energy is stored in advance. Concentrated load shifting-out occurs in the high-pressure window in the evening and at night, indicating that this unit is used to directly reduce the effective load when the power side is tight. In system modeling, this mechanism can provide fast, flexible resources, adjust the net load curve, relieve peak pressure after renewable output decreases, and does not require capacity investment cost, thereby effectively improving the overall economic performance of the system.

The contribution of different types of flexible loads to the equivalent energy storage capacity in this study is presented in Table 5. The flexible loads include thermostatically controlled loads, domestic hot water loads, cooling loads, water pump loads, and deferrable residential electricity loads. According to the scheduling results, the average upward demand response power over the 24-h scheduling horizon is 38 kW, and the total curtailed energy is 912 kWh. Among these loads, thermostatically controlled loads make the largest contribution, with a curtailed energy of 356 kWh, accounting for approximately 39.03% of the total reduction. This indicates that thermostatically controlled loads, such as heating and air-conditioning loads, are the main source of demand-side flexibility. In contrast, deferrable residential electricity loads contribute 91 kWh, accounting for approximately 9.98%. Their flexibility is relatively limited because the start-up time of deferrable appliances is strongly dependent on user behavior, making it difficult to achieve large-scale concentrated reductions during peak periods as in the case of thermostatically controlled loads.

Figure 10 presents the sensitivity analysis of the CVaR confidence level α. Seven groups of comparative experiments are set, and α is uniformly sampled in [0.80,0.95]. In the two-stage solution process, a larger α is equivalent to penalizing the average loss over the worst-case scenario set, meaning that the model focuses more on protecting extreme tail scenarios, and the dispatch strategy becomes more conservative. The bar chart decomposes the total cost of the objective function into three parts: stage-1 cost, stage-2 revenue, and stage-2 adjustment cost. The dashed line is the tail-risk indicator, which measures the CVaR value of losses. The results show that, as α increases, the solution becomes more conservative. When α increases from 0.8 to 0.95, the stage-1 cost increases from 25,218.9 to 25,888.10, with an increase of about 2.65%, indicating that stage 1 tends to invest more resources, such as increasing the gas turbine output trajectory. Stage-2 revenue gradually decreases from 5520.08 to 4909.16, with a decrease of about 11.07%, indicating that the model reduces reliance on market opportunity revenue, such as reducing energy sales, converging sensitive behaviors to energy prices, and instead allocating energy to hedge load pressure. In addition, when α is in a higher interval, such as [0.925,0.95] in the figure, the marginal benefit decreases. In this interval, further increasing α can reduce CVaR, but the economic cost to reduce this risk increases, indicating that when the strategy is already sufficiently conservative, improving tail protection requires disproportionate cost sacrifice.

Figure 11 presents the sensitivity analysis of the relative radius parameter ε in the Wasserstein-DRO framework. ε is uniformly sampled in [0.05,0.35]. A larger ε means a larger ambiguity set around the empirical scenario distribution, i.e., a stronger possible deviation of the true distribution from the empirical distribution. The results show that, with increasing ε, the stage-1 cost and stage-2 adjustment cost gradually increase, stage-2 revenue gradually decreases, and the tail risk decreases from 2745.1 to 715.6, with a decrease of about 73.93%. This indicates that, as the ambiguity set expands, extremely high-loss outcomes are effectively suppressed, and some sacrifice of economic performance is exchanged for faster convergence of the tail-risk upper bound.

Figure 12 shows the three-dimensional surface of time, altitude, and power upper envelope of airborne wind energy. The altitude range of the airborne platform is set to [500,1000]. The upper-envelope power is obtained by the inner-loop model. The outer-loop model embeds this feasible region into the scheduling model, sets climbing and descending power limits of airborne wind energy, and adds the regulation energy consumption in the objective function.

Figure 13 presents the error validation of the AWE power envelope, where the absolute normalized error is 0.84%. The reference power is obtained from the inner-loop dynamic model under different combinations of wind speed and operating altitude, while denotes the envelope-based available power used in the outer-loop optimization. The MAE and RMSE over the altitude range of 500–1000 m are calculated, with values of 2.13 kW and 4.68 kW, respectively. The above three indicators all remain within a reasonable range.

Figure 14 validates the adaptability of the proposed framework under different resource structures by varying the renewable energy penetration level. Five renewable energy penetration levels, namely, 40%, 50%, 60%, 70%, and 80%, are considered, and the proposed method is compared with a conventional two-stage optimization benchmark. The results show that, as the renewable energy penetration increases from 40% to 80%, the system gradually reduces its dependence on gas turbines and purchased electricity, leading to a decrease in the total operating cost. This is because, with a higher level of available renewable generation, the system can still retain a certain amount of usable renewable output even under adverse scenarios, thereby reducing the recourse cost in tail-risk scenarios. Compared with the conventional two-stage optimization method, the proposed method achieves lower operating costs under all penetration levels, indicating that it maintains good economic performance under different renewable energy penetration conditions. In particular, when the renewable energy penetration is high and the associated uncertainty becomes more significant, the advantage of the proposed method is more pronounced.

5 Conclusion

This paper proposes a distributionally robust optimal scheduling method for weak-grid integrated energy systems that incorporates airborne wind energy and virtual battery-based demand response. By embedding CVaR into the objective function, the method achieves a precise trade-off between scheduling costs and extreme operational risks. Numerical results on a township-level energy system in Hami, Xinjiang, China, show that the proposed approach significantly reduces total operating costs while effectively mitigating tail risks to energy supply security.

Firstly, a full-dimensional physical model of the integrated energy system including airborne wind energy is constructed, together with an upper envelope model of available power output that accounts for altitude adjustment and energy consumption. The platform operating height is treated as an optimizable decision variable, thereby fully unlocking its active dispatch potential. With the introduction of airborne wind energy, the total system cost is reduced by 3.67% compared with the conventional wind power scheme.

Secondly, a virtual battery-based demand response model is developed. Adjustable loads are equivalently modeled as virtual battery charging and discharging processes, with SOC used to characterize cumulative peak-shifting capacity and sustainable response margin. Power, capacity, and mutual exclusion constraints enable coordinated optimization with physical storage within the same framework. The inclusion of this model further reduces total system cost by 10.3%.

Finally, a two-stage distributionally robust optimization framework based on Wasserstein distance is established. CVaR is incorporated into the objective function for tail risk quantification, and the model is solved efficiently via the C&CG algorithm, thereby achieving a precise balance between economic performance and risk control.

Although the present study delivers the above improvements, certain limitations remain: the model is developed for a single-region scenario and does not yet incorporate exogenous uncertainties such as market mechanisms or policy factors. Future work will develop an integrated optimization framework that jointly optimizes capacity planning and operational scheduling, thereby forming a more engineering-oriented planning–scheduling co-optimization approach tailored to weak-grid security requirements.

References

• BuY.YuH.JiH.SongG.XuJ.LiJ.et al (2024). Hybrid data-driven operation method for demand response of community integrated energy systems utilizing virtual and physical energy storage. Appl. Energy366, 123295. 10.1016/j.apenergy.2024.123295

• CabotC.VillavicencioM. (2024). The demand-side flexibility in liberalised power market: a review of current market design and objectives. Renew. Sustain. Energy Rev.201, 114643. 10.1016/j.rser.2024.114643

• ChawdaG. S.ShaikA. G.SuW. (2024). Efficient wind energy integration in weak AC grid with a DLMF-Based adaptive approach. Appl. Energy372, 123779. 10.1016/j.apenergy.2024.123779

• ChenH.WangD.LinC.ZhangR.JiangT.LiX. (2025). Risk-averse optimal scheduling of regional electricity-heating integrated energy system considering interface with microgrids. CSEE J. Power Energy Syst.11, 2117–2127. 10.17775/CSEEJPES.2021.07120

• ChengY. (2025). World’s first megawatt-class commercial airborne wind power system completes a Gobi-condition test flight. China News Serv. (Xinjiang). Available online at: https://www.xj.chinanews.com.cn/dizhou/2025-09-25/detail-ihevmzis5432827.shtml (Accessed September 25, 2025).

• DuY.LiX.LiangY.TanZ. (2024). Two-stage multi-objective distributionally robust optimization of the electricity–hydrogen coupling system under multiple markets. Energy303, 131960. 10.1016/j.energy.2024.131960

• FagianoL.QuackM.BauerF.CarnelL.OlandE. (2022). Autonomous airborne wind energy systems: accomplishments and challenges. Annu. Rev. Control Robot. Auton. Syst.5, 603–631. 10.1146/annurev-control-042820-124658

• FanG.PengC.WangX.WuP.YangY.SunH. (2024). Optimal scheduling of integrated energy system considering renewable energy uncertainties based on distributionally robust adaptive MPC. Renew. Energy226, 120457. 10.1016/j.renene.2024.120457

• FuQ.XingZ.ZhangC.XuJ. (2024). A review and prospective study on modeling approaches and applications of virtual energy storage in integrated electric–thermal energy systems. Energies17, 4099. 10.3390/en17164099

• GaoX.WangS.SunY.ZhaiJ.ChenN.ZhangX.-P. (2024). Low-carbon energy scheduling for integrated energy systems considering offshore wind power hydrogen production and dynamic hydrogen doping strategy. Appl. Energy376, 124194. 10.1016/j.apenergy.2024.124194

• GaoY.TahirM.SianoP.BiY.HuS.YangJ. (2025). Optimization of renewable energy-based integrated energy systems: a three-stage stochastic robust model. Appl. Energy377, 124635. 10.1016/j.apenergy.2024.124635

• GeS.LiJ.HeX.LiuH. (2021). Joint energy market design for local integrated energy system service procurement considering demand flexibility. Appl. Energy297, 117060. 10.1016/j.apenergy.2021.117060

• HanF.ZengJ.LiZ.ZhangH.LiY.WangZ. (2023). A novel two-layer nested optimization method for a zero-carbon island integrated energy system, incorporating tidal current power generation. Renew. Energy218, 119296. 10.1016/j.renene.2023.119296

• HaoJ.YangY.XuC.DuX. (2022). A comprehensive review of planning, modeling, optimization, and control of distributed energy systems. Carbon Neutrality1, 28. 10.1007/s43979-022-00029-1

• HassankashiA.DiniA.PirouziS.VeisiM.BahreiniM. (2026). Flexible renewable integrated energy system capabilities to improve voltage stability with power quality and economic environmental operation of smart grid. Sci. Rep.16, 99. 10.1038/s41598-025-29052-9

• HerdingR.RossE.JonesW. R.EndlerE.CharitopoulosV. M.PapageorgiouL. G. (2024). Risk-aware microgrid operation and participation in the day-ahead electricity market. Adv. Appl. Energy15, 100180. 10.1016/j.adapen.2024.100180

• International Renewable Energy Agency (IRENA) (2025). Record-breaking annual growth in renewable power capacity. Available online at: https://www.irena.org/News/pressreleases/2025/Mar/Record-Breaking-Annual-Growth-in-Renewable-Power-Capacity (Accessed July 23, 2026).

• KhurshidY.PaulP.ElhesasyM.AliB.KamraM. M.OkashaM.et al (2025). From inception to commercialization: a systematic review of airborne wind energy systems. Sustain. Energy Technol. Assess.83, 104623. 10.1016/j.seta.2025.104623

• LiH.LiuH.MaJ.ZhangW. (2023). Distributionally robust optimal dispatching method of integrated electricity and heating system based on improved wasserstein metric. Int. J. Electr. Power Energy Syst.151, 109120. 10.1016/j.ijepes.2023.109120

• LiX.TianZ.WuX.FengW.NiuJ. (2024). Optimal planning for hybrid renewable energy systems under limited information based on uncertainty quantification. Renew. Energy237, 121866. 10.1016/j.renene.2024.121866

• LiF.WangD.GuoH.ZhangJ. (2024). Distributionally robust optimization for integrated energy system accounting for refinement utilization of hydrogen and ladder-type carbon trading mechanism. Appl. Energy367, 123391. 10.1016/j.apenergy.2024.123391

• LinB.HuangY.YuD.FuC.ChenC. (2025). Stochastic distributionally robust optimization scheduling of high-proportion new energy distribution network considering detailed modeling of energy storage. Processes13, 2230. 10.3390/pr13072230

• LiuH.FanZ.XieH.WangN. (2022). Distributionally robust joint chance-constrained dispatch for electricity–gas–heat integrated energy system considering wind uncertainty. Energies15, 1796. 10.3390/en15051796

• MaH.ChenQ.HuB.SunQ.LiT.WangS. (2021). A compact model to coordinate flexibility and efficiency for decomposed scheduling of integrated energy system. Appl. Energy285, 116474. 10.1016/j.apenergy.2021.116474

• McKennaR.LilliestamJ.HeinrichsH. U.WeinandJ.SchmidtJ.StaffellI.et al (2025). System impacts of wind energy developments: key research challenges and opportunities. Joule9, 101799. 10.1016/j.joule.2024.11.016

• SchwidtalJ. M.PicciniP.TronciaM.ChitchyanR.MontakhabiM.FrancisC.et al (2023). Emerging business models in local energy markets: a systematic review of peer-to-peer, community self-consumption, and transactive energy models. Renew. Sustain. Energy Rev.179, 113273. 10.1016/j.rser.2023.113273

• SunH.ZhaiH.WuX. (2021). Research and application of multi-energy coordinated control of generation, network, load and storage. Trans. China Electrotech. Soc.36, 3264–3271.10.19595/j.cnki.1000-6753.tces.200188

• SunQ.FuY.LinH.WennerstenR. (2022). A novel integrated stochastic programming–information gap decision theory (IGDT) approach for optimization of integrated energy systems (IESs) with multiple uncertainties. Appl. Energy314, 119002. 10.1016/j.apenergy.2022.119002

• Vahedipour-DahraieM.Anvari-MoghaddamA.GuerreroJ. M. (2018). Evaluation of reliability in risk-constrained scheduling of autonomous microgrids with demand response and renewable resources. IET Renew. Power Gener.12, 657–667. 10.1049/iet-rpg.2017.0770

• WangW.MaH.WangB.ZhengJ.LiuZ.LiT.et al (2026). Scheduling of electricity-hydrogen integrated system under renewable energy sources uncertainty: non-anticipativity and robust feasibility. Appl. Energy404, 126996. 10.1016/j.apenergy.2025.126996

• WenJ.ChenX.JiaR.JiaoY.CaoG.GuoY.et al (2025). Distributed scheduling model of multi-network coupled energy system considering dynamic charging price mechanism and multiple utilization of hydrogen energy. Energy335, 138006. 10.1016/j.energy.2025.138006

• WuG.WangH.WuQ. (2020). Wind power development in the belt and road area of Xinjiang, China: problems and solutions. Util. Policy64, 101024. 10.1016/j.jup.2020.101024

• XuanA.ShenX.GuoQ.SunH. (2021). A conditional value-at-risk based planning model for integrated energy system with energy storage and renewables. Appl. Energy294, 116971. 10.1016/j.apenergy.2021.116971

• YangS.WuH.SongJ.LiH.ChenH. (2024). Two-stage robust optimization scheduling for integrated energy systems considering ammonia energy and waste heat utilization. Energy Convers. Manag.319, 118922. 10.1016/j.enconman.2024.118922

• ZhangM.WuQ.WenJ.PanB.QiS. (2020). Two-stage stochastic optimal operation of integrated electricity and heat system considering reserve of flexible devices and spatial-temporal correlation of wind power. Appl. Energy275, 115357. 10.1016/j.apenergy.2020.115357

• ZhangR.ChenY.LiZ.JiangT.LiX. (2024). Two-stage robust operation of electricity-gas-heat integrated multi-energy microgrids considering heterogeneous uncertainties. Appl. Energy371, 123690. 10.1016/j.apenergy.2024.123690

• ZhangL.SongW.SunE.ZhangQ.WuD.ChenF.et al (2025). A high-altitude wind resource assessment method for decentralized wind power based on improved linear regression. Renew. Energy238, 121968. 10.1016/j.renene.2024.121968

• ZongX.YuanY. (2023). Two-stage robust optimization of regional integrated energy systems considering uncertainties of distributed energy stations. Front. Energy Res.11, 1135056. 10.3389/fenrg.2023.1135056