Distributionally Robust Unit Commitment with N-k Security Criterion and Operational Flexibility of CSP

Abstract

In order to reduce the conservatism of the robust optimization method and the complexity of the stochastic optimization method and to enhance the ability of power systems to deal with occasional line fault disturbance, this paper proposes a distributionally robust unit commitment (DRUC) model with concentrating solar power (CSP) operational flexibility and N-k safety criterion under distributed uncertainty. According to the limited historical sample data, under the condition of satisfying a certain confidence level, based on the imprecise Dirichlet model (IDM), an ambiguity set is constructed to describe the uncertainty of transmission line fault probability. Through the identification of the worst probability distribution in the ambiguity set, the adaptive robust optimal scheduling problem is transformed into a two-stage robust optimization decision model under the condition of deterministic probability distribution. The CSP flexibility column and constraint generation (C&CG) algorithm is used to process the model and the main problem and subproblem are solved by using the Big-M method, linearization technique, and duality principle. Then, a mixed integer linear programming problem (MILP) model is obtained, which effectively reduces the difficulty of solving the model. Finally, case studies on the IEEE 14 bus system and the IEEE 118 bus system demonstrate the efficiency of the proposed method, such as enhancing the ability of power systems to cope with occasional line fault disturbances and reducing the conservatism of the robust optimization method.

1. Introduction

The intermittency, volatility, and uncertainty of new energy generation lead to probabilistic power flow in new power systems [1,2]. It is difficult to accurately estimate the impact of transmission line faults combined with probabilistic power flow fluctuations on the safe operation of the power system [3]. Therefore, the N-k safety criterion considering k line faults should be applied to the economic scheduling, and the optimal scheduling method should be sought on the premise of ensuring the safe operation of the system, taking into account both safety and economy.

Reference [4] proposed a robust optimization model considering N-k faults, but due to the limitations of transmission capacity and other factors, the solution error is large. In reference [5], fault probability information and αcut criterion were combined to solve the model. In order to describe the uncertainty of faults, reference [6] proposes a two-stage robust model, and reference [7] proposes a power system unit combination model based on N-k-ε safety criterion. The probability distribution information is often ignored in the abovementioned pieces of literature, and the number of line faults is assumed to be a deterministic constant, which leads to unreasonable scheduling decisions.

In reality, it is difficult to accurately obtain the probability distribution of line faults, and the probability distribution obtained by mining the law of historical data often has errors, that is, there is a high-order uncertainty problem [8]. In this case, the deterministic line failure decision-making based on optimum scheduling accuracy is easily affected. At present, distributionally robust optimization methods have been employed to solve the optimal scheduling problem of power systems under uncertain conditions. This method combines the advantages of both robust optimization and stochastic optimization, which can be obtained by the statistical information (such as the first moment, the second moment, etc.) In [9], the possible probability distribution functions of random variables are described and constitute the so-called ambiguity set, which is adopted to quantify the statistical laws of uncertain quantification. Then, we can make decisions that are immune to the uncertainty of distribution, that is, we can find the best random decision under the condition of the worst probability distribution in the ambiguity set. Evidently, the distributionally robust optimization method has a better ability to describe realistic decision scenarios.

The construction of ambiguity sets is a key factor affecting the effect of distributionally robust optimization. The construction method of ambiguity set includes the moment information of known uncertain quantity [8,9]. In addition, the distance of the known distribution function [10], data-driven methods [11], and how to accurately construct the ambiguity set and dock it with the traditional robust scheduling model is the key. There have been many studies on the application of the distributionally robust optimization method to solve the uncertain optimal scheduling problem [12,13,14,15], but there are still gaps in solving the uncertain problem of N-k security criteria with the distributionally robust optimization method. Therefore, it is very important to extract the statistical information of transmission line faults from the existing data, incorporate the N-k safety criterion into the economic scheduling, and make economic scheduling decisions that take both economy and security into account.

One of the most important reasons to promote the development of new energy generation is to cope with the climate crisis. However, the uncontrollability of new energy needs the support of flexible regulatory resources. If flexibility is provided by thermal power units, the overall carbon emissions of the power system will be increased. Concentrating solar power (CSP) uses the storage system to store the heat generated by the conversion of solar energy and adjusts the heat release of the heat storage system according to the change in light and load demand to promote the stable output of the turbine, thereby achieving the purpose of controlling the light energy. The generation principle of CSP is similar to that of thermal power plants [16]. In addition, CSP has the ability of rapid adjustment, and the climbing rate can reach 20% per minute, which is higher than that of traditional thermal power units [17]. Giving full play to the energy storage efficiency and flexibility of CSP can improve the economy of the system operation and effectively reduce the carbon emissions of the power system.

Therefore, this paper proposes a distributionally robust unit commitment (DRUC) method based on the imprecise Dirichlet model (IDM) under distribution uncertainty, which takes CSP flexibility and N-k safety criteria into account. An ambiguity set in the form of a probability distribution interval is constructed to describe the worst-case distribution of line fault uncertainty. The Devroye–Wise method is presented to describe the estimation range of the real values of the random variables, and then a two-stage distributionally robust economic scheduling model is constructed taking into account the flexibility of CSP and N-k safety criteria. The Big-M method and column and constraint generation (C&CG) algorithm are proposed to transform and solve the main problem and sub-problem.

The main novelties of this paper are therefore:

• A DRUC with N-k security criterion formwork that can address the line fault distributional uncertainty is proposed for the first time. An imprecise Dirichlet model (IDM)-based non-parametric ambiguity set construction method is proposed. The Devroye–Wise method is presented to describe the estimation range of the real values of the random variables, thereby being able to reflect the available information more objectively.
• The potential flexibility of CSP can be fully tapped into due to the storage system to store the heat generated by the conversion of solar energy and it adjusts the heat release of the heat storage system according to the change of light and load demand to promote the stable output of the turbine, thereby achieving the purpose of controlling the light energy.

2. Construction of the Uncertain Set Based on IDM

2.1. IDM Model

IDM is an extension of the deterministic Dirichlet model [18], which uses a single Dirichlet distribution as a prior distribution to evaluate the probability of occurrence of various states of random variables. The corresponding single prior probability density function can be expressed as:

$$f_1=\frac{\Gamma ({\displaystyle \sum _{n=1}^N{\beta }_n})\underset{n=1}{\overset{N}{\Pi }}{\theta }_n^{{\beta }_n-1}}{\underset{n=1}{\overset{N}{\Pi }}\Gamma ({\beta }_n)}$$

(1)

where ${\theta }_n$ is the exact single-valued probability of the random variable appearing in the NTH state, which satisfies the constraint conditions of ${\theta }_n\ge 0,n=1,2,\dots ,N$ and ${\displaystyle \sum _{n=1}^N{\theta }_n}=1$; N is the total number of possible states of the random variable; ${\beta }_n$ is the positive parameter in the Dirichlet distribution; $\Gamma$ is the Gamma function.

According to the Bayesian statistical theory, the total number of samples observations obtained by Dirichlet distribution is M, and as the prior probability density function is updated by the Bayesian process, the posterior probability density function expression is obtained as follows:

$$f_2=\frac{\Gamma ({\displaystyle \sum _{n=1}^N(m_n+{\beta }_n)})\underset{n=1}{\overset{N}{\Pi }}{\theta }_n^{{}^{m_n+s\mathbf{·}{r}_n-1}}}{\underset{n=1}{\overset{N}{\Pi }}\Gamma ({\beta }_n+m_n)}$$

(2)

Among them:

$${\theta }_n=\frac{({\beta }_n+m_n)}{{\displaystyle \sum _{n=1}^N{\beta }_n+m_n}}$$

(3)

$$\forall r_n\in [0,1],{\displaystyle \sum _{n=1}^N{r}_n=1}$$

(4)

where $m_n$ is the occurrence times of any state of the random variable, and $r_n$ is the prior weight of each state of the random variable, which is the mean of ${\theta }_n$.

Once the sample observation data is missing, then $m_n=0$, at this time, ${\theta }_n={\beta }_n/{\displaystyle \sum _{n=1}^N{\beta }_n}$, the parameter ${\beta }_n$ directly determines the probability ${\theta }_n$; ${\displaystyle \sum _{n=1}^N{\beta }_n}$ in (3) can be replaced by the parameter s, and its value range is usually [19]. Parameter s reflects the relationship between the prior information and the posterior probability [20]. In order to eliminate the influence caused by this relationship, more sample observation data should be used [21].

If $s\mathbf{·}{r}_n$ has the same effect as ${\beta }_n$ in Equation (1), then ${\beta }_n=s\mathbf{·}{r}_n$. When parameter s is determined, the corresponding parameter ${\tilde{\theta }}_1,{\tilde{\theta }}_2,\dots ,{\tilde{\theta }}_n$ of each interval probability can be obtained through Equation (5) [18].

$${\tilde{\theta }}_n=[\frac{m_n}{M+s},\frac{m_n+s}{M+s}]$$

(5)

For comparison, Equation (3) is rewritten as:

$${\theta }_n=\frac{m_n+s\mathbf{·}{r}_n}{M+s}$$

(6)

Compare Equations (5) and (6), in the given interval, and take any parameter $r_n$ value, and the width of the probability interval is $s/M+s$ and the determining factors are parameter s and sample total M. If the total number of samples increases, the width of the probability interval shrinks and gradually converges to the exact single-value probability. Therefore, the size of the total number of samples can be reflected by the width of the probability interval. The upper and lower bounds of the interval probability estimated by IDM do not contain prior information, and the probability interval is used to replace the exact single-valued probability, which can eliminate the error caused by unreasonable prior information on the probability estimation result. Refer to the literature [22] for more descriptions of imprecise probabilities.

2.2. Construction and Transformation of Ambiguity Sets for Line Fault Probability Distribution

Due to the lack of effective samples of line faults, it is difficult to accurately represent the information of line fault probability with accurate single-value probability [23]. As can be seen from Equation (6), the prior weight $s\mathbf{·}{r}_n$ determines the probability of each situation in the random variable. Under the condition of small samples, it is impossible to give the exact $r_n$, which will lead to the deviation of probability estimation. Therefore, the imprecise probability is obtained by data statistics, among which the imprecise Dirichlet model is the most effective and widely used one. Therefore, all sets of Dirichlet model distributions are used for estimation. Under the condition that the parameter s = 1, the set consists of all Dirichlet distributions, and the prior parameters $r_n$ are traversed through the interval [0, 1].

The larger the sample size, the narrower the width of the IDM probability interval, the narrower the ambiguity set, and the higher the accuracy of the probability interval. This method avoids the problem of inaccurate probability estimation of events caused by unreasonable prior setting in the deterministic Dirichlet model.

The advantage of ambiguity set $\rho$ is that it can represent the information of the fault probability distribution interval without assuming any prior knowledge [24]. Figure 1 shows that the richer the number of samples, the narrower the resulting ambiguity set.

At the same time, when the random variable is the number of line faults, the corresponding interval $[k_l,k_u]$ of the number of line faults can be calculated by the probability points $[(1-\gamma )/2,(1+\gamma )/2]$ on the confidence band of the known cumulative probability distribution function of line faults, where $\gamma$ is ${\theta }_n$ represents the minimum value by the confidence interval that can reach a given confidence level (generally $\gamma$ = 0.95), and $k_u$ and $k_l$ are the upper and lower limits of the line fault number interval, respectively.

Figure 1b shows the process of transforming the ambiguity set into the uncertain set. It can be seen that the ambiguity set constructed by the IDM method can be mapped to the boundary of the uncertain set of line faults. By using this method, the uncertainty interval of line faults can be constructed without making assumptions on the distribution of line fault data, which reduces the dependence on the statistics of line fault history data and improves the applicability of the model.

Referring to the historical operation data of a certain area in the past five years as the research sample, the outage time caused by line faults in this area has been about 6 h since the use of the power grid, and the average fault probability of the transmission line in this grid can be calculated as $1.36986\times {10}^{-4}$.Then, Equation (5) can be used to calculate the inexact probability of line failure, so the interval probability is approximately $[1.37\times {10}^{-4},1.6\times {10}^{-4}]$.

Considering the influence of time on the line faults, without considering the association between faulty lines, the multi-period independent uncertain set Z based on IDM can be expressed as follows.

$$Z=\left\{{\displaystyle \sum _{l=1}^L(1-z_{l,t})}\le k,z_{l,t}\in \{0,1\}\right\}$$

(7)

$$k_{\min }\le k\le k_{\max }$$

(8)

where L is the total number of transmission lines; k is the number of faults in the line; $k_{\min }$ and $k_{\max }$ are the minimum value and maximum value of the number of line faults k, respectively; $z_{l,t}$ is the running state of the line during the time period t. The value of 1 indicates the normal operation of the line, and the value of 0 indicates the fault of the line. $t\in N$ and $t\in [0,24]$.

3. Operation Mechanism of Concentrating Solar Power

The photothermal power generation system mainly includes four types of luminescence systems: disc, trough, tower, and linear Fresnel. Its structure is shown in Figure 2, which is mainly composed of three parts: heat collection, energy storage, and a power generation subsystem. The power generation principle of various types of CSP power stations is almost the same, as shown in Figure 3. The operation principle of photothermal power generation is to use a collector to collect solar energy, conduct heat to heat transfer materials, then produce high-temperature water vapor through heat exchanger, and finally use this water vapor to drive turbines to generate electricity.

In this paper, the operation mechanism of a trough CSP power station is analyzed. Under light conditions, the accumulated solar heat heats the hot oil in the tubes, and some of the hot oil is added directly to the hot water to generate high-temperature and high-pressure steam to drive the turbine. The other part is used to heat nitrate in cryogenic tanks and then it is stored in a hot tank. When the light irradiation is insufficient, the heat conduction oil is heated by nitrate and converted into high-temperature and high-pressure steam to drive the turbine to operate. The temperature of the nitrite in the heating tank is reduced and stored in the cold tank. In addition, the turbine can obtain energy directly from the mirror field, from the heat storage area, or from both parts. The energy flow in this process is shown in Figure 4.

The energy balance relationship inside the CSP power station is shown as follows:

$$Q_t^{\mathrm{G}}=k_t^{\mathrm{SF}}(Q_t^{\mathrm{SF}}-Q_{t.\mathrm{LOSS}}^{\text{O-W}})+k_t^{\mathrm{TS}}(Q_t^{\mathrm{TS}}-Q_{t.\mathrm{LOSS}}^{\text{BN-O}}-Q_{t.\mathrm{LOSS}}^{\text{O-W}})$$

(9)

$$Q_t^{\mathrm{SF}}={\eta }_{\mathrm{SF}}{S}_{\mathrm{SF}}{R}_t{Q}_t^{\mathrm{D}}$$

(10)

$$\Delta Q_t^{\mathrm{TS}}=(1-k_t^{\mathrm{SF}})(Q_t^{\mathrm{SF}}-Q_{t.\mathrm{LOSS}}^{\text{O-BN}})-k_t^{\mathrm{TS}}{Q}_t^{\mathrm{TS}}-Q_{t.\mathrm{LOSS}}^{\mathrm{TS}}$$

(11)

In the CSP power station, the relationship between the input and output energy of the turbine is shown in (12):

$$P_t^{\mathrm{CSP}}=\frac{Q_t^{\mathrm{G}}{\eta }_e{\eta }_m{\eta }_g}{3600}$$

(12)

4. IDM-Based Distributionally Robust Economic Scheduling Considering CSP Flexibility and N-k Safety Criteria

4.1. Objective Function

The objective function is to minimize the unit start–stop cost, output cost, and loss of load penalty cost in the worst case of line fault, which can be expressed as:

$$\underset{u_{g,t},v_{g,t}}{\min }{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{g=1}^G(C_{g,\mathrm{U}}{u}_{g,t}+C_{g,\mathrm{D}}{v}_{g,t})+}}\underset{z\in \widehat{Z}}{\max }Q$$

(13)

$$Q=\min {\displaystyle \sum _{t=1}^T({\displaystyle \sum _{g=1}^G{\lambda }_g{p}_{g,t}}}+{\displaystyle \sum _{i=1}^I{C}_{\mathrm{VOLL}}{d}_{i,t}})$$

(14)

where T is the number of time periods; G is the total number of units; I is the total number of nodes; $C_{\mathrm{VOLL}}^{}$ is the loss of load penalty cost; $C_{g,\mathrm{U}},C_{g,\mathrm{D}}$ is the opening and closing cost of unit G; $u_{g,t}$ is the startup state of unit G in time period t, where 1 indicates unit startup and 0 indicates unit shutdown. $v_{g,t}$ is the closed state of unit G in time period t, where the value 1 means the unit is out of service, and the value 0 means the unit is running. ${\lambda }_g$ is the linear cost coefficient of unit output; $p_{g,t}$ is the output of unit g at time t; $d_{i,t}$ is the unbalance power of node i in time period t, $i\in I$.

4.2. Unit Combination Constraints

(1) Constraints on start–stop status and minimum start–stop time of unit

$$y_{g,t}-y_{g,(t-1)}-u_{g,t}\le 0$$

(15)

$$y_{g,(t-1)}-y_{g,t}-v_{g,t}\le 0$$

(16)

$$(y_{g,(t+1)}-y_{g,t})H_{g,\mathrm{on}}-{\displaystyle \sum _{k=t+2}^{\min \left\{t+H_{g,\mathrm{on}},T\right\}}{y}_{g,t}}\le \max \left\{1,H_{g,\mathrm{on}}-T+t-1\right\}$$

(17)

$$(y_{g,t}-y_{g,(t+1)})H_{g,\mathrm{off}}-{\displaystyle \sum _{k=t+2}^{\min \left\{t+H_{g,\mathrm{off}},T\right\}}{y}_{g,t}}\le H_{g,\mathrm{off}}$$

(18)

$$y_{g,t},v_{g,t},u_{g,t}\in \{0,1\}$$

(19)

where $H_{g,\mathrm{on}}$ and $H_{g,\mathrm{off}}$ are the minimum start-up time and minimum shutdown time of unit G, respectively, $y_{g,t}$ is the operating state of unit G in time period t, whose value is 1, indicating unit operation, and whose value is 0, indicating unit outage; In Equations (15), (16), and (19), $g\in G,t\in T$. In Equations (17) and (18), $t=1,\dots ,T-2$.

(2) Interpolated linearized generation cost constraint

Since the unit generation cost is a quadratic function, the interpolated linearization method [10] is used here to represent:

$$\{\begin{array}{l}{\displaystyle \sum _{r=1}^R{\lambda }_{r,g,t}^{}}{p}_{r,g}^{}=p_{g,t}\\ {\displaystyle \sum _{r=1}^R{\lambda }_{r,g,t}^{}}=u_{g,t}\end{array}$$

(20)

where r is the number of interpolation points.

4.3. CSP Internal Constraints on Power Plant Operation

(1) Power generation output constraints:

$$\begin{array}{c}{P}_{s,\min }{I}_{t,s}^{}\le P_{t,s}^{}\le P_{s,\max }{I}_{t,s}^{}\\ I_{t,s}\in \{0,1\}\end{array}$$

(21)

where $I_{t,s}^{}$ is the on/off state of the CSP power plant device at time period t, which is represented by binary variables; $P_{t,s}^{}$ is the scheduled output of CSP power station at time t; $P_{s,\min },P_{s,\max }$ is the lower limit and upper limit of CSP power output; ${\Omega }_T$ is a group of periods of a day.

(2) Minimum on/off time period constraints:

$$\begin{array}{l}(I_{t-1,s}^{}-I_{t,s}^{})T_{s,\mathrm{off}}+{\displaystyle \sum _{j=t-T_{s,\mathrm{on}}}^{t-1}(1-I_{j,s}^{})\ge 0;}\\ (I_{t,s}^{}-I_{t-1,s}^{})T_{s,\mathrm{on}}+{\displaystyle \sum _{j=t-T_{s,\mathrm{on}}}^{t-1}{I}_{j,s}^{}\ge 0;}\end{array}$$

(22)

where $T_{s,\mathrm{off}}$ is the period of CSP power plant unit shutdown; $T_{s,\mathrm{on}}$ is the period when the CSP power station unit is turned on.

(3) Instantaneous thermal power constraint:

$$P_{t,s}^{}/{\eta }_{\mathrm{PB},s}^{}+P_{\mathrm{cha},t,s}^{}-P_{\mathrm{dis},t,s}^{}\le P_{t,\mathrm{fore},s}^{}$$

(23)

where ${\eta }_{\mathrm{PB},s}^{}$ is the power efficiency of the CSP power plant device; $P_{\mathrm{cha},t,s}^{}$ is the charging output power of the CSP power plant at time period t; $P_{\mathrm{dis},t,s}^{}$ is the emission of the heat storage system in the CSP power plant at time period t; $P_{t,\mathrm{fore},s}^{}$ is the available solar thermal power in time period t.

(4) Impulse power constraint of heat storage tank:

$$0\le P_{\mathrm{cha},t,s}\le u_{t,\mathrm{c},s}{P}_{s,\mathrm{ch},\max }$$

(24)

$$0\le P_{\mathrm{dis},t,s}\le u_{t,\mathrm{f},s}{P}_{s,\mathrm{d},\max }$$

(25)

$$u_{t,\mathrm{c},s}+u_{t,\mathrm{f},s}\le 1$$

(26)

where $P_{s,\mathrm{ch},\max }$ and $P_{s,\mathrm{d},\max }$ are the maximum punching power and the maximum releasing power of the heat storage tank, respectively. Ensure that the state of flushing and discharging cannot exist at the same time at a certain moment. $u_{t,\mathrm{c},s}$ and $u_{t,\mathrm{f},s}$ are binary variables of charging and discharging state and thermal state of the heat storage system, respectively (1 for action, 0 for no action).

(5) Constraint of charge and heat state:

$$E_{t,s}^{}=E_{t-1,s}^{}+{\eta }_{\mathrm{TES},s}^{\mathrm{ch}}{P}_{\mathrm{cha},t,s}^{}-P_{\mathrm{dis},t,s}^{}/{\eta }_{\mathrm{TES},s}^{\mathrm{d}}$$

(27)

where $E_{t,s}^{}$ is the state of charge of heat storage system in the CSP power station at time period t; ${\eta }_{\mathrm{TES},s}^{\mathrm{ch}}$ and ${\eta }_{\mathrm{TES},s}^{\mathrm{d}}$ are the charging efficiency coefficient and discharge efficiency coefficient of the heat storage system, respectively.

(6) Boundary constraints of the energy storage system:

$$E_{s,\min }\le E_s\le E_{s,\max }$$

(28)

where $E_s$ is the state of charge of the heat storage system in the CSP power plant, and $E_{s,\min },E_{s,\max }$ are the lower and upper limit of the SOC of the heat storage system in the CSP power plant.

(7) Heat storage capacity constraint:

$$E_{s,t}=E_{s,T}$$

(29)

where $E_{s,t}$ is the heat storage capacity of CSP power station in the initial period, and $E_{s,T}$ is the heat storage capacity of the CSP power station in the terminal period.

4.4. Operation Constraints

(1) Upper and lower limit constraints on unit output:

$$y_{g,t}{G}_{g,\min }\le p_{g,t}\le y_{g,t}{G}_{g,\max }$$

(30)

where $G_{g,\max }$ and $G_{g,\min }$ are the upper and lower limits of the output of generator set G, respectively.

(2) Transmission capacity constraints of transmission lines:

$$z_l{f}_{i,j,\min }^{}\le f_{i,j,t}^{}\le z_l{f}_{i,j,\max }^{}$$

(31)

where $f_{i,j,\max }$ and $f_{i,j,\min }$ are the upper and lower limits of the transmission power of the line $(i,j)$, respectively, and $f_{i,j,t}$ is the power transmitted by line $(i,j)$ in time period t; $\forall (i,j)\in L,t\in T$.

(3) Node phase angle constraint:

$${\varphi }_{i,\min }\le {\varphi }_{i,t}\le {\varphi }_{i,\max }$$

(32)

where ${\varphi }_{i,\max }$ and ${\varphi }_{i,\min }$ are the upper and lower limits of the phase angle of node i, respectively, and ${\varphi }_{i,t}$ is the phase angle of node i in time period t.

(4) Line power balance constraint:

$$\frac{{\varphi }_{i,t}-{\varphi }_{j,t}}{x_{i,j}}-f_{i,j,t}+(1-z_l)B_1\ge 0$$

(33)

$$\frac{{\varphi }_{i,t}-{\varphi }_{j,t}}{x_{i,j}}-f_{i,j,t}-(1-z_l)B_2\le 0$$

(34)

where $x_{i,j}$ is the reactance of the line $(i,j)$, and B1 and B2 are large enough numbers so that Equations (33) and (34) only work for lines that have not failed.

(5) Unit climbing rate constraint:

$$p_{g,t}-p_{g,(t-1)}\le R_{g,\mathrm{up}}·y_{g,(t-1)}+U_{g,\mathrm{up}}·(y_{g,t}-y_{g,(t-1)})+G_{g,\max }·(1-y_{g,t})$$

(35)

$$p_{g,(t-1)}-p_{g,t}\le R_{g,\mathrm{dn}}·y_{g,t}+D_{g,\mathrm{dn}}·(y_{g,(t-1)}-y_{g,t})+G_{g,\max }·(1-y_{g,(t-1)})$$

(36)

where $R_{g,\mathrm{up}}$ and $R_{g,\mathrm{dn}}$ are the upward and downward climbing rates of unit g, respectively (in normal operation) and $U_{g,\mathrm{up}}$ and $D_{g,\mathrm{dn}}$ are the climbing rates of unit G at startup and shutdown, respectively.

(6) Node load balance constraint:

$$-d_{i,t}\le {\displaystyle \sum _{\forall j\in L(·,i)}{f}_{j,i,t}}-{\displaystyle \sum _{\forall j\in L(i,·)}{f}_{i,j,t}}+{\displaystyle \sum _{g\in G}{p}_{g,t}}-D_{i,t}\le d_{i,t}$$

(37)

where $D_{i,t}$ is the load of node i in time period t.

5. Model Transformation and Solution

Since the constructed model has a min–max–min structure [25], which cannot be solved directly, the subproblems are transformed here in pairs [24].

The uncertain set is a discrete variable, and the C&CG algorithm is chosen in this paper [26]. This kind of problem can be transformed into a main problem and a sub-problem, and then solved iteratively [27].

5.1. Main Problem

The first stage corresponding to the main problem is the start–stop decision of the unit before the line fault, which is expressed as:

$$\underset{u_{g,t},v_{g,t}}{\min }{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{g=1}^G(C_{g,\mathrm{U}}^{}{u}_{g,t}+C_{g,\mathrm{D}}^{}{v}_{g,t})+\widehat{W}}}$$

(38)

$s.t.$ Constraints (18) to (19)

$$\widehat{W}\ge \underset{z\in Z}{\max }Q$$

(39)

where $\widehat{W}$ is the auxiliary variable needed for solving, which is used to represent the objective function of the second stage.

5.2. Subproblems

The subproblem is an economic scheduling problem after the occurrence of the worst line fault. The unit output is adjusted to reduce the generation cost and the loss of load penalty cost [28], which can be expressed as:

$$\underset{z\in Z}{\max }Q=\underset{z\in Z}{\max }\min {\displaystyle \sum _{t=1}^T({\displaystyle \sum _{g=1}^G{\lambda }_g{p}_{g,t}}}+{\displaystyle \sum _{i=1}^I{C}_{\mathrm{VOLL}}{d}_{i,t}})$$

(40)

$s.t.$ Constraints (7)~(8), CSP constraints, (21)~(29)

Since the subproblem is constructed in the form of max–min, it is difficult to solve it directly. In this paper, according to the duality theory, the inner min problem is transformed, and the corresponding max problem is obtained. At the same time, the outer max problem is combined for calculation. Since the objective function of the subproblem is nonlinear, in order to change the subproblem of the MILP problem [29] to make the model easier to solve, the Big-M [30] method can be used to express the results linearly.

5.3. C&CG Algorithm

General form of two-stage robust adaptive optimization:

$$\begin{array}{l}\underset{\mathit{y}}{\min }{\mathit{c}}^T\mathit{y}+\underset{\mathit{u}\in U}{max}\underset{\mathit{x}\in \mathit{F}(\mathit{y},\mathit{u})}{\min }{\mathit{b}}^{T}x\\ \mathrm{s}.\mathrm{t}.{\mathit{A}}_1{}^T\mathit{y}\ge \mathit{d}\\ \mathit{F}(\mathit{y},u)=\{{\mathit{A}}_2\mathit{x}:\mathit{x}\ge \mathit{h}-{\mathit{A}}_3\mathit{y}-{\mathit{A}}_4\mathit{u}\}\end{array}$$

(41)

where y and x are the matrices composed of decision variables in the first and second stages of the C&CG algorithm, respectively, and u is an uncertain set, which can be a discrete set or a polyhedron. The optimization constraint F(y, u) in the second stage is a linear function of u. The matrix ${\mathit{A}}_1,\mathit{b},\mathit{c},\mathit{d}$ is the constant coefficient matrix corresponding to the model in the first stage. The matrix ${\mathit{A}}_2,{\mathit{A}}_3,{\mathit{A}}_4,\mathit{h}$ is the constant coefficient matrix corresponding to the model of the second stage.

After transforming both the master problem and the subproblem into a mixed integer programming problem, the C&CG algorithm can be used to solve the problem, and the corresponding solver can be used to obtain the results iteratively [28].

The C&CG steps are as follows:

Step 1: Data initialization: set LB to infinite lower bound, UB to infinite upper bound, iteration number j = 0, and set O to empty set.

Step 2: Solve the MP master problem.

MP:

$$\begin{array}{l}\underset{\mathit{y}}{\min }{\mathit{c}}^T\mathit{y}+\mathit{\eta }\\ \mathrm{s}.\mathrm{t}.{\mathit{A}}_1{}^T\mathit{y}\ge \mathit{d}\\ \mathit{\eta }\ge {\mathit{b}}^T{\mathit{x}}_j,\forall j\in \mathit{O}\\ {\mathit{A}}_2{\mathit{x}}_j\ge \mathit{h}-{\mathit{A}}_3\mathit{y}-{\mathit{A}}_4{u}_j^{*}\end{array}$$

(42)

Obtain the optimal solution $({\mathit{y}}_{k+1}^{*},{\mathit{\eta }}_{k+1}^{*},{\mathit{x}}_1^{*},\dots ,{\mathit{x}}_k^{*})$ and update the lower bound LB = ${\mathit{c}}^T{\mathit{y}}_{k+1}^{*}+{\mathit{\eta }}_{k+1}^{*}$;

Step 3: Solve the SP subproblem.

Plug in $\mathit{y}={\mathit{y}}_{k+1}^{*}$

SP:

$$\mathit{Q}({\mathit{y}}_{k+1}^{*})=\underset{\mathit{u}\in U}{max}\underset{\mathit{x}\in \mathit{F}(\mathit{y},\mathit{u})}{\min }\{{\mathit{b}}^T\mathit{x}:{\mathit{A}}_2\mathit{x}\ge \mathit{h}-{\mathit{A}}_3{\mathit{y}}_{k+1}^{*}-{\mathit{A}}_{4}u\}$$

(43)

The optimal solution $(u_{k+1}^{*},{\mathit{x}}_{k+1}^{*})$ is obtained, and the upper bound of the original problem is obtained, UB = ${\mathit{c}}^T{\mathit{y}}_{k+1}^{*}+{\mathit{b}}^T{\mathit{x}}_{k+1}^{*}$;

Step 4: Determine whether the convergence condition is satisfied.

If $UB-LB\le \epsilon$, return and the program terminates $({\mathit{y}}_{k+1}^{*},{\mathit{x}}_{k+1}^{*})$;

Otherwise:

(a) If $\mathit{Q}({\mathit{y}}_{k+1}^{*})<+\infty$, add variables ${\mathit{x}}_{k+1}^{*}$ and add the following constraints

$$\mathit{\eta }\ge {\mathit{b}}^T{\mathit{x}}_{k+1}$$

(44)

$${\mathit{A}}_2{\mathit{x}}_{k+j}\ge \mathit{h}-{\mathit{A}}_3\mathit{y}-A_4{u}_{k+1}^{*}$$

(45)

Return to the main problem (42). The solution to problem (43) can be found by searching the database. $j=j+1$ then, $\mathbf{O}=\mathbf{O}\cup j+1$, go to Step 2.

(b) $\mathit{Q}({\mathit{y}}_{k+1}^{*})=+\infty$ (For some ${\mathit{u}}^{*}\in U$, if the second stage decision $\mathit{Q}({\mathit{y}}_{k+1}^{*})$ is not feasible, $\mathit{Q}({\mathit{y}}_{k+1}^{*})$ denoted as $+\infty$, add variables ${\mathit{x}}_{k+1}$, and add the following constraints:

$${\mathit{A}}_2{\mathit{x}}_{k+j}\ge \mathit{h}-{\mathit{A}}_3\mathit{y}-{\mathit{A}}_4{u}_{k+1}^{*}$$

(46)

Return to the main problem (42). $j=j+1$, go to Step 2.

6. Case Studies

To verify the effectiveness of the proposed method, the modified IEEE test system [31] and an actual 78 nodes distribution system in China are simulated. The simulation is conducted with GAMS and the problem is solved by CPLEX. The simulation environment is Intel Core i7-8700k, 3.7-GHz CPU with 8 GB of RAM.

6.1. Analysis of Computational Results for IEEE 14-Node System

The IEEE 14-node test system involves five conventional thermal power generators and 20 lines, as shown in Figure 5. The time interval considered by the model in this chapter is 1 h, and the duration is 24 h. The cost of load loss penalty is 3000 USD/(MW·h), which is only used to illustrate the effectiveness of the proposed model because the IEEE 14-node system is a small system. The CSP power station replaces one thermal power unit of the system, and the specific parameters of the unit and line are shown in

Table 1. Cost parameters of units.

Unit	Node	Start-Up Costs/M USD	Marginal Generation Cost/(USD/MWh)
1	1	0.14286	84.26
2	2	0.04286	80.00
3	3	0.14286	85.71
4	6	0.02571	79.45
5	8	0.14286	82.56

and

Table 2. Power generation parameters of units.

Unit	1	2	3	4	5
Node	1	2	3	6	8
Upper power limit/MW	232	114	220	210	100
Lower power limit/MW	30	30	12.5	12.5	12.5
Creep speed/(MW·h−1)	15	15	15	15	15
Minimum down time/h	12	12	12	12	12
Minimum startup time/h	12	12	12	12	12

. The predicted output of wind power PV and load are shown in Figure 6 and Figure 7.

Set the following four solutions for testing:

Scheme 1: Do not consider the CSP power station and N-k safety criteria.

Scheme 2: Consider the N-k safety criterion instead of the CSP power station.

Scheme 3: Consider the CSP power station, without considering the N-k safety criterion.

Scheme 4: Consider the CSP power station and the N-k safety criteria.

(1) Analysis and comparison of optimization operation results.

Table 3. Comparison of optimized operation results of four cases.

Index	Plan	Total System Operation Cost/USD
CSP plants are not considered	Plan 1	8402.78
	Plan 2	8415.58
CSP plants are considered	Plan 3	7629.89
	Plan 4	7674.61

shows the decision results of the model under different schemes. Table 3 shows that considering the CSP plant’s plan 3 and 4 compared with and without considering CSP plant’s plan 1 and 2, the total system operating costs have dropped substantially, which proved that the CSP plant can rise to adjust the system and storage and can replace the poor economic performance of the thermal power unit, which reduces the traditional operation cost of the unit. Although the operation cost of scheme 1 and scheme 3 without the N-k safety criterion is lower than that of scheme 2 and scheme 4 with the N-k safety criterion, it does not take into account the expected fault. Once an accident occurs in the power system, it may lead to the phenomenon of the power flow exceeding the limit, thus affecting the normal working of the whole system. The reason for the increase in the total operating cost after considering the N-k safety criterion is that the objective function of the model considers the penalty cost of the loss of load, and the phenomenon of load cutting occurs after the failure, which makes the cost higher. Therefore, the combination of the CSP power station and the N-k safety criterion can complement each other, so as to ensure that the power system can operate safely and further reduce the operating cost of the system.

(2) Replace the conventional thermal power unit of node 6 with a group of CSP power stations and observe the disconnection of system decisions before and after adding CSP power stations. The results are shown in

Table 4. Unit disconnection before and after adding concentrating solar power.

Plan	The CSP Power Station Is Not Added	The CSP Power Station Is Added
k = 1	L7	L14
k = 2	L7, L10	L14, L15
k = 3	L7, L10, L13	L8, L14, L15
k = 4	L7, L10, L11, L13	L1, L8, L14, L15

below.

It can be seen from Table 4 that the disconnection conditions made by the system before adding the CSP power station are lines L7, L10, L11, and L13, while the disconnection conditions made by the system after introducing the CSP power station are lines L8, L12, L14, and L15, which are quite different from each other. This is because lines L10 and L13 are connected to node 6. When the CSP power station is at node 6, the nearby lines are not easy to break, which indicates that the CSP power station plays an important role in the power system and can change the power flow distribution of the transmission line and realize the overall optimization of the system.

(3) In the IEEE14 node system, two conventional generating units are randomly selected to compare the output status of conventional generating units before and after adding the CSP power station. The result is shown in Figure 8.

The results show that after the introduction of the CSP power station, the total output of conventional units has changed and its output power has slowed down, which indicates that when the power grid load is at a low ebb, the CSP power station uses its own energy storage function to store the remaining power in the grid. When the power system load was in the peak state, in order to compensate the shortage of power supply, the CSP power station started to discharge. The regulating function of the CSP power station can relieve the peak regulating pressure of conventional units to a certain extent and reduce the pressure of the scheduling and distribution of thermal power units.

(4) In the IEEE14 node, the wind farm wind abandonment and line currents are compared before and after the addition of the CSP plant, and the results of the runs are shown in Figure 9, Figure 10 and Figure 11.

As can be seen from the operational results, a certain degree of wind abandonment is generated due to high wind speeds and low loads at night. After joining the CSP power station, the amount of wind abandonment in the power system is reduced, especially at night, by nearly 50% compared to when it was not connected to the CSP power station. This is because the heat storage system of the CSP power station can store and discharge heat at different times, which, together with wind power generation, can achieve peak shaving and valley filling, increase the new grid-connected capacity, and improve the optimal scheduling capacity of the power grid. In addition, it is clear from the changes in line currents before and after the addition of the CSP power station that the CSP power station can improve the active tide distribution, significantly smoothing the power generation output and maximizing the system grid. In summary, CSP power stations have a role in the consumption of new energy, i.e., they have a better optimizing effect on the economy of the power system.

6.2. Analysis of Calculation Results for the IEEE 118 Node System

To test the applicability of the proposed model in a larger scale grid, the IEEE 118 node system was tested and analyzed, which involved 54 conventional generating units, 186 lines, and 91 load nodes, and the structure of the system is shown in Figure 12. CSP power stations were substituted for the system’s 15th, 16th, 36th, 46th, 56th, and 66th power stations. The operating results of the different scenarios and the distribution of transmission line currents before and after the addition of the CSP power station are shown in

Table 5. Comparison of optimized operation results of four cases.

Index	Plan	Total System Operation Cost/USD
CSP plants are not considered	Plan 1	77,315.30
	Plan 2	78,025.14
CSP plants are considered	Plan 3	56,035.04
	Plan 4	56,749.17

, Figure 13 and Figure 14 below.

As can be seen from Table 5, the total system operating costs of option 3 and option 4 considering CSP plants are lower in the IEEE 118 node system tested, with a total cost reduction of around 20,000 USD, compared to option 1 and option 2 without CSP plants. This demonstrates the ability of CSP plants to regulate and store heat in this system, which helps to reduce the start-up and shutdown and the output of less economical conventional units, and it also demonstrates its suitability for larger systems. The total operating costs of the system are slightly higher when the N-k safety criterion is taken into account because of the occurrence of load shedding in the event of an N-k scenario, but this provides a safeguard for the safety of the system operation. In addition, as can be seen from Figure 13 and Figure 14, the CSP plant is able to improve the tidal distribution of the system active power, and this device, when used in conjunction with the thermal storage device, can significantly smooth out the power output, reduce fluctuations in power output, and improve the stability of the system. Therefore, the N-k safety criterion enables the system to operate safely in the event of a fault. On this basis, the introduction of CSP power plants can effectively improve the flexibility of the power system and make the system more economical under the condition of optimizing the grid structure.

7. Conclusions

This paper constructs an ambiguity set of power system line faults by using the existing practical sample data. Considering the flexibility provided by CSP, a distributionally robust economic scheduling model is proposed under uncertain conditions considering the flexibility of CSP and the N-k safety criteria. The flexible regulation ability of CSP can be used to support new energy consumption, which can reduce the carbon emissions of the power system and effectively improve the economy of system operation. Based on the IDM method, a distributionally robust economic scheduling model considering the flexibility of CSP and the N-k safety criteria is constructed, which can integrate the robust optimization and stochastic optimization and reduce the conservatism of the robust optimization method.

In subsequent studies, the uncertainty of new energy generation such as through transmission lines and wind power can be comprehensively considered, so that the proposed method can be extended and applied to the unit combination and optimal scheduling problems of power systems with a high proportion of new energy generation.