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.
4.1. 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 (
costgrid), the energy of the available DGs (
costDG), and eventually, the energy arbitrage of ESSs (
costESS). 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
ycha (stands for charging mode) and
ydch (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 (
cet) 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 (
pcha) or discharging (
pdch) 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,
vit, and the active power of ESS. Since the operation modes are separated using binary variables in Equation (6),
pcha and
pdch 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
Cmin and
Cmax, 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 (
vi) 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.
4.2. 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
costAPLC and costs related to the ESSs working as AHFs, shown with
costESS. 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,
ifth. In this equation, the distortion power is assumed to be limited by the maximum available power of the AHF,
Dmax. 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 (
Pmax–
Pet,1) for idle periods and allowable 10% variation in output power (
Pet,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,
Imax in Equation (17). In Equation (18), the maximum allowable limit of AHF is considered as
Vmax.
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 (
pfth) 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 (
cfth) 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 (
Cft) plus harmonic active powers,
pfth. 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 (
iiht) is determined based on the harmonic spectrum of nonlinear load (
Ispec) and current injection of the load in the fundamental frequency (
Iit,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.