A Kriging Model Based Optimization of Active Distribution Networks Considering Loss Reduction and Voltage Profile Improvement

Abstract

Optimal operation of the active distribution networks (ADN) is essential to keep its safety, reliability and economy. With the integration of multiple controllable resources, the distribution networks are facing more challenges in which the optimization strategy is the key. This paper establishes the optimal operation model of the ADN considering a diversity of controllable resources including energy storage devices, distributed generators, voltage regulators and switchable capacitor banks. The objective functions contain reducing the power losses and improving the voltage profiles. To solve the optimization problem, the Kriging model based Improved Surrogate Optimization-Mixed-Integer (ISO-MI) algorithm is proposed in this paper. The Kriging model is applied to approximate the complicated distribution networks, which speeds up the solving process. Finally, the accuracy of the Kriging model is validated and the efficiency among the proposed method, genetic algorithm (GA) and particle swarm optimization (PSO) is compared in an unbalanced IEEE-123 nodes test feeder. The results demonstrate that the proposed method has better performance than GA and PSO.

1. Introduction

Nowadays, the increasingly prominent environmental pollution contributes the development of renewable energy resources, such as wind and photovoltaic. However, it is known that the renewable energy sources are intermittent and volatile. Hence, how to effectively deal with the power and voltage fluctuation caused by large scale integration of distributed generations (DGs) is a very difficult task [1,2]. Although in recent years, with the development of distributed generation technologies, energy storage technologies and power electronics, the problem of DGs accessing to the grid has been resolved to a certain degree [3], the lack of efficient optimization methods, low degree of automation as well as lack of the participation of demand side all limit the further development of the clean and renewable energy resources [4]. Under these circumstances, active distribution network (ADN) technology came into being, the core of which is based on advanced information and communication, power electronics and automation technology, making full use of controllable resources (distributed generation unit, energy storage, controllable load, etc.). Through coordinated control of “source-network-load”, ADN can achieve the goal of large-scale renewable energy access, thus improving the distribution network operation economy, ensuring the quality of electricity users and power supply reliability [5,6,7].

Optimization strategy of ADN is the key, and a hot spot of active distribution network related technology research [8,9]. The traditional optimal operation and schedule of distribution network is often seen as a Volt/Var control problem with the goal of minimizing the operation cost or reducing the power losses. The controllable resources include the on-load tape changer (OLTC), voltage regulator and switchable capacitor banks [10,11,12]. In recent years, with the rapid development of information and communication technology (ICT) as well as power electronics, more controllable resources can be utilized to serve the optimization problem, which brings opportunities for the development of optimization algorithm. For instance, in [10], considering the penetration of electric vehicles (EVs), Advanced Metering Infrastructure (AMI)-based quasi real-time Volt-VAR Optimization (VVO) is introduced, aimed at minimizing the grid loss and Volt-VAR control assets operating costs. In [13], considering operation cost and pollutant treatment cost, an improved particle swarm optimization (PSO) algorithm combined with Monte Carlo simulation is used to solve the objective function which is maximizing the comprehensive benefits. In [14,15], considering switchable capacitor banks, DGs and energy storages, the NSGA-II multi-objective optimization algorithm is used to minimize the power losses, the electricity generation cost and carbon emissions. In [16], the day-ahead scheduling of distribution network is simulated, which considers different kinds of resources including DGs and demand response, etc. In [17], a tractable min–max–min cost model is proposed to find a robust optimal day-ahead scheduling of ADN considering demand response as one of the important resources.

However, the diversity of controllable resources also deepens the difficulties of the optimization problem. In ADN, the control variables include both continuous variables (DGs’ output, storage charge, discharge power, etc.), and discrete variables (regulator tap positions, switching status of switchable capacitor banks, etc.). In terms of constraints, not only linear constraints (upper and lower bounds of power output, etc.), but also nonlinear constraints (power flow equality constraints, node voltage constraints, etc.) are considered. Therefore, the optimal operation of the ADN is a complex mixed integer nonlinear programming (MINP) problem. Because the solution based on the traditional interior point method is not ideal, the intelligent algorithms with wider adaptability are widely used [18,19,20,21]. In the literature [18], the optimal scheduling of OLTC and capacitor banks is achieved by using genetic algorithm (GA). In [19], the optimal scheduling model of ADN is established, and the intelligent single particle optimization algorithm (IPSO) is used to reduce the running cost of the system. In [20], the hybridized genetic algorithm and the ant colony algorithm (ACO) are used to realize the economic operation of the distribution system. In [22], a new optimization framework is presented to optimize the bidding strategy of a distribution company in a day-ahead energy market and Benders decomposition technique (BDT) is employed to simplify the optimization procedure. A receding horizon optimization strategy is presented in [23] to solve the objective functions of meeting the electricity demand while minimizing the overall operating and environmental costs considering both generation side and demand side. In [24], the PSO combined with fuzzification technique is applied as Volt-VAR Optimization (VVO) algorithm, aiming at minimizing distribution network loss and the operating cost of reactive power injection, improving voltage profile of the system. A chaotic improved honey bee mating optimization is proposed to solve the optimal dispatch of ADN in [25].

Although the intelligent algorithm has been widely researched and applied in the optimization of ADN, there are still the following problems: (1) the lack of effective constraint processing mechanism, resulting in reduced efficiency of the solution; and (2) the speed of the solution is reduced due to the large number of data needed for intelligent algorithms and the simulation of the power flow, especially for three-phase modeling of the complex distribution network. Recently, using surrogate based optimization techniques to solve optimization problems with computationally black-box objective functions have attracted attention of researchers, considering the interrupt load, the economic operation model of ADN is established with the goal of minimizing the operation cost. A modified fuzzy adaptive PSO assisted by Kriging model (KMA-MFAPSO) is developed to solve the optimization problem. Based on the above research background, considering the stability and reliability of ADN, this paper establishes the optimal operation model of AND with two objectives including minimizing the power losses and improving the voltage profile. The resources contain distributed power generation units, energy storage equipment, voltage regulators and switchable capacitor banks. A hybrid algorithm based on Kriging model named ISO-MI (Improved Surrogate Optimization-Mixed-Integer) is proposed. Finally, the validity of the Kriging model and the efficiency of the proposed algorithm are verified by comparison and analysis of the traditional algorithm.

The main contributions of this paper are as follows.

(1)

The optimal operation and schedule model of ADN is proposed considering multiple controllable resources such as battery storage, DGs, etc. The objectives include reducing the power losses and improving the voltage profile.

(2)

The Kriging model is used to approximate the complex active distribution network, speeding up the solving process.

(3)

The Kriging model based optimization method named ISO-MI is proposed to solve the optimization problem, which improves the solving efficiency.

This paper is organized as follows. Section 2 formulates the optimization problem of distribution system. The principle of Kriging model and the basis of the proposed algorithm ISO-MI are discussed in Section 3. The simulation results of the proposed method are included in Section 4. Concluding remarks are given in Section 5.

2. Problem Formulation

This paper focuses on the optimal scheduling and operation of ADN. The diagram of optimal dispatch for ADN is shown in Figure 1. As can be observed, there is a diversity of controllable resources in both supply and demand side in ADN, such as regulators, battery storages, DGs, and thermostat controlled loads (TCLs), etc. The distribution network operator (DNO) [26], which is the control center of ADN, gathers the information from both supply and demand side as well as predicts the scheduling information containing customers’ load profiles and DGs’ outputs. In recent years, with the development of communication and measurement technologies such as advanced metering infrastructure (AMI), the DNO has had the access to acquiring all kinds of necessary information in terms of climate data, market data, and load data. A local database therefore can be constructed, providing data support for the follow-up research. Based on the acquired data and physical model of ADN, the power flow analysis can be accomplished, whose function is to coordinate with the optimization program. Moreover, the scheduling information including the users’ load profiles, output of DGs and so on can be also predicted by using forecast algorithms. Then, based on these information, the optimal dispatch can be achieved by applying the efficient solving algorithm. Finally, the optimal results of ADN for the following scheduling times containing charge/discharge of the battery storages, tap position of regulators and so on can be determined and return to the ADN.

Generally, there are two types of decision variables in distribution network: continuous variables and discrete variables. Continuous variables include power of electric components and node voltage, while discrete variables are characterized by states of devices, such as the operation state of capacitor banks and tap position of voltage regulator or OLTC, etc. Hence, the optimization problem should be formulated as a mixed-integer nonlinear programing (MINP) problem. The objective in this paper is to minimize the power loss of ADN while reducing the fluctuation of nodal voltage profile around the nominal value. The details about the distribution network will be discussed in the fifth part.

2.1. Objective Function

One of the advantages of ADN is that a lot of controllable resources, including DGs, battery storages, voltage regulators as well as switchable capacitor banks can be utilized to achieve the goal of optimal dispatch. In this paper, considering the above resources, the first objective is to minimize the total transmission lines energy losses during the whole schedule time. The first objective function is formulated as follows:

$$f_{\mathrm{I}}={\displaystyle \sum _{t=1}^T{P}_{\mathrm{losses}}\left(t\right)\Delta t}$$

(1)

where $\Delta t$ is the time interval and in this paper, its value is 1 h; the variable T represents the scheduling period, which is set to 24 h. $P_{\mathrm{losses}}\left(t\right)$ is the total transmission lines active power losses in kW at time t of ADN. Taking into account the mutual impedance between phases, $P_{\mathrm{losses}}\left(t\right)$ can be precisely calculated by Equation (2) [27]:

$$\begin{array}{ll}{P}_{\mathrm{losses}}\left(t\right)& =\frac{1}{2}{\displaystyle \sum _{ij\in N_l}{\displaystyle \sum _{\gamma =\mathrm{a},\mathrm{b},\mathrm{c}}{\displaystyle \sum _{\beta =\mathrm{a},\mathrm{b},\mathrm{c}}\{-G_{ij}^{\gamma \beta }\left[\left(e_i^{\gamma }\left(t\right)-e_j^{\gamma }\left(t\right)\right)\left(e_i^{\beta }\left(t\right)-e_j^{\beta }\left(t\right)\right)+\left(f_i^{\gamma }\left(t\right)-f_j^{\gamma }\left(t\right)\right)\left(f_i^{\beta }\left(t\right)-f_j^{\beta }\left(t\right)\right)\right]}}}\\ & -B_{ij}^{\gamma \beta }\left[\left(e_i^{\gamma }\left(t\right)-e_j^{\gamma }\left(t\right)\right)\left(f_j^{\beta }\left(t\right)-f_i^{\beta }\left(t\right)\right)+\left(f_i^{\gamma }\left(t\right)-f_j^{\gamma }\left(t\right)\right)\left(e_i^{\beta }\left(t\right)-e_j^{\beta }\left(t\right)\right)\right]\}\end{array}$$

(2)

In the above equation, $e_i^{\gamma }\left(t\right)$ and $f_i^{\gamma }\left(t\right)$ are the real and imaginary part of the voltage at node i of phase $\gamma$ at time t, respectively. $G_{ij}^{\gamma \beta }$ and $B_{ij}^{\gamma \beta }$ are the real and imaginary part of the element corresponding to the γ phase of node i and β phase of node j in the node admittance matrix ($\gamma ,\beta \in \left(\mathrm{a},\mathrm{b},\mathrm{c}\right)$). $N_l$ is the collection of all branches.

The second objective function of improving the voltage profiles can be represented by the following equation:

$$f_{\mathrm{II}}={\displaystyle \sum _{i=1}^N{\displaystyle \sum _{p\in (\mathrm{a},\mathrm{b},\mathrm{c})}\frac{\left|V_i^p\left(t\right)-V_{\mathrm{N}}\right|}{V_{\mathrm{N}}}}}$$

(3)

where $V_{\mathrm{N}}$ is the nominal voltage; $V_i^p\left(t\right)$ means the voltage magnitude at node i of phase p ($p\in \left(\mathrm{a},\mathrm{b},\mathrm{c}\right)$) at time t; and N is the total node number.

Because the units between the two objectives are inconsistent, it is necessary to normalize $f_{\mathrm{I}}$ and $f_{\mathrm{II}}$. The normalization of the two objective functions can be calculated by:

$$\{\begin{array}{l}{\stackrel{-}{f}}_1=\frac{{\displaystyle \sum _{t=1}^T{P}_{\mathrm{loss}}\left(t\right)\Delta t}}{f_{1,\mathrm{B}}}=\frac{f_1}{f_{1,\mathrm{B}}}\\ {\stackrel{-}{f}}_{\mathrm{II}}=\frac{{\displaystyle \sum _{i=1}^N{\displaystyle \sum _{p\in (\mathrm{a},\mathrm{b},\mathrm{c})}\frac{\left|V_i^p\left(t\right)-V_{\mathrm{N}}\right|}{V_{\mathrm{N}}}}}}{f_{1\mathrm{I},\mathrm{B}}}=\frac{f_{\mathrm{II}}}{f_{1\mathrm{I},\mathrm{B}}}\end{array}$$

(4)

where $f_{1,\mathrm{B}}$ and $f_{1\mathrm{I},\mathrm{B}}$ is the basic value of normalization and, in this paper, the value of $f_{1,\mathrm{B}}$ is the maximum power loss during the day before optimization, while $f_{1\mathrm{I},\mathrm{B}}$ means the maximum voltage offset during the day before optimization.

Then, a comprehensive objective function F can be generated after the normalization, reflecting the combined effects of the abovementioned objectives, as Equation (5):

$$\mathrm{Min}F={\omega }_1{f}_{\mathrm{I}}+{\omega }_2{f}_{\mathrm{II}},{\omega }_1+{\omega }_2=1$$

(5)

where ${\omega }_1$ and ${\omega }_2$ are the weight coefficients of $f_1$ and $f_2$, respectively. In this paper, the two objectives are considered of equal importance and their value are both 0.5.

2.2. Constraints

The minimization of Problem (5) is subjected to the following constraints (Equations (6)–(20)):

(1) Unbalanced three-phase power flow constraints:

$$\{\begin{array}{l}{P}_{i,t}^{\mathrm{G},p}-P_{i,t}^{\mathrm{L},p}=V_{i,t}^p{\displaystyle \sum _{j=1}^{N_{\mathrm{bus}}}{\displaystyle \sum _{m\in (\mathrm{a},\mathrm{b},\mathrm{c})}{V}_{j,t}^m\left(G_{ij}^{pm}\cos {\theta }_{ij,t}^{pm}+B_{ij}^{pm}\sin {\theta }_{ij,t}^{pm}\right)}}\\ Q_{i,t}^{\mathrm{G},p}-Q_{i,t}^{\mathrm{L},p}=V_{i,t}^p{\displaystyle \sum _{j=1}^{N_{\mathrm{bus}}}{\displaystyle \sum _{m\in (\mathrm{a},\mathrm{b},\mathrm{c})}{V}_{j,t}^m\left(G_{ij}^{pm}\sin {\theta }_{ij,t}^{pm}-B_{ij}^{pm}\cos {\theta }_{ij,t}^{pm}\right)}}\end{array}$$

(6)

In the following equation, $P_{i,t}^{\mathrm{G},p}$ and $Q_{i,t}^{\mathrm{G},p}$ are the active and reactive generation power injected into node i of phase p, respectively, while $P_{i,t}^{\mathrm{L},p}$ and $Q_{i,t}^{\mathrm{L},p}$ are the corresponding load power at node i of phase p at time t; $V_{i,t}^p$ is the voltage magnitude of phase p at node i and $V_{j,t}^m$ is that of phase m at node j; ${\theta }_{ij}^{pm}$ means the phase angle difference of phase p and m between node i and j; and $G_{ij}^{pm}$ and $B_{ij}^{pm}$ represent the real and imaginary parts of admittance matrix, respectively.

(2) Nodal voltage magnitude constraint:

$$V^{\min }\le V_i^p\left(t\right)\le V^{\max }$$

(7)

In this paper, the minimum and maximum voltage magnitude are set to 0.95 p.u and 1.05 p.u, respectively.

(3) Constraints of active and reactive power outputs of DGs:

$$P_{\mathrm{DG},n}^{\min }\le P_{\mathrm{DG},n}\left(t\right)\le P_{\mathrm{DG},n}^{\max }$$

(8)

$$Q_{\mathrm{DG},n}^{\min }\le Q_{\mathrm{DG},n}\left(t\right)\le Q_{\mathrm{DG},n}^{\max }$$

(9)

where $P_{\mathrm{DG},n}^{\min }$ and $Q_{\mathrm{DG},n}^{\min }$ are the minimum active and reactive power output of DG n, respectively; $P_{\mathrm{DG},n}^{\max }$ and $Q_{\mathrm{DG},n}^{\max }$ are the maximum active and reactive power output of DG n, respectively; and $P_{\mathrm{DG},n}\left(t\right)$ and $Q_{\mathrm{DG},n}\left(t\right)$ are the actual active and reactive output of DG n, respectively. The power generated by DGs is supposed to be assigned to three phases equally.

(4) Thermal limits of transmission lines constraints:

$$\begin{array}{l}{S}_l\left(t\right)=\sqrt{P_l{\left(t\right)}_2+Q_l{\left(t\right)}_2}\\ s.t\hspace{1em}\{\begin{array}{l}{P}_l\left(t\right)={\displaystyle \sum _{p\in \left(\mathrm{a},\mathrm{b},\mathrm{c}\right)}^n{P}_{i,t}^{\mathrm{G},p}-P_{i,t}^{\mathrm{L},p}}\\ Q_l\left(t\right)={\displaystyle \sum _{p\in \left(\mathrm{a},\mathrm{b},\mathrm{c}\right)}^n{Q}_{i,t}^{\mathrm{G},p}-Q_{i,t}^{\mathrm{L},p}}\end{array}\end{array}$$

(10)

$$S_l\left(t\right)\le S_l^{\max }$$

(11)

where $S_l\left(t\right)$ is the apparent power of line l at time interval t and $S_l^{\max }$ represents the maximum transmission power of line l. $P_l\left(t\right)$ and $Q_l\left(t\right)$ are the active and reactive power of line l at time t.

(5) Tap change constraints of voltage regulator or OLTC:

$$ta{p}_k^{\min }\le ta{p}_k\left(t\right)\le ta{p}_k^{\max },ta{p}_k\left(t\right)\in Z$$

(12)

$${\displaystyle \sum _{t=2}^{T}TapChang{e}_k\left(t\right)}\le TapChang{e}_k^{\max }$$

(13)

$$s.t\{\begin{array}{l}Ta{p}_k\left(t\right)-Ta{p}_k\left(t-1\right)=0,TapChang{e}_k\left(t\right)=0\\ Ta{p}_k\left(t\right)-Ta{p}_k\left(t-1\right)\ne 0,TapChang{e}_k\left(t\right)=1\end{array}$$

(14)

where $ta{p}_k^{\min }$ and $ta{p}_k^{\max }$ are the lower and upper limit of tap position of regulator k; $ta{p}_k\left(t\right)$ is the actual position at time t. Equation (12) constraints the actual tap position between $ta{p}_k^{\min }$ and $ta{p}_k^{\max }$. To avoid increasing cost of maintenance caused by changing the tap position frequently, the maximum tap change number constraint, expressed as Equation (13), is introduced during the whole scheduling period, and $TapChang{e}_k\left(t\right)$ means the tap change status at time interval t.

(6) Switchable capacitor bank switch constrains:

$${\displaystyle \sum _{t=2}^{T}CapChang{e}_n\left(t\right)}\le CapChang{e}_n^{\max }$$

(15)

$$s.t\{\begin{array}{l}CapChang{e}_n\left(t\right)=0,Ca{p}_n\left(t\right)-Ca{p}_n\left(t-1\right)=0\\ CapChang{e}_n\left(t\right)=1,Ca{p}_n\left(t\right)-Ca{p}_n\left(t-1\right)\ne 0\end{array}$$

(16)

where $Ca{p}_n\left(t\right)$ is the switch status of capacitor bank n at time t, for which 0 represents open and 1 means closed. In addition, $CapChang{e}_n\left(t\right)$ is a logic variable reflecting whether the switch state is changed.

(7) Battery storage constraints:

$$P_{b,n,\mathrm{discharge}}^{\max }\le P_{b,n}\left(t\right)\le P_{b,n,\mathrm{charge}}^{\max }$$

(17)

$$SO{C}_n^{\min }\le SO{C}_n\left(t\right)\le SO{C}_n^{\max }$$

(18)

$$s.t\{\begin{array}{l}SO{C}_n\left(t\right)=SO{C}_n\left(t-1\right)+P_{b,n}\left(t-1\right)\cdot {\eta }_n/B_n,P_{b,n}\left(t-1\right)\ge 0\\ SO{C}_n\left(t\right)=SO{C}_n\left(t-1\right)+P_{b,n}\left(t-1\right)/{\eta }_n{B}_n,P_{b,n}\left(t-1\right)<0\end{array}$$

(19)

$$0\le times\left(t\right)\le time{s}^{\max }$$

(20)

where $P_{b,n}\left(t\right)$ is the charging/discharging power of battery storage n at time interval t, for which $P_{b,n}\left(t\right)>0$ means charging and $P_{b,n}\left(t\right)<0$ means discharging; $P_{b,n,\mathrm{discharge}}^{\max }$ and $P_{b,n,\mathrm{charge}}^{\max }$ represent the maximum rate for discharge and charge, respectively; $SO{C}_n^{\min }$ and $SO{C}_n^{\max }$ mean the allowable range of state-of-charge (SOC), which are 0.2 and 0.95 in this paper; and ${\eta }_n$ represents the charging/discharging efficiency of battery n. Considering the times of charging/discharging affect the battery life, the charge and discharge times constraint is set in Equation (20), where $time{s}^{\max }$ means the limit of charge and discharge times and its value is 50 in this paper [28,29].

3. Model Solution

In this section, a Kriging model based improved Surrogate Optimization-Mixture-Integer (ISO-MI) algorithm is proposed to solve the optimization problem. The brief introduction of Kriging model is illustrated first. Then, the procedure of the proposed algorithm is explained in detail and, finally, the flowchart of the proposed method is presented.

3.1. Kriging Model

Kriging model is one type of surrogate model; a surrogate model is essentially an approximate model of the complicated physical system, which expresses the relationship between the inputs and the corresponding outputs with simple equation [30]. Nowadays, there are four kinds of commonly used surrogate models including Response Surface Model (RSM), Kriging models, Radial Basis Function (RBF), and Artificial Neural Network (ANN) [31]. Every model has its own characteristics and there is no specific criteria to measure which is the best. Many comparative studies have been done and insights have been gained through a number of experiments [32].

Among these abovementioned models, the Kriging model, combined with a random process, is more accurate for highly nonlinear problems. Besides, Kriging model is also flexible in either interpolating the sample points or filtering noisy data. On the contrary, a polynomial model is easy to construct, clear on parameter sensitivity, and cheap to work with but is less accurate than the Kriging model [33]. In practical engineering application, Kriging model has been successfully applied in many practical projects, such as structural optimization, airfoil aerodynamic design, missile aerodynamic design, multidisciplinary design optimization and so on. In this paper, considering approximation accuracy and the software tools, the Kriging model is selected to approximate the ADN.

Assume that the sampled points can be expressed by $\mathit{X}={\left[{\mathit{x}}_1,{\mathit{x}}_2\dots ,{\mathit{x}}_n\right]}^T$, and the corresponding response values are $\mathit{Y}={\left[{\mathit{y}}_1,{\mathit{y}}_2\dots ,{\mathit{y}}_n\right]}^T$, where $\mathit{X}$ and $\mathit{Y}$ can be obtained through history data of ADN. Kriging function combines linear regression model and random process model to predict the real function of the relationship between input and output [34]. The general form of the Kriging model can be expressed as Equation (21):

$$\hat{y}\left(\mathit{x}\right)={\displaystyle \sum _{j=1}^p{\beta }_j{f}_j\left(\mathit{x}\right)}+z\left(\mathit{x}\right)$$

(21)

where $\hat{y}\left(\mathit{x}\right)$ represents the predictive value of function at point $\mathit{x}$; $f_j\left(\mathit{x}\right)$ is the polynomial regression function of the sampled points; $\hat{\mathit{\beta }}={\left[{\beta }_1,\dots ,{\beta }_p\right]}^T$ denotes the least square estimator of regression coefficients; and $z\left(\mathit{x}\right)$ is a Gaussian random process. The detailed calculation process of these variables can be found in [35] and, in this paper, a MATLAB tool box called DACE [36] is used to construct the Kriging model.

3.2. Kriging Model Based Optimization

As mentioned above, Kriging model is an approximate model for complex black-box system such as complicated distribution system, it is easy to calculate without losing the approximation accuracy.

The general framework of Kriging model based optimization algorithm is shown in Figure 2. To construct Kriging, first, some test points in the decision space should be generated. It is essential to select the design points properly to capture the characteristics of a complicated system. There are two categories of sampling techniques which are the classical techniques and the space filling techniques [37]. Between the two categories, the latter is widely used to construct Kriging models. The space filling designs include three methods namely Orthogonal arrays, uniform designs and various Latin Hypercube Designs (LHD). Among all of these, LHS is selected because of its good properties of uniformly and flexibility on the size of sampling [38]. Then the simulation software can be called to compute the actual response of the inputs and outputs. After that, a Kriging model can be constructed by using these experimental data. Finally, the Kriging models can be applied to the solving process of the optimization algorithm under the premise of meeting the approximation accuracy.

In the process of optimization, the Kriging model can be kept constant or dynamically updated. If the model remains unchanged, the initial constructed model is required to meet certain approximation accuracy requirements. However, in the optimization process, the approximation error of the Kriging model is likely to result in an erroneous solution. Therefore, an effective method is to update the constructed model in the process of optimization, so as to improve the approximation accuracy with the model. In this work, the developed optimization algorithm is inspired by the concept of dynamic Kriging model. As shown in Figure 2, the dynamic Kriging model adopts dynamic adding design point algorithm (DADPA), in which new sample points will be iteratively selected to do the expensive function simulation and update the model until the stopping criterion is met while the whole search space will not be changed. According to [39], DADPA has good performance in optimization problem.

3.3. Improved Surrogate Optimization-Mixed-Integer Algorithm

Inspired by the optimization algorithm SO-MI (Surrogate Optimization-Mixed-Integer) [40] developed to solve mixed-integer, no-convex and nonlinear programming problem and DYCORS (Dynamic Coordinate Search Using Response Surface models) framework developed for bound constrained optimization in black-box system [41], a modified SO-MI named ISO-MI is proposed in this paper to solve the optimal operation and scheduling problem in ADN. Based on the SO-MI, the method of coordinate perturbation is added to improve the efficiency of local search. At the same time, adaptive adjustment strategy of disturbance range is introduced, which balances local and global disturbances and improves the probability of finding better solutions.

The main procedure of ISO-MI can be divided into four steps and the flowchart of the algorithm is shown in Figure 3.

Step 1—Construct Initial Kriging model. In this step, it is easy to construct the Kriging model with high accuracy using the measured history operation data of distribution system. $S_n^{\mathrm{ini}}\left(x\right)$ denotes the initial surrogate model built by the set of sampled points ${\mathit{B}}_n^{\mathrm{ini}}=\left\{x_1,\dots ,x_n\right\}$;

Step 2—Initialize optimization parameters. In this step, first, a series of parameters need to be set for the beginning of the optimization program, such as maximum evolution number M of expensive simulation, initial coordinate perturbation range $r_0$ and the minimum range $r_{\min }$, etc. Second, the LHD method is used to generate initial feasible design points. To construct a new Kriging model, it is necessary to do the expensive simulations and get the responses from them. Finally, find the best point $x_k^{\ast }$, which represents the control variables of the distribution network, and the corresponding objective function value $f\left(x_k^{\ast }\right)$ in the current design points.

Step 3—Iterate until the evolution number, m > M, maximum evolution number. In this step, first, create four groups of candidate points by randomly perturbing coordinates around $x_k^{\ast }$. The purpose is to improve the efficiency of local search by adding the coordinate perturbation. The generation of the four groups in the global range are shown as follows:

(1)

Group 1: Uniformly and randomly perturb the continuous coordinates of $x_k^{\ast }$ at the range of $r_k$, $x_{\mathrm{G}1,i}=x_k^{\ast }+\alpha \cdot r_k$, where i is the index of points and $\alpha ~N\left(0,1\right)$;

(2)

Group 2: Uniformly and randomly perturb the discrete coordinates of $x_k^{\ast }$ at the range of $r_k$, $x_{\mathrm{G}2,i}=x_k^{\ast }+round\left(\alpha \cdot r_k\right)$;

(3)

Group 3: Uniformly and randomly perturb all coordinates of $x_k^{\ast }$ at the range of $r_k$, $x_{\mathrm{G}3,i}=x_k^{\ast }+\alpha \cdot r_k$ and round the discrete coordinates to the closet integers;

(4)

Group 4: Select candidate points $x_{\mathrm{G}4}$ in the whole design space using LHS.

Second, check the candidate points to ensure that they are in feasible domain and discard the points that violate the constraints. Third, select the best candidate points from four groups of candidate points. The specific method is illustrated as follows:

(1)

Calculate the objective function using the initial Kriging model $S_n^{\mathrm{ini}}\left(x\right)$ and current new Kriging model $S_{n,k}\left(x\right)$ of candidate points in four groups and compute the objective function score $s_1\left(x_i\right)$ and $s_2\left(x_i\right)$ of all points, where $s_1$ and $s_2$ are the normalized objective functions. Their value can be calculated by Euclidean distance in n dimensional space.

(2)

Compute the distance score of all design points.

(3)

Compute the weighted score $S\left(x_i\right)={\omega }_1{s}_1\left(x_i\right)+{\omega }_2{s}_2\left(x_i\right)$, where $S$ denotes the objective function F in this paper, and select the point with minimum score $S\left(x_i\right)$ to add it into design points set $B_{n,k+1}$ and do the expensive function evaluation, again.

Fourth, update the best point $x_{k+1}^{\ast }$ and corresponding objective function value $f\left(x_{k+1}^{\ast }\right)$ in the current design points set $B_{n,k}$, where $f\left(x_{k+1}^{\ast }\right)$ is the sub-objective function. Fifth, check if the objective function value has an improvement ($x_{k+1}^{\ast }<x_k^{\ast }$). If this is true, the perturb range should be enlarged to two times. Otherwise, it will be shrink for the purpose of balancing the local and global search. Finally, update the current Kriging model.

Step 4—Output the value of best solution $x^{\ast }$ and the corresponding objective functions found so far.

4. Simulation and Case Studies

To test the effectiveness of the proposed method, we establish a co-simulation platform combining the two software packages GridLAB-D and MALAB. GridLAB-D is a new type of power system simulation tool, which uses an agent-based approach to simulating smart grids [42]. It contains different modules with specific functions, such as power flow module, residential module and so on. The interaction between the two software is shown in Figure 4. In each time step, first, the GridLAB-D software is called to do the expensive power flow simulations of distribution system. Then part of the simulation data needed for optimization is recorded by GridLAB-D and sent to MATLAB through an external application link called matlab.link [43]. After that, the optimal scheduling of ADN is achieved in MATLAB based on the data provided by GridLAB-D. Finally, some control variables in GridLAB-D are changed according to the optimization results. The simulation is implemented on a PC with 3.7 GHz CPU and 16 G RAM.

4.1. Test System Specification

In this paper, a modified three phase unbalanced IEEE-123 distribution system [44], with detailed data given in [45], is utilized to test the proposed method. The topology of the distribution system is shown in Figure 5.

Four regulators and four capacitor banks are installed in this system and their parameters are listed in

Table 1. Parameters of voltage regulators.

Name	Installed Location	Phases	Tap Range	Voltage Regulation Range	Maximum Operating Times
Reg1	150–149	A-B-C	[−16, +16]	[0.95, 1.05]	10
Reg1	25–26	A-C	[−16, +16]	[0.95, 1.05]	10
Reg1	9–14	A	[−16, +16]	[0.95, 1.05]	10
Reg1	160–67	A-B-C	[−16, +16]	[0.95, 1.05]	10

and

Table 2. Parameters of capacitor banks.

Name	Installed Location	Installed Capacity (kVar)			Maximum Operating Times
		Phase A	Phase B	Phase C
Cap1	83	100	100	100	10
Cap2	88	50	0	0	10
Cap3	90	0	50	0	10
Cap4	92	0	0	50	10

. All the control of taps and switches are on the banked level. Besides, seven DGs, whose characteristics are given in

Table 3. Characteristics of DGs.

Name	Installed Location	Type	Rated Power (kW)	Power Factor
DG1	66	WT	150	0.9
DG1	51	PV	100	0.9
DG1	30	MT	150	0.9
DG1	18	MT	200	0.9–1.0
DG1	60	MT	150	0.9–1.0
DG1	108	MT	200	0.9–1.0
DG1	77	MT	150	0.9–1.0

, are located at different nodes of the system in Figure 5. For the purpose of maximizing the use of renewable energy, the outputs of PVs and wind turbines are assumed to be uncontrolled with generating electricity at their maximum power. The forecast outputs of PVs and wind turbines are shown in Figure 6. In terms of loads, they are modeled as ZIP loads, which are comprised of constant impedance Z, constant current I and constant power P loads. The ZIP coefficients is calculated according to [46], which are given in

Table 4. ZIP coefficients of the loads [46].

ZIP Coefficients	Z	I	P
Active load	0.418	0.135	0.447
Reactive load	0.515	0.023	0.462

. If all nodal voltages are maintained at the rated values, the total load profiles of three phases is shown in Figure 7. The parameters of battery storage are listed in

Table 5. Characteristics of battery storage.

Name	Installed Location	Power (kW)	Capacity (kWh)	Efficiency
				Charging	Discharging
BAT1	86	[–150, 150]	750	0.9	0.9

.

4.2. Accuracy of Kriging Model in Distribution System

The accuracy of Kriging model can be estimated by statistic validation methods [47]. In this section, the following indexes are formulated to test the accuracy of the proposed model.

(1)

Root Mean Square Error (RMSE):

$$RMSE=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(y_i-\widehat{y_i}\right)}^2}}$$

(22)

(2)

Relative Maximum Absolute Error (RMAX):

$$RMAX=\max \left|\left(y_i-\widehat{y_i}\right)/y_i\right|$$

(23)

(3)

Relative Average Absolute Error (RAAE):

$$RAAE=\frac{1}{N}{\displaystyle \sum _{i=1}^N\left|\left(y_i-\widehat{y_i}\right)/y_i\right|}$$

(24)

(4)

R-Square:

$$R^2=1-{\displaystyle \sum _{i=1}^N{\left(y_i-\widehat{y_i}\right)}^2}/{\displaystyle \sum _{i=1}^N{\left(y_i-\stackrel{-}{y}\right)}^2}=1-\frac{MSE}{Variance}$$

(25)

where $y_i$ is the actual response of black-box system which means distribution system in this paper; $\widehat{y_i}$ is the predicted value; and $\stackrel{-}{y}$ is the mean of actual response value. In the above indexes, RMSE and RAAE are used to measure the overall accuracy of the model, while RMAX is used to gauge the local accuracy. In addition, MSE (Mean Square Error) represents the departure of the Kriging model from the real simulation model and the Variance captures how irregular the problem is. Lower values of RMSE, RMAX and RAAE or closer value of R-Square to 1 lead to a more accurate Kriging model [48].

The control variables at each scheduled time in this paper contain tap positions of the four regulators, switching status of the four capacitor banks and the output of the battery storage. The Kriging model is developed with the designed points generated by LHS. One hundred randomly sampled test points are used to validate the accuracy of the model and the test results are given in the following table in which N represents the number of designed points.

As shown in

Table 6. Test results of accuracy of Kriging model.

N		50	100	200
Index
Voltage fluctuation	RMSE	9.4 × 10−5	9.0 × 10−5	9.1 × 10−5
	RMAE	7.5 × 10−4	7.2 × 10−4	6.2 × 10−4
	RAAE	6.3 × 10−5	5.9 × 10−5	5.2 × 10−5
	R2	1.0000	1.0000	1.0000
Power loss	RMSE	2.3 × 10−3	1.9 × 10−3	1.3 × 10−3
	RMAE	3.4 × 10−3	2.5 × 10−3	1.9 × 10−3
	RAAE	1.1 × 10−3	0.9 × 10−3	0.7 × 10−3
	R2	0.9954	0.9973	0.9988

, the Kriging model has high accuracy to approximate the characteristics of the active distribution system and the accuracy is improved with the increase in the number of designed points. In the case of 50 designed points, the average relative absolute error of voltage fluctuation and power loss is less than 1%. Trading off between the approximation accuracy and the time required to construct the model, this paper selects 50 initial points to build the initial Kriging model.

4.3. Solving Efficiency of ISO-MI

In this section, the solving efficiency of the proposed method is illustrated compared with GA and PSO. The maximum expensive function evaluation number is set to 200. Besides, the initial SOC of the battery storage is 0.5. The upper and lower limits of SOC are 0.95 and 0.2, respectively. The schedule period is 24 h, while the time interval is Δt = 1 h. The weights, ω₁ and ω₂, are both 0.5. The simulation results are displayed as follows. The comparisons among ISO-MI, GA and PSO are shown in Figure 8, Figure 9 and Figure 10. The optimization results of tap positions of Reg1–Reg4 are displayed in Figure 11 and the charge power and SOC of BAT1 are shown in Figure 12.

As shown in Figure 8, ISO-MI has a better convergence performance compared with GA and PSO. Using ISO-MI, the smallest objective function value of 0.35 is obtained. In contrast, the objective function values of GA and PSO are more than that of ISO-MI with the same evaluation number of 200, about 0.5 and 0.55. Moreover, the simulation time of ISO-MI is about 360 s, shorter than that of GA and PSO, about 720 s and 540 s, on the condition that the objective function value is around 0.35. As shown in Figure 9, the smallest power loss is obtained when ISO-MI is used, indicating that ISO-MI has better effect on reducing the overall power loss than GA and PSO. It can be calculated that the total power loss is about 1077 kW per day in the case of ISO-MI while the value is 1441 kW and 1550 kW per day when using GA and PSO. As shown in Figure 10, all the nodal voltages are within the allowable range of 0.95–1.05. It can also be observed that the optimal result obtained by ISO-MI has smallest voltage fluctuation, demonstrating that the proposed method is more efficient than the other two. The phenomenon above all demonstrate that the proposed method has better convergence characteristics and is more efficient in the optimization problem.

As shown in Figure 11, it is found out that the taps of voltage regulators, Reg1, Reg2 and Reg3, change more frequently than Reg4. The tap position of Reg4 is kept in the lower positions all the time. This is because the downstream area of Reg4 is small and its effect on the whole active distribution system is little. Considering the ZIP load, the tap position should not be set too high to reduce the power consumptions.

The charging and discharging power as well as the SOC profile of the battery storage BAT1 is demonstrated in Figure 12. It can be observed in Figure 12b that the SOC of the battery storage is kept within the range of 0.2–0.85, satisfying the upper and lower limits. As shown in Figure 12a, to reduce the fluctuation of the voltage and power loss, the battery storage tends to charge when the renewable energy outputs are in their peak and discharge during the time that the renewable energy outputs are low. For example, during 06:00–07:00 and 21:00–24:00, the battery storage is in discharging state, while during 11:00–15:00, the BAT1 is in the state of charging. In addition, similar to Figure 12a, the SOC of the battery storage is also large when the renewable energy outputs are high in Figure 12b.

5. Conclusions

In this paper, the optimal operation model of the ADN is established. The model considers the controllable resources, such as distributed generation units, energy storages, voltage regulators and switchable capacitor banks. Aiming at solving the optimal mathematical model, a Kriging model based optimization algorithm named ISO-MI is proposed to improve the efficiency of the optimization algorithm. Finally, this paper validates the model and the solution method by IEEE-123 test system. The simulation results indicate that ISO-MI yields a better solution than GA and PSO within the predefined expensive evaluation time limit. However, the optimization results of this paper are based on the prediction accuracy of 100% and, due to the existence of errors, the prediction error will inevitably have an adverse impact on the optimization results. Therefore, the optimal operation model of the ADN with multiple time scales should be established, which further improves the economy and reliability of the active distribution network and will be the focus of our future research.