Front. Energy Res. Frontiers in Energy Research Front. Energy Res. 2296-598X Frontiers Media S.A. 843241 10.3389/fenrg.2022.843241 Energy Research Original Research A GPU-Based Resilience Enhanced Voltage Optimization Model for Distribution Networks Liang and Luo GPU Parallelized Voltage Optimization Model Liang Liang 1 * Luo Chuhang 2 1 Harbin Institute of Technology, Shenzhen, China 2 Guangzhou Power Supply Bureau of Guangdong Power Grid Corporation, Guangzhou, China

Edited by: Shunjiang Lin, South China University of Technology, China

Reviewed by: Bin Zhou, Hunan University, China

Haoran Liu, The Chinese University of Hong Kong, Hong Kong SAR, China

 *Correspondence: Liang Liang, liangl@hit.edu.cn

This article was submitted to Smart Grids, a section of the journal Frontiers in Energy Research

 14 03 2022 2022 10 843241 25 12 2021 21 02 2022 Copyright © 2022 Liang and Luo. 2022 Liang and Luo

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Improving the survivability of critical loads after extreme events is essential to enhance the resilience of power systems, especially for distribution networks. A distribution network with various operational resources can be separated into several sub-distribution networks without electrical connections. Maintaining the power supply with acceptable power quality to critical loads in such separated distribution networks is a challenging task for the operators of power systems. In this paper, an optimization model is proposed to maximize the ability to supply power to critical loads in distribution networks. Moreover, a GPU was employed to accelerate the proposed model using genetic algorithm. With the acceleration of the GPU platform, the solving time was reduced and the population size can be enlarged to enhance the convergence rate and convergence quality of the algorithm. Finally, case studies were carried out in IEEE 33-bus and 118-bus systems, and the effectiveness of the method was validated by comparing the solution results on GPU and CPU platforms.

 resilience voltage regulation distribution network gpu genetic algorithm 52077045 National Natural Science Foundation of China10.13039/501100001809 Introduction

In recent years, a number of extreme natural events have brought huge economic losses to the society (Panteli and Mancarella, 2015a; Wang et al., 2016; Mohamed et al., 2019). Critical loads, such as government buildings, communication systems, transportation systems, and the sensors and controllers of power systems, are essential for the operation of power systems during extreme natural events, as well as the restoration of power systems after such extreme natural events. Improving the survivability of critical loads after these extreme natural events can enhance the resilience of the power system. Resilience is defined as the survivability and the ability to defend against extreme natural events. It is defined as the ability of the power system to survive, defend, and recover from extreme natural disasters (Panteli and Mancarella, 2015b; Gao et al., 2016). Nowadays, more and more distributed generation resources, such as photovoltaic (PV) systems, small scale wind turbines, and gas turbines are integrated into distribution networks. Moreover, most loads of power systems, especially critical loads, are supplied by distribution networks. The operation of distribution network during and after extreme natural events is essential to enhance the resilience of power systems.

A huge challenge for the distribution network operation during and after extreme natural events is the insufficiency of controllable operational resources. Conventional fossil fueled power generators with large capacities may not be able to connect to distribution networks because of the damaged transmission networks. In this case, all kinds of operational resources in distribution networks, such as distributed PV resources, distributed wind farms, small scale energy storage systems, diesel generators, demand response devices, and voltage support devices, should be well coordinated to provide sufficient operation resources.

Currently, a large number of distributed PV stations, other renewable resources, and energy storage systems with declining costs are connected to the distribution network (Mahmud et al., 2014). Comparing with traditional distribution networks, these distributed resources make the operating environment of the distribution network more complex and challenging (Wang et al., 2020). As indicated in (Wang et al., 2020), the PV systems could take effects on the voltage quality of distribution network in different locations. On the other hand, these resources can be survived after extreme natural events in most cases and bring more flexible resources to the distribution networks and provide new technical devices and methods to enhance power supply to critical loads. All types of flexible operational resources should be effectively involved in the operation of the distribution networks to ensure a stable power supply of critical loads and support restoration after extreme natural events (Utkarsh et al., 2021). An operation framework is proposed in (Utkarsh et al., 2021) to coordinate the demand response resources, such as the home energy systems.

Regulating bus voltages within acceptable area is an important operational target in distribution networks. Considering the critical loads are sensitive to the voltage quality, the voltage optimization model can improve the survival ability of critical loads. The distribution networks should dispatch various types of reactive power and voltage support devices, such as distribution static synchronous compensators (D-STATCOMs), capacitor banks, and line voltage regulators (LVRs). These devices have different operational characteristics. For instance, the capacitor banks (CBs) are ON/OFF-type devices and several discrete control states are available for LVRs. While D-STATCOMs can adjust its output power continuously. Thus, coordinate these devices with different types of mathematic models is a challenge for the operation of distribution networks. An optimal operational model is proposed in this paper to coordinate the operation of different types of devices with a multi objective optimization model. The proposed model is a mixed integer nonlinear programming (MINLP) model. These models can be solved by traditional solution methods and commercial software, such as the gradient based optimization algorithm (Nicholson and H. Sterling, 1973; Mamandur and Chenoweth, 1981; Grudinin, 1998), and decomposition approach (Deeb and Shahidehpour, 1990). However, the nonlinear characteristics and binary variables of the optimization model make the global optimization results can not be achieved guaranteed. The computing time is varied a lot depends on the proposed mathematic models. The Genetic algorithms (GAs) are metaheuristic algorithms and the first GA was proposed by Professor J Holland in 1975. GAs are widely used to solve the nonlinear programming (NLP) optimization models of power systems and can achieve better robustness for various kinds of mathematic models (Delfanti et al., 2000; Swarup and Yamashiro, 2002; Enacheanu et al., 2008; Queiroz and Lyra, 2009). However, in solving MINLP problems with strong nonlinear objective functions, the calculation time and convergence quality requires a large population number and a long solving time. In this way, the GA is not acceptable for practical applications.

Population size is an essential parameter for the performance of the GA. By choosing a large initial population size, the convergence rate and convergence quality of the algorithm can be enhanced to obtain better results (Harik et al., 1999). However, increasing the population size will greatly increase the computing time, which significantly reduces the value of the algorithm for practical projects. The genetic operation in the evolution process has strong independence, and there is no data dependence between individuals, so the algorithm can be computed parallelly (Xing Chen et al., 2005). Based on this feature, Xing Chen, F. Gonzalez Bulnes et al. used a distributed computing method to calculate the population size generated on multiple computers, and finally mixed the population to obtain the result (Gonzalez Bulnes et al., 2013; de la Calle et al., 2015). However, the distributed computing method needs to switch data results among multiple computers, which requires powerful communication capabilities and increase the total computing time. The cost of setting up powerful control and communication systems limits the application of distributed and parallel GAs.

In this study, the proposed optimization model was solved by GA and accelerated on a graphics processing unit (GPU) platform by implementing a parallel computing structure. With the rapid development and iteration of AI technologies, such as deep learning and reinforcement learning, the hardware performance of GPU platforms have increased rapidly. For instance, the GPU of Nvidia 2080TI has Turing TU102 architecture, 11 GB GDDR6 memory, and the single precision floating point computing performance can exceed 10 Tera floating-point operations per second (TFLOPS). The GPU can process large amounts of data in parallel in the form of single instruction multiple data (SIMD) through its unique hardware architecture, which greatly increases the computing speed of artificial intelligence algorithms. This increased computing speed, in turn, makes artificial intelligence algorithms widely used (Topa et al., 2011a; Topa et al., 2011b). Thanks to the powerful computing performance of GPU platform, a GPU and GA combined structure was applied to solve the GA in this paper. The GPU parallelization technology was used to improve the convergence rate and convergence quality of the algorithm while reducing the calculation time to ensure that the proposed method can potentially be applied to practical projects. Currently, there are some similar works about the parallelized GA algorithms on GPU platforms (Tsutsui and Fujimoto, 2009; Jaros, 2012; Luo et al., 2017). In (Tsutsui and Fujimoto, 2009), a standard quadratic assignment problem was solved by GA on a GPU platform parallelized. The results showed that the GPU can accelerate the quadratic assignment problem with GA effectively. Moreover, multiple GPUs are applied to improve the performance for Gas on a knapsack problem in (Jaros, 2012). In (Luo et al., 2017), the GPU accelerated GA/MLP algorithms are applied to extract and classify the features of electrical signals of neural systems. These indicates the voltage optimization problem is potential to be solved by the GA on GPU platform. In this paper, an improved GA was used to solve the proposed optimization model of the distribution network. This work is a try to apply the parallelized GPU solved GA algorithms for the voltage optimization problem in distribution network.

The major contributions of this paper are listed as follows:

1) An optimization model was proposed to coordinate the different types of operational resources to maximize the resilience of the distribution network. Power flow constraints were also considered in the proposed model.

2) A parallel GA was proposed to solve the proposed model on a GPU platform. The solving performance was improved by the proposed method and has the potential to be applied in practical projects.

3) The critical loads in the distribution network can be reliably supplied with power using the proposed model in this study. Multiple operational targets can be coordinated in the proposed model.

The remainder of this paper is organized as follows: Modeling Typical Reactive Power Compensation Devices presents the models for various kinds of resources. Voltage Optimization Model in Distribution Network presents an optimization model of the distribution network, and describes the large population genetic algorithm and the parallel computing architecture of the GPU in detail. Then, in Case Studies, IEEE 33-and 118-bus systems are tested, and the results demonstrate the effectiveness of the proposed methods. Finally, Conclusion draws the conclusions in this study.

 Modeling Typical Reactive Power Compensation Devices

As shown in Figure 1, after extreme natural events, such as typhoons, snow disasters and storms, distribution networks may lose power supply from bulk power systems. Critical loads such as hospitals, government agencies, and automatic devices associated with the operation of power systems are directly connected to the distribution networks. Once the power supply to these critical loads is lost, the emergency repair and restoration speed of the power systems will be affected and reduced. Moreover, cascading failures may be occurred due to the lost of these critical loads. To ensure a stable power supply of critical loads after extreme natural events, it is necessary to dispatch all types of controllable operational resources in the distribution networks to regulate bus voltages and maintain the system frequency. Modeling these operational resources is an essential task for the proposed optimization model in this paper.

Condition of distribution network after extreme natural events.

 Capacitor Banks

CBs are the cheapest reactive power compensation devices in power systems. The CB is consisted of several group of capacitors in both serial and parallel connection as shown in Figure 2.

Structure of the CBs.

Each group of capacitors can be connected and disconnected from the power systems by the status of the switch. In this way, the CB is a kind of ON/OFF controlled device and can be modeled as a 0/1 variable in mathematic models. Q i C , o u t = Q C ( U i U N ) 2 N C i 0 ≤ N C i ≤ N C max i , ∀ i ∈ Ω C where Q i C , o u t is the output reactive power of the CB at bus i. The U i represents the voltage at bus i and U N represents the nominal voltage. Ω C is the set of nodes where the capacitor banks are arranged, Q C is the rated output reactive power of each group capacitors under the nominal voltage. N C i is the number of switches to be closed to active the capacitors. N C max i is the upper limit of the number of closed switches and is determined by the group number of CBs.

 Line Voltage Regulator

The LVR is consisted of a transformer and several solid-state switches as shown in Figure 3. The LVR could regulate the bus voltage to several setpoints with a fixed step, such as 2.5%, around the nominal voltage. It is achieved by changing the status of the solid-state switches to change the ratio of the transformer connected with the bus.

Structure of the LVR.

In this way, the voltage of the bus connected with LVR can be modeled as a positive discrete variable in mathematic models. U i = U 0 + N L V R i U S t i − N L V R max i ≤ N L V R i ≤ N L V R max i , ∀ i ∈ Ω L V R where Ω L V R is the set of nodes where the LVR is arranged, U i represents the voltage of bus i, U 0 is the nominal bus voltage, U S t i is the voltage step for the LVR installed at bus i, N L V R i is the setpoint for the LVR at bus i and is a discrete variable. N L V R max i is the limit for the setpoint of LVR.

 D-STATCOM

D-STATCOM could compensate the reactive power with fast response and high accuracy. It is consisted of a full bridge power converter and could adjust its output reactive power according to the operation commands as shown in Figure 4. Different from the CB and LVR, D-STATCOM can be modeled as a continuous variable in mathematic models. Q i S , o u t = Q S i P i S , o u t = 0 Where Q i S , o u t represents the output reactive power of D-STATCOM installed at node i. Q S i represents the reactive output power of D-STATCOM for bus i. P i S , o u t represents the output active power of D-STATCOM installed at node i.

Structure of the D-STATCOM.

All of these resources need to be coordinated to regulate the voltages of the distribution power systems.

 Voltage Optimization Model in Distribution Network Optimization Model

Determining the power flows in a distribution network is an important task for the optimization of distribution networks. The radial structure is the most common topology for practical distribution networks, as shown in Figure 5. The Distflow equation can be applied to determine the power flows in the radial distribution networks as follows (Baran and Wu, 1989), P l , j = P b , i j − ∑ k ∈ S ( j ) ( P b , j k + P b , j k 2 + Q b , j k 2 U k 2 R j k ), Q l , j = Q b , i j − ∑ k ∈ S ( j ) ( Q b , j k + P b , j k 2 + Q b , j k 2 U k 2 X j k ), U j 2 = U i 2 − 2 ( P b , i j R i j + Q b , i j X i j ) + ( R i j 2 + X i j 2 ) P b , i j 2 + Q b , i j 2 U i 2 , where k ∈ S ( j ) , the set S(j) refers to the set of child nodes of node j.

Radial network with active power flow.

The objectives of the proposed optimization model network are to minimize the active power loss, voltage deviations of buses, and operational costs by dispatching controllable voltage regulating devices. It includes three main objectives as follows,

1. Bus voltage: This is one of the most important targets for ensuring a stable power supply to critical loads. Due to the fluctuating output power of renewable resources, controllable devices need to respond to fluctuating output power in a short time to reduce voltage deviations.

2. Active power loss: The energy capacity of energy storage systems is limited, and the distributed power generation is random. To maintain power supply for as long as possible, it is necessary to reduce the active power loss of the distribution networks.

3. Cost of devices: Compared to D-STATCOM, capacitor banks have smaller operational costs. To reduce the total cost of the distribution network economy, capacitor banks will be dispatched preferentially.

The optimization model of the distribution network is proposed as follows, min Q S 1 , ... , Q S m , N C 1 , ... , N C m γ v C V Q + γ p C p l o s s + γ c C cos ⁡ t . such that C p l o s s = ∑ ( i , j ) ∈ N d P b , i j 2 + Q b , i j 2 U j 2 R i j C V Q = ∑ i ∈ M ( U i − U r e f ) 2 C cos ⁡ t = λ c ∑ i ∈ Ω ∩ Ω S ( Q S 2 ) 2 P l , j = P b , i j − ∑ k ∈ S ( j ) ( P b , j k + P b , j k 2 + Q b , j k 2 U k 2 R j k ) , ∀ j ∈ N d Q l , j = Q b , i j − ∑ k ∈ S ( j ) ( Q b , j k + P b , j k 2 + Q b , j k 2 U k 2 X j k ) , ∀ j ∈ N d U j 2 = U i 2 − 2 ( P b , i j R i j + Q b , i j X i j ) + ( R i j 2 + X i j 2 ) P b , i j 2 + Q b , i j 2 U i 2 ∀ ( i , j ) ∈ N d Q l , i = Q l , i c − Q l , i c o m , ∀ i ∈ Ω Q l , i c o m = Q S i + Q C ( U i U N ) 2 N C i , ∀ i ∈ Ω U r e f − ξ ≤ U i ≤ U r e f + ξ , ∀ i ∈ N d − Q S max i ≤ Q S i ≤ Q S max i , ∀ i ∈ Ω S , 0 ≤ N C i ≤ N C max i , ∀ i ∈ Ω C , and the constrains as in (Eqs. 1–7).

Where N d is the set of nodes of the distribution network, M is the set of critical nodes, Ω is the set of nodes where the voltage regulation devices are arranged, Ω S is the set of nodes where the D-STATCOM is arranged, Ω C is the set of nodes where the capacitor banks are arranged, C p l o s s is the active power loss of the distribution network, C V Q is the bus voltage offset of critical loads, C cos ⁡ t is the D-STATCOM and CBs operation cost, and γ p , γ v , γ c are the weight coefficients. The LVR is installed at the root node of the distribution network and the voltage of the root bus can be controlled to several discrete values as constrained by Eqs. 3, 4. Constraint (11) represents the active power loss of the distribution network, which is equal to the summation of the active power loss in all branches of the distribution network. Constraint (12) represents the bus voltage deviations of the critical loads, which represent the voltage quality of critical loads. The bus voltage deviations can be calculated based on the bus voltage and reference voltage. The reference voltage is related to the demands of critical loads and is set to 1 p.u. in most cases. Constraint (13) represents the operational costs associated with the operations of the D-STATCOMs, which are related to the reactive power output of the D-STATCOMs. Constraints (14–16) are the Distflow equations. Constraint (17) represents the actual reactive power load of the node i , which is related to the reactive power generated by the voltage regulation devices. Constraint (18) represents the reactive power generated by the voltage regulation devices. The reactive power output of the capacitor banks is related to the number of groups and the bus voltage of the corresponding node. Constraint (19) guarantees that the voltage of each node meets the operating standards of distribution networks, U r e f represents the nominal voltage for all buses and ξ is commonly set to 0.1 p.u. This constraint ensures the voltage deviations for all buses are within secured region. The voltage deviations on the critical loads can be minimized by the objective function described in Eq. 12. Constraint (20) represents the upper and lower limits of the reactive power output of the D-STATCOM at node i. Constraint (21) represents the upper and lower limits of the number of capacitor banks that can be switched at node i.

 Improved Genetic Algorithm

As shown in Figure 6, the standard GA consists of the following parts: population initialization, code, fitness function, selection operator, crossover operator, and mutation operator.

Calculation procedure of GA.

The current mainstream coding methods include binary coding and real number coding. Although binary encoding is simple to operate and can prove the convergence of GA through the mode theorem, The hamming cliff problem is still existed and the binary variable takes up a lot of storage. To avoid these problems, the real number encoding method was chosen in this study.

Population Initialization

Although the GA is widely used to solve various types of optimization model problems, the performance of the standard GA is not good enough to be applied in practical projects that have a large number of variables and strong nonlinear characteristics as an MINLP problem. The convergence quality of the standard GA is poor and it easily converges to local optimal solutions. Population size is the main factor that determines the convergence rate and convergence quality of the GA. When the population size is small, the diversity of individuals is insufficient, and the search space is small. In particular, for the optimization model with a large number of variables, the solution space of the initial individuals covered is limited. Once an individual obtains a local optimal solution, majority of individuals will converge to the solution after several iterations of evolution, and it will be difficult to jump out of the local optimal solution by the operation of the crossover or mutation operator, which results in the prematurity phenomenon. In addition, the penalty function leads to the degradation phenomenon in the early stage of the iteration, which significantly reduces the evolution speed.

Increasing the population size can increase the diversity of the individuals, making the initial solution cover a larger area. This initial solution not only restrains the degradation of the individuals in the early iteration, but also improves the convergence rate, and expands the solution space of the individuals covered. This, in turn, significantly improves the search ability of the GA. Consequently, there is a high probability of obtaining a global optimal solution.

 Fitness Function

The fitness function is essential to the performance of the GA, which evaluates the excellence of an individual. The fitness function is closely related to the objective function of the optimization model, whereas the constraints in the optimization model cannot be neglected. Hence, the fitness function will be constituted by two parts: the first part is the objective function written in Eq. 10, and the other part is the penalty function established to satisfy the constraints. The fitness function can be described using the following equation. f = 1 γ p C p l o s s + γ v C V Q + γ c C cos ⁡ t + λ u ∑ i ∈ N d Δ U i 2 ( U max − U min ) 2 Where λ u is the coefficient of penalty function, is described as: Δ U i = { U max - U i , 0 , U min − U i , i f i f i f U i > U max U min ≤ U i ≤ U max U i < U min

The most common selection operators are the roulette, tournament, and elite methods. In this study, the selection operator is composed of the tournament method combined with the elite method to ensure convergence of the algorithm. In addition to good time complexity, the tournament method has a greater probability of selecting better individuals than the roulette method.

 Crossover and Mutation Operation

The crossover operation is a search process. Individuals can gradually approach the global optimal solution with the crossover operation. Because the coding method is based on real number coding, the simulated binary crossover (SBX) was selected as the crossover operator. This can be described by the following equations. c 1 = 0.5 [ ( 1 + β ) p 1 + ( 1 − β ) p 2 ], c 2 = 0.5 [ ( 1 − β ) p 1 + ( 1 + β ) p 2 ], where is calculated as: β = { ( 2 u c ) 1 / ( η c + 1 )            u c ≤ 0.5 ( 1 2 ( 1 − u c ) ) 1 / ( η c + 1 )    u c > 0.5 , p 1 , p 2 are the parents to be crossed; c 1 , c 2 are the offspring generated by the cross operation; β is the spread factor; u C is a random number between 0 and 1; η c is the cross-distribution index. From Eqs. 18–20, it can be seen that the difference between the values of the offspring and the parent is inversely proportional to η c , which affects the convergence rate and the convergence precision of the GA. The standard GA is usually set η c = 2 to maintain the balance of the convergence precision and convergence rate. However, in the early stage of evolution, a smaller η c would lead to larger individual change steps to improve the convergence rate. In the later stage of evolution, when the population has gathered near the global optimal solution, more attention is paid to the convergence accuracy, and a larger η c can make the search more accurate. In addition, due to the differences in variable types, such as 0/1 variables and continuous variables, the value of the cross-distribution index needs to be set separately. Hence, the value of the cross-distribution index before i i t e r is η c , 0 / 1 , η c , c o n t i , and the value of the cross-distribution index after i i t e r is η c , 0 / 1 ' , η c , c o n t i ' .

The mutation operation aims to avoid premature convergence of solutions by expanding the existing search space. A polynomial mutation was selected as the mutation operator. This can be described by the following equation. ν k ' = ν k + δ ( u k − l k ) , Where δ is calculated as δ = { [ 2 u m + ( 1 − 2 u m ) ( 1 − δ 1 ) η m + 1 ] 1 / ( η m + 1 ) , u m ≤ 0.5 1 − [ 2 ( 1 − u m ) + 2 ( u m − 0.5 ) ( 1 − δ 2 ) η m + 1 ] 1 / ( η m + 1 ) , u m > 0.5 . The δ 1 , δ 2 is described as: δ 1 = ( ν k − l k ) / ( u k − l k ), δ 2 = ( u k − ν k ) / ( u k − l k ), where ν k ' is the offspring generated by the mutation operation, ν k is the parent to be mutated, u k is the upper limit of the variable k, l k is the lower limit of the variable k, δ is the coefficient of variation, u m is a random number between 0 and 1, and η m is the variation distribution index. Similar to the crossover operation, the variation distribution index also needs to be set separately. Before i i t e r , the value of the cross-distribution index is η m , 0 / 1 , η m , c o n t i and after i iter , the value of the cross-distribution index is η m , 0 / 1 ' , η m , c o n t i ' .

The calculation procedure for the improved GA is illustrated in Figure 7. Modifying the individuals that do not meet the constraints, preserving the optimal solution of each iteration and setting the distribution index by stages are added.

The calculation procedure of improved GA.

 Parallel Computing Structure of GPU

Although increasing the population will effectively improve the performance of the algorithm, the population size is generally set from 50 to 100 due to the linear increase in computation time. The typical time intervals are several seconds or minutes depending on the control systems of the distribution networks, which require that the solving time of the optimization model for the distribution networks cannot be too long to update the operation commands. Considering that the calculation time of the fitness function takes up to 90% of the total calculation time for the GA algorithm, which is not a logical and temporal correlation between individuals, it is essential to reduce the calculation time of fitness functions in parallel computing to make the GA algorithm available for practical projects.

At present, artificial intelligence algorithms, such as deep learning and reinforcement learning, are widely used in various fields. To take advantage of AI algorithms, a large amount of data is required to train the model. Hence, new structures of computing systems, such as GPUs and field programmable gate arrays (FPGAs), that can quickly process massive data parallelly, have been developed in recent years.

The GPU has a unique hardware structure. Compared with the CPU, the GPU abandons a large number of storage and control units to configure additional computing units. Hence, compared to the CPU, which is good at processing complex logic operations, GPUs are better at processing large amounts of data with the same instructions. These computing units are called streaming processors (SPs). Multiple SPs form a streaming multiple processor (SM). Each SM has an independent share memory, and all SMs share the global memory and constant memory. Thus, the data in different SMs can be shared. When solving a GA on a GPU platform, each individual can be assigned to different computing cores for parallel computing, and the population size can be set to a larger value than a CPU platform. As a result, increasing the population size of the GA has less effect on the computing time than the CPU platform.

Figure 8 shows the process of using a GPU to accelerate the GA. The data of individuals are transmitted from the host (CPU side) to the device (GPU side), and then a large number of threads are opened on the device to be prepared for calculation. Each thread is responsible for calculating the fitness function of an individual. All the fitness functions of individuals are calculated in parallel mode on the GPU. Finally, the results were transmitted from the device back to the host.

Calculation procedure of fitness function on GPU platform.

 Case Studies

To verify the effectiveness of the proposed optimization model and solving method, two distribution network systems, IEEE 33-bus and IEEE 118-bus as shown in Figure 9 and 10 respectively, were selected for case studies. The computing system includes an Intel i9-9900K CPU, Nvidia 2080 Ti GPU, 32GB RAM, 12GB GPU RAM, and CUDA version 10.2.

IEEE 33-bus distributed network topology.

IEEE 118-bus distributed network topology.

Table 1 shows the configuration of the test systems, including the nodes of the critical load and compensators, as well as the type and capacity of the compensators. The parameters of the GA are selected as follows: n _ g e n e r a t i o n = 500 , crossover probability P c = 0.85 , mutation probability P m = 0.1 , i i t e r = 20 , η c , 0 / 1 = 3 , η c , c o n t i = 6 , η m , 0 / 1 = 1 , η m , c o n t i = 3 , η c , 0 / 1 ' = 5 , η c , c o n t i ' = 9 , η m , 0 / 1 ' = 2 , η m , c o n t i ' = 5 .

Configuration of test systems.

Test system Critical load node Node of compensator Compensator Volume of compensator
IEEE 33-bus 11 7 Capacitor Banks 12 Groups/0.15 MVAR
15 10
21 18
27 22 D-STATCOM −2∼2 MVAR
33 25
33
IEEE 118-bus 35 54 Capacitor Banks 15 Groups/0.15 MVAR
62
51 76
105
71 113
25 D-STATCOM −2∼2 MVAR
95 75
85
111 91
111

When the population size is 50, the solution results are poorly convergent and the maximum standard deviation of the variables is 0.4. Although the standard deviation of the objective function is 0.00064, the binary (discrete) type variables cannot converge to a unique solution. Then, the population size was increased from 50 to 5,000, while the other parameters remained the same. The optimization model is solved again, and new solution results are obtained with different population sizes.

Figure 11 shows a comparison of the convergence speed of the algorithm with different population sizes. This indicates that as the population size increases, the initial solution and convergence speed of the algorithm become better. When the population size is 5,000, the convergence of the solution is completed at the 18th generation, while the population size of 50 achieves convergence of the solution at the 35th generation. In addition, the values of the objective function will be smaller when the population size is larger. This means that better optimized results can be achieved. When the population size increases to 5,000, the convergence quality of the algorithm is improved. The binary (discrete) type variables can converge to a unique solution, while the maximum standard deviation of continuous variables is 0.00088, and the objective function can converge to the same value. The results of control variables of CBs installed at bus 7, 10, 18 are 5, 2, and 1, respectively. The output reactive power of D-STATCOMs installed at bus 22, 25 and 33 are 0.08MVar, 0.283MVar, and 0.543MVar, respectively.

Comparison of consequence speeds with different population size in IEEE 33-bus system.

However, increasing the population size will greatly increase the calculation time of the algorithm, as mentioned before. When the population size is increased from 50 to 5,000, the running time is also changed from 2.18 to 180.69 s, which is an increase of 90 times.

Table 2 shows the comparison of computing time on both the GPU and CPU platforms with different population sizes in the IEEE 33-bus and IEEE 118-bus system. The GA algorithm computed on the GPU platform has been revised to parallelized to enhance the computing performance. For the IEEE 33-bus system, on the GPU platform, the growth in the population size had almost no effect on the running time. Compared with the population size of 50 and the population size of 5,000, the difference in calculation time is only 1.15 s. This means that within a certain population size range, the larger the population size, the better the acceleration effect of the GPU. The algorithm is accelerated by parallel computing on the GPU platform, which greatly reduces the running time and effectively improves the practicality of the algorithm. In testing the IEEE 33-bus system, increasing the population size of the GA can effectively enhance the convergence ability of the algorithm. The long computing time caused by the large population size was effectively reduced through parallel computing on the GPU platform.

Computing time of CPU platform and GPU platform with different populations (IEEE 33-bus system and IEEE-118-bus system).

Test system Population size CPU(s) GPU(s) Speed-up ratio
IEEE 33-bus 50 2.18 2.22 0.98
200 8.99 2.41 3.73
500 17.93 2.37 7.57
2000 63.28 2.53 25.01
5,000 180.69 3.37 53.62
IEEE 118-bus 50 30.29 31.72 0.95
200 110.65 30.93 3.58
500 310.47 33.01 9.41
2000 1341.32 32.77 40.93
5,000 2877.40 37.84 76.04
20,000 11,437.81 46.26 247.25

Furthermore, the IEEE 118-bus system was emulated to test the effect of the method and the speedup effect of calculating on the GPU platform in solving the optimization model of the larger scale distribution network. As shown in Figure 12, in the optimization model of the larger scale distribution network, increasing the population size can also improve the convergence speed of the algorithm. Consider the IEEE 118-bus system contains more buses than the IEEE 33-bus system, the maximal size of population is set to 20,000 as shown in Figure 12 to improve the convergence performance. The results of the population size under 20,000 is marked as the line with rhombus. When the population size is 20,000, binary (discrete) type variables of the result can converge to a unique solution, and the maximum standard deviation of the continuous variables is only 0.011. When the population size is 50, the binary (discrete) type variables cannot complete convergence, and the maximum standard deviation of continuous variables is 0.137. The results of control variables of CBs installed at bus 54, 62, 76, 105, and 113 are 11, 13, 4, 10, and 10, respectively. The output reactive power of D-STATCOMs installed at bus 25, 75, 85, 91, and 111 are 0.32MVar, 0.23MVar, 0.17MVar, 0.36MVar 0.18MVar, and 0.44MVar, respectively. Table 2 also shows in the IEEE 118-bus system, the calculation time for different population sizes on the CPU and GPU platforms. It can be seen that in the IEEE 118-bus system, the speedup effect was further improved, and the speedup ratio reached 247 times when the population size was 20,000. The results show that the effect of speedup in solving the optimization model of a larger scale distribution network will be more powerful.

Comparison of consequence speeds with different population size in IEEE 118-bus system.

Finally, Table 3 shows the bus voltage comparison of the critical loads before and after optimization. By dispatching the voltage regulation devices, the voltages of the critical loads are increased, which guarantees the stable and secure operation of the critical loads. Meanwhile, the active power loss of the system was reduced. The active power loss of the IEEE 33-bus system was reduced from 0.2056 to 0.1353 MW, which is a decrease of 34.14%. The active power loss of the IEEE 118-bus system was reduced from 1.289 to 0.9236 MW, which is a decrease of 28.34%.

The voltage comparison of the critical loads before and after optimization.

Test system Critical load node Node voltage before optimization (p.u.) Node voltage after optimization (p.u.)
IEEE 33-bus 11 0.928 0.954
15 0.917 0.944
21 0.992 0.994
27 0.945 0.962
33 0.917 0.948
IEEE 118-bus 35 0.934 0.955
51 0.915 0.944
71 0.884 0.911
95 0.946 0.954
111 0.905 0.932
Conclusion

This paper proposes an optimization model for distribution networks to dispatch controllable operational resources, which aim to improve the resilience of the distribution networks after extreme natural events. A GPU-based parallelized GA method was proposed to solve the proposed MINLP type model for the operations. A large population size of the GA, and setting the value of the distributed index by two stages can improve the convergence quality and consequence speed of the GA. Meanwhile, the GPU can also accelerate solving of the fitness function of the GA method with large populations, which effectively reduces the excessive calculation time. Finally, case studies with IEEE 33-bus and 118-bus were carried out to show that the proposed GA method can be solved effectively and accelerated to enhance the resilience of distribution networks. The GPU platform can improve the calculation performance and facilitate GA algorithms to be applied in practical projects. Further work will consider adding the electrical vehicles to the optimization model, whose limited energy needs to be allocated on multiple time scales to achieve the best optimization results.

 Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

 Author Contributions

LL-Original paper draft, PI of the support funding CL-Case studies, GA implementation.

 Funding

This research was funded by National Natural Science Foundation of China, grant number 52077045 and Open Fund of State Key Laboratory of Operation and Control of Renewable Energy and Storage Systems (China Electric Power Research Institute, DG80-20-001).

 Conflict of Interest

Author CL was employed by the company Guangzhou Power Supply Bureau of Guangdong Power Grid Corporation.

The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Thanks to the discussions with Dr. Xuewei Pan and his suggestions on this paper.

 References Baran M. Wu F. F. (1989). Optimal Sizing of Capacitors Placed on a Radial Distribution System. IEEE Trans. Power Deliv. 4 (1), 735–743. 10.1109/61.19266 de la Calle F. J. Bulnes F. G. Garcia D. F. Usamentiaga R. Molleda J. (2015). A Parallel Genetic Algorithm for Configuring Defect Detection Methods. IEEE Latin Am. Trans. 13 (5), 1462–1468. 10.1109/tla.2015.7112003 Deeb N. Shahidehpour S. M. (1990). Linear Reactive Power Optimization in a Large Power Network Using the Decomposition Approach. IEEE Trans. Power Syst. 5 (2), 428–438. 10.1109/59.54549 Delfanti M. Granelli G. P. Marannino P. Montagna M. (2000). Optimal Capacitor Placement Using Deterministic and Genetic Algorithms. IEEE Trans. Power Syst. 15 (3), 1041–1046. 10.1109/59.871731 Enacheanu B. Raison B. Caire R. Devaux O. Bienia W. HadjSaid N. (2008). Radial Network Reconfiguration Using Genetic Algorithm Based on the Matroid Theory. IEEE Trans. Power Syst. 23 (1), 186–195. 10.1109/TPWRS.2007.913303 Gao H. Chen Y. Xu Y. Liu C.-C. (2016). Resilience-Oriented Critical Load Restoration Using Microgrids in Distribution Systems. IEEE Trans. Smart Grid 7 (6), 2837–2848. 10.1109/TSG.2016.2550625 Gonzalez Bulnes F. Usamentiaga R. Fernando Garcia D. Molleda J. (2013). A Parallel Genetic Algorithm for Optimizing an Industrial Inspection System. IEEE Latin Am. Trans. 11 (6), 1338–1343. 10.1109/tla.2013.6710381 Grudinin N. (1998). Reactive Power Optimization Using Successive Quadratic Programming Method. IEEE Trans. Power Syst. 13 (4), 1219–1225. 10.1109/59.736232 Harik G. Cantú-Paz E. Goldberg D. E. Miller B. L. (1999). The Gambler's Ruin Problem, Genetic Algorithms, and the Sizing of Populations. Evol. Comput. 7 (3), 231–253. 10.1162/evco.1999.7.3.231 Jaros J. (2012). Multi-GPU Island-Based Genetic Algorithm for Solving the Knapsack Problem. In 2012 IEEE Congress on Evolutionary Computation. IEEE, 1–8. 10.1109/cec.2012.6256131 Luo W. Zhang Z. Wen T. Li C. Luo Z. (2017). Features Extraction and Multi-Classification of sEMG Using a GPU-Accelerated GA/MLP Hybrid Algorithm. Xst 25 (2), 273–286. 10.3233/xst-17259 Mahmud M. A. Hossain M. J. Pota H. R. (2014). Voltage Variation on Distribution Networks with Distributed Generation: Worst Case Scenario. IEEE Syst. J. 8 (4), 1096–1103. 10.1109/jsyst.2013.2265176 Mamandur K. R. C. Chenoweth R. D. (1981). Optimal Control of Reactive Power Flow for Improvements in Voltage Profiles and for Real Power Loss Minimization. IEEE Power Eng. Rev.PER- PER-1 (7), 29–30. 10.1109/mper.1981.5511679 Mohamed M. A. Chen T. Su W. Jin T. (2019). Proactive Resilience of Power Systems against Natural Disasters: A Literature Review. IEEE Access 7, 163778–163795. 10.1109/access.2019.2952362 Nicholson H. H. Sterling M. (1973). Optimum Dispatch of Active and Reactive Generation by Quadratic Programming. IEEE Trans. Power Apparatus Syst. PAS-92 (2), 644–654. 10.1109/tpas.1973.293768 Panteli M. Mancarella P. (2015). Influence of Extreme Weather and Climate Change on the Resilience of Power Systems: Impacts and Possible Mitigation Strategies. Electric Power Syst. Res. 127, 259–270. 10.1016/j.epsr.2015.06.012 Panteli M. Mancarella P. (2015). The Grid: Stronger, Bigger, Smarter?: Presenting a Conceptual Framework of Power System Resilience. IEEE Power Energ. Mag. 13 (3), 58–66. 10.1109/MPE.2015.2397334 Queiroz L. M. O. Lyra C. (2009). Adaptive Hybrid Genetic Algorithm for Technical Loss Reduction in Distribution Networks under Variable Demands. IEEE Trans. Power Syst. 24 (1), 445–453. 10.1109/TPWRS.2008.2009488 Swarup K. S. Yamashiro S. (2002). Unit Commitment Solution Methodology Using Genetic Algorithm. IEEE Trans. Power Syst. 17 (1), 87–91. 10.1109/59.982197 Topa T. Karwowski A. Noga A. (2011). Using GPU with CUDA to Accelerate MoM-Based Electromagnetic Simulation of Wire-Grid Models. Antennas Wirel. Propag. Lett. 10, 342–345. 10.1109/LAWP.2011.2144557 Topa T. Noga A. Karwowski A. (2011). Adapting MoM with RWG Basis Functions to GPU Technology Using CUDA. Antennas Wirel. Propag. Lett. 10, 480–483. 10.1109/LAWP.2011.2154373 Tsutsui S. Fujimoto N. (2009). Solving Quadratic Assignment Problems by Genetic Algorithms with GPU Computation: a Case Study. In Proceedings of the 11th Annual Conference Companion on Genetic and Evolutionary Computation Conference: Late Breaking Papers, 2523–2530. Utkarsh K. Ding F. Jin X. Blonsky M. Padullaparti H. Balamurugan S. P. (2021). A Network-Aware Distributed Energy Resource Aggregation Framework for Flexible, Cost-Optimal, and Resilient Operation. IEEE Trans. Smart Grid, 1. (early acess). Wang P. Liang F. Song J. Jiang N. Zhang X. Guo L. (2020). Impact of the PV Location in Distribution Networks on Network Power Losses and Voltage Swings with PSO Analysis. CSEE J. Power Energ. Syst., 1. (early acess). Wang Y. Chen C. Wang J. Baldick R. (2016). Research on Resilience of Power Systems under Natural Disasters-A Review. IEEE Trans. Power Syst. 31 (2), 1604–1613. 10.1109/tpwrs.2015.2429656 Xing Chen Xing. Xiao-Bang Xu Kama. Xu Xiao-Bang. (2005). Automated Design of a Three-Dimensional Fishbone Antenna Using Parallel Genetic Algorithm and NEC. Antennas Wirel. Propag. Lett. 4, 425–428. 10.1109/lawp.2005.859384