z 1 = ∑ t = 1 T R t , (7) where T is the total number of time intervals in the dispatch period. Taking 1 day as the dispatch period and each time interval as 15 min, then T = 96.

z 2 = ∑ t = 1 T ∑ l = 1 N L P Los , l , t 0 , (8) where N _L is the total number of branches in the system; and P Los , l , t 0 is the network loss of the branch l at time t, which can be expressed as follows (Lin et al., 2017): P Los , l , t 0 = g k m [ ∑ i ∈ ψ ( X k i − X m i ) P i , t 0 − ∑ r ∈ L ( X k r − X m r ) P ld r , t ] 2 , (9) where k and m are the buses at both ends of the branch l; g _km is the conductance of the branch l; X _ki , X _mi , X _kr , and X _mr are the elements of the node impedance matrix in the DC power flow model; ψ is the set of generator buses and injected power buses connecting external power grids; L is the set of load buses; and P _ldr,t is the active power of the load bus r at time t.

z 3 = ∑ t = 1 T [ ∑ i = 1 N 1 ( F i , t 0 + C i U , t + C i D , t ) + ∑ b = 1 N 2 F b , t 0 + ∑ s = 1 N 3 ( C s U , t + C s D , t ) ] , (10) where N ₃ is the total number of PSH units. C _iU,t /C _iD,t and C _sU,t /C _sD,t are the startup/shutdown cost of the TG unit i and PSH unit s at time t, which can be expressed as follows: { C i U , t ≥ K i ( I i , t − I i , t − 1 ) , C i U , t ≥ 0 C s U , t ≥ K s ( Z s , t − Z s , t − 1 ) , C s U , t ≥ 0 C i D , t ≥ J i ( I i , t − 1 − I i , t ) , C i D , t ≥ 0 C s D , t ≥ J s ( Z s , t − 1 − Z s , t ) , C s D , t ≥ 0 , (11) where K _i /J _i and K _s /J _s are the startup/shutdown cost of the TG unit i and PSH unit s. Z _s,t is the on/off state of the PSH unit s at time t, where 1/0 represents the on/off state.

z 4 = ∑ t = 1 T ( ∑ i = 1 N 1 C i , t P i , t 0 + ∑ w = 1 N w C w , t P w , t 0 ) , (12) where C _i,t and C _w,t are the power purchase prices of the TG unit i and the WF w at time t.

z 5 = ∑ t = 1 T ∑ i = 1 N 1 ( B i , 1 P i , t 0 + B i , 0 I i , t ) , (13) where B _i,1 and B _i,0 are the polluted gas emission coefficients of the TG unit i. Because of the small amount of polluted gas emissions of gas-fired units, their coefficients are B _i,1 = B _i,0 = 0.

The power balance constraint is expressed as Eq. 14: ∑ i = 1 N 1 P i , t s c + ∑ b = 1 N 2 ( P dis , b , t s c + P ch , b , t s c ) + ∑ s = 1 N 3 ( P pg , s , t s c + P pp , s , t s c ) = P nld , t s c + ∑ l = 1 N L P Los , l , t s c , (14) where P pg , s , t s c / P pp , s , t s c is the generating/pumping power of the PSH unit s at time t; P Los , l , t s c can be calculated similar to P Los , l , t 0 as in Eq. 9; and P nld , t s c is the system net load at time t, which is equal to the difference between the total load forecast value and the total WF output value, expressed as Eq. 15: P nld , t s c = ∑ r = 1 N LD P ld r , t − ∑ w = 1 N w P w , t s c . (15) The constraints of generators’ output are expressed as follows: { I i , t P i , min ≤ P i , t s c ≤ I i , t P i , max 0 ≤ P w , t s c ≤ P ^ w , t s c P i , t s c − P i , t − 1 s c ≤ r u i T 15 I i , t − 1 + P i , min ( I i , t − I i , t − 1 ) P i , t − 1 s c − P i , t s c ≤ r ⅆ i T 15 I i , t + P i , min ( I i , t − 1 − I i , t ) , (16) { P i , t 0 − P i , t s c ≤ r u i T 10 I i , t P i , t s c − P i , t 0 ≤ r ⅆ i T 10 I i , t   s c = 1 , ... , N s c , (17) where P _i,min/P _i,max is the minimum/maximum output of the TG unit i, and r _ui /r _di is its ramp up/down rate. T ₁₅ is the length of each time interval, which is 15 min. In Eq. 16, the unit maintains the minimum output in the first time interval after starting up and the last time interval before shutting down. T ₁₀ is the response time of the regulating power of the TG units, which is 10 min.

The constraints of minimum continual on/off time of TG units can be expressed as follows (Carrion and Arroyo, 2006): { I i , t = 1 , t ∈ [ 1 , U i ] , U i = min { T , ( T on i − X on i , 0 ) I i , 0 } ∑ n = t t + T on i − 1 I i , n ≥ T on i ( I i , t − I i , t − 1 ) , t ∈ [ U i + 1 , T − T on i + 1 ] ∑ n = t T [ I i , n − ( I i , t − I i , t − 1 ) ] ≥ 0 , t ∈ [ T − T on i + 2 , T ] I i , t = 0 , t ∈ [ 1 , D i ] , D i = min { T , ( T off i − X off i , 0 ) ( 1 − I i , 0 ) } ∑ n = t t + T off i − 1 ( 1 − I i , n ) ≥ T off i ( I i , t − 1 − I i , t ) , t ∈ [ D i + 1 , T − T off i + 1 ] ∑ n = t T [ 1 − I i , n − ( I i , t − 1 − I i , t ) ] ≥ 0 , t ∈ [ T − T off i + 2 , T ] , (18) where U _i /D _i is the number of time intervals during which the unit i must stay on/stay off at the beginning of the dispatch period, which is related to the unit states at the end of the last dispatch period; T _oni /T _offi is the minimum number of continual on/off time intervals of the unit i; and X _oni,0/X _offi,0 is the number of continual on/off time intervals of the unit i at the beginning of the dispatch period.

The network security constraints represented by the DC power flow model are as follows: { P l , k m , t s c = ∑ i ∈ ψ ( X k i − X m i x k m ) P i , t s c − ∑ r ∈ L ( X k r − X m r x k m ) P ld r , t − P ¯ l , k m ≤ P l , k m , t s c ≤ P ¯ l , k m , (19) where P l , k m , t s c is the transmission power of the branch l at time interval t; x _km is the reactance of the branch l; and P ¯ l , k m is the maximum transmission power of the branch l.

The constraints of the power output are as follows: { 0 ≤ P pg , s , t s c ≤ P pg , s , max ⋅ Z pg , s , t P pp , s , max ⋅ Z pp , s , t ≤ P pp , s , t s c ≤ 0 , (20) where P _pg,s,max/P _pp,s,max is the maximum generating/pumping power of the PSH unit s, and Z _pg,s,t /Z _pp,s,t is the generating/pumping state of the PSH unit s at time interval t.

The complementary constraint of operating states is as follows: Z s , t = Z pg , s , t + Z pp , s , t ≤ 1. (21)

The constraint of energy balance within a day is as follows: ∑ t = 1 T P pg , s , t s c + ξ ⋅ ∑ t = 1 T P pp , s , t s c = 0 , (22) where ξ is the conversion efficiency of PSH units, which is 75% in this study.

The constraints of switching time of operating states are as follows: { Z pg , s , t + Z pp , s , t + 1 ≤ 1 t = 1,2 , ⋯ , ( T − 1 ) Z pg , s , t + Z pp , s , t + 2 ≤ 1 t = 1,2 , ⋯ , ( T − 2 ) Z pp , s , t + Z pg , s , t + 1 ≤ 1 t = 1,2 , ⋯ , ( T − 1 ) Z pp , s , t + Z pg , s , t + 2 ≤ 1 t = 1,2 , ⋯ , ( T − 2 ) . (23)

The constraints of discharging/charging power are as follows: { 0 ≤ P dis , b , t s c ≤ P dis , b , max ⋅ Z dis , b , t P ch , b , max ⋅ Z ch , b , t ≤ P ch , b , t s c ≤ 0 , (24) where P _dis,b,max/P _ch,b,max is the maximum discharging/charging power of the BES station b, and Z _dis,s,t /Z _ch,s,t is the discharging/charging state of the BES station b at time interval t.

The complementary constraint of operating state is as follows: Z b , t = Z dis , b , t + Z ch , b , t ≤ 1. (25)

The constraints of the remaining energy are as follows: { E b , t s c = E b , t − 1 s c + ( − P dis , b , t s c / η dis , b − P ch , b , t s c η ch , b ) T 15 0 ≤ E b , t s c ≤ E b , max , (26) where η _dis,b /η _ch,b is the discharging/charging efficiency of the BES station b, and E b , t s c is the stored remaining energy of the BES station b at time interval t.

The system flexibility demand needs to consider not only the uncertainties of the load power and WF output in each time interval but also the change of the net load power in adjacent time intervals. Thus, the upward and downward system flexibility demand constraints are expressed as follows: { P nld , t + 1 0 − P nld , t 0 + L u % P ld , t + 1 + w u % ∑ w = 1 N w P w , t + 1 0 ≤ Δ P fu , t 0 + Δ P lku , t 0 P nld , t 0 − P nld , t + 1 0 + L d % P ld , t + 1 + w d % ∑ w = 1 N w ( P ^ w , t + 1 − P w , t + 1 0 ) ≤ Δ P fd , t 0 + Δ P lkd , t 0 , (27) where L _u%/L _d% is the percentage of the upward/downward system flexibility demand due to the load forecast deviation; w _u%/w _d% is the percentage of the upward/downward system flexibility demand due to the wind power forecast deviation; and Δ P fu , t 0 / Δ P fd , t 0 is the upward/downward system flexibility supply capacity at time interval t.

The system flexibility supply capacity is the sum of the flexibility supply capacity provided by all units, as follows: { Δ P fu , t 0 = ∑ i = 1 N 1 Δ P u , i , t 0 + ∑ b = 1 N 2 Δ P u , b , t 0 + ∑ s = 1 N 3 Δ P u , s , t 0 Δ P fd , t 0 = ∑ i = 1 N 1 Δ P d , i , t 0 + ∑ b = 1 N 2 Δ P d , b , t 0 + ∑ s = 1 N 3 Δ P d , s , t 0 , (28) where Δ P u , i , t 0 / Δ P u , s , t 0 / Δ P u , b , t 0 is the upward flexibility supply capacity provided by the TG unit i/PSH unit s/BES station b at time interval t, and Δ P d , i , t 0 / Δ P d , s , t 0 / Δ P d , b , t 0 is the corresponding downward flexibility supply capacity.

The flexibility supply capacities provided by different units can be expressed as follows: { 0 ≤ Δ P u , i , t 0 ≤ min ( P i , max I i , t − P i , t 0 , r u i T 10 I i , t ) 0 ≤ Δ P d , i , t 0 ≤ min ( P i , t 0 − P i , min I i , t , r d i T 10 I i , t ) 0 ≤ Δ P u , s , t 0 ≤ P pg , s , max Z pg , s , t − P pg , s , t 0 − P pp , s , t 0 0 ≤ Δ P d , s , t 0 ≤ P pg , s , t 0 + P pp , s , t 0 − P pp , s , max Z pp , s , t 0 ≤ Δ P u , b , t 0 ≤ P dis , b , max − ( P dis , b , t 0 + P ch , b , t 0 ) 0 ≤ Δ P d , b , t 0 ≤ P dis , b , t 0 + P ch , b , t 0 − P ch , b , max . (29)

Taking z ₁, z ₂, z ₃, z ₄, and z ₅ as the optimization objectives of the SCUC model, we used the GUROBI solver to solve each single-objective optimization problem (SOP) after convex relaxation of quadratic terms of network loss. We compared the obtained network loss after optimization with the real network loss calculated by substituting the obtained optimal power outputs of units into Eq. 9, and the results are given in Table 1. It is evident that the difference between the network loss obtained by the convex relaxation method and the real network loss was very small, and the maximum relative error was only 0.224%, which demonstrated that the proposed convex hull relaxation method for the quadratic term of network loss had high computational accuracy.

We used the SBB solver to solve the original MINNP of the SCUC model. We compared the CPU time of solving the MINNP model with that of solving the MICP model, as shown in Table 2. We observed that the CPU time of the MICP model decreased more than 75%. Moreover, for the SOP with z ₄ as the objective, the optimal solution could not be obtained by solving the MINNP model, but it could be obtained by solving the MICP model. Thus, the proposed convex hull relaxation method could significantly improve the computational efficiency, and the optimal solution of the SCUC model could be obtained more reliably.

We analyzed the correlation between the objective functions. We solved each of the five SOPs, and the objective function values of the five optimal solutions are shown in Table 3. We used these five optimal solutions to calculate Spearman correlation coefficients between the optimization objectives and the sum of the negative correlation coefficients of each optimization objective. The results are given in Supplementray Table A2. From Supplementray Table A2, we selected the objective with the largest absolute value of the sum of negative correlation coefficients (i.e., z ₅) as the first-level conflicting objective and selected the objective with the next largest absolute value of the sum of negative correlation coefficients (i.e., z ₄) as the second-level conflicting objective. Thus, we obtained the most conflicting objective set S _c = [z ₄, z ₅] by the objective selection method, and z ₅ and z ₄ were the upper layer optimization objectives.

Then, we solved the Pareto frontier of the five-objective SCUC model by the two proposed ε-constraint methods and the NSGA-II algorithm. Considering the value differences of different objectives, we defined the color column coordinates k ^* = z ₄ ^*/z ₅ ^* = (z ₄/z _4min)/(z ₅/z _5min); thus, as k ^* increases, z ₄ increases and z ₅ decreases.

1) Grid-based five-objective ε-constraint method.

We calculated z _5max and z _5min and obtained the four points of z ₅ by setting q ₁ = 3. Then, for each z _5k , according to Eq. 38 and Eq. 38, we obtained z _4,5kmax and z _4,5kmin and obtained the four points of z ₄ by setting q ₂ = 3. Thus, the total number of grid points of the upper layer two objectives was 16. The Pareto frontier surface of the lower layer three-objective optimization model (Eq. 40) corresponding to each grid point was solved. After gathering all the obtained optimal solutions and removing the dominated solutions, we obtained the Pareto frontier composed of 427 non-dominated solutions, as shown in Figure 2. The COS is also obtained and is shown as the red-colored point in Figure 2.

2) Curve-based five-objective ε-constraint method.

We used the NNC method to solve the upper layer two-objective optimization model (Eq. 41). Letting m = 11, the obtained Pareto frontier curve is shown in Supplementary Figure A2. We added the kth point (z _4,k , z _5,k ) values to the lower layer three-objective optimization model (Eq. 42) as constraints and solved the Pareto frontier surface of the model (Eq. 42) corresponding to each point (z _4,k , z _5,k ). Finally, we obtained the Pareto frontier composed of 595 non-dominated points, as shown in Figure 3. We also obtained the COS, which is shown as the red-colored point in Figure 3.

By comparing the computational process and results of the aforementioned two methods, it can be seen that the curve-based ε-constraint method obtained more non-dominated solutions and more compact Pareto frontier surfaces. This is because when the curve-based ε-constraint method was used to solve the five-objective SCUC model, the selected points were within the Pareto frontier curve of the upper layer two objectives and had a good non-dominated relationship. However, for the grid-based ε-constraint method, although the hierarchical division of grid points seemed to be uniform, it ignored the conflicting relationship between these two objectives, resulting in more dominated points being computed in the lower layer three-objective optimization. Thus, the curve-based ε-constraint method could obtain better results because it could obtain more Pareto optimal solutions and supply more trade-off relationship information of different objectives for the final decision of the five-objective SCUC problem.

3) NSGA-II algorithm.

We used the gamultiobj (.) function of the NSGA-II algorithm in MATLAB to solve the Pareto frontier of the five-objective SCUC model, which included 100 non-dominated solutions, as shown in Figure 4. The obtained Pareto optimal solutions were concentrated and supplied with little trade-off relationship information for the final decision of the five-objective SCUC problem.

We compared the CPU time and COSs of the three methods. The CPU times of the grid-based ε-constraint method, curve-based ε-constraint method, and NSGA-II algorithm were 56,703.958 s, 33,561.423 s, and 117,618.003 s, respectively. The two proposed methods had higher computational efficiency, among which the curve-based ε-constraint method had the highest computational efficiency. Since different Pareto optimal solutions can be calculated independently in the two proposed methods, thus, the parallel calculation technique for opening eight threads in the eight cores of the used computer was used to improve the computational efficiency. By the parallel calculation, the CPU times of the grid-based ε-constraint method and the curve-based ε-constraint method were reduced to 7,260.439 s and 4570.952 s, and the acceleration ratios were 7.810 and 7.342. These demonstrated that the two proposed methods had good computational performance for solving the five-objective SCUC model, and they were more suitable for solving MaOPs. The COSs of the three methods are shown in Table 4. The λ value corresponding to the curve-based ε-constraint method was the largest, which indicated that the obtained COS was the best in the coordinate optimization of the five objectives. Thus, this analysis showed that the curve based ε-constraint method could obtain the Pareto frontier and the COS more efficiently, and the obtained COS had a better optimization degree.

The power output schedule of the units corresponding to the COS obtained by the curve-based ε-constraint method is shown in Figure 5. The PSH and BES stations pumped water and charged energy, respectively, during the low-load time intervals, and they generated electricity and discharged energy, respectively, during the high-load time intervals, which has high flexibility for regulating power rapidly.

Comparative results of the network loss of the five SOPs are shown in Table 5, demonstrating the high computational accuracy of the proposed convex hull relaxation method for the quadratic term of network loss in large-scale SCUC problems. Comparison of the CPU time of the original MINNP model with the MICP model is shown in Table 6. The CPU time of the MICP model decreased by more than 90%, and the convergence problem of the MINNP model was also addressed, which demonstrated the high-computational efficiency and convergence of the proposed convex hull relaxation method.

We calculated Spearman correlation coefficients among the five objective functions and the sum of the negative correlation coefficients of each objective, which are shown in Supplementary Table B2. The absolute sums of the negative correlation coefficients of z ₄ and z ₅ were the largest. Thus, we selected S _c = [z ₄, z ₅] as the first-level conflicting objective and z ₄ as the second-level conflicting objective.

The NSGA-II algorithm cannot solve this large-scale five-objective SCUC problem because of the huge size of the optimization model. Thus, we compared only the results of the two proposed ε-constraint methods. For the grid-based five-objective ε-constraint method, by setting q ₁ = q ₂ = 2, the total number of grid points of the upper layer two objectives was nine. The Pareto frontier surface of the lower layer three-objective optimization model corresponding to each grid point was obtained by solving Eq. 40. The final Pareto frontier of the five-objective SCUC model is shown in Figure 6. For the curve-based ε-constraint method, the NNC method was used to solve the Pareto frontier curve of the upper layer two-objective optimization model (Eq. 41) by setting m = 6. Then, the Pareto frontier surface of the lower layer three-objective optimization model corresponding to each point (z _4,k , z _5,k ) was obtained by solving Eq. 42. The final Pareto frontier of the five-objective SCUC model is shown in Figure 7. From Figures 6, 7, it can be seen that more Pareto optimal solutions were obtained by the curve-based ε-constraint method than by the grid-based ε-constraint method. We compared the CPU time of the two methods, the CPU times of the grid-based and curve-based ε-constraint methods by the parallel calculation technique were 24765.4 and 16,792.8 s, respectively. This demonstrated that the curve-based ε-constraint method had higher computational efficiency than the grid-based ε-constraint method. Obviously, if the used computer has more cores, the computational efficiency of the proposed curve-based ε-constraint method will be further increased to satisfy the requirement of actual application. Although the grid-based ε-constraint method solved a higher number of lower layer three-objective optimization models and consumed more CPU time, many dominated solutions in the obtained results had to be deleted in the final Pareto frontier of the five-objective SCUC model. Thus, the curve-based ε-constraint method had better results and higher computational efficiency than the grid-based ε-constraint method.

The objective functions of the COSs obtained by the two proposed methods and the five single-objective optimal solutions are shown in Table 7. The λ value of each solution is also listed in the table. The COS obtained by the curve-based ε-constraint method was the best because it had the largest λ value. The results demonstrated that the proposed curve-based ε-constraint method could obtain a COS with a better comprehensive optimization degree in the coordinate optimization of the five objectives.

The power output schedules of several units and the BES station corresponding to the COS of the curve-based ε-constraint method are shown in Figure 8. The schedules reflected the high flexibility of PSH and BES stations to rapidly regulate their power output.