Solution to dynamic economic dispatch with prohibited operating zones via MILP

1.   Introduction

Dynamic economic dispatch (DED) is a fundamental problem for the optimal economic operation in power system which aims at allocating the customers' load demands among the available thermal power generating units in an economic, secure and reliable way ^[1]. In practice, some thermal or hydro generating units may have prohibited operating zones (POZ) due to the physical limitations of power plant components, e.g., vibrations in a shaft bearing are amplified in a certain operating region ^[2,3]. The resulting disjoint operating regions lead to a discontinuous generation cost function, which complicates the DED problem.

Comparing to the dynamic economic dispatch with prohibited operating zones (DED-POZ), the static economic dispatch with prohibited operating zones (SED-POZ) which handles only a single load level at a particular time instant is more achievable. Thus, in the past few decades, a myriad of optimization methods have been presented to deal with SED-POZ ^[4,5]. They include genetic algorithm (GA) ^[3,6], evolutionary programming (EP)^[7,8], particle swarm optimization (PSO) ^[9,10], $\lambda$-iterative technique ^[11,12], nonlinear programming (NLP) ^[13], semi-definite programming (SDP) ^[14], mixed integer quadratic programming (MIQP) ^[15,16,17], branch and bound (B & B) ^[18], etc. Although various optimization methods have been used to overcome the discontinuous solution space of SED-POZ, it still can not be easily implemented for the solution of DED-POZ. Because when POZ is taken into account, complexity of the DED problem will increase significantly since DED dispatches over a scheduled time horizon instead of one period. Moreover, additional constraints such as ramp rate constraints and spinning reserve constraints make the problem more complicated and can not be tackled easily.

More recent works for DED have been around heuristic methods, such as particle swarm optimization (PSO) ^[19,20], imperialist competitive algorithm (ICA) ^[21], harmony search (HS) ^[22], improved bees algorithm (IBA) ^[23], etc. In order to make the problem more close to practice, some influencing factors such as pollutant emission, renewable energy, demand response are integrated in DED model. In ^[24], the emission is included in the DED problem, and an adaptive tunicate swarm algorithm (TSA) is introduced for solving DED problem. In ^[25], emission, wind and photovoltaic power generation are considered in DED, and then a novel MFO\_PDU algorithm is proposed for such a DED problem to reduce the generation cost and emissions and improve the utilization efficiency of renewable energy. In ^[26], the short-term load forecasting is combined in the DED problem, and two artificial intelligence techniques were selected for solving them. In ^[27], demand response and renewable energy both developed for DED problem, and the effectiveness of the proposed model with respect to the operating cost, peak load reduction, deferred load have been demonstrated. Generally, heuristic techniques are quite sensitive to various parameter settings and solution may be different at each trial due to the intrinsic stochastic characteristic of heuristics. They do not provide an optimality gap so you have no clue how well of a solution you have obtained. Hybrid methods which combine several heuristic techniques or deterministic approach such as bacterial foraging PSO-DE algorithm (BPSO-DE) ^[28] and hybrid EP-PSO-SQP algorithm ^[29] tend to be more efficient than the individual method. However, they still have the intrinsic drawback of the stochastic search method. Unlike heuristic methods based on the stochastic search techniques, deterministic mathematical programming methods based on its solid mathematical foundations and the availability of the powerful software tools have more advantages. 1) They are deterministic for avoiding any stochastic process and various parameter settings, so that a deterministic and identical solution will be obtained at each trial. 2) They can be applied to any iterative solution method, stochastic or deterministic. 3) They can provide an optimal solution under some optimality criterions. 4) They are stable and scalable. So in ^[30], an efficient real-time approach based on optimality condition decomposition (OCD) technique is proposed to solve the DED-POZ. By using the reformulation and OCD technique, the problem could be decomposed to several simpler sub-problems and then the CPU-time can be reduced significantly, but the spinning reserve constraints are not considered. In ^[31], the DED-POZ formulates an MIQP model which can be solved by a mixed integer programming (MIP) solver immediately. Nevertheless, the complicated transmission losses are not included.

As the significant progresses of mathematical programming theory and the improvements of the MIP solvers, the MIP technique has become a popular alternative in the optimal operation of electric power system. So various MIQP formulations which are capable to solve to global optimality directly are presented for SED-POZ and DED-POZ, while the transmission losses are neglected. Although the solution of MIQP with commercial software such as CPLEX or GUROBI have significantly improved in recent years, it still does not a very good choice for DED-POZ on account of its non-linear characteristic. Comparing with the MIQP, mixed integer linear programming (MILP) is more developed because of the vastly superior warm-start capabilities of the simplex method ^[32]. Besides, from the perspective of engineering, an exact optimal solution of the problem is not always necessary and faster approximations have more value. Consequently, a natural alternative for solving DED-POZ is to reformulate the problem, yielding a reliable MILP formulation which can be solved efficiently by an MIP solver.

In this paper, a developed MILP method, which is based on two new MILP formulations, is introduced to solve the DED problem considering POZ, ramp rate constraints, transmission losses and spinning reserve constraints. From the perspective of mathematics, the original DED model considered in this paper is a nonconvex mixed integer nonlinear programming problem wich is very tricky. The characteristics of nonconvex and nonlinear accompany the vast majority models for power system operation. So in the future research, the MILP techniques proposed in this paper may be expanded to tackle more influence factors with nonconvex and nonlinear characteristics. The main contributions are summarized as follows:

$\bullet$ By using the perspective cut reformulation technique, a tight MILP formulation, denoted by MILP-1, which can be solved to global optimality within a preset tolerance via an MIP solver is obtained for DED.

$\bullet$ When transmission losses are considered, an MILP formulation for DED, denoted by MILP-2, is formed to tackle the complicated quadratic terms of transmission losses.

$\bullet$ We propose an MILP-iteration algorithm to solve the DED problem based on MILP-1 and MILP-2. Simulation results show that the proposed MILP formulations and MILP-iteration algorithm can solve to competitive solutions in a short time.

The rest of this paper is organized as follows. The mathematical formulation for the DED-POZ is described in Section 2. The reformulation and MILP-iteration algorithm are presented in Section 3. Simulation results on several test systems are given in Section 4. Finally, conclusions are drawn in Section 5.

2.   Mathematical formulation for the DED-POZ

The objective of DED-POZ is to minimize the total generation cost over the scheduled time horizon

The minimized DED-POZ should be subjected to the constraints as follows.

1) Power balance equations

where the transmission loss $P_{t}^{loss}$ can be calculated based on Kron's loss formula as follow ^[33]:

2) Power generation limits

3) Ramp rate limits

4) Prohibited operating zones limits

The prohibited operating zones of each unit can be characterised by the disjoint operating regions as shown below

where $P_{i1}^{\min} = P_i^{\min}$, $P_{i{n_i+1}}^{\max} = P_i^{\max}$.

5) Spinning reserve constraints

3.   Reformulation and solution for the DED-POZ

Owing to the disjoint operating zones, the classical mathematical programming methods are not suitable for DED-POZ any more. By introducing some auxiliary variables $P_{i, t}^{j}$ and binary variables $u_{i, t}^{j}$, the POZ constraint (2.6) can be equivalent to the following expressions

Consequently, the DED-POZ can be formulated as an MIQP formulation when the transmission losses are not included:

The DED-POZ, as formulated in (3.4), can be solved via MIQP directly. However, solution via MILP tends to be more efficient since the warm start capabilities of the simplex method available in MILP solver are vastly superior in comparison with the interior-point method in MIQP solver. Moreover, from the perspective of engineering, an exact optimal solution of the problem is not always necessary and faster approximations have more value. As a result, we present two reliable MILP formulations for solving the DED-POZ in the next subsections.

3.1.   Reformulation for the DED-POZ

According to $u_{i, t}^{j}P_{ij}^{\min}\le{P_{i, t}^{j}}\le{u_{i, t}^{j}P_{ij}^{\max}}$ and $u_{i, t}^{j}\in\{0, 1\}$, we have $P_{i, t}^{j}\in\{0\}\cup[P_{ij}^{\min}, P_{ij}^{\max}]$, which indicates that $P_{i, t}^{j}$ is a semi-continuous variable. With the help of (3.2), the objective function in (3.4) can be converted into the sum of quadratic functions on the semi-continuous variables space

Perspective cut reformulation technique ^[34,35] proposed by Frangioni et al. can thus be used for constructing a tight MILP formulation of the problem. Introducing some auxiliary variables $z_{i, t}^{j}(j = 1, ..., n_i+1)$, then (3.5) can be approximated by the following perspective cuts

where $P_{i}^{j(l)} = P_{ij}^{\min}+l(P_{ij}^{\max}-P_{ij}^{\min})/L$$(l = 0, 1, ..., L)$ and $L$ is a given parameter.

Then a tight MILP approximation, denoted as MILP-1, for DED-POZ without considering transmission losses can be formed:

3.2.   Linearization for the transmission losses

When the transmission losses are taken into account for the DED-POZ, it makes the optimization more difficult because of its complicated nonlinearity. Note that, the nonlinearity arises from the third term of (2.3), i.e., $P_{t}^{T}BP_{t}$. Replacing $P_{t}^{T}BP_{t}$ with auxiliary variable $c_t$, the transmission loss (2.3) can be rewritten as

As we can see, (3.9) is a linear constraint which can be addressed easily while (3.10) is a complicated quadratic constraint which is hard to tackle. A natural way to conquer this difficulty is to find an efficient linear approximation instead of the quadratic one. Then the first order Taylor expansion is employed for the $P_{t}^{T}BP_{t}$. As a result, the (3.10) can be replaced by

where $P_{t}^{(k)}$ is taken to be a constant vector corresponding to $P_{t}$.

When the transmission losses are considered, based on MILP-1 (3.8) and the above approximations, the DED-POZ can be formulated as the following MILP formulation, denoted as MILP-2,

It is well know that, when the vector $P_{t}^{(k)}$ we employ in (3.11) is sufficiently close to the optimal solution, a high efficiency MILP approximation for the primal problem can be obtained. Then the MILP-2 formulation can provide a near-optimal solution. To exploit such an efficient MILP approximation and solve to a reliable solution, based on MILP-1 and MILP-2 formulation, a straightforward MILP iteration algorithm, is proposed for DED-POZ when transmission losses are considered. In the following, the details of the MILP iteration algorithm are given.

Initialization Step: Choose a scalar $\epsilon > 0$ and a maximum number of iteration $iter_{max}$ to be used for terminating the algorithm and let $iter = 1$, $k = 3$.

Step 1: Solve MILP-1 (3.8) to obtain an optimal solution denoted by $P^{(1)}$, where $P^{(1)} = [P_{1}^{(1)}; P_{2}^{(1)}; ...; P_{T}^{(1)}]$.

Step 2: Solve MILP-2 (3.12) where the linear approximation (3.11) is taken at $P_{t}^{(1)}$, to obtain an optimal solution denoted by $P^{(2)}$.

Step 3: Solve MILP-2 (3.12) where the linear approximations (3.11) is taken at $(P_{t}^{(k-2)}+P_{t}^{(k-1)})/2$, to obtain an optimal solution denoted by $P^{(k)}$.

Step 4: Calculate the violation of power balance $E_t$, where

and the $P_{t}^{loss(k)}$ is calculated according to (2.3).

Step 5: When

or when $iter = iter_{max}$, the procedure is terminated.

Step 6: Let $iter = iter+1$ and $k = k+1$. Go to Step 3.

The procedure of the MILP iteration algorithm can be summarized as follows.

Because the transmission loss at each period in a DED problem is small compared to the corresponding load demand, when transmission losses are ignored, after optimizing MILP-1 (3.8), an initial approximate optimal solution can be obtained in Step 1. Since transmission losses are ignored, the solution obtained in step 1 is "relative small" in some ways. When transmission losses are included, by solving MILP-2 (3.12) where the linear approximation (3.11) is taken at such a "relative small" solution, a "relative large" solution can be got in Step 2. Then average of the former two solutions is used in Step 3 to balance the power balance equation. Repeat from Steps 3 to 6 until all the violations of power balances are smaller than the preset value or the maximum number of iteration is reached. Generally, several iterations are enough in our algorithm since the transmission loss is small and $P_{t}^{T}BP_{t}$ is only a portion of the transmission loss.

4.   Simulation results

To assess the efficiency of the proposed MILP-1 formulation and MILP iteration algorithm for DED-POZ, two cases, ignoring the transmission losses and considering the transmission losses, are simulated in our numerical experiments. The formulation and algorithm are coded with Matlab and optimized by using CPLEX 12.6.2 ^[36]. The machine for all runs is an Intel Core 2.3 GHz HUAWEI-notebook with 16 GB of RAM.

4.1.   Ignoring the transmission losses

In this subsection, a set of different sizes test systems with units ranging from 6 to 180 over a scheduled time horizon of 24 h are adopted for testing the effectiveness of MILP-1 formulation, where the transmission losses are not considered. The 6-unit system data is taken from ^[19]. The 1-h spinning reserve requirement is 5% of the load demand. Other test systems with 30, 60, 120, and 180 units are obtained by duplicating the 6-unit system 5, 10, 20, and 30 times, respectively.

We directly solve the MIQP (3.4) and MILP-1 (3.8) formulations using CPLEX to 0.01% optimality in a time limit of 300 s. In MILP-1 formulation, $L$ is set to 4. The results are compared in Table 1. As we can see in Table 1, CPLEX can solve the MILP-1 faster than the MIQP. For MILP-1, solutions for all systems can be gained in several seconds, while for MIQP, only the 6-, 30- and 60- unit systems can be solved to the preset tolerance within 300 s. And as the size of the problem increases, the computation for MIQP will be very time-consuming. On the other hand, although the total generation costs for the 30- and 60- unit systems obtained by solving MILP-1 are slightly larger than the costs obtained by solving MIQP, but the differences are very small which can be negligible. Actually, in electrical engineering, getting a high quality solution quickly is become more significant than a very time-consuming solution with a slight improvement. Furthermore, for the 120- and 180- unit systems, MILP-1 can exploit more lower total generation costs within several seconds. Based on the above numerical results and analysis, we can see that, when transmission losses are not considered, solving MIQP formulation may get a greater total generation cost than solving MILP-1. It means that, although MIQP formulation can get a global optimal solution in theory. However, in the actual calculation, it is very hard to achieve a truly global optimal solution when the scale of the problem increases. For the approximate model MILP-1, since it benefits from the perspective cut reformulation technique and the vastly superior warm-start capabilities of the simplex method, it can be solved to a high quality solution in a very short time, which confirms the effectiveness of MILP-1. The outputs for the 6-unit system are given in Table 2 for verification.

4.2.   Considering the transmission losses

In this subsection, a 15-unit system considering the transmission losses is adopted for simulation. Although the system has only 15-units, it is a dynamic problem with 750 variables coded in a solution. The characteristics of the thermal units and load demands are taken from ^[19]. Owing to the limits of space, the loss coefficients with the 100-MVA base capacity are not listed here. One can refer to ^[9]. For fair comparison, the spinning reserve requirement is 5% of the load demand. In our MILP-iteration algorithm, parameters $\epsilon$, $iter_{max}$ and $L$ are set to 0.1, 10 and 4, respectively. Meanwhile, in MILP-iteration algorithm, we directly solve the MILP-1 and MILP-2 formulations using CPLEX to 0.01% optimality. The results obtained by MILP-iteration algorithm and other methods ^[19,20,37] for the 15-unit system are shown in Table 3.

In Table 3, PSO ^[19], Hybrid HNN^[37] and EPSO ^[20] are heuristic methods for DED. The total generation costs obtained by these methods are ＄774, 131, ＄759, 796 and ＄759, 410, respectively. However, the total generation cost obtained by our MILP-iteration algorithm is ＄759, 176, which is much less than the costs obtained by other three methods. Particularly, our method reduces the total generation cost of ＄14, 955 comparing with the PSO method. Furthermore, comparing with other three methods, the time consumed for our algorithm is only 0.74 s, which is shorter than other three methods. More specifically, the computational times of PSO, hybrid HNN and EPSO are 4.47, 3 and 13.51 times than that of our method. It means that, our MILP-iteration algorithm can solve to a lower total generation cost in a faster speed than other listed methods. As stated before, lower cost does not mean a good solution. The violations of power balances must be checked. In Table 4, the violations for PSO, EPSO and MILP-iteration algorithm are calculated. The violations for Hybrid HNN is not calculated since the outputs are not given in ^[37].

We can see from Table 4, all the violations for MILP-iteration algorithm are smaller than PSO and EPSO, and most of them are much smaller than them. In PSO, the maximum violation is 1.0418 MW, and in EPSO, the maximum violation is 0.5416 MW. But in our algorithm, the maximum violation is 0.0949 MW and the total violation is only 0.5175 MW. From the simulation results we can conclude that our MILP-iteration algorithm can solve to a better solution in a short time for DED-POZ.

The outputs and transmission losses for the 15-unit system obtained by MILP-iteration algorithm are given in Table 5 for verification.

5.   Conclusions

Difference from various conventional MIQP methods and heuristic methods, in this paper, a developed MILP method has been successfully introduced to solve the DED problem considering POZ, ramp rate constraints, transmission losses and spinning reserve constraints. 1) By using the perspective cut reformulation and first order Taylor expansion techniques, two new MILP formulations, i.e., MILP-1 and MILP-2, are formed for the two cases of the DED-POZ. 2) Based on these two MILP formulations, a straightforward MILP-iteration algorithm which can benefit from the vastly superior warm-start capabilities of the simplex method, is proposed to optimize the DED-POZ. 3) In order to assess the quality of the solutions, not only the optimality but also the feasibility are discussed. The simulation results show that both MILP formulation and MILP-iteration algorithm can solve to competitive solutions in a short time. In other words, the proposed MILP method can solve the DED-POZ problem efficiently.

The impact of pollutant emission, renewable energy and demand response are not taken into consideration in this paper. Therefore, in the future research, these affecting factors will be incorporated into the DED model while the idea and techniques proposed in this paper may be expanded to tackle some of these limitations.

Acknowledgments

This work was supported by Natural Science Foundation of Guangxi (2021GXNSFBA075012, 2020GXNSFDA238017); Middle-aged and Young Teachers' Basic Ability Promotion Project of Guangxi(2020KY08005); the Doctoral Foundation of Guangxi University of Science and Technology (19Z40).

Conflict of interest

The authors declare there is no conflict of interest.