The dispatch optimization problem at the regional grid level studied in this paper involves two types of units: conventional generation units and pumped-storage power stations. The main constraints can be categorized as follows:

Controllable Conventional Unit Constraints:

Ramping Constraint Equation 1: − R U g ≤ P g , p , t + 1 − P g , p , t ≤ R U g (1) In the formula, t denotes the time period index, and p represents the province index. P g , p , t + 1 , P g , p , t are the generation power at time t+1 and t respectively. R U g represents the ramp rate limit of the thermal power unit, indicating the maximum increase or decrease in output per unit time.

Unit Output Boundary Constraint Equation 2: 0 ≤ P g , p , t ≤ P g max (2) In the formula, P g max denotes the maximum output of the unit, and the constraint is defined with a zero lower bound, which applies to the generation stage of all units.

Pumped-Storage Unit Constraints Equations 3–7:

Pumped-Storage Unit Equality Constraints, including water level conversion power constraint and power balance constraint: η c h ⁡ arg ⁡ e Δ E = P c h ⁡ arg ⁡ e t = η c h ⁡ arg ⁡ e m g Δ h (3) Δ E = η d i s c h ⁡ arg ⁡ e P d i s c h ⁡ arg ⁡ e t = η d i s c h ⁡ arg ⁡ e × m g + 1 2 ρ V 2 S × Δ h (4) E t + 1 = E t + P t c h ⁡ arg ⁡ e η c h ⁡ arg ⁡ e + P t d i s c h ⁡ arg ⁡ e η d i s c h ⁡ arg ⁡ e Δ t (5) The boundary constraints of pumped -storage units include energy boundary constraints as well as pumping and generating power constraints: E min ≤ E t ≤ E max (6) 0 ≤ P t c h ⁡ arg ⁡ e ≤ x t c h ⁡ arg ⁡ e P max c h ⁡ arg ⁡ e − x t d i s c h ⁡ arg ⁡ e P max d i s c h ⁡ arg ⁡ e ≤ P t d i s c h ⁡ arg ⁡ e ≤ 0 x t c h ⁡ arg ⁡ e + x t d i s c h ⁡ arg ⁡ e ≤ 1 (7) In the formula, E t represents the equivalent stored energy of the pumped-storage unit at time period t; Δ E denotes the equivalent energy change corresponding to the variation in water level; m , g , h , ρ , V , S respectively refer to the mass of water in the reservoir, gravitational acceleration, water level height, water density, volume, and bottom area; P t c h ⁡ arg ⁡ e , P t d i s c h ⁡ arg ⁡ e are the pumping and generating power of the pumped-storage unit at time period t; η c h ⁡ arg ⁡ e , η d i s c h ⁡ arg ⁡ e are the efficiency of pumping and power generation for the pumped-storage station; E min , E max denote the minimum and maximum stored energy of the pumped-storage station; P max c h ⁡ arg ⁡ e , P max d i s c h ⁡ arg ⁡ e are the maximum output for both pumping and generating; and x t c h ⁡ arg ⁡ e , x t d i s c h ⁡ arg ⁡ e are the pumping status of the pumped-storage station at time period t, and they are Boolean variables.

Unit Output Proportion Constraint Equation 8:

The output of all types of units needs to be allocated to each provincial grid according to the agreed proportions, and during the pumping period, the corresponding proportion of grid power is also consumed accordingly. P g , p , t = ∑ p p r o v i n c e s P g , p , t × R g , p (8) In the formula, P g , p , t represents the output power of the gth unit to province p at time period t; R g , p denotes the power transmission ratio of the gth unit to the province p.

The provincial two-stage robust optimization model presented in this paper is based on the unit dispatch data obtained from the regional grid-level dispatch strategy after load peak shaving. It addresses the source-load uncertainty problem under the integration of renewable energy.

Controllable Generation Units:

Controllable generation units include adjustable gas turbines, diesel generators, and others. After linearization, their cost function can be expressed as a linear function of their generation power Equation 9: C g = ∑ t T a P g , t + b Δ t (9) In the formula, P g , t represents the generation power of the unit at time period t; a, b are the cost coefficients; and Δ t denotes the duration of each time period.

Output Power Constraint Equation 10: P g min ≤ P g , t ≤ P g max (10) In the formula, P g min , P g max denote the minimum and maximum generation power of the unit, respectively.

Energy Storage Components.

The cost of energy storage components is the cost associated with the charging and discharging processes, which can be expressed as Equation 11: C s = ∑ t T K s P t d i s η + P t c h η Δ t (11)

In the formula, P t d i s , P t c h represent the charging and discharging power of the energy storage at time period t; η is the charging and discharging efficiency; and K s denotes the unit charging and discharging cost.

The constraints that the energy storage components must satisfy are as follows Equation 12: 0 ≤ P t d i s ≤ u t s P max 0 ≤ P t c h ≤ 1 − u t s P max η ∑ t T P t c h Δ t − ∑ t T P t d i s Δ t η = 0 E min ≤ E 0 + η ∑ i t P i c h Δ t − ∑ i t P i d i s Δ t η ≤ E max (12) In the formula, P max represents the maximum allowable output of the energy storage; u t s is a Boolean variable indicating the charging or discharging state of the energy storage; and E min , E max denote the lower and upper bounds of the energy storage capacity.

Demand Response Load:

The flexible scheduling process of demand response load is considered in Zhang et al. (2025a), Zhang et al. (2025b), Wang et al. (2025); Li et al. (2025); Shao et al. (2024). Under the condition that the electricity usage characteristics meet the requirements for providing demand response services, the grid can adjust users’ electricity consumption plans while providing appropriate compensation to the users. Thus, the adjustment cost of the demand response load is expressed as Equation 13: C D R = ∑ t T K D R P t D R − P t D R 0 Δ t (13)

The constraints that need to be satisfied are as follows Equation 14: ∑ t T P t D R Δ t = D D R D t min ≤ P t D R Δ t ≤ D t max (14) In the formula, P t D R represents the actual dispatch power of the demand response load at time period t; D D R is the total electricity demand across all time periods; D t min , D t max denote the minimum and maximum demand response load at time period t, respectively.

Since Equation 13 is a nonlinear function, auxiliary variables are introduced to linearize it. The final cost function and constraints are expressed as Equations 15, 16: C D R = ∑ t T K D R P t D R 1 + P t D R 2 Δ t (15) P t D R − P t D R 0 + P t D R 1 − P t D R 2 = 0 P t D R 1 ≥ 0 P t D R 2 ≥ 0 (16) In the formula, P t D R 1 , P t D R 2 are the linearized auxiliary variables.

Power Trading Interaction:

When the generation output cannot meet the actual load demand, power needs to be purchased from other sectors, and when there is excess power, it can be sold to others.

The power trading process needs to satisfy the basic power balance Equation 17: P t b u y − P t s e l l = P g , t + P t c h − P t d i s + P t D R + P t L o a d − P t P V (17) In the formula, P t b u y , P t s e l l represent the power purchased or sold at time period t; P t L o a d , P t p v denote the load data and photovoltaic generation data at time period t, respectively.

The cost function and constraints of power trading are expressed as follows Equations 18, 19: C m = ∑ t T p t P t b u y − P t s e l l Δ t (18) 0 ≤ P t b u y ≤ u t m P max 0 ≤ P t s e l l ≤ 1 − u t m P max (19) In the formula, p t represents the day-ahead transaction price; u t m is a Boolean variable indicating the purchase or sale state; and P max is the maximum exchange power.

In summary, the objective function of this model is to minimize the dispatch cost, expressed as Equation 20: min ⁡ o b j = C g + C s + C D R + C m (20)

Uncertain Parameters:

To address the conservativeness of robust optimization results, uncertain parameter variables and the maximum values of these parameters are incorporated into the model to limit the frequency of photovoltaic and load data reaching the worst-case scenarios, ensuring that the overall model result is not overly conservative Equation 21: B t p v , B t L o a d (21) The uncertain parameter is a Boolean variable indicating whether the worst-case scenario is reached, used to control the actual values of photovoltaic and load data in the optimization model. Typically, the worst-case scenario occurs when photovoltaic generation is insufficient while load demand is high. Thus, the model considers the lower bound deviation of photovoltaic generation and the upper bound deviation of load demand. P t p v = P t p v 0 − B t p v Δ P t p v P t L o a d = P t L o a d 0 + B t p v Δ P t L o a d (22) In the formula Equation 22, P t p v 0 , P t L o a d 0 represent the predicted values of photovoltaic and load data, respectively, while Δ P t p v , Δ P t L o a d denote the deviation values of photovoltaic and load power at time period t obtained through fuzzy clustering.

Additionally, the uncertain parameters need to satisfy boundary constraints Equation 23: B t p v ≤ B max p v B t L o a d ≤ B max L o a d (23) In the formula, B max p v , B max L o a d denote the upper limit of the uncertain parameter. The larger the upper limit, the more conservative the model results will be.

Based on the aforementioned variables and constraints, the optimization model is established. The primary objective at the regional grid-level is to use the output of generation units and pumped-storage power stations to smooth the unmet load curves of each province, thereby reducing the impact of peak loads on provincial-level dispatch. To achieve this goal, referring to the generation dispatch function of the Chinese power system, the minimization of the variance of unmet loads is selected as the objective function (Meng et al., 2023; Wang et al., 2023).

Thus, the objective function of this optimization problem is Equation 24: min ⁡ o b j p = ∑ t T P p , t ′ − ∑ i T P p , i ′ / T 2 (24) Since there is a corresponding load variance for each province, this model is a multi-objective optimization problem. The weighted sum method is used to combine multiple objectives into an equivalent scalar objective.

To avoid unreasonable results caused by large differences in load magnitudes among provincial grids, the total accepted generation output is normalized. Then, the objective functions of different provinces are integrated into a single objective using appropriate weighting coefficients.

Therefore, the composite objective function can be expressed as: min ⁡ o b j = ∑ p p r o v i n c e s w e i g h t p × o b j p ′ (25) In the formula Equation 25, o b j p ′ represents the normalized form of the objective function for province p. o b j p ′ = ∑ t T P p , t ′ − ∑ i T P p , i ′ / T max L o a d p 2 (26) In the formula Equation 26, w e i g h t p is the weight of province p in the objective function, and max L o a d p is the maximum value in the load curve of province p.

By substituting the sum of generation outputs into Equation 25, we obtain Equation 27: min ⁡ o b j = ∑ p p r o v i n c e s w e i g h t p max L o a d p 2 × ∑ t T ∑ g G P g , p , t − ∑ t T ∑ g G P g , p , t T 2 (27) Considering the separate optimization of peak and valley period variances, the two optimization objectives are weighted to transform the multi-objective problem into a single-objective problem again, as expressed by Equation 28: min ⁡ o b j = w e i g h t p e a k × o b j p e a k + w e i g h t v a l l e y × o b j v a l l e y (28) In the formula, w e i g h t p e a k , w e i g h t v a l l e y are the weights for peak and valley periods set in this study, and o b j p e a k , o b j v a l l e y correspond to the objective functions for peak and valley periods, respectively.

Thus, through the above formula transformations, the multi-objective problem is converted into a single-objective problem. In this model, the objective function is a quadratic function of the decision variables.

The power balance constraints in this model cannot be satisfied. However, research tests show that merely removing the power balance constraints leads to passive output from conventional generation units and non-unique solutions. The primary reason is that the relaxation of constraints is too excessive, requiring additional constraints or changes to the objective function.

This study adopts the penalty function method, adding penalty terms to the objective function in Equation 28, as shown in Equation 29: ϕ = ∑ p p r o v i n c e s ∑ t T max 0 , L o a d p , t − ∑ g G P g max × R g , p × L o a d p , t − ∑ g G P g , p , t (29) The construction of the penalty function is derived from the weighting method. In the formula, the max function is used to determine whether the maximum generation output of units at time period ttt can meet the load demand. If it cannot be met, unmet load is generated. The difference is used as a weight and multiplied by the difference between the actual unit output and the load. By adding this penalty term to the original objective function in Equation 28, the issue of passive output from conventional generation units can be resolved.

The objective function after adding the penalty term is shown in Equation 30: min ⁡ o b j = w e i g h t p e a k × o b j p e a k + w e i g h t v a l l e y × o b j v a l l e y + ϕ (30) where: P p , t ′ = ∑ g G P g , p , t (31) In the formula Equation 31, P p , t ′ denotes the total received generation output of grid p at time period t; P g , p , t is the power transmission of unit g to grid p at time period t; and T is the total number of time periods.

The highest degree term of the objective function in this model is quadratic. Therefore, this optimization problem is a quadratic programming problem. To prove that the model can converge to the global optimum, the compact form of the model is presented as follows Equation 32: min ⁡ x T a x + bx s . t . K eq x = k Dx ≥ d (32) In the compact form, K eq , D represent the constraint coefficient matrix in the model, and k , d are the constant matrixes corresponding to the constraints. Since the coefficient matrix a of the quadratic term in the objective function is a symmetric matrix, its eigenvalues are real numbers. By solving the eigenvalue equation, It can be concluded that all eigenvalues are non-negative, indicating that a is a positive semi-definite matrix. Thus, this problem is a convex quadratic optimization problem. For convex optimization, under the conditions that are satisfied, any local optimal solution is a global optimal solution, and this problem is solvable within polynomial time.

The provincial two-stage robust dispatch model constructed in this paper is mainly solved in two stages: the master problem and the sub-problem.

The compact form of the provincial-level dispatch problem is as follows:   min   x max u ∈ U min y ∈ F x , u cy s . t . Gy ≥ h − Ex − Mu G eq y = h eq − E eq x − M eq u (33) In the formula Equation 33, the variables represented by x Equation 34, are: x = u t s , u t m T (34) The variables represented by y Equation 35 are: y = P g , t , P t c h , P t d i s , P t D R , P t D R 1 , P t D R 2 , P t b u y , P t s e l l , P t p v , P t L o a d T (35) Equation 36 indicates the frequency of reaching the worst-case scenarios. u = B t p v , B t L o a d T (36) G , h , E , M and G eq , h eq , E eq , M eq respectively denote the coefficient matrices corresponding to inequality and equality constraints.

For the objective function, it is divided into two layers, the inner and the outer. The outer layer min x * represents finding the minimum value of the inner function and making decisions on the variables x , so that x can be used as a constant term in the inner layer’s solution. The inner function max u ∈ U min y ∈ F x , u cy represents the cost minimization result under the worst-case scenario. The feasible region y ∈ F x , u is defined by the feasible values of variables x and uncertain parameters u . The function min y ∈ F x , u cy represents the optimization result within this feasible region, while max u ∈ U * identifies the scenario with the worst-case photovoltaic and load values among all cost minimization results.

Based on the above analysis, the master problem and sub-problem can be organized as follows:

The master problem focuses on making decisions for the state variables of energy storage and power trading in the day-ahead stage, using the predicted photovoltaic and load data to optimize and obtain the lower bound of the final solution, as shown in Equation 37: min ⁡ a s . t . a ≥ cy Gy ≥ h − Mu − Ex G eq y = h eq − M eq u − E eq x (37) The sub-problem is solved in the intraday stage, where, based on the already determined values of x , the goal is to find the cost-minimizing decision under the worst-case scenario, as shown in Equation 38. max min cy s . t . Gy ≥ h − Mu − E x value G eq y = h eq − M eq u − E eq x value (38) For the two-stage robust optimization model above, this study uses the Column-and- Constraint Generation (C&CG) algorithm to solve it. The C&CG algorithm achieves the optimal solution of the original problem by decomposing it into a master problem and a sub-problem and solving them iteratively in an alternating manner. The main solution logic is illustrated in Figure 2.

In the process of solving the master problem using the C&CG algorithm, variables and constraints related to the sub-problem are continuously introduced, allowing for a tighter lower bound of the original objective function value, thereby effectively reducing the number of iterations.

The master problem and sub-problem of the provincial-level dispatch are given by Equations 37 and 38 above. During the iterative process of the C&CG algorithm, the uncertain parameters and their corresponding constraints are dynamically updated. Thus, the formal representations of the master problem and sub-problem should be expressed as follows:

Master Problem Equation 39: min ⁡ a s . t . a ≥ cy Gy ≥ h − M u k − Ex G eq y = h eq − M eq u k − E eq x (39) Subproblem Equation 40: max min cy s . t . Gy ≥ h − M u k − E x value G eq y = h eq − M eq u k − E eq x value (40) fundinWhere u k denotes the value of the uncertain parameter selected in the kth iteration. In the sub-problem, the variables x determined by the master problem are treated as constants.

Since the sub-problem is a bilevel problem, this study employs the KKT algorithm to reformulate the inner cost minimization problem into KKT conditions, transforming the bilevel problem into a single-level problem, as shown in Equation 41: max ⁡ c y s . t . Gy ≥ h − M u k − E x value G eq y = h eq − M eq u k − E eq x value G T π 1 + G eq T π 2 = c h − M u k − E x value − Gy π 1 = 0 G eq y − h eq + E eq x value + M eq u π 2 = 0 c − G T π 1 − G eq T π 2 y = 0 y ≥ 0 π 1 ≥ 0 (41) It should be noted that the problem contains nonlinear constraints. By using the Big-M method to linearize the constraints, the sub-problem is ultimately converted into a single-level optimization problem, as shown in Equation 42: max ⁡ c y s . t . Gy ≥ h − M u k − E x value G eq y = h eq − M eq u k − E eq x value G T π 1 + G eq T π 2 = c π 1 ≤ M ∞ v Gy − h + E x value + Mu ≤ M ∞ 1 − v π 2 ≤ M ∞ l G eq y − h eq + E eq x value + M eq u ≤ M ∞ 1 − l y ≤ M ∞ w c − G T π 1 − G eq π 2 ≤ M ∞ 1 − w y ≥ 0 π 1 ≥ 0 (42) In the formula, v , l , w represent auxiliary operators in the Big-M method, all of which are Boolean variables; M ∞ is a constant significantly larger than all parameters; π 1 , π 2 denote the Lagrange multipliers in the KKT transformation.