1 Introduction

Hierarchical architectures are used to solve problems with different time scales in power-frequency and var-voltage control systems. In power-frequency control systems, unit commitment (UC) and economic dispatch (ED) are used to determine the generation schedules over the next day. The schedules are also updated approximately every 15 min according to a tertiary frequency control (TFC) algorithm that solves real-time ED problems for ultra-short-term load forecasts. The secondary frequency control (SFC) and primary frequency control (PFC) methods are, respectively, responsible for making and executing commands over shorter time scales. Var-voltage control is similar to power-frequency control in most situations. For instance, tertiary voltage control (TVC) sets the reference pilot bus voltage curves with the same frequency as TFC. However, as multi-timescale coupling is considered insignificant in var-voltage control, longer time scales, such as those relevant to scheduling an entire day, have been considered in only a few studies. This rarely causes problems because, in general, these systems have an abundance of reactive power, and the state of the system is not dominated by time-coupled behaviors.

In China, the stochastic fluctuations in wind power generated by sources that are centrally integrated into weak networks cause significant interference with the voltage stability of the system. The lack of dynamic reactive-power support from conventional power plants means that the available methods for regulating large voltage perturbations are limited. Meanwhile, classical substation control methods use a coordinated secondary voltage control (CSVC) system to determine when to operate the discrete devices so that the voltage remains within acceptable bounds [1,2,3]. CSVC strategies aim to maintain the real-time pilot bus voltage while complying with the guideline to prioritize the use of discrete devices for voltage regulation. As the CSVC optimization algorithm runs every few minutes, the accuracy of the ultra-short-term forecasts is sufficient for ensuring that the control decisions can adapt to the current state of the system. Future changes in wind-power generation are not taken into account by the classical control algorithm. However, each device can only run for limited times each day and should maintain its state for some time after it is operated. If the choice to dispatch these devices is only based on past and real-time information, they may not be operational at crucial moments. Moreover, this undesirable feature of discrete devices could mean that all of the remaining dynamic reactive power generated by the wind farms is required to return the system to its required state.

Consequently, we considered methods to increase the dynamic reactive power capacity of the system, such as deploying additional static synchronous compensators (STATCOMs). There have already been studies assessing the STATCOM capacity required to avoid cascading after (N−1) faults [4,5,6]. However, there have been few studies on improving the operation schedules of discrete devices, such as shunt capacitors/reactors and on-load tap changers (OLTCs), to maintain a more stable dynamic reactive-power reserve. In fact, the current behavior of discrete devices is far from optimal, so it is not surprising that voltage violations occur and that there is a lack of dynamic reactive-power reserves. The difficulties that arise due to uncertainties in wind-power generation and the limited availability of discrete devices must be overcome before we can properly solve the problem.

In the case of power-frequency control systems, uncertainties are taken into account by UC and ED problems. Four methods for managing uncertainties in system operation are discussed in [7]. Many studies have focused on the fourth method, robust optimization, because it uses a modest amount of knowledge and is universally applicable. In particular, researchers use a two-stage version of this method, where the decision variables are divided into two types and are input into different stages of the optimization algorithm. In the first stage, only the on/off states of units that cannot be corrected frequently are considered. This ensures that when another decision variable modifies the state of the system, the system maintains the best performance, even under the worst conditions [8,9,10,11].

Based the above-mentioned approaches, we propose a two-level hierarchical discrete-device control (THDDC) method for solving the var-voltage problem in areas that are integrated with wind farms. At the day-ahead level, a two-stage optimal robust discrete-device commitment (TORDDC) method is applied to decide the permitted time intervals for device operation. These are also the decision variables for the first stage of the algorithm. In the second stage, we select the exact times to activate the discrete devices in response to various wind scenarios. The within-day stage uses a model-predictive control-based discrete-device dispatch (MPCDDD) method with the latest forecasts to decide which specific moments require action. This approach is more effective in that the dispatch discrete devices can be applied in a variety of wind scenarios.

The main contributions of this paper are summarized below:

  1. 1)

    The robust form of the TORDDC model ensures that most wind scenarios are considered. The cost of safely maintaining the voltage can be calculated from the forecast for the day ahead.

  2. 2)

    The correlations over the course of the entire day are taken into account by the TORDDC model. We can generate reasonable schedules for activating the discrete devices, and sufficient dynamic reactive power can be reserved for critical moments.

  3. 3)

    The TORDDC method only selects permitted operational time intervals. As the within-day control algorithm has some freedom to choose when to activate the discrete devices, it is adaptable to changeable wind environments.

  4. 4)

    We significantly reduced the time window for the MPCDDD problem, and thus the computation time, by identifying potential operational periods for the day ahead.

The remainder of this paper is organized as follows. We outline the proposed method and its role in the var-voltage control system in Section 2. We describe how to model the TORDDC and MPCDDD problems in Section 3. In Section 4, we report the results of the proposed method when tested on a model of a typical system integrated with wind-power sources located in China. Finally, we draw our conclusions in Section 5.

2 Architecture

The classical three-level voltage control architecture has been widely used in large-scale power grids in the decades since it was first proposed by EDF in 1987 [12]. The top level uses TVC to provide reference pilot bus voltages in the case of least power loss, based on optimal power flow throughout the entire grid. This step is repeated every 15 min using inputs in real time or a close single-timescale prediction of the system state to the optimization algorithm. The secondary voltage control (SVC) method, which mitigates the effects of higher frequency perturbations or inaccurate forecasts, is run every few minutes. This acts as an auxiliary to the TVC method and regulates the pilot bus voltage when it deviates from the reference value. The reference values for the control devices are updated periodically. The partition control of the SVC method reflects the regional characteristics of var-voltage control. At the bottom level, the primary voltage-control (PVC) method orchestrates the local systems that perform var-voltage regulation for each device. The load demands in transmission networks that are integrated with few renewable energy sources are supported by conventional power plants. The future loads in these networks vary regularly, as opposed to fluctuating significantly. In this case, the reactive power output is sufficient to guarantee the system voltage safety without frequently operating any discrete device, e.g., switching the shunt capacitors/reactors. Considering the negligible time constant of the generator-excitation system in comparison to the TVC/SVC control period, single-timescale optimization is sufficient in this case.

As more renewable energy sources have been centrally integrated, the system voltage stability has become vulnerable to strong fluctuations in wind power, necessitating the curtailment of wind-power generation. In 2017, five provinces had wind curtailment rates of over 10%. Some of the main reasons for this are the limitations on the availability of the discrete devices and lack of reactive power reserves during critical moments. Consequently, developing reasonable discrete-device operation schedules and making better use of the continuous reactive power output from wind farms within a narrow adjustable range appear to be promising alternatives to deploying more STATCOMs, which is expensive.

The architecture of the proposed THDDC method is shown in Fig. 1. As the time couplings between the discrete devices last all day, they are naturally similar to the UC and ED problems in power-frequency control systems. Hence, we modify the classical var-voltage control system by adding the TORDDC forecast for the day ahead. Two-stage robust optimization is helpful when forecasts are uncertain. Initially, we attempted to compile daily schedules for operating the discrete devices, and then apply additional continuous control methods to compensate for the deviations caused by the inaccuracies in the daily forecasts. However, this approach is not valid in networks containing wind power sources, because they lack reactive power support from conventional power plants, i.e., unnecessary wind curtailments are imposed due to overly conservative daily schedules or the feasible domain does not exist. The MPCDDD procedure is inspired by the model-predictive control (MPC) theory and uses more precise forecasts to decide exactly when to operate each discrete device during the day. We could improve the control performance using information updated in real time and a rolling optimization procedure; however, these approaches are costly and take much longer to evaluate individually in within-day control scenarios. Hence, it is not practical to optimize a day-long time window.

Fig. 1
figure 1

Architecture of THDDC method

The speed of two-stage mixed-integer optimization depends strongly on the scales of the integers. To balance the feasibility of the dispatch schedule and the effectiveness of the second-stage algorithm, we must obtain proper restrictions on the allowed discrete-device operation periods for the entire day ahead. Therefore, the key is to identify a group of day-ahead decision variables that represent the correlations between the discrete-device operations at different moments but can decouple the correlations to fit shorter time scales. Ultimately, we chose to use the permitted discrete-device operating time intervals as the day-ahead decision variables and compromise by deciding precisely when to operate them during the day. Figure 2 shows the rolling operation of this method.

Fig. 2
figure 2

Rolling operation of MPCDDD method

Thus, the MPCDDD procedure synchronously tracks which discrete devices are not permitted to operate during specific SVC periods. Control commands are issued from the control center.

3 Methodology

3.1 Variables for modelling discrete devices

Most types of devices are controlled in terms of variables that have intuitive physical meanings. In the case of discrete devices, they refer to the on/off state of a shunt capacitor/reactor, the tap position of an OLTC, etc. We use a unified integer variable \(n_{i}^{\tau }\) to represent a discrete device \(i\) during period \(\tau\). Each specific discrete device has a time-invariant maximum/minimum \(\bar{n}_{i} /\underline{n}_{i}\).

$$\underline{n}_{i} \le n_{i}^{\tau } \le \bar{n}_{i}$$
(1)

The values of operational variables are relative to the differences between their corresponding state variables. In this paper, as the operating states of the devices have no direction or amplitude, we use a binary variable \(b_{i}^{\tau }\) to represent whether discrete device \(i\) is permitted to operate during period \(\tau\). The device is only operated if the value of \(b_{i}^{\tau }\) is 1, subject to the following constraint:

$$- b_{i}^{\tau } \left( {\bar{n}_{i} - \underline{n}_{i} } \right) \le n_{i}^{\tau } - n_{i}^{\tau - 1} \le b_{i}^{\tau } \left( {\bar{n}_{i} - \underline{n}_{i} } \right)$$
(2)

As mentioned in Section 2, the day-ahead decision variables are the permitted operating-time intervals, i.e., a discrete device cannot be operated more than once during each interval. We represent the times of the earliest and latest potential operations in different time intervals by the pair of binary variables \(d_{i}^{\tau }\) and \(e_{i}^{\tau }\), respectively. However, we cannot find an intuitive relationship between these and \(b_{i}^{\tau }\).

We elucidate this problem by analyzing a simple case. Suppose that different load conditions cause the voltage at a substation to vary along the different curves in Fig. 3a, when the influence of shunt capacitors/reactors is neglected. In response to these conditions, Fig. 3b shows that we could switch on a correlated shunt capacitor once the voltage decreases below 0.95 p.u..

Fig. 3
figure 3

Substation voltage and capacitor state

We obtained the operating-time curves shown in Fig. 4a by summing \(b_{i}^{\tau }\). Although the state of the device may not always change between t1 and t2, we still increase the device operation count each time the substation voltage reaches a new steady state. The upper and lower envelopes of these curves, which are the exact summations of \(d_{i}^{\tau }\) and \(e_{i}^{\tau }\), are shown in Fig. 4b. These restrict the operation of the discrete devices to the intervals with nonzero gaps.

Fig. 4
figure 4

Capacitor operation times and its upper and lower envelopes

From the above, we can express the relationship between \(d_{i}^{\tau } /e_{i}^{\tau }\) and \(b_{i}^{\tau }\) as follows:

$$\sum\limits_{\upsilon = 1}^{\tau } {e_{i}^{\upsilon } } \le \sum\limits_{\upsilon = 1}^{\tau } {b_{i}^{\upsilon } } \le \sum\limits_{\upsilon = 1}^{\tau } {d_{i}^{\upsilon } }$$
(3)

We restrict the height of the shaded rectangle to 1 to ensure that none of the discrete devices are operated multiple times during a specific interval.

$$\sum\limits_{\upsilon = 1}^{\tau } {\left( {d_{i}^{\upsilon } - e_{i}^{\upsilon } } \right)} \le 1$$
(4)

As discussed in Section 2, the time couplings and limited availability must be defined exclusively in terms of variables used in the day-ahead forecasting procedure. Hence, they must be represented by day-ahead decision variables, as follows:

$$\sum\limits_{{\tau \in {\mathcal{T}}_{24} }} {d_{i}^{\tau } } \le \mathcal{A}_{i}$$
(5)

where \({\mathcal{T}}_{24}\) is the set of indices of time periods covering the entire day and \(\mathcal{A}_{i}\) is the maximum number of times we can operate discrete device \(i\).

We prevent too much freedom in the choice of when to operate the devices from adversely affecting the optimization computation time by limiting the daily sum of the shaded area in Fig. 4b for discrete device \(i\) to \(\mathcal{B}_{i}\).

$$\sum\limits_{{\tau \in {\mathcal{T}}_{24} }} {\sum\limits_{\upsilon = 1}^{\tau } {\left( {d_{i}^{\upsilon } - e_{i}^{\upsilon } } \right)} } \le \mathcal{B}_{i}$$
(6)

3.2 Variables representing continuous-control devices

The routine method for controlling the voltage in networks with wind farms uses the reactive power output from wind farms and STATCOMs. In addition, wind-power generation is curtailed when these methods cannot maintain the system voltage within acceptable bounds. We use a unified continuous variable \(q_{j}^{\tau }\) to represent the reactive power output from continuous-control device \(j\) during period \(\tau\). The values of these variables are relative to the active power \(p_{j}^{\tau }\), given that the reactive power generated by a wind farm consists of the power from wind-turbine generators (WTGs) and other compensators. The operation of WTGs is restricted by a minimum power factor. We define parameter \(\lambda_{j}\) as the ratio of this component of the reactive power to the active power under this power factor. Meanwhile, \(\bar{q}_{j} /\underline{q}_{j}\) are, respectively, the maximum/minimum reactive power outputs from device \(i\) when the active power is 0. In this case, we obtain the constraints defined below:

$$\underline{q}_{j} - \lambda_{j} p_{j}^{\tau } \le q_{j}^{\tau } \le \bar{q}_{j} + \lambda_{j} p_{j}^{\tau }$$
(7)

In var-voltage control, which is concerned with the reservation of dynamic reactive power, i.e., maximizing the distance between the current reactive power and the closest bound, this problem is more complicated because it requires a “max-min” model. By considering the opposite problem, we define variable \(s_{i}^{\tau }\) to be the magnitude of the reactive power deviation from its mid-value. The problem can be simplified significantly to the following:

$$\left\{ {\begin{array}{*{20}l} { - s_{j}^{\tau } \le q_{j}^{\tau } - \frac{{\bar{q}_{j} + \underline{q}_{j} }}{2} \le s_{j}^{\tau } } \hfill \\ {s_{j}^{\tau } \ge 0} \hfill \\ \end{array} } \right.$$
(8)

As the variable representing the original active power of device \(j\) at period \(\tau\) is \(r_{j}^{\tau }\) with supremum/infimum \(\bar{r}_{j}^{\tau } /\underline{r}_{j}^{\tau }\), respectively, we have

$$\left\{ {\begin{array}{*{20}l} {\underline{r}_{j}^{\tau } \le r_{j}^{\tau } \le \bar{r}_{j}^{\tau } } \hfill \\ {0 \le p_{j}^{\tau } \le r_{j}^{\tau } } \hfill \\ \end{array} } \right.$$
(9)

3.3 Observed values of variables for different devices

We can assess the safety system by measuring the variables representing other devices, \(u_{k}^{\tau }\), in addition to the control devices. The voltage of bus \(k\) in period \(\tau\) should remain within the respective upper/lower bounds \(\bar{u}_{k} /\underline{u}_{k}\). Sometimes, none of the available means maintain the value of the voltage within the acceptable range, thus causing an excess voltage \(v_{k}^{\tau }\). The system is at risk when \(v_{k}^{\tau }\) is not equal to 0. Excess voltage over the maximum \(\bar{v}_{k}\) is strictly prohibited by the following constraint:

$$\left\{ {\begin{array}{*{20}l} {\underline{u}_{k} - v_{k}^{\tau } \le u_{k}^{\tau } \le \bar{u}_{k} + v_{k}^{\tau } } \hfill \\ {0 \le v_{k}^{\tau } \le \bar{v}_{k} } \hfill \\ \end{array} } \right.$$
(10)

The relationship between the bus voltage and the other control variables can be expressed as a linearized power flow sensitivity model:

$$\begin{aligned} u_{k}^{\tau } - \hat{u}_{k}^{\tau } = & \sum\limits_{{i \in \mathcal{I}}} {\xi_{k,i} \left( {n_{i}^{\tau } - \hat{n}_{i}^{\tau } } \right)} \\ & + \sum\limits_{{j \in \mathcal{J}}} {\psi_{k,j} \left( {p_{j}^{\tau } - \hat{p}_{j}^{\tau } } \right)} + \sum\limits_{{j \in \mathcal{J}}} {\zeta_{k,j} \left( {q_{j}^{\tau } - \hat{q}_{j}^{\tau } } \right)} \\ \end{aligned}$$
(11)

We use the 24-hour power flow forecast to calculate the values predicted by this model. \(\hat{u}_{k}^{\tau }\), \(\hat{n}_{i}^{\tau }\), \(\hat{p}_{j}^{\tau }\), and \(\hat{q}_{j}^{\tau }\) are, respectively, the base-state values of \(u_{k}^{\tau }\), \(n_{i}^{\tau }\), \(p_{j}^{\tau }\), and \(q_{j}^{\tau }\). \(\xi_{k,i}\), \(\psi_{k,j}\), and \(\zeta_{k,j}\) are the sensitivity of bus \(k\) to the control state of discrete device \(i\), the active power output of continuous device \(j\), and the reactive power output of continuous device \(j\), respectively. \(\mathcal{I}\) and \(\mathcal{J}\) are the sets of indices of the discrete- and continuous-control devices.

Decisions to curtail wind-power generation are not based solely on voltage safety; rather, a range of factors are taken into consideration. For example, we must prevent the branch power flow from reaching the transmission limit. Here, \(o_{l}^{\tau }\) is the active power flow of branch \(l\) at period \(\tau\), and \(\bar{o}_{l}\) is the transmission limit.

$$- \bar{o}_{l} \le o_{l}^{\tau } \le \bar{o}_{l}$$
(12)

Since the transmission systems of networks containing wind farms use DC power flow, the branch active power has the following linear representation:

$$o_{l}^{\tau } = \sum\limits_{{j \in \mathcal{J}}} {\kappa_{l,j} p_{j}^{\tau } }$$
(13)

where \(\kappa_{l,j}\) is the factor that shifts continuous-control device \(j\) to branch \(l\).

3.4 Formulation of TORDDC model

The day-ahead discrete-device commitment problem is represented by two separate sets of decision variables, which we determine in two stages. We decide which time intervals the devices can operate in, represented by \(d_{i}^{\tau } /e_{i}^{\tau }\), in the first stage. The vector \(\varvec{y}\) denotes all of the first-stage variables, including \(d_{i}^{\tau }\) and \(e_{i}^{\tau }\), and vector \(\varvec{r}\) denotes the uncertainty, \(r_{j}^{\tau }\). The original problem can then be expressed as (14), which has constraints (4)-(6) and (9) with \(\forall r_j^\tau \in [{\underline r}_j^\tau,{\bar r}_j^\tau], \space i \in \mathcal{I}\), \(\tau \in {\mathcal{T}}_{24}\), \(j \in \mathcal{J}\).

$$\mathop {\hbox{min} }\limits_{\varvec{y}} \sum\limits_{{\tau \in {\mathcal{T}}_{24} }} {\sum\limits_{{i \in {\mathcal{I}}}} {\alpha_{i}^{\tau } d_{i}^{\tau } } } + H\left( {\varvec{y},\varvec{r}} \right)$$
(14)

In the above, the objective function includes a term for the cost of activating each discrete device, where \(\alpha_{i}^{\tau }\) is the cost of activating discrete device \(i\) during period \(\tau\). The function \(H\left( {\varvec{y},\varvec{r}} \right)\) is the optimal value of the second-stage problem. We apply a robust method, and transform the problem into a “min-max” problem, as shown in (15), which has the same constraints as (14).

$$\mathop {\hbox{min} }\limits_{\varvec{y}} \left( {\sum\limits_{{\tau \in {\mathcal{T}}_{24} }} {\sum\limits_{{i \in {\mathcal{I}}}} {\alpha_{i}^{\tau } d_{i}^{\tau } } } + \mathop {\hbox{max} }\limits_{\varvec{r}} H\left( {\varvec{y},\varvec{r}} \right)} \right)$$
(15)

In the second stage, we aim to minimize the composite cost of excess bus voltage, active-power curtailment, and reactive-power deviation. These appear in the objective function in order of descending importance. Their weights are, respectively, \(\beta_{k}^{\tau }\), \(\gamma_{j}^{\tau }\), and \(\delta_{j}^{\tau }\). Vector \(\varvec{x}\) denotes all second-stage decision variables. We then express the second-stage problem as (16), which has constraints (1)-(3) and (7)-(13) with \(\forall \tau \in {\mathcal{T}}_{24}\), \(i \in \mathcal{I}\), \(j \in \mathcal{J}\), \(\forall k \in \mathcal{K}\), \(\forall l \in \mathcal{L}\), where the sets \(\mathcal{K}\) and \(\mathcal{L}\) are the indices of the buses and branches, respectively.

$$\begin{aligned} H\left( {\varvec{y},\varvec{r}} \right) = \mathop {\hbox{min} }\limits_{\varvec{x}} \sum\limits_{{\tau \in {\mathcal{T}}_{24} }} {\left[ {\sum\limits_{{k \in {\mathcal{K}}}} {\beta_{k}^{\tau } v_{k}^{\tau } } + \sum\limits_{{j \in {\mathcal{J}}}} {\gamma_{j}^{\tau } \left( {r_{j}^{\tau } - p_{j}^{\tau } } \right)} + \sum\limits_{{j \in {\mathcal{J}}}} {\delta_{j}^{\tau } s_{j}^{\tau } } } \right]} \\ \end{aligned}$$
(16)

An exact mixed-integer algorithm for solving the two-stage robust linear problem based on a column-and-constraint generation approach is presented in [13]. This algorithm successfully solves the TORDDC model in the required form.

3.5 Formulation of MPCDDD model

The solution of the day-ahead TORDDC problem provides the permitted discrete-device operation intervals \(d_{i}^{\tau *} /e_{i}^{\tau *}\), which are represented by the shaded area in Fig. 5a. The rolling MPCDDD procedure, which runs every 15 min, monitors system performance over the next 4-hour period and makes decisions on when to operate each device, as shown in Fig. 5b. \({\mathcal{T}}_{4}\) is the set of indices of periods that cover this 4-hour interval.

Fig. 5
figure 5

Typical results calculated using TORDDC algorithm and rolling within-day MPCDDD algorithm

The latest active power forecast is \(r_{j}^{\tau *}\). As the control systems are similar to those in the second stage of the TORDDC model, we can write the model as (17), which has constraints (1)-(3) and (7)-(13) with \(\forall \tau \in {\mathcal{T}}_{4}\), \(i \in \mathcal{I}\), \(j \in \mathcal{J}\), \(\forall k \in \mathcal{K}\), \(\forall l \in \mathcal{L}\). The main difference between these two models is that the time index \(\tau\) belongs to different sets, \({\mathcal{T}}_{24}\) and \({\mathcal{T}}_{4}\).

$$\mathop {\hbox{min} }\limits_{\varvec{x}} \sum\limits_{{\tau \in {\mathcal{T}}_{4} }} {\left[ {\sum\limits_{{k \in {\mathcal{K}}}} {\beta_{k}^{\tau } v_{k}^{\tau } } + \sum\limits_{{j \in {\mathcal{J}}}} {\gamma_{j}^{\tau } \left( {r_{j}^{\tau *} - p_{j}^{\tau } } \right)} + \sum\limits_{{j \in {\mathcal{J}}}} {\delta_{j}^{\tau } s_{j}^{\tau } } } \right]}$$
(17)

3.6 Extension of THDDC method to other transmission networks

We designed the THDDC method presented in this paper to improve the control of discrete devices in power networks that contain wind farms. The method is intended to compensate for the lack of conventional reactive power sources in these networks. Similar problems could arise in stress transmission power systems if accurate day-ahead load forecasts are difficult to obtain. Under these conditions, the system could get so close to the voltage stability limit that load shedding is required. As better discrete device control would facilitate the maintenance of the voltage stability, we could also apply the THDDC method to this system. In this case, the terms for curtailing the wind power generation should be removed from the objective functions of the second-stage of the TORDDC problem and the MPCDDD problem; terms representing the costs of factors that affect these systems should be included.

In this case, we consider the load and conventional power plants to be continuous control devices. We linearize the nonlinear cost function in a piecewise manner, and thus express the active power cost of continuous control device \(j\) during period \(\tau\) by:

$$C_{j}^{\tau } = \rho_{0,j}^{\tau } + \sum\limits_{n = 1}^{\mathcal{N}} {\rho_{n,j}^{\tau } z_{n,j}^{\tau } }$$
(18)

where \(\rho_{n,j}^{\tau }\) and \(z_{n,j}^{\tau }\) are the cost factor and the active power generation of component \(n\) during period \(\tau\), respectively; \(\mathcal{N}\) is the dimensionality of the components. The total active power generation in terms of binary variable \(m_{n,j}^{\tau }\) is:

$$\left\{ \begin{aligned} p_{j}^{\tau } = \sum\limits_{n = 1}^{\mathcal{N}} {z_{n,j}^{\tau } } \hfill \\ m_{n + 1,j}^{\tau } \bar{z}_{n,j}^{\tau } \le z_{n,j}^{\tau } \le m_{n,j}^{\tau } \bar{z}_{n,j}^{\tau } \hfill \\ \end{aligned} \right.$$
(19)

We will describe this application of the THDDC method in detail in future papers.

4 Simulations

We simulated the Guyuan grid in north Hebei, China. Figure 6 shows that this network is connected to 18 wind farms; 8 × 12 Mvar shunt capacitors and 4 × 12 Mvar shunt reactors are deployed in the CB substation. None of these devices are supposed to turn on/off more than eight times per day. Meanwhile, the ±2 × 2% transformer in the GY substation is limited to no more than four tap changes per day.

Fig. 6
figure 6

Testing system connected to a number of wind farms

As the transient processes are not relevant to all sections of the problem, we simulated the states of the system after activating control procedures using MATPOWER. The mixed-integer linear programs for solving the TORDDC and MPCDDD problems were implemented using the CPLEX kernel. We used a sampling interval of 15 min in the multi-timescale optimization and validated the effectiveness of the proposed method by contrasting its performance with the performance of systems controlled using other methods. This validation is discussed in the later subsections. The first alternative method is the classical discrete-device control (CDDC) method, which only takes into account the present performance of the system. We also analyzed a non-robust method called SDDC, which bases its operation decisions of discrete devices on the day-ahead operation schedule [14].

4.1 Wind conditions and solution of TORDDC problem

The location of the power system in the Guyuan area is far from both the load center and conventional power resources. As a result, the system voltages are strongly affected by wind power. In this study, we used actual wind-power data from a typical day as base values and estimated the prediction error to be 20%. The wind-power forecast for the day ahead is shown in Fig. 7a. It takes 20-30 min to solve the TORDDC problem under these conditions, depending on the convergence gap. The number of potential discrete-device operations can be calculated from the optimal result, as shown in Fig. 7b.

Fig. 7
figure 7

Total wind power generation forecast for day ahead and number of permitted discrete device operations

From this figure, we can see that the potential discrete-device operations are distributed over 24 hours. This guarantees that the discrete devices will not reach their operational limits too quickly when the wind conditions fluctuate violently. In addition, most of the potential operations are clustered in time periods when the wind-power generation may vary significantly, e.g., before 10:00. This avoids compounding the large effect that fluctuations in wind power generation have on the system voltage. The OLTC influences the system voltage much more than the other discrete devices because it is located close to the root of the network. Shunt capacitors/reactors are required because many devices must be available from 14:00, when the OLTC operation limit is reached.

4.2 THDDC/CDDC/SDDC performance under typical wind condition

The within-day wind-power curves are randomly generated from values in the uncertainty set. The total real-time wind power generated is shown in Fig. 8. The curves of the total shunt capacitor/reactor reactive power and OLTC tap positions obtained using the different methods are plotted in Fig. 9.

Fig. 8
figure 8

Actual total wind power generated

Fig. 9
figure 9

Total shunt capacitor/reactor reactive power and OLTC tap positions

From Fig. 9, we can see that the CDDC method frequently uses the discrete devices’ response to fluctuations in wind-power generation between 00:00 and 03:30. This means that there are few shunt capacitors/reactors available after 05:45, and the OLTC tap position cannot change after 01:00. The SDDC and THDDC methods both take the day-long time coupling between the discrete devices into account and avoid using up the permitted number of operations. In addition, the THDDC method decides to operate the devices based on the actual wind power curve instead of the mid-value of the wind forecast for the day ahead, which is used by the SDDC method. We can see from Fig. 9a that the variations in the total shunt capacitor/reactor reactive power reserves are closer to those of the actual total wind power shown in Fig. 8 when they are controlled by the THDDC method.

In addition to the above, we are mostly concerned with the day-long 96-point comparisons among the performance of the three methods. Figure 10a shows the comparison between the total voltage excess:

$$I_{1}^{T/C/S} = \sum\limits_{{k \in {\mathcal{K}}}} {v_{k}^{\tau } }$$
(20)

where superscript T/C/S represents the results of the THDDC/CDDC/SDDC methods. The CDDC method exhausted the permitted discrete-device regulation operations, so the system voltage could barely be maintained within the acceptable boundaries during the period between 17:00 and 20:00. We also compared the total wind-power generation curtailment (21) for each method in Fig. 10b.

Fig. 10
figure 10

Day-long 96-point comparison between index \(I_{1}^{T/C/S}\), index \(I_{2}^{T/C/S}\) and index \(I_{3}^{S}\)-\(I_{3}^{T}\)

$$I_{2}^{T/C/S} = \sum\limits_{{j \in {\mathcal{J}}}} {\left( {r_{j}^{\tau } - p_{j}^{\tau } } \right)}$$
(21)

The THDDC method avoided curtailment in this scenario, while the other two methods performed worse. Although handicapped by the mismatch between the real wind power generation and the forecast, the advantage of the SDDC method over the CDDC method is not affected because the prediction error is relatively small. We then considered the total deviation in the wind farm reactive power output:

$$I_{3}^{T/C/S} = \sum\limits_{{j \in {\mathcal{J}}}} {s_{j}^{\tau } }$$
(22)

We only discuss the performance of the THDDC and SDDC methods because the CDDC method was no longer considered competitive due to its poor performance with respect to the previously discussed objectives, which are relatively important. The majority of data points in Fig. 10c are above the \(I_{3}^{S}\) = \(I_{3}^{T}\) curve, thus proving that the performance of the THDDC method is superior.

4.3 Comparisons among Monte Carlo simulations of each method

Simulation results over a single timescale are not sufficient to reveal the advantages and disadvantages of the three methods. Monte Carlo simulations of a variety of random scenarios are more informative. Owing to the centralized distribution of the wind-farm power output, there are strong correlations between the values of the original active power outputs during each period. Hence, we suppose that the total wind power output during period \(\tau\) is distributed evenly within its boundaries:

$$\chi^{\tau } \sim U\left( {\sum\limits_{{j \in \mathcal{J}}} {\underline{r}_{j}^{\tau } } ,\sum\limits_{{j \in \mathcal{J}}} {\bar{r}_{j}^{\tau } } } \right)\quad \forall \tau \in {\mathcal{T}}_{24}$$
(23)

In each sample, the original wind power output of the wind farm is:

$$r_{j}^{\tau } = \pi_{j}^{\tau } \chi^{\tau }$$
(24)

where \(\pi_{j}^{\tau }\) is the percentage factor of continuous control device \(j\) during period \(\tau\). We selected three typical indices to sketch the control performance.

  1. 1)

    The daily average bus voltage qualification rate:

    $$J_{1}^{T/C/S} = \sum\limits_{{\tau \in T_{24} }} {\sum\limits_{{k \in {\mathcal{K}}}} {\frac{{w_{k}^{\tau } }}{{T_{24} K}}} }$$
    (25)
    $$w_{k}^{\tau } = \left\{ {\begin{array}{*{20}l} 1 \hfill & {v_{k}^{\tau } = 0} \hfill \\ 0 \hfill & {v_{k}^{\tau } > 0} \hfill \\ \end{array} } \right.$$
    (26)

where \(w_{k}^{\tau }\) is a binary variable that represents whether the voltage of bus \(k\) during period \(\tau\) is within the acceptable range; \(T_{24}\) and \(K\) are the numbers of elements in sets \({\mathcal{T}}_{24}\) and \({\mathcal{K}}\), respectively.

  1. 2)

    The daily total wind curtailment:

    $$J_{2}^{T/C/S} = \sum\limits_{{\tau \in {\mathcal{T}}_{24} }} {\sum\limits_{{j \in {\mathcal{J}}}} {\left( {r_{j}^{\tau } - p_{j}^{\tau } } \right)\Delta t} }$$
    (27)

where \(\Delta t\) is the length of each period.

  1. 3)

    The total daily average deviation in the reactive power output from the wind farms:

    $$J_{3}^{T/C/S} = \sum\limits_{{\tau \in {\mathcal{T}}_{24} }} {\sum\limits_{{j \in {\mathcal{J}}}} {\frac{{s_{j}^{\tau } }}{{T_{24} }}} }$$
    (28)

We simulated 1000 random samples, with the prediction error set to 20%. Figure 11 shows comparisons between the performance of the different control methods with respect to these indices, where E(·) is the expectation function and P(·) is the probability function.

Fig. 11
figure 11

Comparison of performances of three methods with respect to indices \(J_{1}^{T/C/S}\), \(J_{2}^{T/C/S}\), and \(J_{3}^{S}\)-\(J_{3}^{T}\) with a prediction error of 20%

It is clear from Fig. 11a that both the THDDC and SDDC methods prevented voltage excess. On the contrary, index \(J_{1}^{C}\) is below 0.98 in the scenarios with prediction errors of 20%. In Fig. 11b, the CDDC method performs the worst, while the SDDC method always initiates a small amount of wind-power generation curtailment due to the prediction errors. Meanwhile, the THDDC method prevents the curtailment of wind-power generation in most scenarios. We only discuss the results for index \(J_{3}\) from the THDDC and SDDC methods. The overwhelming majority of data points in Fig. 11c are above the \(J_{3}^{S}\) = \(J_{3}^{T}\) curve, illustrating that the THDDC method maintains a larger reserve of continuous reactive power than the other methods.

4.4 Effect of prediction accuracy on performance of THDDC/CDDC/SDDC methods

The sample space used in the previous section was not large enough to reach the voltage safety boundary, so none of the three methods necessitated much wind-power generation curtailment. We need to assess how well the system operates when the voltage approaches the stability limit, and whether our proposed method works with larger uncertainties. Thus, we simulated 1000 samples with the prediction error set to 30% and plotted the results in Fig. 12.

Fig. 12
figure 12

Comparison of performances of three methods with respect to indices \(J_{1}^{T/C/S}\), \(J_{2}^{T/C/S}\) and \(J_{3}^{S}\)-\(J_{3}^{T}\) with a prediction error of 30%

It is clear from Fig. 12a that increasing the sample space causes voltage excess in all three methods. The THDDC method is best at maintaining the voltage within the acceptable limits because it estimates the variations in wind-power generation in advance. However, the larger prediction error caused the actual wind-power generation to deviate far from the forecast mid-value in more cases. This caused the voltage maintenance performance of the SDDC method to deteriorate further when the prediction errors were 30% rather than 20%. Figure 12b shows similar results. The THDDC method is again shown to produce the best results in Fig. 11c, as indicated by the majority of data points being above the \(J_{3}^{S}\) = \(J_{3}^{T}\) curve. We can conclude that the decrease in the prediction accuracy does less harm to the performance of the THDDC method than to the other methods.

5 Conclusion

In this paper, we have presented a THDDC method. Our method consists of two stages: ① two-stage robust optimization of the discrete device resource allocation for the day ahead; ② within-day MPC-based discrete-device dispatch. The optimal robust model used for forecasting the day ahead takes into account the time coupling between the discrete devices over the entire day and the uncertainties in wind power generation. We avoided overly conservative results using the time intervals when the discrete devices are permitted to operate as decision variables. The algorithm that runs throughout the day selects the specific moments to operate each device by solving a rolling multi-timescale optimization problem.

We explained the advantages of the THDDC method over other control algorithms under typical wind conditions with respect to control performance and computation speed. We validated the method using Monte Carlo simulations and confirmed its robustness by varying the prediction accuracy.