1. Introduction

Distributed Generation (DG) and Energy Storage (ES) are two technologies that fall under the concept of Distributed Energy Resources (DER) ^[1]. Both technologies can support voltage, reduce losses, provide backup power and ancillary services, improve local power quality and reliability, and defer distribution system upgrade ^[2,3,4,5]. Furthermore, the combined operation of DG and ES can enhance their potential benefits to the grid; for example, ES units can be charged when there is a generation excess, and provide support to the grid when generation decreases ^[6].

Modelling of renewable generation raises several challenges to distribution load flow calculations since capabilities for representing intermittent generators, voltage-control equipment, or multi-phase unbalanced systems are required. In addition, the study of systems with intermittent non-dispatchable resources will usually require a probabilistic approach and calculations performed over an arbitrary time period that may range from minutes to years. The evaluation of ES presents additional challenges: ES operation can be conducted from the perspective of a utility (i.e., aiming to improve network conditions) or operated by a private operator (i.e., with the objective of generating profits) ^[7]. In both cases an optimal operation strategy must be defined in order to maximize the benefits provided by ES; see ^{[8,9,10,11,12,13,14,15,16,17,18]}.

This paper uses a General Parallel Genetic Algorithm (GPGA) to estimate the optimal ratings of multiple fixed-location ES units and photovoltaic (PV) generators in order to minimize distribution system energy losses, annual energy and peak power supplied by the substation transformer. The goal is to explore the potential improvement on the overall operating conditions of the distribution system when ES units and PV generators are operated from the perspective of a utility. The developed methodology should be adequate for performing a technical evaluation of DG and ES as an alternative within the utility's planning scheme.

The optimal sizing of ES units has been conducted from different perspectives; e.g., improving network conditions ^[19,20], minimizing power and operation costs ^[21,22,23], meeting a predefined operation schedule ^[24,25,26], and multi-objective optimization ^[27,28]. Furthermore, various optimization techniques, such as particle swarm optimization ^[20,23], linear, dynamic, and cone programming ^[22,24,27], stochastic optimization ^[21,25], and genetic algorithms ^[26,28] have been employed.

The optimization of a combination of generation and energy storage has also been the subject of several recent works; see ^{[29,30,31,32,33]}. In ^[29,30], and ^[31] the optimization is conducted to minimize operation costs in a microgrid, while ^[32] presents a techno-economical procedure for sizing DG and ES in smart households. Reference ^[33] explores the optimization of pumped hydro storage that is operated in conjunction with wind generation to minimize operation costs or maximize generation and storage penetration. As with optimal sizing of ES, different optimization algorithms have been used; e.g., differential evolutionary algorithms ^[30,31], mix-integer linear programming ^[32], and genetic algorithms ^[33]. Furthermore, references ^[22,26,27] make use of clustering techniques in order to reduce the number of scenarios considered for optimization.

A strategy for operating ES units based on the variation of load and photovoltaic generation over the evaluation period (i.e., one year) is proposed here. This strategy is simultaneously optimized along with the ES units and PV generators ratings and it assumes that the operation curve is the same for all ES units while the days of the evaluated period are grouped according to different classification methods in order to reduce the dimension of the optimization problem. ES operation curves are validated through the simulation of the system under study for the selected evaluation period; only those curves that produce grid compliant solutions are considered as valid operation strategies.

The simulation tools used for this purpose are MATLAB and OpenDSS, a distribution system simulator whose capabilities allow users to represent the most important distribution components and perform multiphase calculations. This environment takes advantage of OpenDSS capabilities, which can be expanded with an external link to other tools, in this case MATLAB ^[34,35]. The use of OpenDSS as a simulation tool can provide an accurate estimation of the network conditions and avoid the use of simplification techniques (such as DC power flow or linearization of the power flow equations).

System behavior over the evaluation period is defined by node load and PV generation yearly curves which have been synthetically generated with different applications previously developed by the authors; see ^[36]. New applications based on clustering techniques have been implemented in this work for classifying and grouping the days of the evaluation period. This work was carried out using a modified version of the built-in OpenDSS storage model ^[37].

The paper has been organized as follows. The main characteristics of the ES model are summarized in Section 2. The procedure followed to obtain the operation curves of ES units using different day-classification methods is presented in Section 3. Section 4 details the GPGA procedure implemented for this work. The application to an introductory case study is presented in Section 5, while the optimization of multiple ES units and PV generators is analyzed in Section 6. A discussion regarding the procedure's main points and future work is presented in Section 7. Finally, the main conclusions drawn from this study are summarized in Section 8.

---

2. Modeling of energy storage units

OpenDSS provides a technology-independent ES model that allows users to perform planning studies under different operation modes ^[37]: Peak Shaving, Follow, Support, Load Level, and Time. Additionally, different charge and discharge strategies can be combined in order to produce the desired effect. The use of this ES model allows the user to set limits for charge/discharge power and minimum/maximum stored energy; moreover, it automatically calculates ES losses and updates the stored energy value.

The OpenDSS ES model can be used in the different solution modes available. Furthermore, the StorageController object allows the user to apply the same control strategy to a fleet of ES units.

To adequate the OpenDSS storage model to the optimization process proposed in this work, new features have been added to its native capabilities:

• ES efficiency is not assumed constant but defined by a function that depends on ES loading with respect to its rated power (in kVA).

• A minimum charge/discharge power as a percentage of the ES rated power (in kVA) has been defined. Regardless of the dispatched power, the ES unit will not absorb/inject power unless this value is greater than the minimum charge/discharge power.

• The model defines a maximum charge/discharge power based on the stored energy and simulation time-step. The objective is aimed at ensuring that the minimum/maximum stored energy values are not violated. For example, assume the present stored energy is 100 kWh and the minimum acceptable value is 80 kWh; given a 1-hour time-step and a discharge efficiency of 100%, the ES unit should only be able to inject a power of 20 kW for a period of 1-hour. Regardless of the predefined dispatched power, the model will enforce the 20 kW limit.

---

3. Operation of energy storage units

---

3.1. Determination of energy storage operation curves

The operation of multiple ES units in a power distribution system can be optimized by establishing coordination among all units; that is, it is advisable for all ES units to follow the same control strategy and dispatch in order to ensure total cooperation since individual operation strategies based on local measurements may not accurately reflect the global system requirements ^[34,37].

An ES dispatch strategy based on the principles of the Peak Shaving and Load Level operation modes has been developed for this work; it is aimed at minimizing daily power peaks and flattening the active power profile served by the substation transformer. This work assumes that the same operation curve is used for all ES units and it is derived from the power curve that would be measured at the substation terminals with DG but without ES.

Two limits are defined for charge and discharge operations: ES discharge can take place if the served power is above the discharge limit, whereas ES charge will occur if the served power is below the charge limit. In both cases the behavior will mirror the daily power curve pattern (see Figure 1). Additionally, charge and discharge curves will be modified by applying pertinent correction factors.

The operation curve of ES units for a specific day is built as follows:

(1) Obtain the daily power curve served by the substation transformer, including PV generation (see Figure 1a).

(2) Calculate the mean value of this daily power curve.

(3) Subtract the mean value from the daily power curve (see Figure 1b).

(4) Set power limits for ES charge/discharge (see Figure 1b).

(5) Find the power values below and above the charge/discharge limits.

(6) Subtract average charge/discharge power limit from the power values found in the previous step (see Figure 1c).

(7) Apply correction factors to the new power values found in step 6.

(8) The corrected power curve is normalized using the total rated ES power connected to the system. Negative values will correspond to ES charge, while ES discharge will be performed following the positive values in the normalized curve (see Figure 1d).

With this approach four parameters are required to determine the operation curve to be followed by all ES units: charge power limit, discharge power limit, charge curve correction factor, and discharge curve correction factor. Charge and discharge power limits are defined as a multiplier of the daily power's standard deviation; for the example presented in Figure 1 charge and discharge limits were set as –0.5 and 0.8 times the standard deviation, respectively.

---

3.2. Classification of days using hierarchical clustering techniques

The operation of ES units is specified in this study for an evaluation period of one year. Since using only one set of parameters for constructing the yearly operation curve is not enough to guarantee an optimal operation and defining a set of parameters for every day of the year leads to a very high number of parameters, different methods for classifying days have been implemented.

A hierarchical clustering algorithm groups a collection of individuals by forming a cluster tree ^[38]. The cluster tree does not separate the individuals into clusters but rather shows the connections among them at several hierarchical levels.

Hierarchical clustering is conducted using the following steps:

(1) Define a similarity measurement (e.g., Euclidean distance).

(2) Calculate similarity for every pair of individuals.

(3) Build the cluster tree.

(4) Form clusters by cutting the cluster tree. The number of clusters can be defined as a fixed parameter or as a height h at which the cluster tree is cut.

The similarity for every pair of individuals is expressed as a row vector of length N × (N–1)/2, where N is the number of individuals (points) in the collection. In this manner, the duplication of similarities is avoided. The cluster tree is built using the similarity vector as input. The algorithm links similar individuals into binary clusters (i.e., with two elements); moreover, these new binary clusters are linked to other binary clusters and individuals. The result is a multilevel tree, where all individuals are connected among them and also to the new binary clusters, see Figure 2a. Once the cluster tree has been built, it must be cut in order to group individuals into clusters. Depending on how the tree is separated, different clusters will be generated. In Figure 2b the cluster tree is cut at a height of 3000 producing 5 clusters, while in Figure 2c the total number of clusters has been previously chosen and is equal to 3.

Classification methods can be used in different applications. Reference ^[39] presented a two-stage pattern recognition of load curves, in which the first stage requires grouping daily load curves into clusters, using different clustering algorithms, and determining the optimal number of typical days (i.e., optimal number of clusters); clusters are constructed based on the Euclidean distance among curves. In ^[40] different similarity indices (including Dynamic Time Warping, DTW) were used as measurements for time series clustering in order to discover building energy patterns. Classification methods based on field indices (e.g., ratio of average to maximum daily load) are also used for categorizing the behavior of distribution customers ^[41].

Three different methods aimed at classifying the days of the evaluation period have been implemented in MATLAB using the Statistics and Machine Learning Toolbox.

(1) Daily served energy and PV generation: In this scenario the information related to the energy daily served by the substation transformer (without PV generation) and energy daily produced by the PV generator is gathered for every day of the year. Next, all days are assigned to a quartile according to the daily served and daily produced energy values. Each quantity (served and produced energy) produces four sub-groups; the combination of these sub-groups creates 16 groups. Finally, each day of the year is placed in one of the 16 groups based on the sub-groups to which they were originally assigned.

(2) Time-series clustering: This method is based on the similarity among daily power curves (including PV generation). Similarity among daily curves is calculated according to the Euclidean distance between two different curves. However, since a simple Euclidean distance may yield an incorrect evaluation of similarity; the minimum distance between two daily curves is calculated using the DTW algorithm ^[42], which grants a better assessment of similarity among the evaluated curves. After obtaining the distances among all daily curves, the day classification is conducted using MATLAB's linkage and cluster functions; the distances calculated with the DTW algorithm are used as inputs for the hierarchical clustering algorithm. The number of groups or clusters can be defined as a fixed parameter or as a height h at which the cluster tree is cut. In this work the number of clusters is defined as a fixed parameter.

Taking into account that the total number of operation parameters depends on the number of clusters, but a too low number of clusters may lead to a poor ES performance, the optimal number of clusters is determined by means of the so-called CH index as follows ^[43]:

$CH\left( K \right) = \frac{{B\left( K \right)}}{{W\left( K \right)}} \cdot \frac{{N - K}}{{K - 1}}$

(1)

$W\left( K \right) = \sum\nolimits_{k = 1}^K {\sum\nolimits_{i = 1}^{{N_k}} {{{\left\| {{X_i} - \overline {{X_k}} } \right\|}^2}} }$

(2)

$B\left( K \right) = \sum\nolimits_{k = 1}^K {{N_k}{{\left\| {\overline {{X_k}} - \bar X} \right\|}^2}}$

(3)

where K is the total number of clusters, N_k is the number of points that belong to cluster k, ${\overline {{X_k}} }$ is the average of points in cluster k, ${\bar X}$ is the overall average, and ${\left\| {\; \; } \right\|^2}$ represents the square of the Euclidean distance between two points. In order to make use of this index, hierarchical clustering must be evaluated for different numbers of clusters; afterwards the CH index is calculated for every evaluation. The optimal number of clusters will be that with the largest CH index.

(3) Daily-values clustering: For this scenario four measurements are calculated: the energy daily served by the substation transformer and the standard deviation of the daily power curve, with and without PV generation. Each day of the year is then defined by these four measurements and hierarchical clustering is performed based on the Euclidean distance among days. As in the Time-series clustering, the optimal number of clusters can be determined by means of the CH index.

The three methods rely on different system quantities: the first method uses a fixed number of groups, whereas the other two can be executed considering different numbers of classification groups. The Time-series approach was implemented in a similar manner to the approaches presented in ^[39] and ^[40], while the other two approaches could be seen as classification methods based on field indices. With these classifications, all days belonging to the same group will share the same operation parameters so the total number of required parameters for defining ES operation will consequently be reduced. To obtain ES operation curves during the evaluation period (i.e., one year), the day-classification is first conducted (according to the selected method); next, a set of operation parameters is specified for every classification group: each day is first defined by its day-type (i.e., the group/cluster number to which it belongs) and then its daily operation curve is constructed according to the procedure presented above. Figure 3 presents a flow chart with the complete process implemented to obtain the operation curve of each ES unit.

---

4. General parallel genetic algorithm

A genetic algorithm (GA) is an optimization technique that evaluates more than one area of the search space and can discover more than one solution to a problem ^[44]. The application of a GA starts with the evaluation of a randomly-generated initial population. After the first evaluation, the fittest individuals are selected to become parents for a new generation, where different methods can be applied for the parent-selection process (e.g., fitness proportionate selection, tournament selection, truncation selection). The crossover genetic operator is used to mate the selected parents in order to create the members of a new generation; finally, the mutation operator introduces a random variation to the members of the new generation.

A complete GA procedure must include the following characteristics ^[45]: (ⅰ) a criterion to create the initial population; (ⅱ) an evaluation function to obtain the fitness of all possible solutions; (ⅲ) a genetic operator to obtain the next generation; (ⅳ) a stopping criterion.

GAs are prone to parallelization, depending on the approach used to evaluate the fitness function and how the genetic operators are applied. Several types of General Parallel Genetic Algorithms (GPGA) have been proposed: global single-population master-slave GAs, single-population fine-grained, multiple-population coarse-grained GAs ^[46]. In global single-population GAs only one population is created and the evaluation of each member of a given generation is independent from the rest of the members; therefore, the evaluation of a generation can be seen as an embarrassingly parallel problem. This approach yields a reduction of the total execution times by distributing the evaluation of individuals in a generation among available cores in a multi-core installation.

GAs have been applied to a wide range of optimization problems related power systems (e.g., ^{[26,28,33,46,47,48]}).

A global single-population GPGA has been implemented in this work to determine the optimal values of the following inputs:

• Nominal PV generator power rating (in kVA).

• Nominal ES unit power rating (in kVA).

• Nominal ES unit energy rating (in kWh).

• Security factor for interconnection transformer (in pu).

• ES unit operation strategy parameters.

The number of inputs in the GPGA will be equal to:

${N_{IN}} = 4 \times {N_U} + 4 \times K$

(4)

where N_IN is the number of inputs to be optimized by the GPGA, K is the number of day groups (or clusters), and N_U is the number of PV generators and ES units under consideration.

The parameters used in the optimization procedure are shown in Table 1.

The GPGA seeks to minimize the following operational variables:

• distribution system energy losses, without including losses in substation transformer;

• annual energy supplied from the MV side of the substation transformer;

• peak power supplied by the substation transformer and measured at MV-side terminals;

• standard deviation of the power values measured at the MV-side of the substation.

Two stopping criteria have been established:

(1) The generation count (i.e., the number of executed generations) is equal to the established maximum number of generations.

(2) The best fitness changes less than ε over δ generations.

The GPGA has been implemented in MATLAB and the evaluation of each member in a population is conducted by simulating the test system during one year using a 1-hour time-step. Individual executions are distributed among 60 cores using the library developed by M. Buehren ^[49]. For more details about the multi-core installation, see ^[50]. The GPGA complete general procedure is depicted in Figure 4.

The fitness function calculated for every evaluation will be equal to the Euclidean distance obtained from the reduction of the operational variables:

$Fitness = \sqrt {{{\left( {1 + \Delta {P_L}} \right)}^2} + {{\left( {1 + \Delta {P_{MAX}}} \right)}^2} + {{\left( {1 + \Delta {P_{STD}}} \right)}^2} + {{\left( {1 + \Delta E} \right)}^2} + {{\left( {\Delta {R_{PF}}} \right)}^2}}$

(5)

$\Delta {R_{PF}} = \left\{ {\begin{array}{*{20}{c}} {1,{P_{N\_MIN}} > 0}\\ {1 + \frac{{{P_{N\_MIN}}}}{{{P_{MAX}}}},{P_{N\_MIN}} \le 0} \end{array}} \right.$

(6)

where ΔP_L is the reduction of energy losses, ΔP_MAX is the reduction in the peak power supplied by the substation transformer, ΔP_STD is the reduction of the standard deviation of the yearly power measured at the MV terminals of the substation, ΔE is the reduction of the annual active energy supplied by the substation transformer, P_{N_MIN} is the minimum yearly power delivered by the substation transformer, and P_MAX is the maximum yearly power in the base case (i.e., without PV generation or energy storage). ΔR_PF has been introduced in the fitness function as penalty for those cases in which reverse power flow occurs. All reductions are calculated in pu with respect to the case without generation and energy storage.

The procedure adds "1" to the reduction values in order to avoid any inconsistencies due to the appearance of negative reductions. For every evaluation, the procedure checks that all voltages are within accepted limits (between 0.95 and 1.05 pu), no network element experiences overload, and the iterative algorithm converges for all solution steps. If one or more of these conditions are not fulfilled, the procedure will assign a fitness value close to zero, which will difficult its participation in the crossover stage of the GPGA.

The following aspects must be taken into account when executing the optimization procedure:

• The resulting ratings for PV generators and ES units will be rounded in steps of 5 kVA, while the rating of the interconnection transformers will be rounded in steps of 10 kVA.

• PV generators and ES units only inject active power, and ES units only absorb active power.

• PV generators are allowed to produce as much energy as possible. Moreover, the PV generator model implemented in OpenDSS assumes that the generator is capable of reaching the Maximum Power Point (MPP) for every simulation step.

• ES units absorb/inject power only when the dispatched power is above 10% of the rated power; in the same manner, the PV generators only inject power if the generated power is above 10% of the rated power.

• The minimum stored energy for ES units is equal to 20% of the ES rated energy.

• Although the procedure relies on a technology-independent ES model, the present application and characteristics of the ES unit (i.e., rated power, rated energy, and energy to power ratio) indicate that Battery Systems are the technology that would best fit a real-life implementation.

• The ES operation curve (see Section 3) is constructed from the active power profile served by the substation transformer when PV generation is present. This curve depends on the actual PV generation and cannot be estimated before the optimization procedure is executed. Therefore, in order to construct the operation curve, for every member of a generation the test system must be first evaluated taking only into account the presence of PV generation; the resulting active power profile will be used by the procedure to construct the yearly ES operation curve.

• The day-classification methods also rely on the active power curve served by the substation transformer, including PV generation. Although a similar approach to that presented in the previous bullet could be used, this could lead to discrepancies in the procedure. The system quantities used as criterion for day-classification depend on the actual value of PV generation and different ratings for PV generators could generate different classifications. If each generation member was to produce different day-classifications, the optimization of the operation parameters would be incorrect because the set of parameters assigned to a specific cluster would not be applied to the same groups of days. Furthermore, the CH index could not be used to predefine the optimal number of clusters to be used in the optimization procedure, since each day-classification method could have a different optimal number of clusters. In order to circumvent this inconvenient, a preliminary study will be conducted; in this study, only the sizing of PV generators will be optimized, while the fitness function is defined by the following equation:

$Fitnes{s_{PV}} = \sqrt {{{\left( {1 + \Delta {P_L}} \right)}^2} + {{\left( {1 + \Delta {P_{MAX}}} \right)}^2} + {{\left( {1 + \Delta E} \right)}^2} + {{\left( {\Delta {R_{PF}}} \right)}^2}}$

(7)

The GPGA will use the parameters presented in Table 1; however, the security factor of the interconnection transformer will be constant and equal to 1.01. The day-classification will be conducted using the optimized PV generator ratings found by this study.

• The maximum rated power for a connection node and the initial values for PV generator and ES unit power ratings are generated according to rules presented in ^[51].

• The term "with/without PV generation" makes reference to the PV generation to be defined by the procedure. Should the system under study have previously installed DG, it would be considered as part of the original design and not included in the optimization. Moreover, the active power profile "with and without PV generation" is obtained through a yearly simulation using the synthetically generated system curves.

• The locations of the PV generators and ES units have been arbitrarily chosen. Although the locations could have been derived from a procedure aimed at optimizing the impact on the grid, there are several non-technical factors that can have a greater impact on selecting the DER locations than the results obtained from such procedure. Land price, environmental permits, availability of space, and the so-called "not-in-my-backyard" stance are some of the aspects that can limit the options for the location of PV generators and ES units.

---

5. Case study 1: optimization of a single energy storage unit

---

5.1. Test system configuration

The test system is a three-phase 60-Hz overhead system serving two spot loads from a High-Voltage/Medium-Voltage (HV/MV) substation transformer; see Figure 5. The PV generator and ES unit are connected to the MV network through a distribution Medium-Voltage/Low-Voltage (MV/LV) transformer.

The main characteristics and operating conditions are as follows:

• HV rated voltage: 230 kV

• HV actual voltage: 1.05 pu

• MV rated voltage: 4.16 kV

• LV rated voltage: 0.4 kV

• Substation transformer rated power: 1000 kVA

• Total rated active load: 550 kW.

System loads have been modelled as voltage-independent and characterized by curve shapes derived with the procedure presented in ^[36]. The efficiency curve used for the ES unit is shown in Figure 6. The same curve has been assumed for both charge and discharge cycles.

---

5.2. Simulation results

The GPGA procedure was executed for three different scenarios based on the day-classification methods proposed in Section 3 and using the values shown in Table 1. Operational variables without ES and PV generators are as follows:

• Annual energy: 1877440.1 kWh

• Yearly energy losses: 40209.7 kWh

• Maximum yearly power: 427.0 kW

• Yearly power standard deviation: 51.2 kW.

Day-classification was conducted using a 60-kW PV generator (obtained during the initial study). A sample day-classification is shown in Figure 7, where a bar graph is used to show the cluster/group number to which each day of the year is assigned; that is, all days whose bars present the same height have been placed in the same cluster/group and thus, share the same parameters used in the construction of the ES operation curve. As expected, the three classification methods produce different graphs. For the Time-series and Daily-values classification methods the optimization procedure was carried out considering different numbers of clusters. In all scenarios the GPGA seeks to optimize all 4 operational variables (see Section 4).

The results are presented in Table 2. It can be observed that all scenarios converge due to the δ/ε criterion and that best fitness values are consistent throughout the executions. Although operational variables exhibit some variations (based on the scenario under evaluation), they are sufficiently small to consider them negligible and regard the results as acceptable. Except for yearly energy losses (on which the optimization procedure seems to have a very limited effect), there are significant reductions of all operational variables with respect to the case without generation and storage.

The optimal values for ratings show larger variations than those exhibited by the operational variables; however, in general these values appear to be bound to a specific range, so they can be used as a reference for choosing PV generator and ES unit ratings in real-life applications. The executions based on the Time-series and Daily-values classification methods show a dependency on the number of clusters used for the day-classification; better results are obtained with a smaller number of clusters for both methods. Under these circumstances the CH index could be used to calculate the optimal number of clusters; this approach could prevent the repeated execution of the optimization procedure with different number of clusters for the day-classification.

The calculation of the CH index has determined that the optimal number of clusters for the Time-series and Daily-values Clustering methods are 5 and 16 clusters, respectively (the CH index was calculated for a number of clusters ranging from 2 to 24). The results obtained when executing the optimization procedure with the optimal number of clusters for the Time-series and Daily-values classification methods are shown in Table 3. It can be observed that the Time-series method produced a higher value for best fitness with respect to those presented in Table 2 (although energy losses exhibit a slight increment), whereas the results obtained with the Daily-values method show no improvement; in fact, the day-classification with 6 clusters produced better results.

Overall, the highest fitness value was obtained using the Daily-values method with 6 clusters; however, the difference with respect to the highest valued found with the Time-series method is negligible; moreover, the Time-series method required only half of the execution time.

This behavior shows that the Time-series classification method can be used in combination with an optimal number of clusters derived from the calculation of the CH index in order to obtain adequate results from the GPGA. This approach requires the lowest number of clusters, which implies a minimization of the dimension of the optimization problem and a reduction in total execution times.

The optimal operation parameters for the Time-series classification method using 5 clusters are shown in Table 4. Figure 8 depicts the behavior of the yearly active power curve obtained following these parameters. It can be observed that PV generation alone is not capable of reducing the yearly maximum power but due to its effect the power curve presents lower minimum values and also exhibits larger variations, see Figure 8b. The introduction of ES manages to reduce the yearly power peak and produce a flatter profile than the curves depicted in Figure 8a and 8b; however, the yearly curve still presents power spikes, see Figure 8c.

---

5.3. Discussion

One important aspect that should not be ignored is the effect that the fitness function can have on the GPGA performance and how the devices affect the operational variables. PV generation can help to reduce both the annual energy provided by the substation transformer and the yearly energy losses; the peak power occurs at a time when PV generation is zero and the yearly power standard deviation is actually increased due to the presence of PV generation.

As for ES, operation parameters will determine how the operational variables will be affected by ES operation. The operation band defined by the charge and discharge limits will delimit the expected standard deviation by constraining the active power profile to this band; this action will also help to reduce the maximum peak power delivered by the substation transformer. On the other hand, total energy losses will be a result of the interaction between the ES charge (loss increment) and discharge (loss reduction) states. The ES operation curve mirrors the behavior of the active power profile; through this mirrored behavior the ES unit will adjust its power injection/absorption to the system's loading. However, the actual power value that needs to be generated/absorbed in order to produce an overall improvement in system conditions cannot be pre-assumed. Therefore, the initial operation curve defined by the charge and discharge limits must corrected in order to produce the greatest impact on the system: the charge correction factor will ensure that the ES unit has enough energy to meet the pre-dispatched discharge without producing an unacceptable increment in energy losses, whereas the discharge correction factor will help to improve the overall system conditions (standard deviation, peak power, and energy losses) while balancing the use of energy so it is not depleted too quickly. In general, the interactions among the different operation parameters will produce much larger reductions for the yearly maximum power and yearly power standard deviation than for energy losses; therefore, the procedure will drive the optimal operation towards benefiting these operational variables.

Defining a fitness function that accurately reflects the objectives behind the optimization procedure is a critical task. Furthermore, the interconnection transformers represent a new source of losses (i.e., no-load and load losses), which in some cases cannot be compensated by the reduction in conduction losses (i.e., losses in distribution lines) and can limit the connection of larger PV generators and ES units.

The resulting ES operation curves correspond to a system behavior under known load and generator curves; therefore, they cannot be used for actual day-to-day operation. However, the optimal operation parameters found by the GPGA can be used to define a preliminary operation curve based on day-ahead load and generation forecast. For this purpose, the distribution system must be simulated using the forecasted load and generation curves; then the resulting active power profile will be used to construct a preliminary ES operation curve by applying the corresponding operation parameters. The preliminary operation curve can be used as a base to create a set of operation curves (by introducing random variations) that can serve as the initial population of a GPGA aimed at determining the optimal daily operation curve by optimizing every single ES operation point. Using a set of curves based on the preliminary operation curve should speed-up algorithm convergence when compared to a randomly generated initial population.

---

6. Case study 2: optimization of multiple energy storage units

---

6.1. Test system configuration

A new test system has been studied in order to test the optimization of multiple ES units; the optimization procedure presented in the previous sections was implemented to be independent of test system size as well as of the number of PV generators and ES units under evaluation.

The new test system is a three-phase 60-Hz overhead system with a simplified representation of the HV transmission system (see Figure 9). ES units and PV generators have the same configuration and characteristics that in the previous test system; see Figures 5 and 6. By default, all distribution loads are connected to the system through MV/LV transformers.

Some important information related to the new test system follows:

• HV rated voltage: 230 kV

• HV actual voltage: 1.03 pu

• MV rated voltage: 4.16 kV

• LV rated voltage: 0.4 kV

• Substation transformer rated power: 2500 kVA

• Total rated active load: 1700 kW.

Operational variables for the new test system (without ES and PV generators) are as follows:

• Annual energy: 7459245.8 kWh

• Yearly energy losses: 199519.6 kWh

• Maximum yearly power: 1473.8 kW

• Yearly standard deviation: 183.9 kW.

As in Test System 1, loads have been modelled as voltage-independent and curve shapes have been derived with the procedure presented in ^[36]. Day-classification was conducted with the following rated powers for PV generators 1–4: 115, 40,115, and 270 kW.

---

6.2. Simulation results

Simulation results for all three scenarios are presented in Table 5; for the Time-series and Daily-values methods, the optimal number of clusters derived from the CH index (2 and 3 clusters, respectively) were used for the ES operation, see Figure 10.

Note that optimal numbers of clusters are much lower than those obtained with the previous test system. This behavior is due to the smoothening effect produced by the aggregation of loads: Test System 1 only served two loads, which caused a more varying power profile, whereas the larger number of loads in the new test system generates a smoother profile. It can be observed that the Time-series classification produces the highest best fitness value and requires the shortest simulation times; these results are in line with those presented in the previous section. With the exception of energy losses, all classification methods were able to produce reductions in the operational variables when compared to the values without ES and PV generation; moreover, the Time-series and Daily-values methods produced better results than the Daily served energy and PV generation method. Consequently, the optimization procedure can be regarded as satisfactory, since the increment in energy losses is clearly outweighed by the improvement in the rest of the operational variables.

The optimal values for PV generators and ES units present large variations for the ratings of individual units; however, the Time-series and Daily-values classification methods present similar penetration factors (i.e., total installed power) for both generation and storage.

The results presented in ^[51] show that different combinations of rated powers and locations of PV generators can produce similar values of energy losses; therefore, it can also be expected that similar values for the operational variables can be obtained with different ratings of PV generators and ES units.

All GPGA executions converged due to the δ/ε criterion; however, since GA is a stochastic optimization procedure, it is possible for it to converge before reaching the optimal solution. If the test system size and number of clusters used in the ES operation mode result in affordable simulation times, the GPGA can be executed without the δ/ε criterion, that is, the procedure will stop when it has reached the maximum number of generations.

The optimization procedure was executed for 200,250, and 300 generations using the Time-series classification method with 2 clusters. Results in Table 6 show an improvement with respect to the values without ES and PV generation; moreover, the best fitness value seems to indicate that 200 generations are enough for obtaining optimal results. Note that operational variables present consistent values among executions. The optimal ratings of PV generators and ES units show clear variations among executions, whereas penetration factors exhibit a more constant trend, a behavior that has remained constant throughout all conducted studies. Consecutive executions of the procedure can be used to produce several combinations of rating values that can be selected in a real-life application or to explore the optimal penetration factor for PV generation and energy storage.

---

7. Discussion and future work

The main focus of most previous works on the main topics covered in this paper is the optimum sizing and operation of ES units within microgrids or small-size power systems, being the typical optimization period for ES operation 24 hours; larger evaluation periods are generally divided into a reduced number of representative days (usually through clustering techniques) and ES operation is optimized for each representative day. For this work, an evaluation period of one year with a time step of one hour has been considered. Since only through the optimization of every single operation point (i.e., every hour of the year) the global optimum can be reached, the problem dimension becomes very large, and even with the aid of a multi-core installation simulation times are prohibitive. Furthermore, a stochastic optimization method, such as GA, cannot guarantee the optimal solution for an optimization problem with such a large number of input parameters. Therefore, there must exist a tradeoff between optimality and the achievement of a feasible solution.

A comparison between the results derived from the developed methodology and other previous works has not been possible due to the difficulty (or impossibility) to adequate other procedures to the conditions established in the present work. The long term evaluation (i.e., one year) and system solution by means of a power flow simulator appear to be incompatible with the formulation used in other methodologies. On the other hand, assuming the restrictions imposed in other works (such as the evaluation period represented as a group of representative days or the introduction of simplifications into the solution of the electrical system) would defeat the purpose of a procedure that seeks to solve the problem through a detailed representation of the system under study and its behavior during the entire evaluation period.

Future work could be aimed at improving the criteria used in the classification of the evaluation period days, expanding the evaluation period (i.e., several years) and taking into account yearly load variation, refining the proposed operation mode, and defining a fitness function that best balances the improvements in each of the operational variables. Further work could also be based on an optimization procedure that does not rely on dimension reduction techniques and could cope with time steps shorter than one hour (e.g., 15 minutes).

---

8. Conclusions

This paper has presented the application of a General Parallel Genetic Algorithm to determine the optimal ratings of multiple fixed-location PV generators and ES units when they are operated from the utility's perspective. The evaluation is performed with OpenDSS, which allows for a detailed modeling of the distribution system and an accurate estimation of network conditions; that is, no simplifications are introduced into the solution of the electrical system, the evaluated system can be of any size, and detailed models and advanced control strategies algorithms for system components can be included.

The optimization of ES ratings cannot be done without considering an optimal dispatch strategy; the operation mode proposed in this work attempts to simplify the arduous task of optimizing the complete yearly ES dispatch. The implemented strategy reduces the number of variables to be optimized by using clustering techniques that group the days of the year into different clusters. However, unlike other works, each group of days is not independently optimized. Each possible solution is evaluated by simulating the test system during the entire evaluation period. This approach ensures that the evaluation period is assessed as a whole (not as a decoupled group of representative days) and the ES current state is the result of the complete ES operation history.

The introduction of parallel computing into the optimization procedure has allowed for affordable execution times. Although each individual simulation does not require large simulation times, the repeated evaluation of the system in a sequential manner would result in prohibitive execution times. Therefore, only optimization techniques that are prone to parallelization could be used in conjunction with a power flow simulator and a long term evaluation, since sequential techniques would be unfeasible due to the resulting execution times.

It is important to emphasize that the locations of PV generators and ES units were arbitrarily chosen in this work; however, one must keep in mind that the location of a power source (generation or consumption) can have a non-negligible impact on energy losses and system voltages.

The main features of the implemented methodology may be listed as follows: (ⅰ) system conditions are estimated by means of a power flow simulator, (ⅱ) the methodology is independent of system size and the number of units to be optimized, (ⅲ) ES operation is optimized for an evaluation period of one year, (ⅳ) the operation strategy of ES units is based on the daily variation of load and generation curves, (ⅴ) hierarchical clustering is used to group daily power curves in order to reduce the dimension of the optimization problem, (ⅵ) time-driven simulation ensures that the evaluation period is assessed as a whole and the ES current state is the result of the ES operation history, and (ⅶ) the parallel approach used in the GA implementation allows for affordable execution times.

---

Conflict of interest

All authors declare no conflict of interest in this paper.

---