Coordinated Operation of Energy Storage Systems for Distributed Harmonic Compensation in Microgrids

Energy storage systems (ESSs) bring various opportunities for a more reliable and flexible operation of microgrids (MGs). Among them, energy arbitrage and ancillary services are the most investigated application of ESSs. Furthermore, it has been shown that some other services could also be provided by ESSs such as power quality (PQ) improvements. This issue could be more challenging in MGs with widespread nonlinear loads injecting harmonic currents to the MG. In this paper, the feasibility of ESSs to act as coordinated active harmonic filters (AHF) for distributed compensation was investigated. An optimization model was proposed for the coordination of the harmonic compensation activities of ESSs. The model takes into account the various technical and systematic constraints to economically determine the required reference currents of various AHFs. Simulation cases showed the performance of the proposed model for enhancing the harmonic filtering capability of the MG, reduction in the compensation cost, and more flexibility of the distributed harmonic compensation schemes. It was also shown that ESS activities in harmonic compensation do not have much of an effect on the ESSs revenue from energy arbitrage. Hence, it could make ESSs more justifiable for use in MGs.


Introduction
Energy storage systems (ESSs) have received special attention due to their great flexibility and applicability in power systems. The smart charging and discharging behavior of ESSs could produce many advantages for various power system applications. The energy arbitrage for the operation of ESSs as consumers of low-price energy at off-peak periods and suppliers of high-price energy at peak hours is the most favorable application of ESSs. Furthermore, they could provide some advantages such as providing ancillary services and deferring investment costs in various sections of the power and distribution systems (DSs) [1][2][3][4]. All of these applications are related to the supply and security aspects of the energy. However, utilization of ESSs for the quality aspect of the energy for power quality (PQ) improvements has not gained much attention until now. With respect to the major costs that could be imposed by PQ related problems [5,6], it seems that investigations should be carried out regarding possible the PQ advantages of ESSs. Since the penetration of nonlinear power electronic (PE) based loads is increasing, it seems that ESSs could play an important role in enhancing PQ levels affected by these new loads in future DSs.
It should be noted that the distributed compensation scheme has better applicability to cope with the harmonics of the dispersed nonlinear loads of the MG. Furthermore, for reliable operation and unwanted stability problems, the coordinated operation of the AHFs should be ensured. Therefore, the proposed ESS-AHF model could successfully meet the requirements for distributed compensation and stability considerations.
The rest of this paper is organized as follows. In Section 2, active harmonic filters are introduced. Various ESS applications and characteristics are investigated in Section 3. A harmonic model of ESSs in harmonic compensation activities is also represented in this section. The proposed ESS-AHF model is introduced in Section 4. In Section 5, the results of the implementation of the introduced models are provided and analyzed. The paper then presents our concluding remarks in Section 6.

Active Harmonic Filters
Conventional mitigation methods of harmonics can be categorized as passive filters, static compensators, and active power line conditioners (APLCs). Currently, passive filters are the most prevalent technology for harmonic compensation. They consist of a single tune, second tune, and multi tune topologies. The passive filters are economic solutions and there is a great body of literature regarding their sizing and application for DSs [17]. Static compensators follow a principle similar to that of passive filters. However, these solutions have some limitations. Insufficient flexibility is the biggest disadvantage of passive solutions. Passive filters are usually designed for a specific harmonic condition, and hence change in the harmonic distortion or working conditions could degrade the applicability of the passive filters. Furthermore, resonance problems could affect the performance of passive filters and could increase harmonic distortion in some points of the network. Due to higher flexibility, there is an increasing approach in PE-based active harmonic filters (AHFs) such as APLCs. They could be optimal choices to cope with PQ related problems. In addition to the conventional APLCs, almost all energy resources with PE interfaces may be controlled to act, also, as harmonic active filters. Among them, ESSs equipped with PE converters could also be accounted for a possible AHF resource [18][19][20][21].
APLCs are considered as one of the most promising solutions for harmonic compensation activities. These devices do not suffer from conventional limitations of passive filters such as limited tuning, and especially resonance problems. They have the flexibility to work at a broad range of harmonics at any working point. Hence, APLCs are going to be more employed in electric networks [22][23][24]. An APLC is, in fact, a voltage or current source inverter that is controlled to produce harmonic compensation currents or voltages at the output. The APLC injects a harmonic compensation voltage or current to the network and eliminates harmonics in the source side. Unlike passive filters, APLC produces no resonance in the network. Hence, they are more prevalent, despite their higher investment costs and more complex control methods [22][23][24].
The schematics of two common types of APLCs called shunt and series topologies are represented in Figure 1. In both types, the voltage and current measurements of the network are processed by the APLC controller. The controller determines the required reference harmonic currents that should be injected by the filter. In [25,26], common control strategies for producing the desired reference current have been investigated. The calculated reference currents are then injected into the grid by a proper switching method. The required power for producing compensation currents can be supplied by an external source or taken by the grid using the same power converter. It has been shown that distributed compensation schemes could bring more flexibility and filtering capability for harmonic compensation [5,16]. Since the harmonic compensation is provided all over the grid, it is preferable to compensate current harmonics as close as possible to the disturbing source. In this way, harmonic losses can be reduced in the network. Harmonic loads are usually positioned at various locations with different nonlinear characteristics. Therefore, the distributed compensation scheme could have better applicability to cope with the harmonics of the dispersed nonlinear loads. Furthermore, for reliable operation of the MG, the coordinated operation of the AHFs should be ensured. It has been highlighted that instability problems can arise in the case of uncoordinated operation [5,16].
It should be noted that the penetration of nonlinear loads is increasing in MGs and compensation methods solely based on APLCs could be expensive and insufficient. Investment and operation costs are usually high for APLC devices [5]. Hence, other possible AHFs based on other PE converters should be utilized for more economic and flexible solutions. There is a great body of literature about the possible utilization of ESSs for PQ improvement [2,4,6,7,11,13,27]. However, to the best of the author's knowledge, the coordinated operation of ESSs for PQ problems has not been paid attention thus far. In this paper, the possibility of storage systems for the coordinated mitigation of harmonics is addressed. The coordinated operation of ESSs was achieved by an optimization problem including various constraints that could guarantee the economic and reliable operation of PQ improvement activities. In the next section, the possible utilization of ESSs for PQ improvement was investigated.

Energy Storage Systems for Power Quality Improvement
Energy storage systems have received special attention due to their great flexibility and applicability to the various power system sections. In MGs, an ESS converts electrical energy via a power interface, which is usually based on a PE converter. The general characteristics of ESSs are typically determined by their PE interfaces. Aside from PE interfaces, various management and control components are included in ESSs that enable them to have various functionalities in the network.
The main purpose of incorporating ESSs is to participate in energy markets. The term energy arbitrage refers to the charging of ESSs in the off-peak periods and discharging in the peak periods. The ESS revenues come from the difference between the off-peak and peak energy prices. Transmission and distribution services could also benefit from ESSs. The extension of network lines could be deferred using the flexibility of ESSs. This may avoid congestion problems in both the transmission and distribution grids, which result in more reliable operations of the system as well as reduced prices of energy [2]. Much work has been done on the integration of DERs in MGs. In MGs, the fluctuating output power of DERs in both grid-connected and islanded modes can cause voltage instability. This can also result in intermittent energy resources such as wind and solar farms. The ESSs could play a role for energy buffering and could smooth fluctuations of intermittent DER by a smart charging and discharging procedure. ESSs could then participate in voltage regulation and compensate reactive power to enhance the stability of the system [2].
One of the main applications of ESSs is in the area of ancillary services. Among various services, frequency regulation is the main service that could be supplied due to the fast response speed of PE interfaces. ESSs with high energy capacity could provide fast-spinning reserves. Furthermore, ESSs could act as AHFs and compensate for various PQ problems. Integrating ESSs in various DERs and flexible alternating current transmission system (FACTS) devices could enable this equipment for the compensation of extra reactive and harmonic currents. They could also avoid the PQ problems related to the application of DER in the MG.
ESSs vary in a broad range of technologies and categories. A common classification method is based on the form of stored energy. Mechanical (such as pumped hydro energy storage, compressed air energy storage (CAES) and flywheel), electrochemical (such as secondary batteries as the most useful technology), chemical (such as fuel cell), thermal (such as cryogenic energy storage and high-temperature thermal energy storage), and electrical (such as ultra-capacitors and superconducting energy storage systems (SMES)) are the main categories. Review papers about various ESSs technologies are available in [7,10]. In [10], the current practical implementations of ESSs all around the world were investigated.
Typical applications of ESSs are usually specified based on their properties in terms of power and energy capabilities. In [4], a comprehensive comparison was made among various ESSs in which output powers could be represented based on output duration, which is shown in Figure 2. Various battery types and mechanical energy storage systems such as pumped hydro and CAES with high energy density could be used in applications such a load-leveling and emergency power. Some other ESS types such as SMES and flywheels with conventional bearings have high power capacity. Hence, they are suitable for some PQ problems such as instantaneous voltage drop, voltage flicker, and short duration interruptions. Other ESSs such as flywheels with levitation bearing and super-capacitors can be employed for high power requests for a short time [4]. However, the duration depends on the requested power from the ESS. The general connection of some ESSs suitable for PQ mitigation activities is represented in Figure 3. The energy train may consist of a back-to-back converter for power conditioning between the MG and the main energy source. For example, flywheels work with AC power with their permanent magnet motors. Hence, they need an AC/DC converter for charging the DC link. Sometimes, a DC/DC converter can be employed before the DC link, as shown in Figure 1. After the DC link, a grid side converter (GSC) is used to deliver the energy from the DC link to the PCC. In order to enable the ESSs to participate in PQ mitigation activities, the control methods of the GSC should be modified.

Application of ESSs for Harmonic Compensation
Aside from AHFs and passive filters, ESSs could play an important role in the mitigation of PQ and harmonic problems [2,4,7,10]. The ability to control the PE converters of ESSs could enable them to act as harmonic active filters. In the literature, there are some examples of the modification of the control algorithms of PE converters for harmonic compensation. The battery storage system is currently the most promising ESS due to economic considerations and mature technology [27]. A smart battery controller was proposed in [11] for enhancing the PQ and adjusting steady-state voltage and frequency of MG. An investigation was made in [13] into the PQ of wind turbines with and without ESSs and the results showed the positive effects of ESSs as a solution for the mitigation of wind turbine harmonics.
Optimal harmonic compensation and the economic employment of ESSs could be achieved by the coordinated control of PE-based inverters [5]. It seems that the optimal operation of ESSs could efficiently consider various technical and systematic constraints in the calculation harmonic sharing of each ESS. Hence, the coordinated control of ESSs was proposed in this paper for harmonic filtering activities as an extension for various applications of ESSs. Therefore, the utilization of various ESSs for harmonic compensation and the coordinated control of ESSs for a global distributed harmonic compensation scheme were proposed and formulated. To this end, a harmonic model of ESSs is required in the coordination algorithm. Based on Figure 3, the equivalent circuit model of ESSs in harmonic mitigation can be represented as in Figure 4. The internal current source of the ESS is controlled based on the internal parameters of the ESS such as the voltage of the DC link, control, and switching method of the inverter. The GSC is controlled in the current control mode for harmonic compensation [27]. Internal impedance Z f shows the Thevenin impedance of the ESS seen in the output of the GSC. The internal impedance depends on the internal parameters, the control blocks, and the topology of the GSC. Some calculation methodologies of the GSC internal impedance as a VSI was investigated in [28][29][30]. A simple model for the Z f can be represented as Equation (1), which includes two resistive (R fh ) and inductive parts (X fh ). Both parts are considered as functions of various internal parameters such as the voltage of the DC bus (V DC ), internal inductance (L), and also as gain factors of the internal control blocks (G fh ) of the GSC. The current gain, the voltage gain, the gain of current control compensator, the feedforward gain, and the modulator gain are internal control parameters affecting G fh , which should be considered in the calculation of Z fh [28].
Finally, it should be noted that the reactance of X fh was employed for the output filter of the ESSs, as shown in Figure 3. The proposed model in Figure 4 could be employed for the coordination of the distributed ESSs for harmonic compensation. In the next section, the proposed model for harmonic compensation with ESSs is introduced.

Coordination Model of ESSs for Harmonic Filtering
In this section, the proposed model for the coordination of ESSs for harmonic mitigation is introduced. The model was developed based on the concept of the harmonic power flow [5,31] and active power filter sizing and placement [32]. Since the main power transactions of ESSs are in the fundamental frequency, the side operation for harmonic compensation could be implemented as a sequential operation model. Hence, the harmonic spectrum approach could be employed to calculate the harmonic filtering sharing of each ESS after running the day-ahead energy market [33]. In the following, the optimization model for incorporating ESSs in the operation of the MG at the fundamental frequency is first presented. After that, the proposed model for the coordination of ESSs working as AHFs (named the ESS-AHF model) and related working constraints are introduced.

Operation in Fundamental Frequency
In the fundamental frequency, ESSs can be used in the daily operation of the MG for energy arbitrage. This problem can be formulated as an operation problem including ESSs [34,35]. The goal of the MG operator is to economically supply the load demand using energy purchased from the main grid (cost grid ), the energy of the available DGs (cost DG ), and eventually, the energy arbitrage of ESSs (cost ESS ). Hence, the operation cost (OC) can be mathematically formulated as the minimization problem shown in Equation (2).
The objective function in Equation (2) consists of three main costs. The energy purchased from the main grid can be calculated with respect to the price of energy and power imported from the main grid to the MG as shown in Equation (3). Furthermore, DGs could also inject power to the MG whose related costs can be calculated as per Equation (4) in each operation period. For ESSs, the costs should be calculated with respect to the injected power in discharging mode and consumed power in charging mode. This can be mathematically formulated as in Equation (5). In Equations (3)-(5), the terms of π stand for prices of energy in the main grid for the injected power of DGs and prices for the charging and discharging of ESSs, respectively.
There are some technical and systematic constraints for the minimization of the OC, which should be taken into account. For ESSs, the working constraints are shown in Equations (6)-(10). Three possible working modes of ESSs can be considered: charging, discharging, and idleness modes.
These modes can be modeled as Equation (6) by using the binary variables of y cha (stands for charging mode) and y dch (stands for discharging mode). If these two variables are simultaneously adjusted to zero, then the ESS is idle in that period. The state of charge (SOC) of ESS (c et ) was calculated in Equation (7) in each operation period. It was determined based on the SOC at the beginning of the period in addition to power transactions in charging (p cha ) or discharging (p dch ) modes. In Equation (7), charging and discharging efficiencies are also considered as η cha and η dch , respectively. Constraint (8) is proposed for the calculation of the injection current (i) of the ESS to the connection bus in each mode using the conjugate of bus voltage, v it , and the active power of ESS. Since the operation modes are separated using binary variables in Equation (6), p cha and p dch do no occur simultaneously. Hence, in Equation (8), either charging or discharging of the ESS could take place. Consequently, the injection current could be positive or negative, with the injection current considered to be positive in charging mode. In Equation (9), the working limits of the SOC are considered using C min and C max , respectively, for the minimum and maximum limits. The maximum and minimum limits of active power are also considered in Equation (10).
Power flow constraints are the main system constraints in the optimization model. Current injections of loads, DGs, and nonlinear loads could be calculated in Equation (11) using the apparent power (S) and the voltage of the connecting bus (v i ) at each operation period. The injection shown in Equation (11) is valid for the set of load buses (Ω PQ ), the set of DG buses (Ω DG ), and finally the set of nonlinear buses (Ω NL ). Without losing generality, it is assumed that the DGs work in constant power mode. Anyway, other operation modes could simply be applied in the model [31]. Furthermore, the efficient linear power flow constraints of [31] were employed for modeling the network. In this power flow model, the MG was modeled with some incidence matrices using graph theory. Power flow constraints are modeled in Equations (12) and (13) using A and B incidence matrices employed for modeling the MG. In these equations, the matrices of the line currents, bus injections, and voltage drops of buses are shown with [U], [I], and [∆V], respectively. It has been shown in [31] that this power flow model has no convergence problem in MGs and can easily handle meshed topologies and unbalanced networks.
Nevertheless, modeling the energy arbitrage mechanism of ESSs was not the main concern of this paper. The output of the fundamental frequency optimization model is the SOC of each ESS and the working points in the operation horizon were used for calculations in harmonic frequencies.

Coordinated Operation of ESSs as AHF
In this subsection, the proposed model for employing ESSs as AHFs (named the ESS-AHF model) is introduced. The goal of the MG operator is to minimize the cost of procurement of PQ improvement actions. Since the harmonics were only investigated in this paper, this goal could be reached using AHF resources, APLCs, and ESSs. Hence, the objective function can be mathematically shown as Equation (14).
The cost of PQ improvement actions includes two main costs: costs related to APLCs denoted by the cost APLC and costs related to the ESSs working as AHFs, shown with cost ESS . The calculation of these two cost terms is also represented in Equation (14) using the amount of the provided compensation distortion power (d) and the offered compensation price (π) for each AHF (f ) at each period. In harmonic frequencies, due to the increase in power losses of the AHF, some amount of distortion power is imposed on the AHF. The payments are required for covering such costs for the owners of AHFs. Hence, offered prices of harmonic compensation, currents should be capable to refund extra costs imposed on the filter. On the other hand, since the required compensation actions are supplied in a competitive mechanism, the offered prices should ensure sufficient participation of the ESS owner in PQ improvement activities. Moreover, the offered prices of ESSs are generally lower than the prices of APLCs, since the harmonic compensation action is taken as a side function from ESSs [5]. Hence, some operation costs related to the investment cost can be decreased for ESSs. Anyway, the cost of all AHFs can be calculated using the offered price of each AHF and the related provided compensation distortion power as shown in Equation (14).
The distortion power generally includes the current distortion power [Marini, 2019 #445]. Hence, the provided compensation distortion power could be shown as Equation (15) using the provided harmonic compensation currents in each harmonic order, i fth . In this equation, the distortion power is assumed to be limited by the maximum available power of the AHF, D max . For APLCs that do not contribute at the fundamental frequency, this limitation is simply considered as the rating volt-ampere of the filter. However, for AHFs based on ESSs, considerations should be made for power transactions in the fundamental frequency. This limitation is set as the maximum value of available power of the ESS (P max -P et,1 ) for idle periods and allowable 10% variation in output power (P et,1 ) for charging or discharging periods. This constraint is mathematically shown as the maximum value of these two values in Equation (16). Furthermore, with respect to the general model of AHF as shown in Figure 4, working constraints could be formulated as Equations (17) and (18). In harmonic studies, harmonic voltages and currents are superimposed to the fundamental voltage and current. Hence, voltage and currents should be replaced with harmonic root mean square equivalents, as shown in Equations (17) and (18). The maximum current of AHF is limited to the maximum value, I max in Equation (17). In Equation (18), the maximum allowable limit of AHF is considered as V max .
The SOC of the ESS should also be taken into account in harmonic conditions. Although harmonic power transactions have low values, their effects on the SOC of the ESS should be considered. Harmonic active power (p fth ) can be calculated using the real value of the product of harmonic voltage and current as shown in Equation (19). The SOC of the ESS should be modified with respect to the provided active harmonic power. Since the PE converter of the ESS is generally a four-quadrant converter, harmonic active powers could be positive and negative. This could lead the ESS to be charged or discharged to provide the harmonic compensation current. Variations in the SOC in harmonic condition (c fth ) with respect to the provided active harming powers is mathematically shown in Equation (20). In Equation (20), the SOC in the harmonic condition is supposed to be equal to the SOC for fundamental frequency (C ft ) plus harmonic active powers, p fth . Harmonic reactive power is usually achieved by changing the phase angle and modulation method of the converter. It can be considered as extra losses of active power and is usually considered as a percentage of active harmonic power [36]. This modification could be applied in Equation (20), which could be a subject for future research. Finally, the minimum and maximum limitations of the SOC should also be met in harmonic activities, which is mathematically shown in Equation (21).
Harmonic standards are usually reported as total harmonic distortion (THD) and individual harmonic distortion (IHD). Typical levels of harmonic standards for grid voltage are 5% and 3% for THD and IHD, respectively [37]. The important goal of the MG operator is the economic satisfaction of the desired harmonic levels across the entirety of the MG at each period. These constraints are shown in Equations (22)-(23), respectively, for THD and IHD. In each operation period, the THD could be met for each bus and the IHD should be met for each bus in each harmonic order.
The harmonic power flow constraints should also be considered in the model. Since the provided harmonic current compensation is injected to the MG in different nodes, these constraints are essential for calculations of the flow of harmonic powers. The constraints of Equations (24) and (25) are the harmonic counterparts of fundamental power flow previously shown in Equations (12) and (13). Furthermore, the harmonic spectrum approach is used for the harmonic modeling of nonlinear loads [33]. It is required to determine the magnitude and phase angle of harmonic injections of nonlinear loads in this approach. In Equation (26), the magnitude of harmonic injection (i iht ) is determined based on the harmonic spectrum of nonlinear load (I spec ) and current injection of the load in the fundamental frequency (I it,1 ) for all buses of the network with nonlinear load (Ω NL ). The phase angle of harmonic injection (θ ith ) was also determined in Equation (27) using the phase spectrum of the nonlinear load (θ spec ) and the phase angle of its current injection in the fundamental frequency (θ it,1 ).
The developed model in Equations (14)- (27) is the complete model proposed for the coordinated operation of ESSs working as AHFs. It should be noted that the operation model of Equations (2)- (13) and the ESS-AHF model of Equations (14)- (27) are related and the ESS-AHF model will be implemented after the operation model. Hence, the ESS-AHF model could be accounted for as a sequential model, as shown in Figure 5. The output nonlinear load current of the operation model is the setpoints for the calculation of harmonic injections. These dependencies are mathematically shown in Equations (26) and (27). Harmonic injections are calculated in Equations (26) and (27) based on the harmonic spectrum of the nonlinear load. Hence, these dependencies could affect the ESS-AHF model since the harmonics will be related to harmonic injections of nonlinear loads. Furthermore, since ESSs are employed as AHFs in addition to APLCs, the working conditions of ESSs could be affected by the output of the operation model. The maximum available capacity of the distortion power for ESSs is set as the maximum value of available power of the ESS and allowable 10% variation in output power for charging or discharging periods, as shown in Equation (16). Hence, these setpoints could change the ESS activities as AHFs. The output of the ESS-AHF model is the calculated harmonic reference currents of the ESSs, which should be injected into the MG. These reference currents are then transmitted to the ESSs by the available infrastructures of the MG, and the AHF works in the current control mode. Various technical constraints, harmonic conditions of the network, and the status of AHFs should be sent to the MG control center for consideration in the optimization problem. All of these require a proper communication infrastructure that could be available in future MGs. The availability of communication links in AHF buses, its bandwidth, and time delay are important challenges of the communication system that could affect the applicability of the proposed optimization method, which should be taken into account in future studies. In the next section, the results of the implementation of the proposed model will be represented.

Simulations
A modified version of the IEEE 33-bus test network was used to demonstrate the performance of the ESS-AHF model. The test grid consists of 32 branches and five tie-switches, as shown in Figure 6. Line and bus data were taken from [38] and are provided in Appendix A. In order to increase the effects of harmonics, line impedances were increased by a factor of three. Data for the daily load demand were taken from [39].
Four nonlinear loads were connected to the test network. The locations and active and reactive powers of nonlinear loads are shown in Table 1. A general non-linear load model as a six-pulse diode bridge rectifier was adopted for all nonlinear loads [33]. The represented harmonic loads are the main sources of harmonic distortions in the MG. The harmonic spectrum of nonlinear loads is reported in Table 2, where the data were taken from [33]. In the harmonic spectrum approach, the magnitude and phase of the harmonic injections of nonlinear loads are used to calculate the flow of harmonic power in the network individually at each harmonic frequency.
Five DERs were added to the test network in selected buses, as shown in Figure 6. A parking lot of electric vehicles was considered at bus #12, a wind turbine was installed at bus #21, a photovoltaic (PV) generator was installed at bus #24, a battery energy storage was connected at bus #28, and finally, a PV array was connected at bus #31 equipped with a battery energy storage. The parking lot at bus #12 included a battery ESS. Hence, three ESSs have been placed in the test network which could be used as possible AHFs. Furthermore, two APLCs were connected to the network in buses #2 and #9, as shown in Figure 6. Hence, the total number of AHFs was equal to five filters. Data of the DER and APLCs are presented in Table 3. The maximum and minimum of the SOC and maximum power of each ESS are also reported in this table. In the last column, the offered price of distortion power is shown based on $/kVARh. The offered prices for ESSs were lower than the APLCs since harmonic filtering is considered as a side function from ESSs, which leads to a reduced distortion power cost.  The lower and upper bounds of voltage were taken as 0.95 pu and 1.05 pu, respectively. The operation horizon for unit commitment problem is considered the next day, which is divided into hourly periods. The proposed optimization models (in the form of a mixed-integer nonlinear programming model) were implemented in GAMS software (GAMS Development Corporation, Washington, D.C., USA) [40] in a personal computer equipped with a 2.93 GHz Intel processor and 8 GB memory. Without loss of generality, it was assumed that the MG works in grid-connected mode. Hence, the voltage at the PCC is assumed to be sinusoidal. The simulation results are reported in the following three subsections.

Case I: The Base Case
In the base case, the unit commitment model of 4.1 was implemented on the test network. The various cost terms after running the model are summarized in Table 4. The total OC was equal to 2410 k$, 87% of which was due to the energy supplied by the main grid. The other 13% of the OC was due to the energy produced or exchanged by DGs and ESSs. The cost of the energy supplied by DGs was equal to 248 k$ and represented the revenue of DG owners for participating in the energy market. In the revenue of DG #2, the PV at bus #24 was higher compared to other DERs for its proper location. Moreover, smart charging and discharging of ESSs implied a revenue of 58 k$ from energy arbitrage. This revenue is, of course, a cost for the MG operator. In particular, ESS #2, the battery storage at bus #28, earnt about 50% of the total ESSs revenues. The charging profile as well as the SOC of ESSs in the operation horizon is shown in Figure 7. All ESSs were charged from 3 am to 6 am, which are the off-peak periods of the operation horizon. Similarly, ESSs are discharged in peak periods from 5 pm to 7 pm. In Figure 7b, the SOCs of ESSs are represented. Allowable limits in Table 4 for the SOC and maximum power are met by the optimization model. Variation of the SOC for ESS #2 with more energy capacity at off-peak and peak periods can be seen in this figure.
However, the harmonic conditions should also be considered for PQ monitoring of the MG. In Figure 8, the voltage THD of the MG was analyzed. The voltage THD of all buses at 5 am (as the lowest) and 6 pm (as the highest load demand), are shown in Figure 8a. Obviously, either for off-peak periods or peak periods, the THD was above the standard level for nearly all system buses. All other periods may have a THD trend between these two selected periods. Bus #17 had the maximum THD compared to other buses. In Figure 8b, the voltage THD of this bus was represented in the operation horizon. The standard THD level is never met and hence, compensation actions are necessary.  In order to mitigate the harmonics effects from the grid, two simulation cases are provided to show the performance of the ESS-AHF model.

Case II: Harmonic Compensation Using APLCs
In Case II, only APLCs were employed for harmonic mitigation. APLCs were placed in the network based on filtering requirements and economic conditions by optimization models. The voltage THD of the system buses and the maximum THD after compensation of the APLCs are shown in Figure 9 for this case. Although the harmonic pollution was improved when compared to Case I, it was still above the standard level of 5%. This situation is worthwhile in peak periods, as shown in Figure 9a, since harmonic injections of nonlinear loads are increased in these periods. The maximum THD was for bus #32, which was about 11%. In Figure 9b, the maximum THD is shown for this bus. For off-peak periods, APLCs could adjust the voltage THD at the standard level. However, most of the time, the voltage THD is high, which could be the result of insufficient harmonic filtering capacity in the MG. In Figure 10 the loading of each APLC is shown for Case II. In almost all periods (except off-peak ones), the APLCs were loaded at their maximum rating. However, they could not satisfy the harmonic standard in the MG. Hence, more harmonic filters are required to be employed for more filtering capacity in the MG. This could be achieved using more APLCs. However, it seems that this is not an economic solution. With respect to the high investment costs of the APLCs, it could significantly increase the cost of PQ compensation actions. In the next case, the performance of ESSs for improving the harmonic situation of the MG as an economical alternative to APLCs is presented.

Case III: Harmonic Compensation Using ESSs
In order to increase the harmonic filtering capability of the MG as well as reduce the costs of PQ improvement actions, ESSs are also employed in this case for harmonic compensation. To this end, the unloaded power capacity of ESSs is used for harmonic filtering activities. Even if the ESS does not contribute to energy arbitrage, it can provide harmonic filtering capacity for the network. The results of the implementation of the ESS-AHF model are represented in this subsection.
The voltage THD and the maximum THD at each hour are represented in Figure 11. As it could be seen, the maximum THD in all periods of the operation horizon was adjusted at the standard level. For either the off-peak periods or peak periods, the voltage THD was in the acceptable range. Hence, using ESSs as AHFs has provided more harmonic filtering capacity for the MG. Furthermore, since the ESSs are dispersed in the MG, the provided compensation currents could be considered as distributed harmonic compensation actions. In this scheme, harmonic compensation currents are injected in the grid, and hence harmonic propagation and the resulting losses are strongly reduced. In Figure 13b, the maximum THD of voltage buses is shown, which was adjusted to the 5% standard level. This shows the effects of ESSs for the adjustment of harmonic standards for MGs. The loadings of various AHFs are shown in Figure 12 for Case III. AHFs based on ESSs are almost loaded at their maximum capacity. This is the result of their lower prices for compensation distortion power as well as their location. Since the APLCs have expensive distortion power, their compensation loadings are performed after ESSs and generally have lower values compared to ESSs. APLC #2 had more compensation loading despite its higher cost. This shows the effects of filter location on its provided harmonic compensation to the MG. The harmonic active power and the SOC of ESSs after harmonic compensation activities are reported in Figure 13. The provided harmonic active power was negative for all ESSs, which shows that it worked in discharging mode for harmonic compensation. However, the harmonic injection powers had low values compared to the fundamental active power previously shown in Figure 7. This means that participating in harmonic compensation activities does not have a significant effect on the SOC of ESSs, which can be seen Figure 13b. The trend of the SOC was close to the trend of that in Figure 7 and all allowable limits were met in the harmonic conditions. Hence, the harmonic utilization of ESSs could bring various advantages for the MG, and meanwhile, does not affect the ESS parameters for energy arbitrage. The ESSs earn from both energy arbitrage and harmonic compensation. The required reference currents of the ESSs at each period are determined in the ESS-AHF model. Fundamental harmonic and reference currents for ESS #2 at the peak period (i.e., 6 pm) are shown in Figure 14. The bold curve is the required reference current calculated by the ESS-AHF model, which should be injected by the ESS. The reference current should then be sent to the ESS by the available infrastructures of the MG. The ESS is responsible to inject the required current by the MG operator with a suitable control and switching method. For economic comparisons, the total injections and compensation costs for Case III are presented in Table 5. The total required distortion power was 8.95 pu and 6.29 for Case II and Case III, respectively. Although the compensation distortion power in Case III was about 70% of that in Case II, it could successfully meet the desired level of harmonics in the MG, as shown in Table 5. Furthermore, the loading of compensation currents was more achieved by the ESSs in Case III as they were more economic AHFs compared to the APLCs. This reduces the total compensation cost from 7868 k$ in Case II to 3424 k$ in Case II. It showed a 44% decrease in the compensation cost. Hence, the MG operator could keep the desired harmonic level with only 44% of the cost of the APLCs. This could save remarkable costs for the MG operator, either in the operation of the APLCs or as requirements for the installation of new APLCs. The simulations showed the performance of the ESS-AHF model in the coordination of distributed ESSs for harmonic compensation in the MG. More filtering capacity and reduced compensation costs are the main advantages of ESSs in harmonic compensation. Since harmonic compensation is utilized as a side function, the provided distortion power has a lower price than APLCs. This implies that compensation loading is first demanded by the ESSs. Furthermore, ESSs are dispersed in the MG, and hence harmonic compensation will be distributed in the MG. The provided compensation scheme will be accounted for by distributed harmonic compensation with more flexibility and applicability for PQ improvements. The ESS owners could also benefit from both energy arbitrage and harmonic compensation. This could decrease the investment costs of the ESSs and make them more justified for use in the MG. The simulations showed that ESS compensation activities do not have much of an effect on the SOC and ESS revenues from energy arbitrage. Finally, it should be noted that, unlike energy and ancillary markets, PQ improvement actions are not security obligations of the network. Therefore, any time that the working of ESSs in PQ improvement actions affects their application in ancillary services and the security of the network, it could be temporarily ignored.

Conclusions
In this paper, the applicability of ESSs for harmonic filtering activities was investigated. The common procedure for this was to modify the control methods of the PE converters of ESSs to enable them for the injection of compensation currents. However, in this paper, the coordinated operation of ESSs for a distributed harmonic compensation scheme was addressed. The proposed coordination model employed conventional control algorithms for the injection of a reference current that was determined in the proposed ESS-AHF optimization model. Hence, the provided compensation scheme could be considered as a distribution compensation scheme. The simulation cases showed the performance of the proposed model in enhancing the filtering capability of the MG, the reduction in compensation cost, and the greater flexibility of the distributed compensation schemes. Hence, requirements for new filters could be deferred and both the MG operator and the ESS owner could benefit from reduced compensation costs.
ESSs can be used to perform several ancillary services where most of them need the energy stored to be accomplished. What is demonstrated in this paper is that the ancillary service for harmonic power filtering can be performed without significantly affecting the battery SOC. For this reason, this service can be achieved for free, without compromising the others, and therefore the revenues from energy arbitrage and other ancillary services can be maintained. This could make ESSs more economical for use in future power systems. Investigations should be carried out for various technical and systematic challenges of the proposed ESS-AHF model in future research. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflicts of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Nomenclature
Note that single parameters, single variables, and matrices are shown with capital letters, small letters, and capital letters in brackets, respectively, in the paper. When each notation comes with a specific index and superscript, it stands for the corresponding variable or parameter of related indices. Repetitive definitions can be ignored using this method. For example, the letter 'v' always stands for voltage. When it comes in brackets, it stands for the matrix of all voltage variables; when it comes individually with index i, it stands for a single variable of voltage for bus i, and V it,1 stands for the voltage parameter of bus i at period t for the fundamental frequency, h = 1. The notations are as follows: Abbreviations