A New Approach to Optimal Allocation of Reactive Power Ancillary Service in Distribution Systems in the Presence of Distributed Energy Resources

Abstract

One of the most important Distribution System Operators (DSO) schemes addresses the Volt/Var control (VVC) problem. Developing a cost-based reactive power dispatch model for distribution systems, in which the reactive powers are appropriately priced, can motivate Distributed Energy Resources (DERs) to participate actively in VVC. In this paper, new reactive power cost models for DERs, including synchronous machine-based DGs and wind turbines (WTs), are formulated based on their capability curves. To address VVC in the context of competitive electricity markets in distribution systems, first, in a day-ahead active power market, the initial active power dispatch of generation units is estimated considering environmental and economic aspects. Based on the results of the initial active power dispatch, the proposed VVC model is executed to optimally allocate reactive power support among all providers. Another novelty of this paper lies in the pricing scheme that rewards transformers and capacitors for tap and step changing, respectively, while incorporating the reactive power dispatch model. A Benders decomposition algorithm is employed as a solution method to solve the proposed reactive power dispatch, which is a mixed integer non-linear programming (MINLP) problem. Finally, a typical 22-bus distribution network is used to verify the efficiency of the proposed method.

1. Introduction

The conventional Volt/Var control aims to find appropriate coordination between the on-load tap changer (OLTC) and all of the switched shunt capacitors (Sh.Cs) in the distribution networks. The main goal of a VVC system is to achieve an optimum voltage profile over the distribution feeders and optimum reactive power flows in the system [1,2,3]. Recently, because of the integration of various types of Distributed Generations (DGs), some challenging issues have emerged in distribution system operation. Generally, DGs are expected to increase the number of switching operations of conventional voltage control devices, such as OLTCs and Sh.Cs [4]. At present, inverters coupled with DG units can provide both active and reactive power based on the DSO request. Hence, by growing DG penetration into distribution systems, DG units could incorporate daily VVC.

Nowadays, all studies of VVC can be classified into two main frameworks: centralized offline control and real-time control. Studies on centralized offline control aim to determine dispatching schedules of VVC devices according to the day-ahead load forecast. In these studies, various objective functions such as total energy cost offered by generation units, electrical energy losses, voltage deviations, and total emission of generation units have been adopted as the methodologies for managing VVC. In [5] an Ant Colony Optimization (ACO) algorithm has been adopted to optimize the total cost of electrical energy generated by Distribution Companies (Discos) and DGs in the daily VVC problem. In [6] a fuzzy price-based compensation methodology has been proposed to solve the daily VVC problem in distribution systems in the presence of DGs. In [7] a new optimization algorithm based on a Chaotic Improved Honey Bee Mating Optimization (CIHBMO) has been implemented to determine control variables including the active and reactive power of DG units, reactive power values of capacitors, and tap positions of transformers for the next day. Also, in [8] a Fuzzy Adaptive Chaotic Particle Swarm Optimization (FACPSO) has been introduced to solve the multi-objective optimal operation management of distribution networks including fuel-cell power plants. In [9], minimization of active power losses and micro-generation shedding have been suggested as a methodology for optimized and coordinated voltage support in distribution networks with large integration of DGs and micro-grids. In [10], an analytic hierarchy process (AHP) strategy and Binary Ant Colony Optimization (BACO) algorithm have been employed to solve the multi-objective daily VVC in distribution systems. In [11] a multi-objective θ-Smart Bacterial Foraging Algorithm (Mθ-SBFA) has been used for daily VVC, considering environmental and economic aspects as well as technical issues of distribution networks. According to progress in the Wind Turbine (WT) technology, Discos have paid attention to WTs more than any other Renewable Energy Sources (RESs). The stochastic nature of the wind speed may cause a fluctuation of electrical power in the distribution systems. Thus, a probabilistic analysis of distribution systems is required to cope with all uncertainties caused by the wind speed variations and load fluctuations [12,13]. The ability of doubly fed induction generator (DFIG)-based wind farms to deliver multiple reactive power objectives considering variable wind conditions is examined in [14].

On the other hand, studies on real-time control methods have helped to control the VVC equipment based on real-time and local measurements and experiences. The second framework of VVC requires a higher level of distribution system automation and more hardware and software supports [15]. This control methodology generally provides no coordination between devices, and is often limited to a unidirectional power flow. Furthermore, it is very difficult for a real-time controller to take into account the overall load change as well as the constraints of the maximum allowable daily operating times of switchable equipment. Based on SCADA capabilities and communication infrastructure, in [16,17], a real-time reactive power control has been implemented to optimally control the switched capacitors in distribution systems in order to minimize system losses and maintain admissible voltage profile. In [18], a new real-time voltage control method has been discussed, using load curtailment as a part of demand response programs to regulate voltage of the distribution feeders within their allowable ranges.

However, few studies have been done about the daily VVC problem in distribution systems, in which the cost of reactive power support by various types of DERs has been considered. Hence, the establishment of a fair payment method for reactive power ancillary service of DER is necessary. Also, the reactive power capability (Q-capability) of DERs, especially the Q-capability of WTs considering wind speed fluctuations, has not been taken into account in previous studies of the VVC problem. These deficiencies motivated us to formulate a new pricing model for reactive power service from DERs, including synchronous machine-based DGs and WTs. This article also presents a new day-ahead active power market to minimize the electrical energy costs and the gas emissions of generation units. Along with this active power market, a novel reactive power dispatch framework is introduced to minimize the total cost of the following components: adjustment of the initially scheduled active powers, total active power losses, reactive power provided by DERs and Disco, and depreciation cost of the switchable facilities such as OLTC and Sh.Cs in order to achieve an economic plan for the daily VVC problem. Due to the presence of control devices such as DERs, OLTCs, Sh.Cs, etc., the daily VVC of a distribution system is a Mixed Integer Nonlinear Programming (MINLP) optimization problem. The complexity of solving this nonlinear optimization problem comprising integer variables and a great number of continuous variables forced us to employ decomposition techniques, such as the Benders decomposition algorithm.

The innovative contributions of this paper are summarized as follows:

• A new methodology is presented to determine the cost of reactive power support provided by DERs.
• The Q-capabilities of DGs and renewable energy sources are involved in the daily VVC problem.
• A novel reactive power dispatch framework is developed for daily VVC.

The remainder of this paper is organized as follows: Section 2 describes the cost of reactive power production from DERs. In Section 3, the day-ahead active power market is introduced. The proposed daily VVC model is presented in Section 4. Section 5 presents the Benders decomposition algorithm, which is used as a solution methodology to solve the proposed reactive power dispatch model. The simulation results are discussed in Section 6, and the conclusions of this paper are reported in Section 7.

2. Cost of Reactive Power Production from DERs

2.1. Synchronous Machine-Based DGs

An offer price framework used by a synchronous generator was developed in [19] for a reactive power market. In this section, a new framework is formulated to determine the cost of reactive power support for a synchronous machine-based DG (briefly speaking, DG). The capability curve of a synchronous generator is depicted in Figure 1. This diagram demonstrates the relationship between active and reactive power generated by this generator. In this study, it is assumed that there is a mandatory reactive power requirement approved by grid code for connectivity of a DG to the grid. We consider the grid code requirement to be such that the DG units should operate between a mandatory leading power factor ($p{f}_{mand} )$ and lagging $p{f}_{mand}$ at any operating point.

Figure 1. Capability diagram of a synchronous generator.

Figure 1. Capability diagram of a synchronous generator.

Three operating regions for a DG on the reactive power generation (at an active power $P_A$) can be identified as follows:

Region I: $-Q_{mand}^A$to $Q_{mand}^A$: Reactive power produced in this region is based on the grid code requirement. It is suggested that the DG is paid an availability cost (with a price ${\mathsf{\rho }}_0$ in $), which is a fixed component implying the portion of a supplier’s capital cost that can be contributed to reactive power production.

Region II: ($Q_{mand}^A$to$Q_A$ ) and ( $Q_{DG}^{min}$to$-Q_{mand}^A$): In this area, the DG is demanded by the DSO for additional reactive power provision beyond the amount of mandatory reactive power support, without requiring reduction of its generated active power. Due to increased losses in the windings, the DG can expect to receive a payment for its services. This payment consists of two components, referred to as the availability component and the cost of loss component. In this region two loss payment components are defined as price ${\mathsf{\rho }}_1$. in $/MVArh for operating in the under excitation mode and price ${\mathsf{\rho }}_2$ in $/MVArh for operating in over excitation mode [19].

Region III: ( $Q_A$to$Q_{DG}^{max}$): In this region the DG is requested by DSO to reduce its active power production so that the system reactive power requirements are fulfilled. Thus, together with the two other components, it deserves to receive an additional payment according to its opportunity cost of reduced active power generation. In this study, an adjustment bid is used for opportunity cost payment [20]. The adjustment bids are employed by the generation units to indicate information about the prices received for the reduction (change) in its initial active power scheduled by the Market Operator (MO). Also, along with this price, the maximum changes regarding the initial active power accepted are specified and sent to the DSO. In region III, if the corresponding change in the active power, which has a negative value, is represented by the $\Delta P^{Adj}$ variable, the opportunity cost payment can be calculated by multiplying $\left|\Delta P^{Adj}\right|$ by the adjustment price (${\mathsf{\rho }}_{\text{Adj}}$).

Moreover, a change in the initial generation schedule of a DG may be due to the enforcement of operation constraints or to allow a rescheduling of initial active power output to compensate for active power reductions of other DGs operated in region III. This change can take both positive and negative values, which in both cases are compensated for by the adjustment price.

Accordingly, the cost of reactive power production for each DG unit can be mathematically represented by the following equation, assuming that $P_A=P_0$ represents the initial active power schedule:

$$COST\left(Q_{DG}\right)=\{\begin{array}{c}{\rho }_0\\ {\rho }_0-{\rho }_1(Q_{1DG}+Q_{mand}^A)\\ {\rho }_0+{\rho }_2(Q_{2DG}-Q_{mand}^A)\\ {\rho }_0+{\rho }_2(Q_B-Q_{mand}^B)+{\rho }_{Adj}\left|\Delta P^{Adj}\right|\end{array}\begin{array}{c}-Q_{mand}^A\le Q_{0DG}\le Q{}_{mand}^A\\ Q_{DG}^{\min }\le Q_{1DG}\le -Q_{mand}^A\\ Q_{mand}^A\le Q_{2DG}\le Q_A\\ Q_A\le Q_{3DG}=Q_B\le Q_{DG}^{\max }\end{array}$$

(1)

where $\Delta P^{Adj}$ is:

$$\Delta P^{Adj}=P_B-P_A$$

(2)

In order to model the constraints corresponding to the capability diagram of a DG, the curve between $\left(0,Q_{DG}^{max}\right)$ and $\left(P_{DG}^{max},Q_{DG}^x\right)$ in the capability diagram can be approximated by linear expression, as will be used in Equation (6). Therefore, the DSO will need to know the values of $P_{DG}^{max}, Q_{DG}^{max}, Q_{DG}^{min}, \text{and} Q_{DG}^x$for each DG unit. From a mathematical viewpoint, Equation (1) can be modeled using the binary variables $W_0, W_1$, $W_2$, and a set of algebraic relations as follows:

$$COST\left(Q_{DG}\right)=W_0·{\rho }_0+W_1·\left({\rho }_0-{\rho }_1\left(Q_{1DG}-Q_{mand}\right)\right)+W_2·\left({\rho }_0+{\rho }_2\left(Q_{2DG}-Q_{mand}\right)\right)+{\rho }_{Adj}·\left|\Delta P^{Adj}\right|$$

(3)

$$-W_0·Q_{mand}\le Q_{0DG}\le W_0·Q_{mand}$$

(4)

$$W_1·Q_{min}\le Q_{1DG}\le -W_1·Q_{mand}$$

(5)

$$W_2·Q_{mand}\le Q_{2DG}\le W_2·\left(Q_{DG}^{max}-\frac{Q_{DG}^{max}-Q_{DG}^x}{P_{DG}^{max}}(P_0+\Delta P^{Adj})\right)$$

(6)

$$Q_{mand}=(P_0+\Delta P^{Adj})·tan\left(co{s}^{-1}\left(p{f}_{mand}\right)\right)$$

(7)

$$Q_{DG}=Q_{0DG}+Q_{1DG}+Q_{2DG}$$

(8)

$$W_0+W_1+W_2\le 1$$

(9)

In Equation (7), when the DG is operated in region III, the $\Delta P^{Adj}$ variable is activated; otherwise it will be zero. Equation (9) ensures that the DG operates in only one of the three defined regions. Finally, Equation (3) can be rewritten as:

$$COST\left(Q_{DG}\right)=C_{DG}\left(Q_{DG}\right)+{\rho }_{Adj}·\left|\Delta P^{Adj}\right|$$

(10)

where $C_{DG}\left(Q_{DG}\right)$ represents the cost of reactive power support regardless of LOC as follows:

$$C_{DG}\left(Q_{DG}\right)=W_0{\rho }_0+W_1·\left({\rho }_0-{\rho }_1\left(Q_{1DG}-Q_{mand}\right)\right)+W_2·\left({\rho }_0+{\rho }_2\left(Q_{2DG}-Q_{mand}\right)\right)$$

(11)

2.2. Wind Turbines

In order to accomplish a cost-based reactive power dispatch, the different cost components associated with the reactive power of WTs should be determined. Determination of these cost components will help the DSO in managing appropriate financial compensation mechanisms in the daily VVC problem.

2.2.1. Reactive Power Capability of a WT

In [21], the capability curve for a WT with a full power back-to-back converter together with its reactive power cost model have been well addressed; the presented results are utilized and modified in this paper. The active power generation of a WT is a function of the wind speed, as given in [22]. The maximum reactive power provision capability of the WT is given as:

$$Q_{WT}=min\left\{Q_c , Q_v\right\}$$

(12)

where:

$$Q_c=\sqrt{{\left(V_g{I}_{c,max}\right)}^2-P_{WT}^2}$$

(13)

$$Q_v=\sqrt{{\left(\frac{V_g{V}_{c,max}}{X}\right)}^2-P_{WT}^2}-\frac{V_g^2}{X}$$

(14)

where X represents the total reactance of the WT transformer, the grid filters, and the reactance of the transformer adapting the WT’s voltage to the grid voltage. Due to the deviation of the actual wind power from the forecasted value, delivery of the reactive power given in Equation (12) based on the predicted active power value ($P_{WT}^h)$ cannot be uniformly ensured by the WT over the entire hour of operation. It is assumed that the maximum hourly variation of wind power from the forecasted value ($\Delta P_{max}$) is estimated based on the previous meteorological data. Hence, the available reactive power of the WT at the hth hour ($Q_{WT}^{h,av}$) can be determined as follows [21]:

$$Q_{WT}^{h,av}=\min \left\{Q_c^{h,av}, Q_v^{h,av}\right\}$$

(15)

where

$$Q_c^{h,av}=\sqrt{{\left(V_g{I}_{c,max}\right)}^2-{\left(P_{WT}^{h,max}\right)}^2}$$

(16)

$$Q_v^{h,av}=\sqrt{{\left(\frac{V_g{V}_{c,max}}{X}\right)}^2-{\left(P_{WT}^{h,max}\right)}^2}-\frac{V_g^2}{X}$$

(17)

$$P_{WT}^{h,max}=min\left\{\left(P_{WT}^h+\Delta P_{max}\right),P_{WT}^{rated}\right\}$$

(18)

2.2.2. Cost Components of Reactive Power of a WT

Fixed Cost Component

This component implies the additional cost imposed by the modifications of the converter to meet the reactive power support of the WT in accordance with network code requirements. Therefore, operating within the leading power factor $p{f}_{mand}$to lagging power factor $p{f}_{mand}$ at any active power generation can be compensated for by a fixed cost ($a_0)$.

Cost of Losses Component

Due to the increased reactive power demanded by DSO between the mandatory and available values, the WT will suffer an extra active power loss that should be determined for the purpose of financial compensation. The increased active power loss resulting from the increased reactive power supply from $Q_{mand}$ (mandatory reactive power production specified by grid codes) to $Q_1$ at active power production of $P_0$is expressed as:

$$\Delta P_{loss}=a_l\left(I_1-I_0\right)+b_l\left({I_1}^2-{I_0}^2\right)$$

(19)

$$I_0=\frac{\sqrt{P_0^2+Q_{mand}^2}}{V_g}$$

(20)

$$I_1=\frac{\sqrt{P_0^2+Q_1^2}}{V_g}$$

(21)

where $a_l$ and $b_l$ represent the active power loss constants. Using the MCP, the cost of losses component to cover the increased active power losses in the converter is calculated as [21]:

$$CLC=MCP·\Delta P_{loss}$$

(22)

Opportunity Cost Component

If the WT is requested to reduce its active power generation to provide the required reactive power, the WT forgoes the revenue cost due to the lost opportunity to sell its active power in the energy market. The available reactive power capability from a wind generator is$Q_{WT}^{h,av}$ as illustrated in Figure 2.

If the WT is requested by the DSO to produce $Q_{WT}^{h,req}$, which is more than$Q_{WT}^{h,av}$, then its active power generation is restricted to$P_{WT}^{h,lim}$. Thus, the WT has to forgo $\left|\Delta P_{WT}^{Adj,h}\right|$ active power production in the energy market, where

$$\Delta P_{WT}^{Adj,h}= P_{WT}^{h,lim}-P_{WT}^{h,max}$$

(23)

Therefore, the LOC can be calculated using the adjustment price of WT as follows:

$$LOC={\rho }_{Adj}·\left|\Delta P_{WT}^{Adj,h}\right|$$

(24)

Cost of reactive power production from a WT can be expressed as the following equation:

$$COST\left(Q_{WT}^h\right)=\{\begin{array}{c}{a}_0\\ a_0+MC{P}^h·\left(\Delta P_{loss}^h\left(Q_{1WT}^h\right)\right)\\ a_0+MC{P}^h·\left(\Delta P_{loss}^h\left(Q_{WT}^{h,av}\right)\right)+{\rho }_{Adj}.\left|\Delta P_{WT}^{Adj,h}\right|\end{array}\begin{array}{c}-Q_{mand}^h\le Q_{0WT}^h\le Q_{mand}^h\\ Q_{mand}^h\le Q_{1WT}^h\le Q_{WT}^{h,av} or Q_{WT}^{min}\le Q_{1WT}^h\le -Q_{mand}^h\\ Q_{WT}^{h,av}\le Q_{2WT}^h=Q_{WT}^{h,req}\le Q_{WT}^{max}\end{array}$$

(25)

where $\Delta P_{loss}^h\left(Q_{1WT}^h\right)$ and $\Delta P_{loss}^h\left(Q_{WT}^{h,av}\right)$represent the increased active power losses caused by increased reactive power production from$Q_{mand}^h$to$Q_{1WT}^h$and $Q_{WT}^{h,av}$, respectively. Similar to Equation (10), we can rewrite Equation (25) as follows:

$$COST\left(Q_{WT}^h\right)=C_{WT}\left(Q_{WT}^h\right)+{\rho }_{Adj}\eta \left|\Delta P_{WT}^{Adj,h}\right|$$

(26)

Figure 2. Reactive power capability of a WT considering the hourly wind power fluctuations.

Figure 2. Reactive power capability of a WT considering the hourly wind power fluctuations.

3. Day-Ahead Active Power Market

In the day-ahead active power market, an Initial Active Power Dispatch (IAPD) will be obtained by the MO for the forecasted load demand. This issue represents a bi-objective optimization problem in order to minimize the electrical energy costs and the gas emissions related to DERs and Disco. The generation units send their hourly selling bids, which consist of coupled quantity and price, to the MO. The total electrical energy costs generated by generation units are defined as:

$$COS{T}_{Energy}={\displaystyle \sum }_{h=1}^{Nh}{P}_{Disco}^h·{\pi }_{Disco}^h+{\displaystyle \sum }_{h=1}^{Nh}{\displaystyle \sum }_{i=1}^{NDG}{P}_{DG,i}^h·{\pi }_{DG,i}^h+{\displaystyle \sum }_{h=1}^{Nh}{\displaystyle \sum }_{i=1}^{NWT}{P}_{WT,i}^h·{\pi }_{WT,i}^h$$

(27)

One of the most important emissions from the electricity sector is CO₂, which is represented by released pollution in terms of tons per MW. In order to economically illustrate the harmful effects of the emissions from the electricity sector activities on the environment, different techniques, such as penalties through carbon taxes or cap-and-trade technique, have been adopted [23]. In this study, the penalty cost function of the CO₂ emissions related to Disco, DGs, and WTs is calculated as follows:

$$\begin{array}{c}COS{T}_{Disco}^{Emission}=\{\begin{array}{cc}{\displaystyle {\displaystyle \sum }_{h=1}^{Nh}}{P}_{Disco}^h·\left(C{O}_{2,Disco}-C{O}_{2,cap}\right)·EP{C}_{CO2}& if C{O}_{2,Disco}>C{O}_{2,cap}\\ 0& if C{O}_{2,Disco}\le C{O}_{2,cap}\end{array}\\ COS{T}_{DG,i}^{Emission}=\{\begin{array}{cc}{\displaystyle {\displaystyle \sum }_{h=1}^{Nh}}{P}_{DG,i}^h·\left(C{O}_{2,DG,i}-C{O}_{2,cap}\right)·EP{C}_{CO2}& if C{O}_{2,DG,i}>C{O}_{2,cap}\\ 0& if C{O}_{2,DG,i}\le C{O}_{2,cap}\end{array}\\ COS{T}_{WT,i}^{Emission}=\{\begin{array}{cc}{\displaystyle {\displaystyle \sum }_{h=1}^{Nh}}{P}_{WT,i}^h·\left(C{O}_{2,WT,i}-C{O}_{2,cap}\right)·EP{C}_{CO2}& if C{O}_{2,WT,i}>C{O}_{2,cap}\\ 0& if C{O}_{2,WT,i}\le C{O}_{2,cap}\end{array}\end{array}$$

(28)

Hence, the total penalty cost of CO₂ emissions produced by generation units is expressed as:

$$COS{T}_{Emissions}=COS{T}_{Disco}^{Emissions}+COS{T}_{DG}^{Emissions}+COS{T}_{WT}^{Emissions}$$

(29)

The MO runs the IAPD problem, which can be modeled by Equations (30)–(34), as a uniform price auction and determines the accepted selling bids and active power market schedule for the next day. The hourly market clearing price (MCP) is determined as the maximum selling bid price accepted for each hour.

Objective function:

$$\min Z=COS{T}_{Energy}+COS{T}_{Emissions}$$

(30)

Constraints:

$$0\le P_{Disco}^h\le PBi{d}_{Disco}^h$$

(31)

$$0\le P_{DG,i}^h\le PBi{d}_{DG,i}^h$$

(32)

$$0\le P_{WT,i}^h\le PBi{d}_{WT,i}^h$$

(33)

$$P_{Disco}^h+{\displaystyle \sum }_{i=1}^{NDG}{P}_{DG,i}^h+{\displaystyle \sum }_{i=1}^{NWT}{P}_{WT,i}^h={\displaystyle \sum }_{i=1}^{Nbus}{P}_{D,i}^h$$

(34)

Equations (31)–(34) represent the limits on generation and Equation (34) represents the constraint of demand/supply balance. The result of the initial schedule is submitted to the DSO to examine it from a technical viewpoint.

4. Proposed Daily VVC Model

To perform the proposed daily VVC model, it is required that all DERs submit their offers to the DSO according to the components of their reactive power cost function discussed in Section 2. Moreover, the results of IAPD obtained by MO are sent to the DSO. After receiving the offers, the DSO executes the reactive power dispatch using an optimization problem as follows:

4.1. Objective Function

The objective function to be minimized through the proposed algorithm is the total payment of the DSO, consisting of five parts, as follows:

4.1.1. Cost of Total Active Power Losses

In order to balance the total active power losses, changing the IAPD of generation units is inevitable. Thus, the non-negative variables corresponding to the change in the IAPD results are allocated to generation units to represent the contribution of each source to balance active power losses. The cost of total active power losses is achieved multiplying these variables by the MCP as follows:

$$f_1={\displaystyle \sum }_{h=1}^{Nh}\Delta P_{Disco}^{L,h}·MC{P}^h+{\displaystyle \sum }_{h=1}^{Nh}{\displaystyle \sum }_{i=1}^{NDG}\Delta P_{DG,i}^{L,h}·MC{P}^h+{\displaystyle \sum }_{h=1}^{Nh}{\displaystyle \sum }_{i=1}^{NWT}\Delta P_{WT,i}^{L,h}·MC{P}^h$$

(35)

4.1.2. Cost of Adjustment of the IAPD

Adjustment of the IAPD of generation units may be due to the enforcement of operation constraints or to enable a certain level of reactive power provision (lost opportunity cost). The corresponding cost is calculated by multiplication of the absolute value of generation adjustment by the respective adjustment prices as follows:

$$f_2={\displaystyle \sum }_{h=1}^{Nh}\left|\Delta P_{Disco}^{Adj,h}\right|·{\rho }_{Ad{j}_{Disco}}^h+{\displaystyle \sum }_{h=1}^{Nh}{\displaystyle \sum }_{i=1}^{NDG}\left|\Delta P_{DG,i}^{Adj,h}\right|·{\rho }_{Ad{j}_{DG,i}}^h+{\displaystyle \sum }_{h=1}^{Nh}{\displaystyle \sum }_{i=1}^{NWT}\left|\Delta P_{WT,i}^{Adj,h}\right|·{\rho }_{Ad{j}_{WT,i}}^h$$

(36)

4.1.3. Cost of Reactive Power Support from DGs and WTs Regardless of LOC

Since the cost of lost opportunity is calculated in the second part of the objective function, the rest of the cost of reactive power will be assessed as the third term of the objective function by Equation (37):

$$f_3={\displaystyle \sum }_{h=1}^{Nh}{\displaystyle \sum }_{i=1}^{NDG}{C}_{DG,i}(Q_{DG,i}^h)+{\displaystyle \sum }_{h=1}^{Nh}{\displaystyle \sum }_{i=1}^{NWT}{C}_{WT,i}(Q_{WT,i}^h)$$

(37)

4.1.4. Cost of Reactive Power Support from Disco

$$f_4={\displaystyle \sum }_{h=1}^{Nh}COST(Q_{Disco}^h)$$

(38)

4.1.5. Costs of the Switching Operations of the OLTC and Switched Sh.Cs

The fifth part of the objective function is the total depreciation cost of switchable devices. These devices are allowed to operate for a limited number of switching operations during their entire lifetime. So, devaluation in the capital cost that arises from each switching operation is represented in terms of $/switching operation [24]. Thus, the costs of the switching operations of the OLTC and switched Sh.Cs are expressed as:

$$f_5={\displaystyle \sum }_{h=1}^{Nh-1}{D}_{Tap}\left|Ta{p}^{h+1}-Ta{p}^h\right|+{\displaystyle \sum }_{h=1}^{Nh-1}{\displaystyle \sum }_{i=1}^{NSC}{D}_S\left|CSte{p}_i^{h+1}-CSte{p}_i^h\right|$$

(39)

The objective function of daily VVC problem can be expressed as follows:

$$\min f=f_1+f_2+f_3+f_4+f_5$$

(40)

4.2. Constraints

In order to achieve optimal scheduling for the daily VVC problem, network equality and inequality constraints should be satisfied. The constraints of the daily VVC problem are defined as follows:

−

Power flow constraints:

$$\begin{gathered} P_{G,i}^h-P_{D,i}^h=\left|V_i^h\right|{\displaystyle \sum }_{j=1}^{NBus}\left|V_j^h\right|\left|Y_{ij}\right|cos\left({\theta }_i^h-{\theta }_j^h-{\phi }_{ij}\right) \\ {Q}_{G,i}^h-Q_{D,i}^h=\left|V_i^h\right|{\displaystyle \sum }_{j=1}^{NBus}\left|V_j^h\right|\left|Y_{ij}\right|sin\left({\theta }_i^h-{\theta }_j^h-{\phi }_{ij}\right) \end{gathered}$$

(41)

−

Bus Voltage magnitude:

$$V_{min}\le \left|V_i^h\right|\le V_{max}$$

(42)

−

Limits of the active power adjustments of Disco, DGs and WTs:

$$\begin{gathered} -x_{Disco}^{max}·P_{Disco}^{ini,h}\le \Delta P_{Disco}^{Adj,h}\le x_{Disco}^{max}·P_{Disco}^{ini,h} \\ -x_{DG,i}^{max}·P_{DG,i}^{ini,h}\le \Delta P_{DG,i}^{Adj,h}\le x_{DG,i}^{max}·P_{DG,i}^{ini,h} \\ -x_{WT,i}^{max}·P_{WT,i}^{ini,h}\le \Delta P_{WT,i}^{Adj,h}\le x_{WT,i}^{max}·P_{WT,i}^{ini,h} \end{gathered}$$

(43)

−

Limits of generation capacity:

$$\begin{gathered} 0\le P_{Disco}^{ini,h}+\Delta P_{Disco}^{L,h}+\Delta P_{Disco}^{Adj,h}\le P_{Disco}^{max} \\ {P}_{DG,i}^{min}\le P_{DG,i}^{ini,h}+\Delta P_{DG,i}^{L,h}+\Delta P_{DG,i}^{Adj,h}\le P_{DG,i}^{max} \\ {P}_{WT,i}^{min}\le P_{WT,i}^{ini,h}+\Delta P_{WT,i}^{L,h}+\Delta P_{WT,i}^{Adj,h}\le P_{WT,i}^{max} \end{gathered}$$

(44)

−

Limit of reactive power generation of Disco:

$$-\sqrt{{\left(S_{Disco}\right)}^2-{\left(P_{Disco}^h\right)}^2}\le Q_{Disco}^h\le \sqrt{{\left(S_{Disco}\right)}^2-{\left(P_{Disco}^h\right)}^2}$$

(45)

−

Limit of transformers tap:

$$Ta{p}_{min}\le Ta{p}^h\le Ta{p}_{max}$$

(46)

−

Limit of steps of capacitors:

$$CSte{p}_i^{min}\le CSte{p}_i^h\le CSte{p}_i^{max}$$

(47)

−

Maximum permissible daily number of OLTC operations:

$${\displaystyle \sum }_{h=1}^{Nh-1}\left|Ta{p}^{h+1}-Ta{p}^h\right|\le N_{OLTC}^{max}$$

(48)

−

Maximum permissible daily number of switching operations of Sh.Cs:

$${\displaystyle \sum }_{h=1}^{Nh-1}\left|CSte{p}_i^{h+1}-CSte{p}_i^h\right|\le N_{Sh.C,i}^{max}$$

(49)

−

Limits of the binary variables related to the cost function of reactive power of DGs and WTs:

$$\begin{gathered} {\displaystyle \sum }_{k=0}^2{W}_{k,D{G}_i}^h\le 1 \\ {\displaystyle \sum }_{k=0}^2{W}_{k,W{T}_i}^h\le 1 \end{gathered}$$

(50)

5. Benders Decomposition Algorithm

The proposed daily VVC model addressed in this paper is formulated as an MINLP problem. The complexity of solving nonlinear optimization problems with integer variables and a great number of continuous variables motivates us to implement decomposition techniques, such as the Benders decomposition algorithm. The general MINLP problem has been divided by means of the Benders decomposition algorithm into two structures, master and slave, which provides an iterative procedure between both structures in order to achieve an optimal solution [25]. The proposed solution methodology is modeled in GAMS software using the CPLEX solver for solving the Mixed Integer Programming (MIP) of the master problem and the CONOPT solver for solving the Non-Linear Programming (NLP) of the slave problem.

The Benders decomposition algorithm is described as follows:

Master Problem: The master problem decides on tap setting of OLTC, steps of Sh.Cs, and values of binary variables related to the cost of reactive power provided by DGs and WTs in order to minimize the costs of switching operations. The master problem solution is transferred to the slave problem. The objective function of master problem minimizes:

$$f_{master}=f_5+\alpha$$

(51)

Subject to the constraints Equations (46)–(50) and the Benders linear cuts as:

$$\begin{gathered} \alpha \ge f_{slave}^{\left(m\right)}+{\displaystyle \sum }_{h=1}^{Nh}{\displaystyle \sum }_{i=1}^{NSC}{\lambda }_i^{h^{\left(m\right)}}\left(CSte{p}_i^h-CSte{p}_i^{h^{\left(m\right)}}\right)+{\displaystyle \sum }_{h=1}^{Nh}{\mu }^{h^{\left(m\right)}}\left(Ta{p}^h-Ta{p}^{h^{\left(m\right)}}\right) \\ +{\displaystyle \sum }_{k=0}^2{\displaystyle \sum }_{h=1}^{Nh}{\displaystyle \sum }_{i=1}^{NDG}{\gamma }_{k,D{G}_i}^{h^{\left(m\right)}}\left(W_{k,D{G}_i}^h-W_{k,D{G}_i}^{h^{\left(m\right)}}\right)+{\displaystyle \sum }_{k=0}^2{\displaystyle \sum }_{h=1}^{Nh}{\displaystyle \sum }_{i=1}^{NWT}{\gamma }_{k,W{T}_i}^{h^{\left(m\right)}}\left(W_{k,W{T}_i}^h-W_{k,W{T}_i}^{h^{\left(m\right)}}\right) \\ m=1,2, \dots ,v-1 \end{gathered}$$

(52)

where v is the iteration counter. The only real variable $\alpha$ in Equations (51) and (52), which contain the infeasibility costs, is an underestimation of the slave problem costs. The Benders linear cuts, which are updated for every iteration, join the master and slave problems. An additional cut is added to the master problem at each iteration with information about the objective value of slave problem in the previous iteration and the dual variables associated with the decision variables fixed by the master problem in the previous iteration. This information helps the master problem to make a new decision and reach the optimal solution.

Because of the absolute (ABS) function included in the master problem, this model expresses a NLP model with discontinuous derivatives (DNLP). The only reliable way to solve a DNLP model is to reformulate it as an equivalent smooth model. The standard reformulation approach for the ABS function is to replace the ABS function with the auxiliary positive variables, ${\text{g}}^{+}$ and${\text{ g}}^{-}$, as follows:

$$\left|\text{g}\left(\text{X}\right)\right|={\text{g}}^{+}+{\text{ g}}^{-}$$

(53)

providing that:

$$\text{g}\left(\text{X}\right)={\text{g}}^{+}-{\text{ g}}^{-}$$

(54)

Therefore, the discontinuous derivative from the ABS function has disappeared and the part of the model shown here is smooth. So, the master problem is formulated as a MIP problem.

Slave problem: The slave problem formulation is nearly similar to the main problem in that all the integer variables are fixed to given the value obtained by the master problem. However, there could be some cases where the master problem’s solution makes the NLP slave problem infeasible. To avoid these cases at each iteration, artificial variables are added to some constraints and embedded in the objective function so that the objective function minimizes the technical infeasibilities of operation [26]. Therefore, the slave problem not only verifies the technical feasibility of the master problem solution but also gives the optimal dispatches of generation units. At the last iteration, the final solution of the problem has to be feasible and optimal, that is all of these artificial variables should be equal to zero. The slave problem is formulated below:

$$f_{slave}=f_1+f_2+f_3+f_4+M{\displaystyle \sum }_{h=1}^{Nh}{\displaystyle \sum }_{i\in L}\left(u_i^h+w\right)$$

(55)

It is subject to the constraints of Equations (42)–(45) as well as the following constraints:

$$P_{G,i}^h-P_{D,i}^h=\left|V_i^h\right|{\displaystyle \sum }_{j=1}^{NBus}\left|V_j^h\right|\left|Y_{ij}\right|cos\left({\theta }_i^h-{\theta }_j^h-{\phi }_{ij}\right)\left(Q_{G,i}^h+q_i^h\right)-Q_{D,i}^h=\left|V_i^h\right|{\displaystyle \sum }_{j=1}^{NBus}\left|V_j^h\right|\left|Y_{ij}\right|\sin \left({\theta }_i^h-{\theta }_j^h-{\phi }_{ij}\right)$$

(56)

$$q_i^h=\{\begin{array}{cc}{u}_i^h-w& i\in L\\ 0& otherwise\end{array}$$

(57)

$$CSte{p}_i^h=CSte{p}_i^{h^{\left(v\right)}} : {\lambda }_i^{h^{\left(v\right)}}$$

(58)

$$Ta{p}^h=Ta{p}^{h^{\left(v\right)}} : {\mu }^{h^{\left(v\right)}}$$

(59)

$$W_{k,D{G}_i}^h=W_{k,D{G}_i}^{h^{\left(v\right)}} : {\gamma }_{k,D{G}_i}^{h^{\left(v\right)}}$$

(60)

$$W_{k,W{T}_i}^h=W_{k,W{T}_i}^{h^{\left(v\right)}} : {\gamma }_{k,W{T}_i}^{h^{\left(v\right)}}$$

(61)

where $u_i^h$ and $w$ are the positive artificial variables of optimization problem and M is a large enough positive constant. The constraints of Equations (58)–(61) demonstrate the dual variables (sensitivities) associated with the discrete variables specified previously by the master problem. These dual variables and the objective value computed by the slave problem are applied to create new Benders cuts for the subsequent iteration.

Figure 3. Flowchart of the optimization problem based on the Benders decomposition algorithm.

Figure 3. Flowchart of the optimization problem based on the Benders decomposition algorithm.

Stopping Criterion

The Benders decomposition algorithm ends when (a) the solution created by the master problem is feasible; and (b) the slave problem costs computed through the slave problem (upper bound) and the variable $\alpha$ (lower bound) computed through the master problem are close enough.

The procedure of the Benders decomposition algorithm for the proposed two stage model of daily VVC is shown in Figure 3.

According to the flowchart of the Benders decomposition algorithm, the following steps are implemented:

Step 1: Initialization. Initialize the iteration counter, ν = 1. Solve the initial MIP master problem with objective function Equation (51) and subject to Equations (46)–(50). Its solution provides $CSte{p}_i^{h^{\left(v\right)}}$, $Ta{p}^{h^{\left(v\right)}}$, $W_{k,D{G}_i}^{h^{\left(v\right)}}$, $W_{k,W{T}_i}^{h^{\left(v\right)}}$, and ${\alpha }^{\left(v\right)}$. In this step, Equation (52) is not considered for initialization.

Step 2: Slave problem solution. Solve the NLP slave problem for all hours. The solution of this problem is the optimal dispatches of generation units, i.e. $Q_{DG,i}^h$, $Q_{WT,i}^h$, $Q_{Disco}^h$, $\Delta P_{Disco}^{L,h}$, $\Delta P_{DG,i}^{L,h}$, $\Delta P_{WT,i}^{L,h}$, $\Delta P_{Disco}^{Adj,h}$, $\Delta P_{DG,i}^{Adj,h}$, and $\Delta P_{WT,i}^{Adj,h}$, with dual variable values ${\lambda }_i^{h^{\left(v\right)}}$, ${\mu }^{h^{\left(v\right)}}$, ${\gamma }_{k,D{G}_i}^{h^{\left(v\right)}}$, and ${\gamma }_{k,W{T}_i}^{h^{\left(v\right)}}$.

Step 3: Convergence checking. If the convergence criterion is satisfied, the optimal solution is the last values attained for variables. Otherwise, the algorithm continues with the next step.

Step 4: Master problem solution. Update the iteration counter, ν = v + 1. Solve the MIP master problem with objective function Equation (51) and subject to Equations (46)–(50) and (52).

The algorithm continues with Step 2.

6. Simulation Results

The 22 bus 20-kV radial distribution test system [17] is modified and used for this study. The single line diagram of the test distribution network is shown in Figure 4. The system data are given in the Appendix. The HV/MV tap changing transformer has 11 tap positions and each tap ratio is 0.01 p.u. The maximum permissible daily number of operations of the transformer is 20. Voltage limits are considered to be ±5% of nominal voltage. The capacitor banks 1 and 2 have the capacity of 1000 kVAr with five switching steps of 200 kVAr. Three DERs including two DGs and a WT are connected at the distribution buses, as shown in Figure 4. General parameters such as the capacity, location, and CO₂ emissions of generation units have been reported in Table 1. The values of $C{O}_{2,cap}$ and $EP{C}_{CO2}$ are assumed to be 0.1 (ton/MWh) and 40 ($/ton), respectively [27].

Figure 4. The single line diagram of a 22-bus distribution test system.

Figure 4. The single line diagram of a 22-bus distribution test system.

Table 1. Capacity, location, and emissions of generation units.

neration Unit	Capacity	Location	CO2 (ton/MWh)
Disco	20 MVA	-	0.927
DG1	400 kW	Bus 11 in feeder 2	0.489
DG2	1000 kW	Bus 9 in feeder 1	0.724
WT	500 kW	Bus 8 in feeder 1	0

Table 1. Capacity, location, and emissions of generation units.

neration Unit	Capacity	Location	CO2 (ton/MWh)
Disco	20 MVA	-	0.927
DG1	400 kW	Bus 11 in feeder 2	0.489
DG2	1000 kW	Bus 9 in feeder 1	0.724
WT	500 kW	Bus 8 in feeder 1	0

Table 2 shows the generation unit selling bids, including the blocks of bid power and the generation bid prices.

Table 2. Selling bids offered by generation units.

Generation Unit	Hour	Block Number	Quantity (kW)	Price ($/kWh)
DG1	1–24	1	200	0.045
		2	100	0.05
		3	100	0.06
DG2	1–24	1	400	0.042
		2	400	0.05
		3	200	0.06
WT	1–6	1	100	0.041
	7–10	1	366	0.041
	11–15	1	442	0.041
	16–18	1	500	0.041
	19–22	1	480	0.041
	23–24	1	100	0.041
Disco	1–8	1	8000	0.04
	9–10	1	8000	0.045
	11–13	1	8000	0.05
	14–15	1	8000	0.052
	16	1	8000	0.053
	17	1	8000	0.054
	18–19	1	8000	0.053
	20–22	1	8000	0.052
	23–24	1	8000	0.05

Table 2. Selling bids offered by generation units.

Generation Unit	Hour	Block Number	Quantity (kW)	Price ($/kWh)
DG1	1–24	1	200	0.045
		2	100	0.05
		3	100	0.06
DG2	1–24	1	400	0.042
		2	400	0.05
		3	200	0.06
WT	1–6	1	100	0.041
	7–10	1	366	0.041
	11–15	1	442	0.041
	16–18	1	500	0.041
	19–22	1	480	0.041
	23–24	1	100	0.041
Disco	1–8	1	8000	0.04
	9–10	1	8000	0.045
	11–13	1	8000	0.05
	14–15	1	8000	0.052
	16	1	8000	0.053
	17	1	8000	0.054
	18–19	1	8000	0.053
	20–22	1	8000	0.052
	23–24	1	8000	0.05

6.1. Active Power Market Schedule Obtained by the MO

The linear programming of IAPD is programmed in the GAMS and solved by the solver CPLEX. The output of this program provides the hourly MCP and the accepted bid power of each generation unit for the next day, as indicated in Table 3. The schedule of each generation unit is sent to the DSO for verifying technical validation.

6.2. Daily Optimal Dispatches of VVC Devices

Table 4 and Table 5 present the data of DERs, including information characterizing the capability diagrams of the generator, components of offered prices of reactive power, the generator adjustment price, and the maximum admitted change as a percentage of its IAPD. The cost of reactive power generated by Disco, corresponding adjustment price, and maximum admitted change in its IAPD are 0.022 ($/kVArh), 0.09 ($/kWh), and 40%, respectively. The mandatory power factor ($p{f}_{mand }$) has been taken to be 0.95. The maximum variability of the actual hourly wind power from the forecasted value is considered to be 10% of the forecasted value.

Table 3. Initial active power schedule.

Hour	$P_{Disco}^{ini}$ (kW)	$P_{DG1}^{ini}$ (kW)	$P_{DG2}^{ini}$ (kW)	$P_{WT}^{ini}$ (kW)	$MCP$ ($/kWh)	Hour	$P_{Disco}^{ini}$ (kW)	$P_{DG1}^{ini}$ (kW)	$P_{DG2}^{ini}$ (kW)	$P_{WT}^{ini}$ (kW)	$MCP$ ($/kWh)
1	422.62	300	400	100	0.05	13	4264.6	400	800	442	0.06
2	76.12	300	400	100	0.05	14	2056.3	400	1000	442	0.06
3	468.89	300	400	100	0.05	15	260.66	400	1000	442	0.06
4	626.39	300	400	100	0.05	16	187.89	400	1000	500	0.06
5	659.85	300	400	100	0.05	17	2008.1	400	1000	500	0.06
6	691.35	300	400	100	0.05	18	5069.76	400	1000	500	0.06
7	488.35	300	400	366	0.05	19	5897	400	1000	480	0.06
8	801.39	300	400	366	0.05	20	5756.23	400	1000	480	0.06
9	1375.4	400	800	366	0.06	21	4646.8	400	1000	480	0.06
10	4003.9	400	800	366	0.06	22	2925.01	400	1000	480	0.06
11	4651.49	400	800	442	0.06	23	598.89	400	800	100	0.06
12	4758.8	400	800	442	0.06	24	80.12	400	800	100	0.06

Table 3. Initial active power schedule.

Hour	$P_{Disco}^{ini}$ (kW)	$P_{DG1}^{ini}$ (kW)	$P_{DG2}^{ini}$ (kW)	$P_{WT}^{ini}$ (kW)	$MCP$ ($/kWh)	Hour	$P_{Disco}^{ini}$ (kW)	$P_{DG1}^{ini}$ (kW)	$P_{DG2}^{ini}$ (kW)	$P_{WT}^{ini}$ (kW)	$MCP$ ($/kWh)
1	422.62	300	400	100	0.05	13	4264.6	400	800	442	0.06
2	76.12	300	400	100	0.05	14	2056.3	400	1000	442	0.06
3	468.89	300	400	100	0.05	15	260.66	400	1000	442	0.06
4	626.39	300	400	100	0.05	16	187.89	400	1000	500	0.06
5	659.85	300	400	100	0.05	17	2008.1	400	1000	500	0.06
6	691.35	300	400	100	0.05	18	5069.76	400	1000	500	0.06
7	488.35	300	400	366	0.05	19	5897	400	1000	480	0.06
8	801.39	300	400	366	0.05	20	5756.23	400	1000	480	0.06
9	1375.4	400	800	366	0.06	21	4646.8	400	1000	480	0.06
10	4003.9	400	800	366	0.06	22	2925.01	400	1000	480	0.06
11	4651.49	400	800	442	0.06	23	598.89	400	800	100	0.06
12	4758.8	400	800	442	0.06	24	80.12	400	800	100	0.06

Table 4. Characteristics of the DGs.

R	$P_{DG}^{max}$ (kW)	$Q_{DG}^{max}$ (kVAr)	$Q_{DG}^{min}$ (kVAr)	$Q_{DG}^x$ (kVAr)	${\mathsf{\rho }}_{0 }$ ($)	${\mathsf{\rho }}_1$ ($/MVArh)	${\mathsf{\rho }}_2$ ($/MVArh)	${\mathsf{\rho }}_{Adj}$ ($/kWh)	$x^{max}$
DG1	400	400	−400	400	0.6	23	23	0.08	50%
DG2	1000	1000	−600	800	0.8	30	30	0.075	50%

Table 4. Characteristics of the DGs.

R	$P_{DG}^{max}$ (kW)	$Q_{DG}^{max}$ (kVAr)	$Q_{DG}^{min}$ (kVAr)	$Q_{DG}^x$ (kVAr)	${\mathsf{\rho }}_{0 }$ ($)	${\mathsf{\rho }}_1$ ($/MVArh)	${\mathsf{\rho }}_2$ ($/MVArh)	${\mathsf{\rho }}_{Adj}$ ($/kWh)	$x^{max}$
DG1	400	400	−400	400	0.6	23	23	0.08	50%
DG2	1000	1000	−600	800	0.8	30	30	0.075	50%

Table 5. Characteristics of the WT.

DER	$V_{c,max}$ (p.u.)	$I_{c,max}$ (p.u.)	X (p.u.)	$Q_{WT}^{min}$ (kVAr)	$a_{0 }$ ($)	${\mathsf{\rho }}_{Adj}$ ($/kWh)	$a_l$	$b_l$	$x^{max}$
WT	1.4	1.24	0.2	−250	0.9	0.095	10.14	0.003	50%

Table 5. Characteristics of the WT.

DER	$V_{c,max}$ (p.u.)	$I_{c,max}$ (p.u.)	X (p.u.)	$Q_{WT}^{min}$ (kVAr)	$a_{0 }$ ($)	${\mathsf{\rho }}_{Adj}$ ($/kWh)	$a_l$	$b_l$	$x^{max}$
WT	1.4	1.24	0.2	−250	0.9	0.095	10.14	0.003	50%

In order to calculate the depreciation cost of OLTC, the installation cost of an OLTC is assumed to be $\text{\$}400,000$ $\left(=\text{\$ }20,000/\text{MVA}\times 20 \text{MVA}\right),$ as reported in [26]. The total number of acceptable switching operations of each OLTC can be 143,080 times (=20 steps/day × 365 days/year × 20 years × 0.98 availability factor). Thus, the depreciation cost for each step change $(D_{Tap})$ is $2.79. Also, the installation cost for each Sh.C is $11,600 $\left(=\text{\$ }11,600/\text{MVA}\times 1 \text{MVA}\right)$ as given in [26]. Each Sh.C can be operated 71,540 times (=10 switching operations/day × 365 days/year × 20 years × 0.98 availability factor). Therefore, the depreciation cost for each switching operation of Sh.C $(D_S)$is equal to $0.162.

The daily VVC problem proposed in this paper will be tested on three different cases, as follows:

6.2.1. Case 1 (Base Case) and Implementing the Benders Decomposition Algorithm

By implementing the Benders decomposition algorithm, Figure 5 shows the hourly optimal dispatch results of the Sh.Cs for case 1. The number of switching operations for C1 and C2 is 9 and 7, respectively, which is less than the maximum allowable daily number of operations (10). Moreover, because of the high depreciation cost of OLTC, the tap position of OLTC is fixed at 1.04 pu during the whole day. Table 6 shows the optimal dispatch of generation units. In this case, the change in the IAPD corresponding to operation or security enforcements is not necessary. So, all $\Delta P^{Adj}$ variables are zero. Regarding the $\Delta P^L$ variables, these variables are zero for DG1 and WT for all hours. Hence, the power losses of network are balanced in the two buses using Disco and DG2. The voltage profiles at the buses of feeders 1 and 2 are illustrated in Figure 6. In this figure, the voltages of buses are brought back to the acceptable range of 0.95 to 1.05 pu for 24 h.

Figure 5. Daily optimal dispatches of the Sh.Cs for case 1.

Figure 5. Daily optimal dispatches of the Sh.Cs for case 1.

Figure 6. Voltage profile of feeders 1 and 2 after Volt/Var control for case 1.

Figure 6. Voltage profile of feeders 1 and 2 after Volt/Var control for case 1.

6.2.2. Case 2: Considering Functions ${\text{f}}_1$, ${\text{f}}_2$, ${\text{f}}_3{\text{ and f}}_4$. and Implementing the DICOPT Solver

In order to validate the results obtained by the proposed approach, the same problem has been tested by relaxing all constraints of the daily number of switching operations. In this case, since the optimal dispatch results of OLTC and Sh.Cs for each hour are not correlated with the solutions of the other hours, the daily VVC problem can be solved separately for each hour using the DICOPT solver. The daily optimal dispatch results of OLTC and Sh.Cs corresponding to the objectives of $f_1$ to $f_4$ are illustrated in Figure 7 and Figure 8, respectively. As seen in Figure 8, the number of operations of C1 and C2 are 12 and 16, respectively, which is higher than the acceptable maximum. Table 7 provides the optimal reactive power dispatch of generation units in case 2. In order to confirm the effectiveness of the proposed method based on the Benders decomposition algorithm, Table 8 provides a comparison between the results of cases 1 and 2. From the comparison of results shown in Table 8, it has been observed that despite the limited switching operations included in the Benders decomposition-based proposed method, the total cost decreases to $449.26, compared to $477.761 when no restriction is imposed on the switched devices; this confirms the effectiveness of Benders decomposition. On the other hand, when the limitation on switching operations of control devices is not considered, the active power loss has decreased from 998.134 kW to 954.506 kW. This is due to the fact that an unlimited number of switching operations provides a more flexible control of the power flow to regulate the voltage within its admissible range and to fulfill the operational constraints.

Table 6. Daily optimal dispatch of generation units for case 1.

Hour	Disco			DG1			DG2			WT
	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVAr)	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVAr)	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVAr)	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVAr)
1	2.31	0	0	0	0	52.64	0	0	−23.82	0	0	−32.87
2	2.93	0	0	0	0	−53.89	0	0	−131.48	0	0	−32.87
3	2.30	0	0	0	0	98.61	0	0	−35.10	0	0	−32.87
4	2.88	0	0	0	0	−35.58	0	0	131.48	0	0	32.87
5	0	0	0	0	0	59.91	2.72	0	44.53	0	0	32.87
6	3.07	0	0	0	0	−7.21	0	0	131.48	0	0	32.87
7	0	0	0	0	0	12.60	3.77	0	63.91	0	0	120.30
8	0	0	0	0	0	98.61	4.38	0	132.92	0	0	171.89
9	10.51	0	0	0	0	−30.40	0	0	262.96	0	0	120.30
10	0	0	1043.76	0	0	131.48	63.38	0	283.79	0	0	471.50
11	0	0	1432.18	0	0	395.42	89.31	0	292.32	0	0	384.72
12	0	0	1515.26	0	0	400.00	94.75	0	294.10	0	0	384.72
13	0	0	1223.6	0	0	313.03	73.41	0	287.09	0	0	384.72
14	17.01	0	0	0	0	131.48	0	0	265.07	0	0	275.70
15	12.31	0	0	0	0	131.48	0	0	−223.85	0	0	−145.29
16	12.99	0	0	0	0	131.48	0	0	−219.35	0
17	16.99	0	0	0	0	131.48	0	0	328.70	0	0	157.86
18	111.50	0	1854.32	0	0	400.00	0	0	328.70	0	0	366.61
19	161.90	0	2472.03	0	0	400.00	0	0	328.70	0	0	366.61
20	162.84	0	2648.13	0	0	131.48	0	0	328.70	0	0	366.61
21	99.95	0	1850.04	0	0	131.48	0	0	328.70	0	0	366.61
22	33.50	0	560.77	0	0	131.48	0	0	328.70	0	0	366.61
23	0	0	0	0	0	131.48	6.50	0	265.10	0	0	28.26
24	6.94	0	0	0	0	−28.98	0	0	93.81	0	0	32.87

Table 6. Daily optimal dispatch of generation units for case 1.

Hour	Disco			DG1			DG2			WT
	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVAr)	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVAr)	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVAr)	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVAr)
1	2.31	0	0	0	0	52.64	0	0	−23.82	0	0	−32.87
2	2.93	0	0	0	0	−53.89	0	0	−131.48	0	0	−32.87
3	2.30	0	0	0	0	98.61	0	0	−35.10	0	0	−32.87
4	2.88	0	0	0	0	−35.58	0	0	131.48	0	0	32.87
5	0	0	0	0	0	59.91	2.72	0	44.53	0	0	32.87
6	3.07	0	0	0	0	−7.21	0	0	131.48	0	0	32.87
7	0	0	0	0	0	12.60	3.77	0	63.91	0	0	120.30
8	0	0	0	0	0	98.61	4.38	0	132.92	0	0	171.89
9	10.51	0	0	0	0	−30.40	0	0	262.96	0	0	120.30
10	0	0	1043.76	0	0	131.48	63.38	0	283.79	0	0	471.50
11	0	0	1432.18	0	0	395.42	89.31	0	292.32	0	0	384.72
12	0	0	1515.26	0	0	400.00	94.75	0	294.10	0	0	384.72
13	0	0	1223.6	0	0	313.03	73.41	0	287.09	0	0	384.72
14	17.01	0	0	0	0	131.48	0	0	265.07	0	0	275.70
15	12.31	0	0	0	0	131.48	0	0	−223.85	0	0	−145.29
16	12.99	0	0	0	0	131.48	0	0	−219.35	0
17	16.99	0	0	0	0	131.48	0	0	328.70	0	0	157.86
18	111.50	0	1854.32	0	0	400.00	0	0	328.70	0	0	366.61
19	161.90	0	2472.03	0	0	400.00	0	0	328.70	0	0	366.61
20	162.84	0	2648.13	0	0	131.48	0	0	328.70	0	0	366.61
21	99.95	0	1850.04	0	0	131.48	0	0	328.70	0	0	366.61
22	33.50	0	560.77	0	0	131.48	0	0	328.70	0	0	366.61
23	0	0	0	0	0	131.48	6.50	0	265.10	0	0	28.26
24	6.94	0	0	0	0	−28.98	0	0	93.81	0	0	32.87

Figure 7. Daily optimal tap positions of OLTC for case 2.

Figure 7. Daily optimal tap positions of OLTC for case 2.

Figure 8. Daily optimal dispatches of the Sh.Cs for case 2.

Figure 8. Daily optimal dispatches of the Sh.Cs for case 2.

Table 7. Daily optimal reactive power dispatches of generation units for case 2.

Hour	Disco (kVar)	DG1 (kVar)	DG2 (kVar)	WT (kVar)	Hour	Disco (kVar)	DG1 (kVar)	DG2 (kVar)	WT (kVar)
1	0	87.86	41.02	66.65	13	1228.46	307.42	286.66	384.72
2	0	65.16	77.62	37.94	14	0	131.48	264.96	275.53
3	0	93.53	57.34	79.36	15	0	131.48	112.66	115.30
4	0	98.61	105.25	124.29	16	0	131.48	166.54	246.82
5	0	98.57	110.98	127.49	17	0	131.48	318.69	366.61
6	0	98.61	120.54	137.52	18	1852.32	400	328.70	366.61
7	0	98.61	132.63	165.34	19	2469.08	400	328.70	366.61
8	0	98.61	95.55	209.17	20	2368.11	400	328.70	366.61
9	0	131.48	190.22	229.27	21	1848.24	131.48	328.71	366.61
10	1013.82	160.49	283.26	471.50	22	560.19	131.48	328.70	366.61
11	1436.91	389.77	291.78	384.72	23	0	131.48	109.34	183.21
12	1514.24	400.00	293.51	384.72	24	0	131.48	235.04	130.39

Table 7. Daily optimal reactive power dispatches of generation units for case 2.

Hour	Disco (kVar)	DG1 (kVar)	DG2 (kVar)	WT (kVar)	Hour	Disco (kVar)	DG1 (kVar)	DG2 (kVar)	WT (kVar)
1	0	87.86	41.02	66.65	13	1228.46	307.42	286.66	384.72
2	0	65.16	77.62	37.94	14	0	131.48	264.96	275.53
3	0	93.53	57.34	79.36	15	0	131.48	112.66	115.30
4	0	98.61	105.25	124.29	16	0	131.48	166.54	246.82
5	0	98.57	110.98	127.49	17	0	131.48	318.69	366.61
6	0	98.61	120.54	137.52	18	1852.32	400	328.70	366.61
7	0	98.61	132.63	165.34	19	2469.08	400	328.70	366.61
8	0	98.61	95.55	209.17	20	2368.11	400	328.70	366.61
9	0	131.48	190.22	229.27	21	1848.24	131.48	328.71	366.61
10	1013.82	160.49	283.26	471.50	22	560.19	131.48	328.70	366.61
11	1436.91	389.77	291.78	384.72	23	0	131.48	109.34	183.21
12	1514.24	400.00	293.51	384.72	24	0	131.48	235.04	130.39

Table 8. Comparison of optimal results in cases 1 and 2.

Different Parts of Objective Function	Case 1	Case 2
	Benders Decomposition	DICOPT Solver
Number of switching operations of OLTC	0	4
Number of switching operations of C1	9	12
Number of switching operations of C2	7	16
Total active power losses (kWh)	998.134	954.506
Depreciation cost of OLTC ($)	0	11.16
Depreciation cost of Sh.Cs ($)	2.592	4.536
Total cost of active power losses ($)	59.64	57.07
Total cost of reactive power ($)	387.03	404.995
Total cost ($)	449.26	477.761

Table 8. Comparison of optimal results in cases 1 and 2.

Different Parts of Objective Function	Case 1	Case 2
	Benders Decomposition	DICOPT Solver
Number of switching operations of OLTC	0	4
Number of switching operations of C1	9	12
Number of switching operations of C2	7	16
Total active power losses (kWh)	998.134	954.506
Depreciation cost of OLTC ($)	0	11.16
Depreciation cost of Sh.Cs ($)	2.592	4.536
Total cost of active power losses ($)	59.64	57.07
Total cost of reactive power ($)	387.03	404.995
Total cost ($)	449.26	477.761

6.2.3. Case 3: Stressed Operation for Peak Hours and Implementing the Benders Decomposition Algorithm

In this case, due to the operational issues, the maximum reactive power of Disco was limited to peak hours (18–21) to get a more stressed operation situation. The maximum reactive power of Disco was decreased by 1300 kVAr, 1920 kVAr, 1860 kVAr, and 1050 kVAr at hours 18, 19, 20, and 21, respectively. Table 9 presents the results obtained for reactive power dispatches of generation units for peak hours. As indicated in Table 9, the $\Delta P^{Adj}$ variable is not zero for Disco and DG2. Also, the reactive power of DG1 is increased by 400 kVar (maximum reactive power capability) at hours 20 and 21. Furthermore, since the active and reactive power values are coupled via the capability diagram, a reduction in the active power generation of DG2 takes place to adjust for the reactive power requirements of the system. The other results of daily optimal dispatches of generation units, OLTC, and Sh.Cs are not changed. Concerning the objective function, its final value is $714.99. This value is much greater than that of case 1 because it comprises the adjustment costs of the initially scheduled active powers.

Table 9. Daily optimal reactive power dispatches of generation units for peak hours for case 3.

Hour	Disco			DG1			DG2			WT
	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVar)	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVar)	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVar)	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVar)
18	125.73	478.79	1300	0	0	400	0	−478.79	895.76	0	0	366.61
19	174.714	461.01	1920	0	0	400	0	−461.01	892.20	0	0	366.61
20	147.39	172.38	1860	0	0	400	0	−172.38	834.48	0	0	366.61
21	96.68	286.48	1050	0	0	400	0	−286.48	857.3	0	0	366.61

Table 9. Daily optimal reactive power dispatches of generation units for peak hours for case 3.

Hour	Disco			DG1			DG2			WT
	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVar)	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVar)	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVar)	$\Delta P^L$ (kW)	$\Delta P^{Adj}$ (kW)	Q (kVar)
18	125.73	478.79	1300	0	0	400	0	−478.79	895.76	0	0	366.61
19	174.714	461.01	1920	0	0	400	0	−461.01	892.20	0	0	366.61
20	147.39	172.38	1860	0	0	400	0	−172.38	834.48	0	0	366.61
21	96.68	286.48	1050	0	0	400	0	−286.48	857.3	0	0	366.61

7. Conclusions

This paper presented a new approach based on the energy market and reactive power dispatch considering the reactive power cost of DERs in VVC in distribution networks. For this purpose, a new pricing framework was presented to determine the cost of reactive power produced by DERs including synchronous machine-based DGs and WTs. In the proposed method, the initial scheduled active powers of Disco and DERs were determined in the day-ahead active power market. Then, the results were surrendered to the DSO, who performed daily VVC testing based on the proposed active power dispatch in order to determine the optimal dispatches of VVC devices. Due to the coupling between active and reactive power of DERs, in the proposed model, this issue was addressed by considering the capability diagram of DERs defined in the PQ plane. A Benders decomposition technique was applied to cope with the complexity of solving the daily VVC problem as a MINLP problem. It was concluded from the research results that incorporating the cost of reactive power support of generation units in the developed model would encourage the DERs to actively contribute to the reactive power support as an ancillary service.