Two-Stage Robust Optimization for Large Logistics Parks to Participate in Grid Peak Shaving

: As new energy integration increases, power grid load curves become steeper. Large logistics parks, with their substantial cooling load, show great peak shaving potential. Leveraging this load while maintaining staff comfort, product quality, and operational costs is a major challenge. This paper proposes a two-stage robust optimization method for large logistics parks to participate in grid peak shaving. First, a Cooling Load’s Economic Contribution (CLEC) index is introduced, integrating the Predicted Mean Vote (PMV) and Sales Pressure Index (SPI). Then, an optimization model is established, accounting for renewable energy uncertainties and maximizing large logistics parks’ participation in peak shaving. Results illustrate that the proposed method leads to a reduction in the peak shaving pressure on the distribution network. Specifically, under the scenario tolerating the maximum potential uncertainty in renewable energy output, the absolute peak-to-valley difference and fluctuation variance of the park’s net load are decreased by 45.82% and 54.59%, respectively. Furthermore, the PMV and the SPI indexes are reduced by 39.12% and 26.36%, respectively. In comparison with the determined optimization method, despite a slight cost increase of 20.06%, the proposed method significantly reduces EDR load shedding by 98.1%.


Introduction 1.Background
As the proportion of renewable energy integration increases, the equilibrium and symmetry between the supply side and demand side of the power system become progressively more disrupted.To maximize the utilization of renewable energy resources, it becomes imperative to fully engage the demand-side components.Also, the intensifying urban population aggregation in major cities leads to a steady surge in electricity demand, especially in core urban areas.The distinct counter-peaking attributes associated with the substantial integration of renewable energy sources further widen the peak-to-valley disparities, thus exacerbating localized supply-demand imbalances and posing significant risks to the stable operation of power grids.

Recent Works
The traditional method of managing peak loads, which involves adjustments to power generation, is economically impractical due to its consequences of inefficient plant usage and elevated pollution levels.Specifically, in recent years, the large-scale integration of intermittent and uncertain renewable energy sources, such as wind power and photovoltaics, has complicated the supply-demand balance within the electrical grid, subsequently augmenting the complexity of peak shaving.Under the above conditions, numerous contemporary studies focus on strategies to address the escalating peak shaving demands, particularly examining responses from both the supply and demand sides.
In the electricity market environment, the participation of power plants in peak shaving has significantly increased.Their active engagement in the scheduling and flexible adjustment of power generation has effectively reduced the peak-to-valley difference in the grid, alleviated operational pressures within the power system, and significantly enhanced its safety.Reference [1] proposes a cost analysis model to address the economic and environmental implications of deep peak shaving for coal-fired units.Reference [2] presents an assessment method that reflects the marginal costs of peak and off-peak electricity generation for combined heat and power plants engaged in peak shaving, discovering that the use of electric heat pumps for waste heat recovery achieves a significant depth of peak shaving.Reference [3] examines the complementary characteristics of wind power and nuclear power in peak shaving, proposing a multi-power dispatch model to optimize nuclear power output.Reference [4] addresses the serious issue of wind curtailment, proposing an optimal configuration method for power-to-heat equipment in combined heat and power plants with a focus on static payback period.Reference [5] proposes a solution combining typical peak shaving output curves with water abandonment adjustment strategies to address the challenges of short-term peak shaving operation.In reference [6], the role of carbon capture devices in peak shaving is explored.A virtual energy storage model is constructed, and a joint peak shaving strategy for carbon capture devices and virtual energy storage is proposed.Reference [7] proposes a model for pumped-storage power stations to participate in the peak shaving ancillary service market, helping to share the peak shaving pressure of thermal power units.Reference [8] explores the coordinated operation of hydropower and renewable energy to mitigate fluctuations and facilitate peak shaving, addressing the challenges posed by the intermittency and uncertainty of wind and solar energy to grid dispatch.Reference [9] presents a new peak shaving model that utilizes mixed integer linear programming without presupposing or fixing the peak shaving order of power stations.Reference [10] proposes a multi-objective unit commitment model that combines the concentration of solar power plants and wind farms for peak shaving, taking operational risks into account.
On the other hand, providing compensation to the demand side through contracting can effectively stimulate load-side resources and attract more load-side participation in the peaking process.This mechanism not only significantly reduces the peak-to-valley difference but also further enhances the safety and stability of the power system by optimizing the allocation of power resources.Reference [11] provides a comprehensive overview of peak load reduction strategies, focusing on the impact of three major strategies: demand-side management, energy storage systems, and grid-connected electric vehicles.Reference [12] constructs a demand response model considering dual uncertainty to address carbon emissions and supply-demand balance issues.In reference [13], a demand response energy management strategy adopting a peak rebate program is proposed.Reference [14] presents a joint peak shaving strategy for carbon capture equipment and virtual energy storage.It establishes a two-stage peak shaving model for day-ahead and intra-day operations, aimed at minimizing the peak-to-valley difference in the load and reducing system operating costs to ensure economically efficient peak shaving.Reference [15] explores the application of demand-side management in air conditioning and heat pumps, particularly focusing on the heat pump sector.Reference [16] proposes a community-based home energy management system that utilizes a particle swarm optimization algorithm and user-defined constraints to achieve peak load reduction.Reference [17] addresses peak load issues in district heating systems by employing a differential return water temperature adjustment strategy for peak shaving.The aforementioned research indicates that the linkage between the power generation side and the demand side can effectively coordinate, flattening the load curve.As a significant power load, the fluctuation of chillers significantly affects the safe and stable operation of the grid.Consequently, incorporating chillers into peak shaving practices is of great significance for enhancing the reliability and economy of power supply.Large chillers, due to their inherent controllability, can effectively reduce the load pressure during peak hours through reasonable control strategies while ensuring the satisfaction of daily living and production needs.On the other hand, buildings equipped with chillers have a certain cold storage capacity, allowing for minor changes in operating status without significantly affecting indoor environmental temperature, demonstrating a high level of inertia.In large cities, air conditioning loads during summer peak hours account for a significant proportion of the total load, becoming a major contributor to peak loads, with significant adjustment potential.Reference [18] proposes a state queueing model that utilizes a bidirectional information channel to collect thermostatically controlled load information, and employs tracking curves and state prediction methods for the control of thermostatically controlled load operations.Reference [19] converts the thermodynamic equations of air conditioning loads into a finite-dimensional state-space model through the finite difference method.Additionally, reference [20] develops a stochastic energy storage charging parameter model for aggregated air conditioning loads.Reference [21] presents a method for adjusting the power of aggregated air conditioning systems through broadcasting temperature setpoints.The study in [22] proposes a comfort-constrained heat pump load management method and realizes closed-loop control.

Motivations and Contributions
While the above methods have achieved the comprehensive optimization and regulation of the chiller, most of them targeted distributed and autonomous air conditioning loads.Currently, there is no large-scale practice involving large logistics parks in peak shaving, mainly due to three reasons: (1) There is a concern that power outages or insufficient power supply during peak shaving may lead to the damage of frozen products within the logistics park.Cold storage facilities require continuous refrigeration to maintain the freshness of frozen products, and any disruption to their stable operation caused by peak shaving may result in product losses.(2) Logistics parks house a large number of office staff, whose office buildings are also cold-consuming structures with high power consumption.However, their demand for cooling loads differs significantly from that of cold storage facilities.It is necessary to establish separate mathematical models for the comfort level of office staff and the storage duration of frozen products, and incorporate them into a comprehensive consideration framework under the same dimension.(3) With the increasing penetration of renewable energy into distribution networks, the difficulty of forecasting new energy sources has increased.The traditional "forecast value as the planned value" approach is hardly applicable, and the challenge of peak shaving for cold storage facilities under such uncertainties has significantly increased, potentially exposing logistics parks to additional risks of frozen product damage.Therefore, a robust optimization approach is needed to mitigate the impact of random errors.
To solve the problems above, this paper proposes a novel optimization method for large logistics parks to participate in grid peak shaving in a two-stage robust manner.The following contributions correspond to the motivations for this study: (1) The inertial traits inherent in the substantial cooling load within the park are exploited to engage in peak shaving, while adhering to temperature constraints.The formulated objective function factors in the comfort of the office setting within the park and the sales pressure associated with cold storage, thereby achieving an optimal balance between operational efficiency, temperature tolerance, and peak shaving engagement.
(2) A composite metric, termed Cooling Load's Economic Contribution (CLEC), is devised.This metric integrates both the Predicted Mean Vote (PMV) for human comfort and the Sales Pressure Index (SPI) for frozen goods, translating these factors into economic costs for a unified optimization process.
(3) A two-stage robust optimization model, encompassing both day-ahead and intraday planning, is devised.This model aims to maximize the park's cooling load participation in the peak shaving efforts, while maintaining personnel comfort and ensuring the quality of frozen products.Additionally, it accounts for the uncertainties introduced by renewable energy sources like wind and solar.
The structure of this paper is organized as follows.Section I introduces the research background and primary research focus.Section II focuses on modeling a large logistics park, establishing the CLEC index to ensure that the park participates in peak shaving while maximizing the normal operational demand for cooling load.Section III provides a twostage robust optimization model for a large logistics park's participation in peak shaving, covering both day-ahead and intra-day periods.Section IV introduces the methodology to solve the problem using the Column and Constraint Generation (C&CG) algorithm.Section V outlines a group of case studies, creating robust optimization scenarios that consider various levels of uncertainty.Pros and cons of the proposed method are compared and summarized.Section VI presents the conclusions of this paper.

Total Framework with Diverse Cooling Loads
Figure 1 depicts the park model constructed in this paper.The power supply consists of photovoltaic and wind power, two representative renewable energy sources, integrated into the park in a decentralized manner.The park is connected to the distribution network through a PCC node, which includes a group of gas turbines that compensates for any shortcomings in the park's renewable energy generation capacity.This gas turbine unit also has reserve capacity, providing backup when the power transmission pressure between the park and the distribution network becomes excessive.Meanwhile, the park has a significant amount of cooling load to meet the needs of frozen goods storage and the cooling demands of staff within the buildings.Traditionally, parks control the output of chillers by adjusting the temperature inside buildings to meet cooling load requirements.However, buildings possess a "cold storage" characteristic, and to avoid "wasting" this stored cold inertia, the CLEC index is employed in the park model, shown in Figure 1, to intelligently regulate the output characteristics of the chillers.In this paper, EDR is activated when the peak shaving pressure on the distribution network becomes too great, and the central optimization controller removes some of the park's electrical load.
background and primary research focus.Section II focuses on modeling a large logistics park, establishing the CLEC index to ensure that the park participates in peak shaving while maximizing the normal operational demand for cooling load.Section III provides a two-stage robust optimization model for a large logistics park's participation in peak shaving, covering both day-ahead and intra-day periods.Section IV introduces the methodology to solve the problem using the Column and Constraint Generation (C&CG) algorithm.Section V outlines a group of case studies, creating robust optimization scenarios that consider various levels of uncertainty.Pros and cons of the proposed method are compared and summarized.Section VI presents the conclusions of this paper.

Total Framework with Diverse Cooling Loads
Figure 1 depicts the park model constructed in this paper.The power supply consists of photovoltaic and wind power, two representative renewable energy sources, integrated into the park in a decentralized manner.The park is connected to the distribution network through a PCC node, which includes a group of gas turbines that compensates for any shortcomings in the park's renewable energy generation capacity.This gas turbine unit also has reserve capacity, providing backup when the power transmission pressure between the park and the distribution network becomes excessive.Meanwhile, the park has a significant amount of cooling load to meet the needs of frozen goods storage and the cooling demands of staff within the buildings.Traditionally, parks control the output of chillers by adjusting the temperature inside buildings to meet cooling load requirements.However, buildings possess a "cold storage" characteristic, and to avoid "wasting" this stored cold inertia, the CLEC index is employed in the park model, shown in Figure 1, to intelligently regulate the output characteristics of the chillers.In this paper, EDR is activated when the peak shaving pressure on the distribution network becomes too great, and the central optimization controller removes some of the park's electrical load.
In this paper, the park's central optimization controller takes into account the park's operating costs, the comfort level of cooling load users, the need to ensure the quality of frozen goods in the park, and the requirement to participate in peak shaving.In this paper, the park's central optimization controller takes into account the park's operating costs, the comfort level of cooling load users, the need to ensure the quality of frozen goods in the park, and the requirement to participate in peak shaving.
The ambient temperature of a chiller can be adjusted by modifying its cooling load (CL), thus varying the cooling capacity.Both humans and warehouses have specific acceptable temperature ranges, with warehouses favoring lower temperatures for longer food storage.As chiller-driven buildings have some heat storage capacity, CL changes do not immediately affect indoor temperature [23].Assuming that the heat balance between the CL and the building is maintained for a short period when the outdoor temperature is constant, the electric power consumed by the chiller is where T in represents the indoor temperature, T out is the outdoor temperature, and λ is the energy efficiency ratio.R 1 is the equivalent heat resistance.A cold-storage building model can describe its temperature dynamic characteristics via a mathematical model that gives the relationship between indoor and outdoor cold and heat sources and room temperature changes.Its temperature dynamic characteristics can be described by an equivalent thermal parameter (ETP) model, as shown in the following equation [24]: where Q L,t is the total cooling load in the duration of t.R and C are the indoor equivalent thermal resistance and equivalent thermal capacity of cold-storage buildings, respectively.
T n,t and T W,t are the indoor and outdoor temperatures, respectively.∆t is the dispatching duration.Depending on the needs of the user, the indoor temperature is usually limited to a specific range [25], as shown in the following equation:

Modeling of Diverse Cooling Load's Economic Contribution
Thermal comfort, which encapsulates individuals' subjective assessments and sensations regarding their thermal environment, serves as the primary indicator of users' ambient temperature preferences.To objectively quantify the effect of temperature on user comfort, this study employs the Predicted Mean Vote (PMV) index [26], as detailed in Table 1.The mathematical form is given in Equation (4).
where M and W represent the human energy metabolic rate and mechanical power, respectively.f cl is the ratio between the clothed body surface area and the naked body surface area.h c is the convective heat transfer coefficient.P a is the water vapor pressure of the air around the human body.t a , t r and t cl represent the air temperature around the human body, mean radiant temperature, and clothing surface temperature, respectively.This paper mainly focuses on the cooling capacity, and temperature is the most intuitive perception of indoor thermal comfort for humans.Therefore, except for the air temperature around the human body, t a , it is assumed that other parameters are given values.For the chiller-driven warehouses, there are upper and lower limits for the ambient temperature.For instance, frozen foods need to be stored below −18 • C, while refrigerated medicines require storage in a cold storage facility between 2 • C and 10 • C [27].Additionally, the optimal refrigeration temperature for different frozen products should be considered.To assess the negative impact caused by deviations from the optimal refrigeration temperature, a Sales Pressure Index (SPI) is established.This index primarily focuses on the quality and freshness retention of frozen products.A lower SPI indicates that the actual temperature of the cold storage is closer to the ideal storage temperature, thereby reducing product quality loss and sales pressure.In this paper, we only focus on the value of frozen products at the end of the dispatching period, assuming no transfer of frozen products during this period.The mathematical equation for the SPI at the end of the period is as follows: where λ SPI represents the comprehensive SPI of the cold storage within the microgrid at the end moment; the superscript i denotes the type of frozen product; T i r,t and T i ideal,t are the actual and ideal storage temperatures, respectively, for type i frozen products at time t; and α i is the sales weighting coefficient for type i frozen products.T dev max is the maximum permitted temperature bias for frozen goods.The combination of PMV and SPI modeling gives rise to the definition of CLEC, denoted as I CLEC , which quantifies the economic impact of the microgrid's cooling load.
In Equation ( 6), a t signifies the temperature sensitivity coefficient of the building's cooled occupants at a specific time t.This coefficient varies temporally, peaking during typical rest periods such as 12:00 to 14:00 and 21:00 to 24:00.The term |I PMV,max | represents the absolute maximum PMV index value recorded among the cooled individuals within the premises.N T stands for the predetermined scheduling duration.Additionally, ω pmv and ω spi are the weighted coefficients assigned to the cooling loads associated with the building's occupants and the frozen items in cold storage, respectively, during the park's operational hours.Meanwhile, T out and T opt refer to the current outdoor temperature and the ideal outdoor temperature, respectively.k serves as an adjustment coefficient that regulates the sensitivity of the I CLEC index to fluctuations in external temperature, thereby controlling its responsiveness to temperature variations.
The additional cost for using cooling loads is defined as where F CLEC represents the cost conversion coefficient that comprehensively considers the loss of human comfort and the sales value of frozen goods.Meanwhile, C CLEC denotes the additional economic costs incurred due to the inefficient use of cooling load within the park.

Methodology
To alleviate the peak shaving pressure on the grid side and leverage the significant inertial cooling load present within the park, this study explores the utilization of economic incentives to motivate the park's participation in the grid's peak shaving initiatives.Specifically, the park receives economic compensation for actively engaging in peak shaving, implemented through a two-stage robust optimization framework that encompasses both day-ahead and intra-day timeframes.
During the day-ahead phase, the system computes the minimum operating costs and the associated peak shaving compensation, leveraging predicted values of photovoltaic and wind power generation.Moving into the intra-day phase, the framework incorporates backup resources from the distribution grid to mitigate any uncertainties within predefined uncertainty sets, thereby ensuring consistent, economical, and reliable system operation.The objective is to minimize daily operating costs while simultaneously reducing the peak shaving pressure to its lowest possible level.The model's objective function is articulated as follows: where C 0 and C S denote the optimization objectives for the initial and subsequent stages, respectively, while U signifies the comprehensive uncertainty set, encompassing uncertainties inherent in photovoltaic and wind power generation.
In Equation ( 9), N T represents the time period and C G(t) denotes the power distribution network's output cost during period t.C pls represents the economic compensation received by the park for participating in peak shaving.
(1) Gas Turbine The power supply from the distribution network to the park primarily relies on gas turbines.The generation cost of these turbines can be represented by a linear function, taking into account the reserve capacity.Therefore, the cost function of the gas turbine is where C G (t) indicates the generation cost of the generator during period t; ∆t is the scheduling step size, set to 1 h.N G stands for the number of generators.F G represents the cost coefficient of the gas turbine.The superscript 0 designates the day-ahead stage.P 0 g,t is the output power of the g-th generator during period t; C R+ g and C R− g are the upward and downward reserve cost coefficients of the generator, respectively; and R + g,t and R − g,t are the upward and downward reserve capacities provided by the generator during period t.Reserve capacity constraints and output power limitation constraints are taken into account, as shown in Equations ( 11)-( 14): where R + max and R − max represent the maximum upward and downward reserve capacities that the unit can provide when transitioning from the first stage to the second stage.Meanwhile, R + min and R − min denote the minimum upward and downward reserve capacities required by the system from the unit.P max g and P min g indicate the maximum and minimum output powers of the gas turbine, which are limited by the rated power and minimum load rate of the generator set, respectively.
(2) Evaluation of peak shaving performance The effectiveness of peak shaving can be measured by the economic compensation C fls obtained by reducing the fluctuation in the grid's net load and the peak-valley difference.C fls can be expressed as ∆P pd = P be pd − P pd (16) ∆P var = ∆P be var − ∆P var (17) where ∆p pd represents the reduction in the peak-valley difference (PVD) of the net load; F pd is the compensation coefficient for peak-valley difference; ∆P var denotes the decrease in fluctuation variance of the net load; and F var is the compensation coefficient for variance.
The superscript indicates the state of the grid's net load before the park participates in peak shaving.P pd and P var are the PVD and variance of the microgrid's net load after peak shaving.P pd = P max − P min (18) In Equation ( 18), P max is the maximum value of the net load during the scheduling period and P min is the minimum value.The absolute PVD P pd reflects the difference between the extreme values of the load peaks and valleys.
In Equation ( 19), N T is the total number of samples and P l.t represents the load value of the t-th sample.The fluctuation variance P var of the net load reflects the degree of load dispersion.A smaller fluctuation variance of the net load indicates a lower degree of dispersion.

Constraints
(1) Power-balance constraints In Equation (20), G b , W b , and P b represent the sets of generators, wind farms, and photovoltaic plants connected to node b, respectively.The notations l = o1_b and l = o2_b denote the feeder-injecting and exporting power from node b, respectively.P 0 W,t and P 0 PV,t refer to the actual consumption of wind and photovoltaic power at time t.P 0 l,t represents the transmission power of line l at time t.L b,t signifies the load connected to node b at time t, which comprises both the base electric load and the load from the chillers.
(2) Cooling load constraints Herein, Q AC,t represents the cooling power of the electric chiller in time period t, while T min in and T max in are the lower and upper limits of the indoor temperature, respectively.The indoor cooling constraints also include Equation (2), which describes the ETP model of the thermal inertia of the cooling system.
(3) Chiller output constraints P AC,t represents the electric power consumed by the electric chiller in time period t, η AC stands for the coefficient of performance (COP) of the electric chiller, and Q min AC /Q max AC represent the lower/upper limits of the output of the electric chiller.
(4) Constraints on power flow In the above equation, θ 0 o1_l and θ 0 o2_l represent the voltage phase angles at the ingress and egress nodes of feeder l, respectively; x l is the reactance of transmission feeder l; P max l is the upper limit of power flow in feeder l, with positive and negative values indicating direction; and θ min and θ max are the minimum and maximum limits of the node voltage phase angles, respectively.
(5) Constraints of renewable energy outputs In the above equation, A 0 PV,t and A 0 W,t are the predicted outputs of the photovoltaic and wind power in the time period t.

Uncertainty Set
In this paper, the uncertainty set in robust optimization consists of uncertainties in the photovoltaic and wind power output, which can be expressed by Equations ( 30)- (35).Here, N T represents the scheduling period of this paper, which is 24 h.Z + W,t and Z − W,t are the status flags for upward and downward fluctuations of wind farm W in period t, respectively.When the value is 1, it indicates the existence of fluctuations, and when it is 0, there are no fluctuations.A S W,t represents the second-stage wind power prediction, which is composed of the first-stage wind power forecast superimposed with power fluctuations.u t,pv and u t,w are the robustness indicators of this paper, representing the time uncertainty limits of photovoltaic and wind power output, respectively.
In the above equation, the superscript S represents the second intra-day stage; C awp refers to the economic loss cost coefficient caused by the failure of the actual wind power output to meet the expected value in the second stage; C S,R+ g and C S,R− g represent the upward and downward reserve dispatch rates of the unit in the second stage, respectively; and R S+ g and R S− g represent the actual upward and downward reserve dispatch capacities of the unit in the second stage.
where P S g,t represents the actual output of the gas turbine in the second stage during time period t, which is composed of the actual output of the generator in the first stage during time period t and the actual reserve capacity dispatched for the unit during time period t.This can be expressed as (2) Feeder power flow constraints −P max l ≤ P S l,t ≤ P max l ∀l, t (40) (3) Emergency demand response constraint where L S,EDR max represents the maximum amount of emergency demand response that each node in the system can participate in during a given time period.
(4) Constraints of renewable energy outputs The above indicates that the actual output of photovoltaic and wind power in the second stage must not exceed the available forecast value after considering uncertainties.
The disparities between the model proposed in this paper and the existing models in the literature are outlined as follows: (1) Modeling differences: While Reference [28] also leverages the cold load for peak shaving, its scope extends to an integrated energy system including cold, heat, electricity, and gas.Conversely, the large logistics park examined in this paper solely includes electrical and cold loads, thereby exhibiting a lower level of energy diversity.Consequently, this model is more adept at illustrating and assessing the park's peak shaving capacity in scenarios where only electrical and cold loads are present.
(2) Indicator differences: Existing models that utilize cold load for peak shaving typically consider only the costs of electrical and thermal power consumed for cooling.This paper, however, refines the usage scenarios of cold loads by incorporating the PMV and SPI, which are amalgamated into the CLEC index for peak shaving control.Additionally, the emergency demand response is introduced to evaluate and compare the resilience to uncertainties among different strategies.
(3) Method differences: Reference [29] employs fuzzy chance-constrained programming for optimized peak shaving, adopting a relatively aggressive strategy that diverges from robust optimization in terms of balancing risks and benefits.Conversely, the model in this paper establishes a two-stage robust optimization framework that accounts for uncertainties in actual wind and photovoltaic operations, while also considering the presence of gas turbines (or steam turbines).The results demonstrate that the utilization of robust optimization enhances the system's capability to manage uncertain factors, rendering it advisable to adopt this more conservative strategy when prioritizing grid security.

Model Solution
Similar to the Benders decomposition method, the C&CG algorithm decomposes the original problem into a master problem and a max-min subproblem [30].It then converts the bi-level optimization subproblem into a single-level optimization form through the Karush-Kuhn-Tucker (KKT) conditions or the strong duality theory (SDT).The optimal solution of the original problem is then obtained by alternately solving the master problem and the subproblem.The key difference between the two methods lies in the fact that during the solution process of the C&CG algorithm, the subproblem continuously introduces relevant variables and constraints into the master problem, resulting in a more compact lower bound for the original objective function, improved efficiency, and effectively reduced iteration times.The model constructed from Equations ( 8)-(44) can be written in the following compact form: where c 0 represents the coefficient column vector corresponding to the objective function (9); Equation (45) denotes the objective function of the two-stage robust optimization problem; x 0 and x S are the relevant control variables in the first-stage and second-stage objective functions, respectively; and Equations ( 46) and (47) represent the constraint conditions for the first and second stages, including Equations ( 10)- (29) and Equations ( 37)-(44), respectively.Among them, A, E, F, G, H, I, and j are the constant column vectors corresponding to the objective functions and constraint conditions.In Equation (48), U represents the uncertainty set, where u 0 w and u 0 p are the available power generation from photovoltaic and wind power in the first stage, respectively, while u w and u p represent the available wind power in the second stage for photovoltaic and wind power, respectively.After decomposition by the C&CG algorithm, the master problem includes the firststage model and the worst-case operating constraints identified by the subproblem: where m represents the current iteration number; u * w,k , u * p,k , and Z * k are the worst-case operating conditions solved by the lower-level problem; and η is the solution of the subproblem to be optimized.The subproblem is a two-level max-min problem.Through the SDT, it can be transformed into a max form.The transformed subproblem form in the m-th iteration is where π is the dual variable of the second-stage constraint conditions.Note that there exists a bilinear term Z T π in Equation (50), as Z is a 0-1 integer variable, resulting in a product form of binary and continuous variables.Therefore, the auxiliary variable and related constraints can be introduced to linearize it using the big M method [31], as shown below: where ω is the introduced continuous auxiliary variable, and M is a sufficiently large constant vector.Thus, the solution process is outlined below: (1) Initialization: Set the upper bound UB of the objective function under the final scheduling plan as +∞, the lower bound LB as −∞, the iteration count m as 1, and the convergence threshold as ε.
(2) Use a set of uncertain variables u w , u p , and Z as the initial worst-case scenario.
( The total methodology presented in this paper is illustrated in Figure 2. The total methodology presented in this paper is illustrated in Figure 2.

Case Studies
Utilizing the MATLAB R2018b platform and the YALMIP toolbox, this paper establishes a load model for large-scale logistics parks and optimizes the results by invoking the Cplex solver.The case study employs the mirrored simulation of the IEEE-6 node

Case Studies
Utilizing the MATLAB R2018b platform and the YALMIP toolbox, this paper establishes a load model for large-scale logistics parks and optimizes the results by invoking the Cplex solver.The case study employs the mirrored simulation of the IEEE-6 node system [32], with the system structure illustrated in Figure 3.The system comprises three gas turbine units and eleven transmission lines, with the parameters detailed in Table 2.
Symmetry 2024, 16, x FOR PEER REVIEW 1 system [32], with the system structure illustrated in Figure 3.The system comprise gas turbine units and eleven transmission lines, with the parameters detailed in Ta    [33] is employed to select typical days from a large number of previous scenarios as the predicted wind-power output within a day.With a 24 h period and a 1 h scheduling unit, four scenarios are established: Scenario 1: Considering only the uncertainty of photovoltaic power, with robustness indices u t,p and u t,w set to 8 and 0, respectively.Scenario 2: Taking into account solely the uncertainty of wind power, where the robustness indices u t,p and u t,w are set at 0 and 8, respectively.Scenario 3: Accounting for the uncertainties of both photovoltaic and wind power, the robustness indices u t,p and u t,w are both set to 8. In this scenario, the uncertainty set implies that the total duration of fluctuations in photovoltaic and wind power does not exceed 8 h within a scheduling period.
Scenario 4: Also addressing the uncertainties of photovoltaic and wind power, but with robustness indices u t,p and u t,w both increased to 16.
Table 3 presents the optimization results under various scenarios, as well as the indices and peak shaving performances of each scenario under their worst operating conditions.C CLEC reflects the comprehensive satisfaction level of the two types of cooling loads used within the park, with a lower C CLEC indicating higher satisfaction.C f ls reflects the peak shaving effect of the park's participation in the distribution network, with a higher C f ls indicating better peak shaving performance.As evident from Table 3, there is minimal difference in various indices between Scenario 2 and Scenario 1, with the emergency demand response participation in Scenario 2 being 18.7% less than that in Scenario 1.In comparison to Scenarios 1 and 2, Scenarios 3 and 4 exhibit significant reductions in both the reserve capacity utilization and EDR participation. Figure 4 illustrates the reserve capacity utilization of power generators under Scenario 3. When compared to Scenario 2, Scenario 3 considers the uncertainties of both photovoltaic and wind power, resulting in an increase in total system cost and total reserve capacity.Specifically, the C CLEC cost increased by 5%, while the C fls compensation decreased by 10%.However, the reduction in EDR load shedding reached 61%.
Both Scenarios 3 and 4 take into account the uncertainties in the actual operation of photovoltaic and wind power, with the difference being that Scenario 4 incorporates greater uncertainty and a more severe worst-case scenario.It can be observed that compared to Scenario 3, Scenario 4 exhibits minimal differences in total operating cost, with slight decreases in the CLEC index and peak shaving performance.However, under the worstcase scenario, the emergency load shedding is reduced to 10 MW.Overall, Scenario 4 demonstrates highest robustness in microgrid operation, with no emergency load shedding even under the most severe conditions.decreased by 10%.However, the reduction in EDR load shedding reached 61%.
Both Scenarios 3 and 4 take into account the uncertainties in the actual operation of photovoltaic and wind power, with the difference being that Scenario 4 incorporates greater uncertainty and a more severe worst-case scenario.It can be observed that compared to Scenario 3, Scenario 4 exhibits minimal differences in total operating cost, with slight decreases in the CLEC index and peak shaving performance.However, under the worst-case scenario, the emergency load shedding is reduced to 10 MW.Overall, Scenario 4 demonstrates highest robustness in microgrid operation, with no emergency load shedding even under the most severe conditions.For this study, the robustness index under Scenario 4 is selected for peak shaving, and the iteration count is set to 13 for system scheduling and control.The resulting convergence criterion is illustrated in Figure 5.In this paper, the fluctuation deviations of wind power and photovoltaic output are set based on historical maximum fluctuation deviations.By observing the final data after the convergence of iterations, a comparison curve of the actual and predicted outputs of wind power and photovoltaic is obtained, as shown in Figure 6.The results of optimization are presented in Figure 7 and Figure 8. Figure 7 depicts the effectiveness of utilizing the CLEC index to engage chiller loads in peak shaving under the final scheme.Figure 8, on the other hand, illustrates the electrical balance histogram corresponding to this scheme.For this study, the robustness index under Scenario 4 is selected for peak shaving, and the iteration count is set to 13 for system scheduling and control.The resulting convergence criterion is illustrated in Figure 5.In this paper, the fluctuation deviations of wind power and photovoltaic output are set based on historical maximum fluctuation deviations.By observing the final data after the convergence of iterations, a comparison curve of the actual and predicted outputs of wind power and photovoltaic is obtained, as shown in Figure 6.The results of optimization are presented in Figures 7 and 8. Figure 7 depicts the effectiveness of utilizing the CLEC index to engage chiller loads in peak shaving under the final scheme.As evident in Figure 7, prior to peak shaving by the logistics park, the net loa riences fluctuations in proximity to peak hours from 10:00 to 15:00, influenced by voltaic and wind turbine outputs alongside the inherent load curve.Conversely, the early morning hours from 00:00 to 06:00, the net load remains at its nadir.Ho once engaged in peak shaving activities, the net load curve undergoes smoothing tively achieving peak shaving.A comparative analysis of various optimization before and after the process is detailed in Table 4    As evident in Figure 7, prior to peak shaving by the logistics park, the riences fluctuations in proximity to peak hours from 10:00 to 15:00, influen voltaic and wind turbine outputs alongside the inherent load curve.Conv the early morning hours from 00:00 to 06:00, the net load remains at its na once engaged in peak shaving activities, the net load curve undergoes smo tively achieving peak shaving.A comparative analysis of various optimiz before and after the process is detailed in Table 4.As evident in 7, prior to peak shaving by the logistics park, the net load experiences fluctuations in proximity to peak hours from 10:00 to 15:00, influenced by photovoltaic and wind turbine outputs alongside the inherent load curve.Conversely, during the early morning hours from 00:00 to 06:00, the net load remains at its nadir.However, once engaged in peak shaving activities, the net load curve undergoes smoothing, effectively achieving peak shaving.A comparative analysis of various optimization metrics before and after the process is detailed in Table 4.In Table 4, the PVD and Fluctuation Variance (FV) are employed as key indicators to gauge the net load variations within the park.Ideally, lower values of these indicators signify reduce peak shaving pressure on the distribution network.Furthermore, the extent of reduction in PVD and FV following optimization using the methodology outlined in this study serves as a measure of the park's effectiveness in participating in peak shaving (note that PMV and SPI values have been normalized for ease of comparison).Notably, there is a substantial decrease of 45.82% and 54.59% in PVD and FV, respectively, indicating a significant alleviation of peak shaving pressure on the distribution network.Additionally, the implementation of CLEC-based optimization has led to a 39.12% improvement in human comfort levels and a 26.36% reduction in sales pressure for frozen goods stored in cold storage facilities.

Non-Domination Verification and Comparative Validation of Indicators
While considering the uncertainties in the actual operation of wind and solar power, this paper also verifies the feasibility of the CLEC index in participating in peak shaving in deterministic optimization.
With the objectives of minimizing the target value of fluctuation degree and minimizing the economic loss caused by the inadequate utilization of the cooling load, modeling is conducted using the YALMIP toolbox in the MATLAB environment.To address the nonlinear issues in the model, the Cplex toolbox is utilized for solving the problem, inputting equality and inequality constraints to obtain the optimal boundary of the Pareto solution set.As the final decision support, the optimal compromise solution is selected from the Pareto solution set using the fuzzy membership degree method [34][35][36].The specific expression of the fuzzy membership function is where µ j i represents the membership function of the j-th solution for the i-th objective function F j i , and F max i and F min i indicate the maximum and minimum values of the i-th objective function among all non-dominated solutions, respectively.The optimal compromise solution set is where M represents the number of non-dominated solutions and n represents the number of objective functions.Figure 9 depicts the impact of the CLEC index I CLEC on the indoor cooling system in deterministic optimization.The colored dashed lines in Figure 9 represent the fluctuation range of indoor temperatures in cooling buildings under corresponding CLEC indices.As shown, the smaller the amplitude of the CLEC index, the stricter the indoor temperature requirements, resulting in a more comfortable environment for humans.The case study in this section for deterministic planning shares the same content parameters as the case study described in Section 5.1, except that it does not consider the uncertainties in the actual operation of wind and solar power.The Pareto frontier of the multi-objective optimization considering both CLEC and net load fluctuation is shown in the figure.As can be seen from Figure 10, the obtained Pareto solution set is concentrated and continuous on the frontier, which facilitates the selection of the most satisfactory solution.It can also be observed from Figure 10 that load fluctuations and the CLEC index are mutually independent, with a higher CLEC index when the fluctuation is small, and vice versa.Continuing with the fuzzy membership function method to solve Figure 10, the optimal solution with the highest satisfaction value is obtained.Further solving under this optimal solution yields the peak shaving results shown in Figure 11.Table 5 compares the various optimization indicators under the deterministic optimization scenario with those under robust optimization Scenario 4 described in Section 5.1.The case study in this section for deterministic planning shares the same content parameters as the case study described in Section 5.1, except that it does not consider the uncertainties in the actual operation of wind and solar power.The Pareto frontier of the multi-objective optimization considering both CLEC and net load fluctuation is shown in the figure .As can be seen from Figure 10, the obtained Pareto solution set is concentrated and continuous on the frontier, which facilitates the selection of the most satisfactory solution.It can also be observed from Figure 10 that load fluctuations and the CLEC index are mutually independent, with a higher CLEC index when the fluctuation is small, and vice versa.Continuing with the fuzzy membership function method to solve Figure 10, the optimal solution with the highest satisfaction value is obtained.Further solving under this optimal solution yields the peak shaving results shown in Figure 11.
Table 5 compares the various optimization indicators under the deterministic optimization scenario with those under robust optimization Scenario 4 described in Section 5.1.
the figure.As can be seen from Figure 10, the obtained Pareto solution set is concentrated and continuous on the frontier, which facilitates the selection of the most satisfactory solution.It can also be observed from Figure 10 that load fluctuations and the CLEC index are mutually independent, with a higher CLEC index when the fluctuation is small, and vice versa.Continuing with the fuzzy membership function method to solve Figure 10, the optimal solution with the highest satisfaction value is obtained.Further solving under this optimal solution yields the peak shaving results shown in Figure 11.Table 5 compares the various optimization indicators under the deterministic optimization scenario with those under robust optimization Scenario 4 described in Section 5.1.By integrating the data in Table 5 and comparing Figure 11 with Figure 7, it is observed that after considering the uncertainties, decreased from 5560.1 to 4291.4,indicating a 22.8% reduction in the park's contribution to the peak shaving.Simultaneously, increased from 1516.8 to 1991.2, representing a 31% decrease in the satisfaction level of the park's cooling load utilization.However, in Scenario 4, which incorporates robust optimization, even under the worst operating conditions, the amount of load shed through EDR was only 10 MW, significantly outperforming the 512.9MW in deterministic the figure.As can be seen from Figure 10, the obtained Pareto solution set is concentrated and continuous on the frontier, which facilitates the selection of the most satisfactory solution.It can also be observed from Figure 10 that load fluctuations and the CLEC index are mutually independent, with a higher CLEC index when the fluctuation is small, and vice versa.Continuing with the fuzzy membership function method to solve Figure 10, the optimal solution with the highest satisfaction value is obtained.Further solving under this optimal solution yields the peak shaving results shown in Figure 11.Table 5   By integrating the data in Table 5 and comparing Figure 11 with Figure 7, it is observed that after considering the uncertainties, decreased from 5560.1 to 4291.4,indicating a 22.8% reduction in the park's contribution to the peak shaving.Simultaneously, increased from 1516.8 to 1991.2, representing a 31% decrease in the satisfaction level of the park's cooling load utilization.However, in Scenario 4, which incorporates robust optimization, even under the worst operating conditions, the amount of load shed through EDR was only 10 MW, significantly outperforming the 512.9MW in deterministic  By integrating the data in Table 5 and comparing Figure 11 with Figure 7, it is observed that after considering the uncertainties, C f ls decreased from 5560.1 to 4291.4,indicating a 22.8% reduction in the park's contribution to the peak shaving.Simultaneously, C CLEC increased from 1516.8 to 1991.2, representing a 31% decrease in the satisfaction level of the park's cooling load utilization.However, in Scenario 4, which incorporates robust optimization, even under the worst operating conditions, the amount of load shed through EDR was only 10 MW, significantly outperforming the 512.9MW in deterministic optimization.This indicates that, although a scheduling scheme considering uncertainties may slightly decrease the satisfaction level of the park's cooling load utilization and diminish its role in peak shaving, the park's operation becomes safer and more robust due to the significant decrease in load shedding through EDR to 10 MW and the 98.1% reduction in emergency load shedding.

Comparison Study: Steam Turbine Replacement for Gas Turbine
Gas turbines are recognized for their swift start-up, flexible adjustment, high ramp rates, and rapid response capabilities, enabling significant energy output within short durations.Conversely, steam turbines excel in terms of efficiency, demonstrating higher energy conversion rates and offering advantages in availability and cost.In this case study, the steam turbine is utilized to replace the gas turbine for supplying electrical power to the park, facilitating a comprehensive comparison between the two technologies.
The steam turbine derives its energy from coal-fired and wood-fired boilers.The generated high-temperature and high-pressure steam is then harnessed in a waste heat boiler to drive the steam turbine for power generation.Considering the relatively slower power response rate of steam turbines in comparison to hourly scheduling, it becomes imperative to account for both their output power constraints and ramp rate constraints.When the demand for thermal load utilization is disregarded, the constraints governing the electrical power output of the steam turbine are outlined as follows: where P 0 st,t represents the output electric power of the steam turbine at time t in the first stage; P min st and P max st are the minimum and maximum output powers of the steam turbine during operation; I t st is the operation status indicator of the steam turbine in time period t; F h st,t denotes the amount of steam consumed by the steam turbine in time period t; η se,e represents the power generation efficiency of the steam turbine; η st,loss is the heat loss rate of the steam turbine; Q coal t and Q wood t are the heat generated by coal-fired and wood-fired boilers in time period t; and η coal and η wood are the working thermal efficiencies of the waste heat boilers for coal-fired and wood-fired boilers, respectively.Ramp constraints are as follows: where R D st and R U st represent the downward and upward ramp rates of the steam turbine, respectively.Due to the high ramp rate and weak reserve capacity of the steam turbine, this section's case study in the second stage considers that when the power supply demand of the park cannot be met, the park can purchase electricity from the distribution network to replace the reserve capacity of the gas turbine.To avoid excessive electricity purchases from the distribution network by the park, which would weaken the peak-shaving effect of electric refrigeration air conditioning in the park, power constraints are imposed on the park's electricity purchases from the distribution network: where P s buy,t represents the amount of electricity purchased by the park from the distribution network in time period t of the second stage; R max g,t is the sum of the upper limits of the reserve capacities of all gas turbines in time period t under the gas turbine case study model in this paper.Considering the uncertainties of photovoltaic and wind power simultaneously, robustness indicators u t,p and u t,w are set to 16 and 16, respectively, for the solution.The scheduling optimization results are shown in Figures 12 and 13.
x FOR PEER REVIEW max , , 0 , where , represents the amount of electricity purchased by the park from the bution network in time period t of the second stage; , is the sum of the uppe of the reserve capacities of all gas turbines in time period t under the gas turbin study model in this paper.Considering the uncertainties of photovoltaic and wind simultaneously, robustness indicators ut,p and ut,w are set to 16 and 16, respectively, solution.The scheduling optimization results are shown in Figures 12 and 13  By comparing Figure 12 with Figure 7 and incorporating the optimization data from various scenarios presented in Table 3, the comparative outcomes between gas turbine and steam turbine outputs under identical uncertainty and robustness conditions are presented in Table 6.As Table 6 indicates, employing a steam turbine for power supply to the park results in a 17.89% reduction in the total operating cost, and this is the benefit of the steam turbine.Figure 12 demonstrates that the steam turbine still achieves satisfactory peak clipping and valley filling effects.However, when compared to scenarios utilizing gas turbines, the improvement in the PVD metric decreases by 14.74%, and the enhancement in the FV metric decreases by 12.32%.Furthermore, the indicator under the steam turbine power supply exhibits only a 4% increase compared to the gas turbine power supply, suggesting a 4% rise in economic losses due to the suboptimal utilization of the park's cooling load.It is evident that using a steam turbine for power supply meets the park's normal operational requirements and satisfies certain peak shaving demands.However, in emergency By comparing Figure 12 with Figure 7 and incorporating the optimization data from various scenarios presented in Table 3, the comparative outcomes between gas turbine and steam turbine outputs under identical uncertainty and robustness conditions are presented in Table 6.As Table 6 indicates, employing a steam turbine for power supply to the park results in a 17.89% reduction in the total operating cost, and this is the benefit of the steam turbine.Figure 12 demonstrates that the steam turbine still achieves satisfactory peak clipping and valley filling effects.However, when compared to scenarios utilizing gas turbines, the improvement in the PVD metric decreases by 14.74%, and the enhancement in the FV metric decreases by 12.32%.Furthermore, the C CLEC indicator under the steam turbine power supply exhibits only a 4% increase compared to the gas turbine power supply, suggesting a 4% rise in economic losses due to the suboptimal utilization of the park's cooling load.It is evident that using a steam turbine for power supply meets the park's normal operational requirements and satisfies certain peak shaving demands.However, in emergency situations, the load shedding amount for demand response in the park reaches 191.9 MW, significantly higher than the 10 MW observed with gas turbines.Thus, while the steam turbine may slightly reduce the park's overall operating cost, it compromises the ability to respond effectively to emergencies.

Conclusions
This paper proposes a novel optimization method for large logistics parks to participate in grid peak shaving in a two-stage robust manner.The conclusions are as follows: (1) This paper introduces a methodology for constructing the CLEC index, which comprehensively integrates the cooling comfort level, measured by PMV, of occupants within the park's buildings, and the sales pressure, indicated by SPI, of frozen goods stored in cold storage facilities.Furthermore, it takes into account the peak shaving demands of the interconnected distribution network.Validation studies confirm that the CLEC index and the park's contribution to peak shaving operate independently, affirming the rationality of this index.
(2) To address the inherent uncertainty in wind and solar power generation within the park, a two-stage robust optimization approach is employed to assess the effectiveness of the proposed method across various operational scenarios.The findings reveal that as the uncertainty surrounding renewable energy output escalates, both the satisfaction level associated with cooling load utilization and the park's peak shaving capabilities decline

Figure 1 .
Figure 1.Configuration of the large logistics park.

Figure 1 .
Figure 1.Configuration of the large logistics park.

)
Solve the master problem (Equation (49)) based on the worst-case scenario.The objective function of the master problem serves as the new lower bound LB = K m , and the control variables x 0 m of the first stage are obtained simultaneously.(4) Substitute the solution x 0 m of the master problem into Equation (50) to solve the subproblem, obtaining its objective function N m and the corresponding uncertain variables u * w,k , u * p,k , and Z * k under the worst-case scenario.Update the upper bound of the objective function as UB = U B , C 0 T x 0 m + N m .(5) If UB − LB ≤ ε, stop the iteration and set the objective function value as UB; otherwise, continue the iteration, increment m by 1, and return to step (3).
subproblem, obtaining its objective function Nm and the corresponding uncertain variables , * , , * , and * under the worst-case scenario.Update the upper bound of the objecthe iteration and set the objective function value as UB; otherwise, continue the iteration, increment m by 1, and return to step (3).

StartFigure 2 .
Figure 2. A flowchart for the proposed method.

Figure 2 .
Figure 2. A flowchart for the proposed method.

Figure 7 .
Figure 7. Peak shaving and valley-filling results in Scenario 4.

Figure 10 .Figure 11 .
Figure 10.The Pareto frontiers between CLEC index and power fluctuation degree.

Figure 10 .
Figure 10.The Pareto frontiers between CLEC index and power fluctuation degree.

Figure 11 .
Figure 11.Peak shaving performance using a deterministic environment.

Figure 11 .
Figure 11.Peak shaving performance using a deterministic environment. .

Figure 13 .
Figure 13.The 24 h electric power balance under steam turbine substitution output.

Table 3 .
Performance under different levels of uncertainties. .

Table 4 .
Performance of the optimization.
compares the various optimization indicators under the deterministic optimization scenario with those under robust optimization Scenario 4 described in Section 5.1.The Pareto frontiers between CLEC index and power fluctuation degree.

Table 5 .
Comparison to the deterministic optimization method.

Table 6 .
Comparison results between gas and steam turbines.

Table 6 .
Comparison results between gas and steam turbines.