Real-time low-carbon scheduling for the wind–thermal–hydro-storage resilient power system using linear stochastic robust optimization

With the large-scale wind power integration, power systems have to address not only the conventional power demand fluctuations but also the wind uncertainty. To improve the economical effectiveness, resilience, and environmental protection of power systems in the source-load uncertainty, a real-time low-carbon scheduling for the wind–thermal–hydro-storage integrated system is proposed. The power imbalance caused by the uncertainty is neutralized by the synergetic linear decision of multiple resources. To address the source-load uncertainty, a stochastic robust optimization is introduced, which establishes the system constraints by robust optimization for the resilience operation, while optimizing the expected operation cost in the empirical uncertainty distribution for economic efficiency. Moreover, a multi-point estimation is applied to formulate the expected operation cost precisely and quickly. By using the dual theory, the proposed real-time power scheduling is derived as a mixed integer bilinear constrained programming. A multi-step sequential convexified solution is developed to solve the complex scheduling problem, which linearizes the bilinear constraints with alternate optimization and relaxes the state variables of energy storages with an “estimation–correction” strategy. Finally, case studies show the superiority of the proposed scheduling method and convexified solution.


Introduction
Large-scale wind power integration is of great significance to energy conservation and emission reduction in power systems (Zhu et al., 2021;He et al., 2022). However, due to the fluctuating nature of wind power, the resilience operation of power systems faces the uncertainty challenge of both wind power and power demand. Real-time scheduling, as a part of multiple time-scale power scheduling, plays an important role in restoring the power balance and improving the system resilience (Surender et al., 2015;Zhou et al., 2018). Therefore, how to address the source-load uncertainty with the co-regulation of wind farms, thermal units, hydro units, energy storages, and other regulation resources in real-time scheduling to ensure the system security should be paid more attention.
At present, methods for the uncertainty problems mainly contain the stochastic optimization and robust optimization. The stochastic optimization Li et al., 2023) formulates the source-load uncertainty with discrete scenarios. The method OPEN ACCESS EDITED BY Zhengmao Li, Nanyang Technological University, Singapore samples multiple scenarios according to the probability distribution of the source-load uncertainty and optimizes the expected operation cost formulated with these sampled scenarios (Li and Xu, 2019). However, it is difficult for these sampled scenarios to cover the whole uncertainty distribution, and hence, the system security cannot be guaranteed in some extreme scenarios. Moreover, the huge sample size and complex calculation also make it difficult for the stochastic optimization to meet the timeliness in real-time scheduling. By contrast, robust optimization methods are widely applied since they can cover the extreme scenarios of the source-load uncertainty. The widely applied robust optimization includes conventional robust optimization (CRO) (Zeng and Zhao, 2013), distributionally robust optimization (DRO) (Ding et al., 2022), and chance-constrained robust optimization (Qi et al., 2020;Li et al., 2022). CRO formulates the source-load uncertainty with the uncertainty set, i.e., the maximum possible space of uncertain scenarios. A two-stage robust unit commitment based on it is established (Bertsimas et al., 2013;Zhou et al., 2019), and the day-ahead scheduling is usually solved with the Benders duality method (Yang and Su, 2021) or column and constraint generation (C&CG) methods . Moreover, the regulation of energy storages is incorporated in the unit commitment, and a nested C&CG is developed to address the integer state variables of energy storages in the intraday regulation (Cobos et al., 2018). CRO only takes the boundary information of the uncertainty into account, while ignores the distribution information, and the methods usually lead to a conservative solution. In view of this, DRO formulates the sourceload uncertainty with the ambiguous set (Wei et al., 2016;Zheng et al., 2021), i.e., the maximum possible space of uncertainty distribution, and optimizes the expected operation cost in the worst-case distributions. However, chance-constrained robust optimization (Qi et al., 2022) assumes that the constraints hold with a certain confidence level. The method improves the flexibility of the scheduling, while the system security is affected.
The aforementioned robust methods can be concluded as the adaptive optimization decision-based robust optimization. In the second stage of the scheduling (i.e., the adjustment stage), the adjustment of regulation resources adapting to the immediate source-load uncertainty is determined through an optimization procedure. However, the system power balance in the real-time scheduling is restored through the automatic generation control (AGC) system (Wu et al., 2017). According to the frequency deviation and tie-line power deviation of the regional power system, the AGC system obtains the power fluctuation from the deviation signals and allocates the adjusted power to regulation resources proportionally. The AGC system has a very high requirement for the decision time, and hence, the aforementioned adaptive optimization decision-based robust methods are inapplicable.
In addition to the adaptive optimization decision-based robust methods, the linear decision-based CRO and DRO (LCRO and LDRO) are also proposed. When the source-load uncertainty occurs, the regulation tasks are allocated to regulation resources, according to the affine factors, i.e., allocation coefficients. Therefore, the linear decision-based robust methods adapt to the real-time scheduling. The co-regulation of generation units and energy storages to wind uncertainty is realized based on LCRO (Jabr 2013). To improve the applicability of LCRO, the allowed interval of wind power is optimized in (Li X et al., 2015;Wang et al., 2019). The allowed upper bound is lower than the forecast upper bound, and hence the strategy makes LCRO applicable for systems with inadequate regulation capacity. In view of the conservativeness of LCRO, which optimizes the operation at the expected scenario, the LDRO-based scheduling model is developed for power systems (Xiong et al., 2017;Alismail et al., 2018). Considering that the worst-case distributions of the source-load uncertainty rarely happen, LDRO may lead to a non-optimal economical effectiveness in common distributions. To enhance the economic efficiency of robust methods in common distributions, a linear decision-based stochastic robust optimization (LSRO) is proposed in (Qu et al., 2020). The method establishes the system constraints using the robust optimization to guarantee the system resilience, while optimizes the operation cost in the empirical distribution of uncertainties for the economical effectiveness. The empirical distribution of uncertainties is the fitting of the existing historical scenarios, and it is closest to the real distribution of uncertainties. Therefore, LSRO shows an economical effectiveness in common distributions of uncertainties.
In view of the robustness and economical effectiveness of LSRO, this paper establishes the real-time scheduling for the wind-thermal-hydrostorage power system based on the method. This paper upgrades the traditional three-point estimation method to multi-point estimation and takes both wind uncertainty and power demand uncertainty into account. Moreover, the extreme scenario-based security constraints (Qu et al., 2020) tend to be conservative, and this study converts the system security constraints, containing the source-load uncertainty to deterministic bilinear constraints, according to dual theory. The proposed real-time power scheduling is derived as a mixed integer bilinear constraint programming, and a convexified solution is developed for the problem. The contributions are concluded as three points. 1) A real-time low-carbon scheduling for the wind-thermal-hydrostorage power system with the source-load uncertainty is established. The power imbalance of the system caused by the uncertainty is neutralized by the co-regulation of multiple resources. The source-load uncertainty is split into two uncertainty subsets to improve the flexibility of the linear decision. However, the carbon trading cost is introduced to realize a low-carbon system. 2) An improved LCRO is introduced to address the source-load uncertainty, which establishes the system constraints using the robust optimization for the resilience operation, while optimizes the expected operation cost in the empirical uncertainty distribution for the economical effectiveness. In addition to, a multi-point estimation is upgraded and applied to calculate the expected operation cost. 3) A multi-step sequential convexified solution is developed to solve the mixed integer bilinear constrained scheduling problem. The sign variables in the objective function are approximated with an initial convex model. The bilinear constraints in the problem are linearized with alternate optimization (AO). Also, the state variables of energy storages are addressed by an "estimation-correction" strategy such that the problem is derived as a convex quadratic programming.
2 Real-time power scheduling model with linear stochastic robust optimization 2.1 Source-load uncertainty set The available wind power of wind farms and the real power demand equal the sum of forecast expectations and forecast errors as follows: where P W,E m,t is the forecast available wind power of the wind farm m, u m,t is the forecast error of the wind power, P D,E n,t is the forecast power demand of the node n and, v n,t is the forecast error of the power demand.
Considering that the regulation resources may be insufficient, a certain amount of wind curtailment is allowed to ensure the system security. Therefore, the system will determine an allowed upper bound of total wind fluctuations (Ũ t ).Ũ t is smaller than the forecast upper bound ( U t ). The superfluous wind power exceeding the allowed threshold will be curtailed as follows: where PW m,t is the baseline power output of the wind farm m,P W m,t is the real wind power output,Û t is the total available wind fluctuations, and U t is the lower bound of wind fluctuations.
When considering the wind curtailment, the real wind fluctuation is smaller than the available wind fluctuation, and thus the system security constraints are only associated with the allowed wind fluctuation (i.e., real wind fluctuation). The uncertainty set of the allowed wind fluctuation is established as where U t is the total allowed wind fluctuation andŨ t is the allowed upper bound of the total wind fluctuation.
As for the power demand, load shedding is not allowed, and thus the uncertainty set for the forecast error of the power demand is formulated as where V t is the total power demand fluctuation. Therefore, the source-load uncertainty set is a union set of the wind uncertainty set and power demand uncertainty set, which is given as: Ξ t = ΞU t∪ΞV t.

Overall architecture of real-time scheduling
Real-time power scheduling optimizes the baseline power output of wind farms, thermal units, hydro units, and baseline charging power of energy storages, as well as the affine factors of non-wind regulation resources for the source-load uncertainty. The real-time power scheduling is usually performed in a rolling optimization way to determine the scheduling plan for the following 1 h period with 5 min resolution. When the wind power and power demand deviate from the baseline value, nonwind regulation resources including thermal units, hydro units, and energy storages co-regulate to restore the system power balance with the affine factors.
where P T i,t , P H j,t and P S k,t are the baseline power output of the thermal unit i, hydro unit j, and baseline charging power of the energy storage k;P T i,t , P H j,t , andP S k,t are their real power outputs or charging power under the source-load uncertainty and; α − t /α + t are downward/upward affine factors of non-wind regulation resources for the source-load uncertainty.
To sum up, the strategy for balancing the source-load uncertainty is presented as follows: Power demand: Load shedding is not allowed, and power demands always trace the actual demands.
Wind farms: If the total available wind fluctuation is below the allowed upper boundŨ t , each wind farm operate at the available value, otherwise, the superfluous wind power beyondŨ t is curtailed, as shown in Eq. 2.
Non-wind resources: Thermal units, hydro units, and energy storages co-regulate as Eqs 5-7 to balance the source-load uncertainty.

System security constraints
1) Power balance: 2) Constraints of wind farms: 3) Constraints of thermal units: where R − i,t /R + i,t are the downward/upward reserves produced by the thermal unit i for the source-load uncertainty, Eq. 12 is the lower and upper power output limits of thermal units, Eqs 13-14 are ramping limits of thermal units, and r − i /r + j are downward/upward ramping rates of thermal units.

4) Constraints of hydro units:
Frontiers in Energy Research frontiersin.org 03 where R − j,t /R + j,t are the downward/upward reserves produced by the hydro unit j for the source-load uncertainty, Eq. 16 is the lower and upper power output limits of hydro units, Eqs 17-18 are ramping limits of hydro units, r − i /r + j are downward/upward ramping rates of hydro units, Eqs 19-21 are water consumption limits of hydro units (Liu et al., 2019), g is the hydro-to-power conversion coefficient, η j is the conversion efficiency of hydro units, H j,t is the height of the water head, andQ − j, t/Q j,t are the extreme water consumptions of hydro units in the adjustment process.

5) Constraints of energy storages:
where R − k,t /R + k,t are the downward/upward reserves produced by energy storage k for the source-load uncertainty,P − S k,t andP S k,t are the extreme charging power of the energy storage k in the adjustment process,C S k,t andC S k,t are the extreme stored power in the energy storage, C S k,0 and C S,R k are the initial stored power and required stored power after the scheduling, respectively, η + k and η − k are the charging and discharging efficiency of the energy storage,π S k,t ,π S k,t , and π S k,t are the binary state variables in the extreme scenario and baseline scenario. Ifπ S k,t ,π S k,t , and π S k,t equal to 1, the energy storage charges power, otherwise, the energy storage discharges power.

6) Transmission limits of lines:
The maximum/minimum operation of wind farms, thermal units, hydro units, and energy storages occur at the extreme total source-load fluctuation, and hence their constraints are established using the extreme scenario method. As for transmission lines, their transmitted power is related to the power output of each wind farm and power demand of each node, which cannot be described with the total source-load fluctuation. Therefore, we formulate the transmission limits of lines as that concerning the source-load uncertainty as follows: where L l,t is the baseline transmitted power of line l, h l,i , h l,j , h l,k , h l,m , and h l,n are transfer distribution factors from the thermal unit i, hydro unit j, energy storage k, wind farm m, and power demand n to the line l, and T l / T l are transmitted power limits of line l. Ξ t,1 and Ξ t,2 are two subsets of the sourceloaduncertaintyset.InΞ t,1 ,U t -V t ≤0,thesystemadoptsthedownwardaffine factors to balance the uncertainty, while in Ξ t,2 , U t -V t ≥0, the system adopts the upward affine factors to balance the uncertainty.
The transmission limits of lines can be described as the following compact form: where A, B, C are constant vectors, d is a constant scalar, x is the decision variable vector, α j is the affine factor matrix, concerning the jth uncertainty subset, and π, λ, τ, and ϑ are dual variables, corresponding to constraints in the uncertainty subset Ξ t,j . By adopting the dual theory , the constraint (29) with the source-load uncertainty is converted to a deterministic constraint as follows: where e is a vector whose elements are all equal 1.

Optimized objective
The scheduling optimizes the baseline operation cost (i.e., the first stage) and the incremental cost expectation in the adjustment state (i.e., the second stage) as follows: where F base is the baseline operation cost, P is the empirical distribution of the source-load uncertainty, ΔF is the incremental adjustment cost under the source-load uncertainty (û, v). This paper introduces carbon trading (Wang et al., 2020a) to improve the environmental benefits of the resilient power system. The operation cost in the baseline state includes fuel costs of thermal units, system carbon trading costs, and penalty costs on the wind curtailment.
where C i () is the quadratic fuel cost function of the thermal unit i, f c is the carbon trading costs of the system, and W() is the quadratic penalty cost function of the wind curtailment.

Frontiers in Energy Research frontiersin.org
This study adopts the baseline method (Wang et al., 2020b;Guo et al., 2020) to calculate the carbon quota of the system, and the carbon trading cost of the system is formulated as follows: where p c is the carbon trading price, E R and E Q are the baseline carbon emission and carbon quota of the system, respectively, β i,t is the carbon emission coefficient of the thermal unit i, anf η c is the carbon quota for per unit of power generation. The incremental cost expectation in the adjustment state includes the adjustment cost of thermal units and energy storages, incremental penalty costs on the wind curtailment, and incremental carbon trading costs. The traditional random sampling methods to formulate the output expectation usually requires a huge sampling scale to achieve an accurate estimation. In this paper, the multi-point estimation (Li Z et al., 2015) can uniformly cover the source-load uncertainty with a small scale of estimation points.
The multi-point estimation actually first determines the estimate vectors in the standard Gaussian distribution, then converts these estimate vectors to that of the source-load uncertainty, and finally calculates the output expectation (Li et al., 2008). The optimal estimate point number will be determined through sensitivity analysis in numerical simulations. The source-load uncertainty is assumed to be irrelevant in the time dimension, and thus only two uncertain variables are considered in the formulation of the incremental cost expectation, i.e. the available wind fluctuation and forecast error of the power demand. The incremental cost expectation in the adjustment state is formulated as U t,r,s , V t,r,s F −1 t Ψ z r,s , where r∈N is the index of uncertain variables, s∈M is the index of estimate points, ω s is the weight of the sth estimate point, (r, s) corresponds to the estimate vector concerning the rth uncertain variable and the sth estimate point, and z r,s is the estimate vector in the standard Gaussian distribution. In the estimate vector, the uncertain variable r equals the sth estimate point and other uncertain variables equal 0. Ψ is the cumulative distribution function of the standard Gaussian distribution. F t −1 is the inverse of the cumulative distribution function of the source-load uncertainty.
Assume that the sign variable concerning the adopted affine factors in each estimate vector is σ t,r,s . If σ t,r,s = 1, the source-load power fluctuation is larger than 0, the system adopts the upward affine factors for the adjustment; else if σ t,r,s = 0, the system adopts the downward affine factors for the adjustment. In each estimate vector, the incremental cost produced by thermal units and energy storages are calculated as U t,r,s ≤Û t,r,s , U t,r,s ≤Ũ t , D t,r,s U t,r,s − V t,r,s , 2σ t,r,s − 1 D t,r,s ≥ 0, Δf T i,t,r,s ≥ − 1 − σ t,r,s ρ i,t α − i,t D t,r,s − σ t,r,s ρ i,t α + i,t D t,r,s Δf S k,t,r,s ≥ 1 − σ t,r,s ρ − k,t α − k,t D t,n,s + σ t,r,s ρ + k,t α + k,t D t,r,s , where U t,r,s is the allowed wind fluctuation in the estimate vector (r, s), ρ i,t, and p i,t are the incremental cost coefficient and regulation cost coefficient of the thermal unit i, respectively, and ρ− k, t ,and ρ+ k,t are the downward and upward incremental cost coefficient of the energy storage k. Since hydro units do not produce fuel costs and carbon emissions in the operation, their incremental cost coefficients equal 0. As for energy storages, since they also play a role of peak-load shifting and standby power generation during intraday operation, the stored power in storages after the scheduling shall not deviate too far from the required value. Therefore, the incremental cost coefficients are set as follows to reduce the adjustment of energy storages: ρ− k,t = -1.1 × max i (ρ i,t ), ρ+ k,t = -0.9 × max i (ρ i,t ).
To sum up, the incremental cost in each estimate vector is calculated as where ρ w is the penalty coefficient on the wind curtailment.

Solution methodology
The proposed scheduling model is a mixed integer bilinear constraint programming. The integer variables include the sign variables σ t,r,s in each estimate vector concerning the incremental cost expectation and the charging/discharging state variables of energy storages. The bilinear constraints include Eq. 11 (15) (23) (30). Due to the timeliness required in the real-time scheduling, the integer variables and bilinear constraints shall be addressed, and such that the scheduling model can be converted to a convex quadratic programming.

Approximation for the sign variables
To determine the sign variables σ t,r,s in each estimate vector concerning the incremental cost expectation, an approximated value for the allowed upper bound of wind fluctuation (Ũ 0 t ) should be first obtained. This paper adopts the extreme scenario methodbased initial model to obtainŨ 0 t .

1)
Objective of the initial model: The objective should be convex and can be close to the original objective (31). The objective of the initial model is 2) Constraints of the initial model: Reformulate the original bilinear constraints concerning the reserve constraints of non-wind regulation resources and transmission limits of lines with the extreme scenario method, while other security constraints remain unchanged. Therefore, the original bilinear constraints are converted to slightly conservative but convex linear constraints.
σ l,m,t ≤ 0 ≤ σ l,m,t , σ l,n,t ≤ 0 ≤ σ l,n,t σ l,t ≤ 0 ≤ δ l,t σ l,m,t ≤ h l,m u m,t ≤ σ l,m,t , σ l,m,t ≤ h l,m u m,t ≤ σ l,m,t σ l,n,t ≤ h l,n v n,t ≤ σ l,n,t , σ l,n,t ≤ h l,n v n,t ≤ σ l,n,t δ l,t ≤ i h l,i R Therefore, the initial model does not contain the sign variables σ t,r,s and bilinear constraints. Compare the available wind power in each estimate vector (Û t,r,s ) with the approximated value for the allowed upper bound of wind fluctuation (Ũ 0 t ) and obtain the sign variables σ t,r,s in each estimate vector.

Linearization for bilinear constraints
This paper adopts the AO  to address the bilinear constraints in the real-time scheduling, which fixes partial decision variables to optimize the rest of the decision variables in the bilinear constraints. The real-time scheduling with bilinear constraints can be described as the following compact form: where A i is a constant matrix, g i is a convex function concerning decision variables, y corresponds to the allowed upper bound of wind fluctuation (Ũ t ), and x corresponds to the remaining decision variables. Detailed steps for AO to solve the real-time scheduling model are presented in Algorithm 1.

1) Initialization: Solve the initial model as shown in
Eqs 42-45 and obtain the initial value of the allowed upper bound of wind fluctuation (Ũ 0 t ), i.e., y 0 . Set the iteration step n = 0.
2) Solve the sub-problem A: Fix y = y n and solve the below problem. The obtained x is marked as x n .
min F x, y n s.t. x ⊤ A i y n + g i x, y n ≤ 0, ∀i .
3) Solve the sub-problem B: Fix x = x n and solve the below problem. The obtained y is marked as y n+1 , and the obtained objective is marked as F n+1 .
min F x n , y where the convergence gap ε is set to 10 −6 .
Algorithm 1 AO to solve the bilinear problem.
Note that AO is a sub-optimal method, whose optimization performance is associated with the initialization. To improve its optimization performance, the initial value of the allowed upper bound of wind fluctuation (Ũ 0 t ) can be evenly initialized with multiple values, and the most optimal solution is finally selected. The initial model as shown in Eqs 42-45 will determineŨ 0 t , and then updateŨ 0 where β is the confidence level corresponding to the upper bound of the total available wind fluctuation, i.e. the cumulative distribution function of U t, β . Evenly set the confidence level and then determine multiple initial values of the allowed upper bound of wind fluctuation.

Approximation for state variables of energy storages
This paper adopts an "estimation-correction" strategy (Shi and Guo, 2022) to approximate the state variables of energy storages. The estimation stage relaxes the state variables of energy storages, and the correction stage directly determines the state variables according to the results obtained in the estimation stage.
Therefore, the scheduling problems in these two stages do not contain the state variables of energy storages. The relaxed model of energy storages in the estimation stage is formulated as where P + k,t , P − k,t ,P − + k,t ,P − − k,t ,P + k,t ,P − k,t are the charging power and discharging power of energy storages in the baseline scenario and two extreme scenarios.
The state variables of energy storage are directly fixed in the correction stage. For the sake of convenience, the following takes the baseline scenario as an example, while the other two extreme scenarios can be disposed in a similar way.
where P + k,t.0 , P − k,t,0 are the charging power and discharging power of energy storages obtained in the estimation stage.
In fact, there is generally little difference between the optimal state variables of energy storages in the estimation stage and the correction stage. On the one hand, the charging and discharging efficiency of energy storages are generally higher. In the real-time scheduling, only the energy storages with the rapid adjustment ability could be involved in the AGC system, such as the lithium storages, and the efficiency of lithium storages can reach about 90% at present (Shi and Guo, 2022). On the other hand, the current capacity of energy storages is generally small compared with the total power demand, and thus the relaxed

Numerical simulation
To illustrate the effectiveness of the proposed real-time lowcarbon scheduling and solving method, numerical simulations are carried out on the IEEE 39-node power system (Zimmerman et al., 2011). The system has seven thermal units, one hydro unit, and two energy storages with 50MW capacity. Moreover, four wind farms are deployed. The scheduling cycle is from 6:00 pm to 7:00 pm with 5 min time resolution. The forecast expectation and forecast interval of wind power and power demand are presented in Figure 2 and Figure 3. Finally, 5,000 scenarios of the source-load uncertainty are sampled by Monte Carlo simulation (MCS) to verify the scheduling performance. The optimizations are solved with the CPLEX solver on a computer with Intel(R) Xeon(R) Silver 4216 CPU, 2.1GHz clock speed, and 16GB RAM.

Effects of the allowed upper bound of wind fluctuation
In the proposed scheduling model, the allowed upper bound of wind fluctuation actually has two functions: 1) it reasonably determines the wind power accommodation and hence enhances the system flexibility, 2) when the system regulation capacity is insufficient, discard the excess wind power to ensure the system

FIGURE 1
Overall flowchart for the convexified solution.

FIGURE 2
Forecast expectation and forecast interval of wind power.

FIGURE 3
Forecast expectation and forecast interval of the power demand.

FIGURE 4
Allowed upper bound (UB) of wind fluctuation in the system with sufficient regulation capacity.

FIGURE 5
Allowed upper bound (UB) of wind fluctuation in the system with insufficient regulation capacity.

Frontiers in Energy Research
frontiersin.org security. Figure 4 and Figure 5 present the allowed upper bounds (UB in the figures) of wind fluctuation in the system with sufficient and insufficient regulation capacity, where the maximum capacity of wind power accommodation is the obtained allowed upper bound of wind fluctuation when adding the wind penalty objective, i.e., the second line in Eq. 42, to the original objective. Moreover, in the system with insufficient regulation capacity, change the ramp rates of thermal units and hydro units as well as the maximum charging and discharging power of energy storages to 75% of the original. As shown in Figure 4, although the system with sufficient regulation capacity can completely accommodate wind power, partial wind power is still discarded. This is mainly because the wind fluctuation exceeding the allowed upper bound rarely appears, and hence these wind powers are discarded to improve the system flexibility. In addition to, as shown in Figure 5, the system with insufficient regulation capacity cannot completely accommodate wind power, and hence the allowed upper bound is below the forecast upper bound to ensure the system security. To further validate the effects of the allowed upper bound in improving the system flexibility, Table 1 compares the optimization results and MCS results with two wind power accommodation strategies. The maximum accommodation strategy is to absorb wind power as far as possible, while the flexible accommodation strategy is that in the proposed scheduling model. As shown in Table 1, although the allowed upper bound of wind fluctuation in the flexible strategy is obviously smaller than the maximum strategy, the discarded wind power is very slight. This is because the wind fluctuation exceeding the allowed upper bound rarely appears. Through the flexible accommodation strategy, the system operation cost is reduced, which validates the effects of allowed upper bound in improving the system flexibility. Figure 6 presents the baseline outputs of thermal units and hydro units, and the baseline charging power of energy storages. Hydro units always operate at the maximum output since they do not produce fuel costs and carbon trading costs. Thermal units produce fuel costs and carbon trading costs, and hence, they arrange the outputs according to the cost-optimal principle. In the scheduling cycle of the testing case, the load peak is in the front period, while the load trough is in the latter period. Therefore, the energy storage is discharged according to the required state in the front period to realize the peak load shifting and reduce the operation cost. Figure 7 and Figure 8 present the affine factors of thermal units, hydro units, and energy storages for the downward and upward

FIGURE 6
Baseline states of thermal units, hydro units, and energy storages.

FIGURE 7
Affine factors for the downward source-load uncertainty.
Frontiers in Energy Research frontiersin.org 08 source-load uncertainty, respectively. Affine factors of regulation resources are closely related to their incremental cost coefficients. When the source-load uncertainty fluctuates downward, the system operation cost is increased. The incremental cost coefficients of regulation resources satisfy ρ − k,t > ρ − i,t > ρ − j,t , and hence the regulation priority for the downward source-load uncertainty satisfies: hydro units > thermal units > energy storages. Since hydro units operate at the maximum output and have no room for the upward adjustment, they do not undertake the regulation task for the downward sourceload uncertainty, as shown in Figure 7. During periods t = 1-6, the source-load uncertainty is relatively small, thermal units have sufficient regulation capacity, and hence the regulation tasks for the downward source-load uncertainty are fully undertaken by the thermal units. While during periods t = 7-12, the regulation capacity of the thermal units is insufficient, and hence energy storages participate in the regulation for the downward source-load uncertainty. When the source-load uncertainty fluctuates upward, the system operation cost is decreased. The incremental cost coefficients of regulation resources satisfy ρ + i,t > ρ + k,t > ρ + j,t , and hence the regulation priority for the upward source-load uncertainty satisfies thermal units > energy storages > hydro units. As shown in Figure 8, during periods t = 1-7, thermal units completely undertake the regulation task for the upward source-load uncertainty, while energy storages and hydro units participate in the regulation for the upward source-load uncertainty during next periods t = 8-12. As indicated in the simulation results, the LSRO reasonably arranges the affine factors, according to the incremental cost coefficients of regulation resources, and thus it can reduce the system incremental cost in the adjustment state.
Moreover, the existing AGC system usually utilizes units for the regulation. Due to the rapid adjustment speed, energy storages can participate in the AGC system to regulate the source-load uncertainty. To validate the effects of energy storages in the regulation, Table 2 compares the optimization results and MCS results without and with the regulation of energy storages.

Comparison of different robust methods
In order to further verify the economic efficiency of the LSRO, we compare it with LCRO and LDRO. The robust methods all optimize the allowed upper bound of wind fluctuation. The difference of these three robust methods is that, LSRO optimizes the expected cost under the empirical distribution of the source-load uncertainty, LCRO minimizes the baseline cost and maximizes the allowed upper bound of wind fluctuation, and LDRO minimizes the expected cost under the worst-case distributions of the sourceload uncertainty. As shown in Table 3, LCRO minimizing the baseline cost ignores the impact of the adjustment process on the operation cost, and hence its adjustment is uneconomical. LDRO minimizes the expected cost under the worst-case distributions, however, the worst-case distributions rarely appear, and hence LDRO actually sacrifices the economic efficiency in general distributions for robustness in worst cases. The empirical distribution is the fitting of the historical scenarios, which is closest to the real distribution of the source-load uncertainty, and hence LSRO shows the best economic efficiency.

Impacts of estimate point numbers
This paper upgrades the traditional three-point estimation to multi-point estimation. The estimate point numbers are set at 3,

FIGURE 8
Affine factors for the upward source-load uncertainty.  Figure 9 shows the allowed upper bound of wind fluctuation obtained by multi-point estimation with different point numbers. The maximum cumulative distribution function covered by estimate points in three-point estimation, five-point estimation, and sevenpoint estimation are 0.9584, 0.9979, and 0.9999, respectively. The distribution of uncertain variables covered by three-point estimation is limited, while the uncertainty distribution covered by five-point estimation and seven-point estimation is large enough. Therefore, the allowed upper bound of wind fluctuation obtained by three-point estimation is much smaller compared to that obtained by five-point estimation and seven-point estimation. In addition, Table 4 compares the results of optimization and MCS with different point numbers, and it can be seen that there is little difference between the results. Considering that the seven-point estimation can cover the maximum distribution of uncertain variables, and the approximate degree of the expected value evaluation generally enhances with the increase in the point number (Li X et al., 2015), the estimate point number is set at 7 in this paper to ensure the robustness of the multipoint estimation method.

Validation of solving methods
As indicated in subsection 3.2, AO is a sub-optimal method, whose optimization is related to initialization. This paper evenly sets the confidence level corresponding to initialization and updates the initial value of the allowed upper bound of wind fluctuation as Eq. 51. The results of optimization and MCS associated with the confidence level are presented in Table 5. Generally, the obtained allowed upper bound of wind fluctuation decreases with decrease in the confidence level. When the confidence level is high, the obtained allowed upper bound is large, and hence the wind curtailment is small. However, the system flexibility is affected due to excessive wind power accommodation, and thus the system operation cost increases. As the confidence level decreases, the obtained allowed upper bound decreases, the system flexibility improves, and the system operation cost is reduced. However, when the confidence level reaches 97.5%, the obtained allowed upper bound is too small, and hence the wind power accommodation is inadequate and the system operation cost gradually increases. Considering that the wind power prediction  accuracy has a certain degree of error, the confidence level can be set at 99.7% to absorb wind power as far as possible without compromising the system's flexibility.
In addition, the integer variables of charging/discharging state of energy storages are addressed with the estimation-correction strategy. Table 6 compares the optimization results with different strategies for the integer variables. The precise strategy takes the state variable of energy storage as variables in the optimization, and hence the scheduling problem is formulated as an iterative mixed integer quadratic constraint programming, which takes a long time to solve. The estimation stage relaxes the state variables, and the correction stage directly fixes the state variables as parameters, according to the results obtained in the estimation stage. Therefore, the scheduling problem is formulated as an iterative quadratic constraint programming with the estimation-correction strategy, and the computation time is short. On the other hand, the theoretically optimal objectives of these three strategies satisfy the following relations: F estimation < F precise < F correction . Since the deviation between the optimization results of the estimation stage and the correction stage is small, it can be deduced that the estimation-correction strategy can find a good enough solution.

Impacts of carbon trading price
Carbon trading cost is an important part of the proposed realtime low-carbon scheduling. To explore the influence of the carbon trading price on the scheduling, this subsection changes the carbon trading price from 0 to 60$/t and obtains the variation trend of the operation cost and carbon emission of the system relative to the carbon trading price. As indicated in Figure 10, carbon trading has little impact on the system operation when the carbon trading price is low, and thus the operation cost and carbon emissions change slightly. With the increase in the carbon trading price, the proportion of the carbon trading cost in the operation cost increases, the system gradually utilizes the low-carbon resources for power generation, and the carbon emission show a significant downward trend. When the carbon trading price rises to 26$/t, the system can make profits by the carbon trading mechanism, and the operation cost gradually decreases. When the carbon trading price rises to 40$/t, the carbon reduction space of the system becomes smaller, and the downward trend of carbon emissions slows down.

Simulation results
To further verify the effectiveness of LSRO, supplementary simulations are performed in a large-scale IEEE 118-node power system. The system has seventeen thermal units, two hydro units, and four energy storages with 40MW capacity. Furthermore, six wind farms are deployed. As shown in Table 7, the results in the IEEE 118-node system are similar to those in the IEEE 39-node system. LCRO focuses on the baseline cost while ignores the influence of the adjustment on the operation cost, and the adjustment of the method is uneconomical. LDRO minimizes the expected cost under the worst-case distribution, in view of the extremely low probability of the worst-case distribution, and the method is not effective enough in common distributions. Compared to LCRO and LDRO, LSRO minimizes the expected cost under the empirical distribution of the source-load uncertainty. In view of that, the empirical distribution is closest to the real distribution of the source-load uncertainty, and LSRO shows the best economic efficiency.

Further discussion on the robust method
In fact, the performance of LSRO is related to the accuracy of the empirical distribution of the source-load uncertainty. If the

FIGURE 10
Variation trend in operation cost and carbon emissions relative to the carbon trading price. scale of historical scenarios is not huge enough or the adopted fitting method to determine the empirical distribution has a poor performance, the empirical distribution may deviate far from the real distribution. Therefore, the performance of the robust method in other distributions is affected in these cases. To improve the robustness of LSRO, our future work will combine the empirical distribution and worst-case distribution together in the robust method.

Conclusion
This paper proposes a real-time low-carbon scheduling for the wind-thermal-hydro-storage resilient power system based on the linear decision rule-based stochastic robust optimization. A multistep alternate optimization-based convexified solution is formulated to solve the complex scheduling problem. The allowed upper bound of wind fluctuation provides the system with more flexibility. Monte Carlo simulation shows that very slight redundant wind power is sacrificed to significantly reduce the system operation cost. By the stochastic robust optimization, the regulation task to regulation sources is reasonably arranged according to the incremental cost coefficients of resources, and thus an economically effective adjustment is realized. The comparison of the adopted stochastic robust optimization with conventional robust optimization and distributionally robust optimization in two case studies validates the economic efficiency of the method. Numerical simulations indicate that the performance of alternate optimization is related to initialization, and the confidence level concerning the initialization can be reasonably set to absorb wind power as far as possible without compromising the system flexibility.
The empirical distribution of the source-load uncertainty may not be an accurate estimation of the real distribution, when the number of historical scenarios is inadequate or the fitting method is inaccurate. Our future work will combine the empirical distribution and worstcase distribution to improve the applicability of the linear decisionbased robust method.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.