Solution to dynamic economic dispatch with prohibited operating zones via MILP

Dynamic economic dispatch (DED) problem considering prohibited operating zones (POZ), ramp rate constraints, transmission losses and spinning reserve constraints is a complicated non-linear problem which is difficult to solve efficiently. In this paper, a mixed integer linear programming (MILP) method is proposed to solve such a DED problem. Firstly, a novel MILP formulation for DED problem without considering the transmission losses, denoted by MILP-1, is presented by using perspective cut reformulation technique. When the transmission losses are considered, the quadratic terms in the transmission losses are replaced by their first order Taylor expansions, and then an MILP formulation for DED considering the transmission losses, denoted by MILP-2, is obtained. Based on MILP-1 and MILP-2, an MILP-iteration algorithm (MILP-IA) is proposed to solve the complicated DED problem. The effectiveness of the MILP-1 and MILP-IA are assessed by several cases and the simulation results show that both of them can solve to competitive solutions in a short time.


Introduction
Dynamic economic dispatch (DED) is a fundamental tool 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. They include genetic algorithm (GA) [3,4], evolutionary programming (EP) [5,6,7,8], particle swarm optimization (PSO) [9,10,11,12], λ-iterative technique [13], nonlinear programming (NLP) [14], semi-definite programming (SDP) [15], mixed integer quadratic programming (MIQP) [16,17,18], branch and bound (B&B) [19], 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. 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-POZ have been around heuristic methods, such as particle swarm optimization (PSO) [20,21,22], imperialist competitive algorithm (ICA) [23], harmony search (HS) [24], improved bees algorithm (IBA) [25], etc. However, heuristic techniques are quite sensitive to various parameter settings and solution may be different at each trial due to the intrinsic stochas-tic 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) [26] and hybrid EP-PSO-SQP algorithm [27] tend to be more efficient than the individual method. However, they still have the intrinsic drawback of the stochastic search method we mentioned above. Unlike heuristics, deterministic algorithms can solve to a robust result due to the solid mathematical foundations and the availability of the powerful software tools. In [28], 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 [29], 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. In [16,17,18,29], 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 [30]. 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.
To find such an MILP formulation for DED-POZ, auxiliary variables are introduced not only for the constraints but also the objective function in our work.
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. When transmission losses are consid- 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-IA are presented in Section 3. Simulation results on several test systems are given in Section 4. Finally, conclusions are drawn in Section 5.

Mathematical formulation for the DED-POZ
The objective of DED-POZ is to minimize the total generation cost over the scheduled time horizon where P i,t is the power output of unit i in period t; N is the total number of units; T is the total number of periods; α i , β i and γ i are the cost coefficients of unit i.
The minimized DED-POZ should be subjected to the constraints as follows.

1) Power balance equations
where D t is the load demand in period t; P loss t is the transmission loss in period t, which can be calculated based on Kron's loss formula as follow [31]: where B 00 , B 0 , B are B-coefficients; P t = [P 1,t , P 2,t , ..., P N,t ] T is the power output vector in period t.

2) Power generation limits
where P min i and P max i are the minimum and maximum power outputs of unit i.

3) Ramp rate limits
where RD i and RU i are the ramp-down and ramp-up rates of unit i.

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 min , n i is the number of prohibited operating zones of unit i.

5) Spinning reserve constraints
where SR i,t is the spinning reserve provided by unit i in period t and R t is the system spinning reserve requirement in period t.

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 j i,t and binary variables u j i,t (j = 1, ..., n i + 1), the POZ constraint (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: (2), (4), (5), (7), (8), (9), (10).
The DED-POZ, as formulated in (11), 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.

Reformulation for the DED-POZ
, which indicates that P j i,t is a semi-continuous variable. With the help of (9), the objective function in (11) can be converted into the sum of quadratic functions on the semi-continuous variables space Perspective cut reformulation technique [32,33] proposed by Frangioni et al.
can thus be used for constructing a tight MILP formulation of the problem.
Introducing some auxiliary variables z j i,t (j = 1, ..., n i + 1), then (12) can be approximated by the following perspective cuts where P .., L) and L is a given parameter.

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 (3), i.e., P T t BP t . Replacing P T t BP t with auxiliary variable c t , the transmission loss (3) can be rewritten as As we can see, (16) is a linear constraint which can be addressed easily while (17) 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 (17) can be replaced by where P (k) t is taken to be a constant vector corresponding to P t .
When the transmission losses are considered, based on MILP-1 and the above approximations, the DED-POZ can be formulated as the following MILP for- It is well know that, when the vector P given.

Initialization
Step: Choose a scalar > 0 and a maximum number of iteration iter max to be used for terminating the algorithm and let iter = 1, k = 3.
Step 4: Calculate the violation of power balance E t , where and the P loss(k) t is calculated according to (3).
Step 5: When or when iter = iter max , the procedure is terminated.
Because the transmission loss at each period in a DED problem is small com-

Simulation results
To assess the efficiency of the proposed MILP-1 formulation and MILP-IA 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 [34]. The machine for all runs is an Intel Core 2.5 GHz Dell-notebook with 8 GB of RAM.

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, where the transmission losses are not considered. The 6unit system data is taken from [20]. The 1-h spinning reserve requirement is 5% of the load demand. Other test systems with 30, 60, 120, and 180 units can be obtained by duplicating the 6-unit system 5, 10, 20, and 30 times, respectively.
In our MILP-1 formulation, L is set to 4. We directly solve the MIQP (11) and MILP-1 (15) formulations using CPLEX to 0.01% optimality in a time limit of 300 s, and the results are compared in Table 1. As we can see in Table 1 The outputs obtained by solving MILP-1 for the 6-unit system are given in Table 2 for verification.

Considering the transmission losses
In this subsection, a 6-unit system and a 15-unit system, which have been widely used for testing [20,21,25] [20]. 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-IA, parameters , iter max and L are set to 0.1, 5 and 4, respectively. Meanwhile, in MILP-IA, we directly solve the MILP-1 and MILP-2 formulations using CPLEX to 0.01% optimality.
Since different computers are used, the run time in different methods may not be directly comparable. In order to have a fair comparison regarding the computational effort, the CPU chip frequency from the used computer is used to convert the CPU times obtained from different methods into a common base [35]. The scaled CPU time of an algorithm can be computed as follow: For fair comparison, the scaled CPU time is used in this subsection and the base CPU speed is 2.0 GHz.

6-unit system considering the transmission losses
In this subsection, a 6-unit system considering the transmission losses is adopted for simulation. In this system, all the units have POZ constraints. The results obtained by MILP-IA and IBA [25] are shown in Table 3. We can see that our MILP-IA can save much more execution time than IBA. Although IBA can obtain much lower total generation cost, but we should note that, not only the optimality but also feasibility should be taken into account for evaluating the quality of the solution. In DED, the constraints such as power generation limits, ramp rate limits, POZ limits and spinning reserve constraints are linear constraints which can be easily satisfied in most algorithms. But the power balance equations are quadratic equality constraints which are hard to meet entirely. In fact, thanks to the MIP solver, the feasibility of those linear constraints can be easily guaranteed in MILP-IA. Therefore, our discussion with respect to the feasibility validation mainly focus on the violations of power balances.
Since the losses coefficients are specified in per unit on a 100-MVA base in our test systems, then the real power loss is given by [31]: In the cost function P i,t is expressed in MW. Therefore, the real power loss in terms of MW generation is [31]: Hence, by using (24) and (20), the violations of power balances can be calculated easily.
The comparison results with respect to the violations for the 6-unit system between IBA and MILP-IA are shown in Table 4.  As we can see in Table 4 shows that MILP-IA can solve to an acceptable near-optimal solution for the 6-unit system within 1 s.
The outputs and transmission losses for the 6-unit system obtained by MILP-IA for DED-POZ are given in Table 5 for verification.

15-unit system considering the transmission losses
In this subsection, a 15-unit system considering the transmission losses is adopted for simulation. In this system only 4 units have POZ constraints. The results obtained by MILP-IA and other two PSO methods [21] are shown in Table 6.
As we can see, our MILP-IA can solve to a lower total generation cost in a faster speed than other two PSO methods. As stated before, lower cost does not mean a good solution. The violations of power balances must be checked. In Table 7, only the violations for EPSO and MILP-IA are calculated, since the solutions for PSO are not available in [21]. Figure 1 depicts the details as well.  The outputs and transmission losses for the 15-unit system obtained by MILP-IA are given in Table 8 for verification.

Conclusion
In this paper, an MILP method has been successfully introduced to solve the DED problem considering POZ, ramp rate constraints, transmission losses and spinning reserve constraints. By using the perspective cut reformulation and first order Taylor expansion techniques, two MILP formulations, i.e., MILP-1 and MILP-2, are formed for the two cases of the DED-POZ. Based on these two MILP formulations, a straightforward MILP-IA is proposed to optimize the DED-POZ. 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-1 and MILP-IA can solve to competitive solutions in a short time. In other words, the proposed MILP method can solve the DED-POZ problem efficiently.