**1. Introduction**

Recently, with the increased popularity and application of power electronic devices in modern industries, some large electrical loads such as from an AC traction system (single-phase), electric furnaces and modern technologies based on high-tech microprocessors in renewable energy applications cause increasing power quality (PQ) problems [1,2]. Power quality problems appear in electric utility as increased harmonic distortions, phase imbalances and low power factors. Electric devices based on power electronics present themselves as non-linear loads, which are characterized as sources of harmonics in the power system. The harmonic currents cause an increase in the RMS value of the current and let the neutral current make circulations in the electric distribution system. The capacity of the distribution system and power losses are affected by the presence of harmonics [2,3].

The effects of nonlinear unbalanced loads on the performance of the power system is an important issue for the power system operators and attracts the interest of many researchers. Many research papers have proposed different techniques and optimization algorithms to improve power quality indicators in the presence of nonlinear and/or unbalanced loads [4]. The static synchronous compensator (STATCOM) is one of the FACTS devices that introduces an efficient and flexible solution for power quality improvement and disturbance mitigation [5–7]. The practical circuit of the STATCOM usually includes two-level voltage source converters (VSC), and a line-frequency transformer is adopted for the enhancement of its voltage and current ratings; however, it results in the heavy, unreliable and expensive design of the compensator [5,6].

The multilevel converters-based STATCOM has received considerable attention from researchers due to its higher capacity that makes it reliable for application in medium- or high-voltage high grids in the absence of line-frequency transformers. Since 1975, the concept of the multilevel converter has been introduced [8]. The multilevel term began with the third-level topology [8,9]. Subsequently, many multilevel converter topologies were developed [10–19]. However, the basic concept of the multilevel converter is the use of a series of semiconductor power switches with several low DC voltage sources to conduct energy conversion and obtain higher energy. The topology of the multilevel converters has succeeded in addressing some of the main problems in the traditional converters [10–19]. The most important of these problems are the problems of power quality, especially the voltage and current harmonics [20–23]. As the number of converter levels increases, the output voltage signal will be as good and as close as possible to the desired waveform [9–19]. But of course, advanced control techniques are needed to improve the performance of multilevel converters and optimize the workflow.

Three major multilevel converters have been applied in high- and medium-power applications, such as the neutral-point clamped multilevel converter (NPCMC), the flying capacitor multilevel converter (FCMC) and the cascaded H-bridge multilevel converter (CHBMC) [8,9]. The major disadvantage of the NPCMC is the inherent voltage imbalance problems, which require an extra auxiliary compensation circuit [8,9]. As for FCMC, besides the need for highly expensive flying capacitors, mainly at low carrier frequencies, this type of multilevel converter is not suitable for those applications with 90◦ leading or lagging currents. Moreover, for both topologies, as the voltage level has to be increased in high power applications, an increase in the size of the capacitors or in the climbing diodes is required, which makes the control circuit and converter structure more complex. Compared with the previously explained two former topologies, CHBMC requires fewer components in the circuit design, which makes it convenient for physical layout and packaging. The major disadvantage of CHBMC is that it cannot perform properly under imbalanced conditions without feeding from isolated DC sources supplied from multi-winding transforms, which gives rise to the same problems associated with the use of the line-frequency transformer [8,9].

To overcome the previously mentioned problems, a new modular multilevel converter-based STATCOM (MMC STATCOM) is recommended for medium and high voltage applications. The MMC STATCOM is characterized by its availability, compact and modular construction, generation of least harmonics, etc. Moreover, the MMC STATCOM is able to exchange both active and reactive power via the common DC-bus. Accordingly, the active power, which is redistributed in the internal circuits, is utilized for the purpose of negative sequence balance. Consequently, the MMC STATCOM can operate under three-phase imbalance without any intermission.

In this paper, an optimization problem is introduced for the optimal design of the multilevel converter and for determining the optimal parameters of the DC voltage controller. The minimization of total harmonic distortion (THD) in the power system under consideration is introduced as the objective function while keeping the ripple voltage within the allowable ranges and obtaining minimum circulating current. Novel optimization algorithms, namely Atom Search Optimization (ASO) and Harris Hawk optimization (HHO) that have not been reported in the literature for tackling the aforementioned problem of optimization of the values of the capacitors of the MMC STATCOM have been applied. Moreover, to ensure the effectiveness of the proposed optimization techniques, the

results obtained from the application of these methods have been tested under two different case studies, nonlinear load and three-phase imbalance.

The paper is organized as follows: in Section 2, a brief review of MMC topologies is introduced in detail; Section 3 introduces the Model of the Modular Multilevel Converter (MMC); the Modular Multilevel Converter-based STATCOM (MMC STATCOM) is introduced in Section 4; Section 5 presents the problem formulation; Section 6 presents the Harris hawk's optimization (HHO) and atom algorithm in detail; the result simulation presents in Section 7; and the conclusion is stated in Section 8.

#### **2. A Brief Review of MMC Topologies**

There are many topologies of the MMC. A comprehensive comparison among these topologies is presented in Table 1 [24,25].

**Table 1.** Comprehensive comparison between Modular Multilevel Converter (MMC) topologies.

For more explanation to state the motivation of this paper, the level of harmonics in the output voltages and currents of the MMC is very essential and has to be within the allowable limits. The harmonic content could be decreased without the need for any filters with the help of modern control methods [26,27]. In order to obtain the output voltage with an acceptable level of THD in the AC side cascaded configuration, the switching frequency of the two-level converter (main converter) is adjusted to 400 Hz, while the frequency of carriers is taken as 2 kHz. Due to the possible mismatch of synchronization between the main converter and the active filter (H-bridge cascaded model), the fifth and seventh harmonic components may arise, which must be canceled using a number of filters.

As example, reference [26] presented a comparison between the output voltage and the harmonic spectra generated from the single-phase AC cascaded configuration with two H-bridges on the AC side, which are presented in Figures 1 and 2, and the output voltage and the harmonic spectra generated from the single-phase AC cascaded configuration with 4H-bridges on the AC side, which are presented in Figures 3 and 4. It can be seen from Figure 2 that the third harmonic component has a large value and also that the value of THD is extremely high. That is because the method of third harmonic subtraction allows the multilevel converters to work independently of the load power factor and without the problems of capacitor voltage balancing [26]. The data presented in Figure 4 prove the superiority of the last configuration of 4H-bridges, as the profile of the output voltage harmonic is slightly better, but with an increase in the number of the power electronic devices and the cost of the system. Moreover, reference [26] concludes that the increasing of the series sub-modules (SMs) in the case of using two H-bridges or 4H-bridges will lead to a reduction in the harmonics and over-shoot in the output voltages of the MMC.

**Figure 1.** Voltage for a single AC side cascaded configuration including two H-bridges [26].

**Figure 2.** A representation of the harmonic component in the output voltage of the AC side cascaded configuration using two H-bridges [26].

**Figure 3.** Voltage for a single AC side cascaded configuration including 4H-bridges [26].

**Figure 4.** A representation of the harmonic component in the output voltage of the AC side cascaded configuration using 4H-bridges [26].

The introduced discussion leads the authors to present an optimal design of the control system and the optimal capacitor value of the MMC for the STATCOM application, using the half-bridge SMs to decrease the number of devices and the cost of implementation.

#### **3. Model of the Modular Multilevel Converter (MMC)**

The modular multilevel converter (MMC) has been introduced firstly by Marquardt in 2001, and is presented as the developed configuration including a number of subsystems for achieving the desired voltage level [15,18,24]. This type of converter is designed for use in medium and high transmission voltages. The overall configuration of the modular multilevel converter is illustrated in Figure 5, where

each leg of this converter consists of a lower and an upper arm connected to each other by a DC voltage link.

**Figure 5.** The circuit configuration of the MMC.

Each arm includes a number of SMs, and in turn, every sub-system includes two switches and two reverse diodes connected to each other via a DC capacitor. The switches S1 and S2 are operated in reverse. Each SM could be configured as a half-bridge, full-bridge, H-bridge or clamp-double circuit [18]. Due to its low cost and high efficiency, the MMC based on the half-bridge SMs shown in Figure 6 has been widely used in industrial applications and high voltage DC projects. It is possible to perform a selectable individual control for each SM in the converter leg. Principally, the converter legs provided a controllable Voltage Source Converter (VSC) for each phase. The desired waveform voltage can easily be achieved at the terminal by adjusting the terminal voltage of the SM's leg.

The equivalent circuit presentation of the typical MMC using half-bridge SMs is shown in Figure 7, where *uvj* is the phase converter output voltage and *ivj* denotes the line current; *ipj* and *inj* are the upper and lower leg currents, respectively, that are presented as [8,10,27]:

$$\begin{aligned} i\_{pj} &= \frac{i\_{vj}}{2} + i\_{diffj} \\ i\_{nj} &= \frac{i\_{vj}}{2} - i\_{diffj} \end{aligned} \tag{1}$$

where *idi*ff*j* denotes the phase *j* inter-circulating current that circulates between the upper and lower legs and is expressed as below:

$$i\_{diff\dot{j}} = \frac{i\_{p\dot{j}} - i\_{n\dot{j}}}{2} \tag{2}$$

**Figure 6.** A schematic of the configuration of the MMC based on half-bridge and full-bridge sub-modules (SMs).

**Figure 7.** The one phase equivalent circuit.

As discussed in [7], the MMC operates according to the following equations:

$$u\_{vj} = c\_j - \frac{R\_{arm}}{2} \cdot i\_{vj} - \frac{L\_{arm}}{2} \cdot \frac{di\_{vj}}{dt} \text{ ( $j = a, b, c$ )}\tag{3}$$

$$R\_{arm} \cdot i\_{diff} + L\_{arm} \cdot \frac{di\_{diff}}{dt} = \frac{V\_{dc}}{2} - \frac{u\_{n\bar{j}} + u\_{p\bar{j}}}{2} \ (j = a, b, c) \tag{4}$$

where *ej* denotes the inner *emf* produced in phase *a* and is given as follows:

$$\mathbf{e}\_{j} = \frac{\mathbf{u}\_{nj} - \mathbf{u}\_{pj}}{2} \tag{5}$$

The internal dynamic performance of theMMC is explained in Equation (4) and can be characterized as follows:

$$u\_{diff} = R\_{arm} \cdot i\_{diff} + L\_{arm} \cdot \frac{di\_{diff}}{dt} = \frac{V\_{dc}}{2} - \frac{u\_{nj} + u\_{pj}}{2} \text{ ( $j = a, b, c$ )}\tag{6}$$

where *udi*ff*j* denotes the phase j inner unbalance voltage.

The MMC is characterized by its small size and simplicity in installation, and as the number of SMs increases, there is a decrease in harmonics and a matching of the required power quality specifications. However, the multilevel converter su ffers from converter inter-circulating current, which increases the switches' conduction loss and increases the thermal stress on the capacitors and switches' IGBTs. Therefore, maintaining the balance of the capacitor voltages is one of the main rules of the design of the MMC.

#### **4. Modular Multilevel Converter-Based STATCOM (MMC STATCOM)**

In order to simplify the analysis of the proposed network, the upper and lower arms of the MMC STATCOM are proposed to be equivalent to one power unit, which includes two IGBTs shunted by opposite-connected diodes and a capacitor forming a bidirectional chopper. Figure 8 shows a simplified one-phase equivalent circuit of the proposed MMC STATCOM. Point P in the circuit refers to the positive DC-link, while point N refers to the negative DC-link. In addition, the following assumptions are provided for more analysis,


**Figure 8.** The one phase equivalent circuit of the MMC static synchronous compensator (STATCOM).

The terminal voltage ua and output current ia are assumed to be pure sinusoidal:

$$
\mu\_4 = \mathcal{U}\sin(\omega t)\tag{7}
$$

$$\dot{a}\_a = I \sin(\omega t + q \dot{\imath}\_i) \tag{8}$$

where *U* denotes the phase voltage peak value, *I* denotes the phase current peak value, ω is the angular frequency and ϕ*i* is the phase shift between voltage and current.

For a symmetrical three-phase system, the phase current *ia* is divided equally between the two legs. Thus, the resulting share from the upper and lower legs can be identified as follows [27]:

$$\dot{m}\_{\rm Pl} = \frac{1}{2}I[K\_I + \sin(\omega t + \varphi\_i)]\tag{9}$$

$$\dot{q}\_{\rm Na} = \frac{1}{2}I[\mathbf{K}\_I - \sin(\omega t + q\_i)]\tag{10}$$

where *KI* is the current coefficient that is expressed as

$$K\_I = \frac{2}{3} \frac{I\_d}{I} \tag{11}$$

The differential equations of the upper and lower arms of Figure 8 can be introduced as follows:

$$C\frac{dV\_c}{dt} = s\_1 i\_{\rm Pa} \tag{12}$$

$$\mathcal{C}\frac{dV\_{\mathcal{L}}}{dt} = \mathbf{s}\_2 i\_{\text{Na}}\tag{13}$$

$$L\frac{d\dot{p}\_{\rm fl}}{dt} + R\dot{p}\_{\rm fl} = \frac{V\_d}{2} - lI\sin(\omega t) - s\_1 V\_\varepsilon \tag{14}$$

$$L\frac{di\_{\rm Na}}{dt} + Ri\_{\rm Na} = \frac{V\_d}{2} + \mathcal{U}\sin(\omega t) - s\_2 V\_c \tag{15}$$

where *Vc* is the voltage of the capacitor, *s1* and *s2* are the upper and lower leg switching functions, *L* denotes the AC inductance of the leg inductor and *R* is the equivalent resistance of the leg inductor.

## *4.1. DC-Link Capacitor*

In this section, the calculations are only provided on the upper leg, as there is no difference between the upper and lower legs [28]. For more simplification, the calculations are performed taking into account the following assumptions:


In the steady-state condition, the upper leg switching function is defined as:

$$s\_1(t) = M \sin(at)\tag{16}$$

where M denotes the modulation index. From Equations (9), (12) and (16), the following formula is produced:

$$\mathcal{K}\frac{dV\_c}{dt} = \frac{1}{4}Ml\cos\varphi\_i + \frac{K\_I}{2}Ml\sin(\omega t) - \frac{1}{4}Ml\cos(2\omega t + \varphi\_i) \tag{17}$$

After integration, the following equation exists:

$$V\_{\mathcal{L}} = V\_{\text{DC}} - \frac{M\mathcal{K}\_I \cos(\omega t)}{2a\mathcal{C}} - \frac{M\sin(2\omega t + q\_i)}{8a\mathcal{C}}\tag{18}$$

where *VDC* denotes the average DC voltage. The ripple in the current of the upper leg is described by the following formula:

$$
\Delta V\_{\varepsilon} = -\frac{M I K\_I \cos(\omega t)}{2\omega \mathcal{C}} - \frac{M I \sin(2\omega t + \varphi\_{\bar{l}})}{8\nu \mathcal{C}} \tag{19}
$$

By similarity, the ripple in the current of the lower leg is calculated as follows:

$$
\Delta V\_{\varepsilon} = -\frac{M I K\_I \cos(\omega t)}{2\omega \mathcal{C}} - \frac{M I \sin(2\omega t + \varphi\_i)}{8\nu \mathcal{C}} \tag{20}
$$

When the phase angle between the voltage and current ϕ*i* equals ±90◦ during the static VAr generating condition and by ignoring the resistance of the inductor, Equation (19) will be converted to:

$$
\Delta V\_c = -\frac{MK\eta\cos(\omega t)}{2\omega \text{C}} + \frac{MI\sin(2\omega t)}{8\omega \text{C}}\tag{21}
$$

To obtain the maximum value of the peak-to-peak ripple voltage of the DC capacitor, Equation (21) is differentiated, which results in:

$$
\Delta V\_{\rm cpp} = \frac{Ml(K\_l + 1)^2}{4\omega \mathcal{C}} \tag{22}
$$

Accordingly, the capacitance of the DC capacitor can be defined as:

$$C \geq \frac{Ml(K\_I + 1)^2}{4\alpha\Delta V\_{cpp}} \tag{23}$$

From Equations (14) and (20) it is clearly shown that the voltage ripple contains the fundamental and the second harmonic, which dramatically increases the rating of the capacitor, resulting in an increase in the total cost of the system. An accurate determination of the capacitance value in each sub-model must be taken into account, as a higher value of the capacitance leads to an increase in the total cost of the MMC converter, and on the other hand, a lower value increases the ripple level in the output voltage.

#### *4.2. Design of SM Capacitance*

The value of the capacitor affects the power quality indices such as the level of the THD and the voltage balance of the system. Many methods in literature have been introduced to determine the values of the capacitors. These methods can be divided into two main methods; the first is based on the idea that the average sum capacitor voltage per arm should be kept constant. The second is based on keeping constant the individual average capacitor voltages for the sub-modules for each arm. Both of the two methods are based on determining the required stored energy of the STATCOM-based MMC. The conclusion based on the literature review is that the first method is considered to be a better choice than the second one because using the first method improves the reliability of the system, and it is more convenient for the sub-module capacitor, as well as the semiconductor device, to perform under low voltage [20]. A general description of the first method has been presented.

Based on the energy storage demands of the converter, the capacitance of the SM can be determined. The minimum value of the capacitance of the SM is calculated as follows [29]:

$$C = \frac{2 \ast N \ast E}{V\_{dc}^2} \tag{24}$$

where *E* denotes the lowest limit of energy storage for each arm. Due to the similarity between both arms, E for the upper arm can be expressed as reported in [25]:

$$E = \frac{\Delta E}{k\_{\text{max}}^2 - \max\left(\frac{n\_u^2 - v\_{u,v}k\_{\text{max}}^2}{1 - v\_{u,v}}\right)'}\tag{25}$$

where kmax refers to the capacitor voltage upper boundary. Typically, kmax = 1.1 is commonly used.

Finally, as ΔE and *E* are in direct dependence on the rated power of the converter, the demands for energy storage of the converter can be given as:

$$\mathcal{W} = \frac{\theta}{\mathcal{S}\_{\text{n}}} E\_{\text{nown}} \tag{26}$$

where *W* denotes the energy storage needed for each MVA rating. According to [24,30,31], the capacitance in each sub-module is selected so that the energy stored in all capacitors of the converter sub-modules is about 30–40 kJ/MVA. As the capacitance of the sub-module is not chosen yet, the ripple level of the sub-module capacitor voltage in the operating point is simply determined. There are two possible ways to obtain an anticipatory value of the voltage ripple. The first way is to assume a ripple level of ±10%, whereas the second way is based on assuming a value for the capacitor in each sub-module in the range, which achieved 30–40kJ/MVA stored energy, then performing simulation in the time domain to determine the level of the voltage ripple.

#### *4.3. Design of Arm Inductance*

Typically, the inductors used in the arms of the MMC are dry-type air-cored ones. The value of the inductance of the arm can be determined based on the way that the short circuits on the DC-side will be demonstrated. In addition, the inductance of the converter arm plays an important role in enhancing the characteristics of the inter-circulating current and limiting the fault current. Where bypass thyristors are used in the sub-modules or where they are separated using additional diode-based power modules, the AC sharing of the fault current on the DC-side does not need to be restricted to provide protection for the sub-modules' anti-parallel diodes. The arm inductance helps in the limitation of the high frequency harmonics in the inter-circulating current and to provide a smooth control for the circulating current. A typical value of the arm inductance occurs around 0.05 p.u. [24].

On the other hand, when bypass thyristors are not used in the sub-modules or they are not separated using additional diode-based power modules, protection for the anti-parallel diodes on the DC-side of the converter have to be designed to protect the sub-modules from the high sharing from the AC-side to the fault current in the DC-side. Consequently, as the slow-acting breakers interrupt under this fault current, the arm inductors must restrict the current magnitude to give protection for the anti-parallel diodes, which act as a conduction path for the fault currents. The value of the arm inductance is determined using DC-fault actual time simulation, taking into account the interruption time of the circuit breakers, the time delay until fault detection and the strategy of performing the fault in the control system. Typically, the value of the inductance occurred in a range of 0.10–0.15 p.u.

Moreover, a short circuit between the terminals of the DC-buses is performed in order to present the most critical fault condition. In order to make a limitation on the short circuit current, the value of the arm inductance should be determined according to [29]:

$$\mathbf{L}\_{\text{arm}} = \mathbf{V}\_{\text{dc}} / 2\,\text{\AA} \tag{27}$$

where α (kA/s) denotes the maximum rising rate of the current. Based on (27), when the maximum rising rate is = 0.1(kA/μs) [29], for grid connected converters, the p.u. value of the arm inductance values are restricted at 0.3 pu for Double-Star Chopper Cell (DSCC-MMC), as well as for Single-Delta Bridge Cell (SDBC-MMC) configurations.

Consequently, modern metaheuristic techniques such as Atom Search Optimization (ASO) and Harris hawk's optimization (HHO) have been proposed in this work for the optimization the of the values of capacitors and the coe fficients of the DC voltage regulator of the MMC STATCOM, in order to minimize the total harmonic distortion (THD) in the power system under consideration, keeping the ripple voltage in the allowable ranges and obtaining minimum current circulation within the converter.

## **5. Problem Formulation**

Usually, in electric power systems, shunt compensating devices are used for voltage and reactive power control. A MATLAB model for the system under study is shown in Figure 9. The parameters of each component in the system are summarized in Table 2. The system under study compromises the utilization of the MMC STATCOM, which is increasingly used in modern electric networks. The proposed MMC STATCOM is based on full-bridge MMC to form 22 modules per phase power converter. The STATCOM can produce or take reactive power from the grid. The transfer of both types of reactive power is obtained via phase reactance.

**Figure 9.** A model of the power system under study.



The proposed MMC generates a three-phase voltage, which has the same phase as that of the grid voltage at the point of coupling. When the amplitude of the voltage generated from the MMC is lower than the grid voltage amplitude, the STATCOM performs like inductance and absorbs reactive power from the grid utility. When the amplitude of the MMC voltage is higher than that of the bus of common coupling, the STATCOM works like a capacitor and reactive power will be injected into the grid. The

control circuit of the MMC is shown in Figure 10a. The controller of the high voltage side includes a phase-locked loop (PLL), transformation and measurement blocks, DC-bus voltage regulator and current regulator. Meanwhile, the controller of the low voltage side includes individual and phase balancing controllers of the DC-bus voltage, which are responsible for generating the control signals for producing the pulses applied to the gates of the transistors. Moreover, the optimized controller for phase balancing based on the HHO or ASO optimization techniques is shown in Figure 10b, and the optimized controller of DC voltage based on the HHO or ASO optimization techniques is shown in Figure 10c.

(**c**) 

**Figure 10.** The controller of the MMC converter: (**a**) Simulink configuration; (**b**) An optimized controller for phase balancing based on the Harris hawk's optimization (HHO) or Atom search optimization (ASO) optimization techniques; (**c**) An optimized controller of DC voltage based on the HHO or ASO optimization techniques.

As discussed in Section 2, the impact of the second component harmonic voltage of the module capacitor on the MMC impedance behavior is expressed. Optimal determination of the value of the module capacitor results in an improvement in the second component of harmonic in the voltage of the module capacitor. Moreover, by reducing the capacitance value, in increase in the second component of harmonic voltage will be observed. Furthermore, it is easy to find that with an increase in the magnitude of the module capacitor voltage second harmonics, an e ffective increase is observed in the magnitude of the impedance response at low frequencies, but there is no change noticed in positive-sequence. Starting from the introduced explanation, the optimal determination of the module capacitor can be addressed as the design of the impedance of the MMC. However, if a stable performance of the STATCOM is desired, the proposed design should give the minimum demands for energy storage and filtering ripples arising in the capacitance voltage, which shape the constraints of controlling module capacitance to produce an optimal impedance response. Therefore, this paper proposes that the value of the capacitor should be selected to be optimized to reduce the THD using the HHO and ASO algorithm. The objective function to achieve optimal capacitor values leading to the minimum quadratic error of the THD can be calculated as follows:

$$F\_1 = \int\_0^\infty t \cdot THD^2 dt\tag{28}$$

Moreover, to achieve the optimal values of both *Kp* and *Ki* for each PI controller, which lead to a minimum quadratic error, the objective function can be calculated as follows:

$$F\_2 = \int\_0^\infty t \cdot e^2 dt\tag{29}$$

where *e*. is the error signal between the desired and the typical signals and is fed to the PI controller. The objective function to minimize F can be determined as follows:

$$\text{Minimize } \left| F(\mathbf{C}\_{\prime}(\mathbf{K}\_{\mathbf{p}}, \mathbf{K}\_{\mathrm{i}})\_{\text{VDC}\_{\prime}}(\mathbf{K}\_{\mathbf{p}}, \mathbf{K}\_{\mathrm{i}})\_{\text{phase balance}} \right| \tag{30}$$

where *F* is the summation of the *F*1 and *F*2.

#### **6. Proposed Optimization Methods**

In the past decades, with the increasing development in the fields of society, economy and industry, numerous complicated and hardly-solved optimization problems have arisen in all fields. Recently, metaheuristic techniques of optimization have been widely implemented to solve complicated engineering problems, as they possess high capabilities of exploration and exploitation to reach global solutions in a short time when compared with traditional optimization methods [31–34].

Taking into account the No Free Lunch Theorem of Optimization [35], the best solution of di fferent optimization problems cannot be fulfilled based on one optimization technique. Based on this theorem, the field of creating and developing better optimization algorithms is still active and full of grea<sup>t</sup> achievements. As a result of these achievements, a recent swarm algorithm inspired by physics, called Atom Search Optimization (ASO), is proposed for tackling the optimization problem provided in this study. ASO presented superior performance over conventional and recent optimization techniques when validated on numerous mathematical problems. In addition, the proposed ASO algorithm provided a successful performance when utilized in the optimization problem of hydrogeologic parameter estimation.

In a related context, like the most of existing optimization methods, the searching process in the proposed ASO algorithm is performed in two phases; exploration and exploitation stages [30]. In the first phase, the optimization algorithm should make utilization and promotion of the randomly initialized operators for the deep exploration of the optimal solution within the search space. Thus, the exploratory behavior of a well-developed optimization algorithm should be randomly enriched to generate more random solutions to di fferent areas of the addressed problem formed in the early stages of the research process [36].

#### *6.1. Harris Hawk's Optimization (HHO)*

The HHO algorithm was developed in 2019, which in its operation mimics the behavior of the hunting of the Harris hawks, a group of intelligent birds that live in the USA [32,33]. The hunting manner of these intelligent birds depends on surprise. A group of hawks appear to their prey (rabbits) from di fferent sides at the same time to surprise it, at which point the best-fit hawk of the group (the leader) surrounds it. Mathematically, the technique of catching the prey of the hawks can be divided into three stages: (i) the exploration stage, (ii) the stage of transition from exploration to exploitation and finally (iii) the exploitation stage. The first phase of the attack is the most important as it determines the way that the hawks will act. Firstly, all hawks sit randomly in a high place to give them a good chance of observing the environment and wait for their chase. When the prey appears in the surrounding area, the hawks may attack it in several possible manners, and only the leader will decide the way of attack. All hawks might cooperate and surround the prey, or the leader will give the opportunity to one of them, and it will depend on the behavior of the prey to escape. For each of the two possible tactics, the hawks may take their places depending on the location of the neighboring hawk as expressed in (31) or be distributed randomly as described in (32):

$$X(t+1) = \left(X\_{\text{best}}(t) - X\_{\text{avg}}(t)\right) - \psi\left(LB + \pi\left(LB - LB\right)\right)\alpha < 0.5\tag{31}$$

$$X(t+1) = X\_{num}(t) - \beta \left| X\_{rand}(t) - 2\rho X(t) \right| \text{ a} \ge 0.5 \tag{32}$$

where *X(t)* is the position of the hawks at the present iteration *t*; *X(t*+*1)* is a vector of the new positions of the hawks in the next iteration; *Xbest(t)* is the position of the chase (rabbits); and α*,* β*,* ψ*,* ϕ and τ are randomly distributed numbers in the range from 0 to 1. *UB* and *LB* denote the upper and lower limits of the position variables. *Xrand(t)* denotes a hawk, which is randomly selected from the present population, and *Xavg(t)* denotes the mean position of the hawks. The average position of the hawks is determined as follows:

$$X\_{\text{avg}}(t) = \frac{1}{N} \sum\_{i=1}^{N} X\_i(t) \tag{33}$$

where *Xi*(*t*) denotes the position of each hawk in iteration *t* and *N* indicates the number of hawks. Depending on the energy of the prey during the escaping process, HHO can make a transition from exploration to exploitation. The energy of the prey (rabbit) is determined as follows:

$$E = 2 \times E\_0 \left( 1 - \frac{t}{T} \right) \tag{34}$$

where, *E0* denotes the initial energy in each iteration, which is randomly taken from [−1, 1], and *T* presents the maximum iterations.

As the hawks have detected the prey as explained in the previous stages, they have to make a sudden attack, which depends on the probability of escape. In the mathematical modelling, the probability of escape *r* will be <*0.5* for a successful escape and *r* ≥ *0.5* if the prey unsuccessfully gets away. Depending on the way that the prey will take during escape, the hawks will carry out a soft or hard siege. Regardless of the energy of escape, a hard or soft siege will occur. When |*E*| ≥ *0.5* and *r* ≥ *0.5*, the prey has enough energy to make random jumps but fails to escape. During these miserable

attempts, the hawks surround the prey until it is exhausted by its force and then suddenly attack it. Thus, the soft siege is described as:

$$X(t+1) = \Delta X(t) - E\left| (I \times X\_{\text{best}}(t)) - X(t) \right| \tag{35}$$

where Δ*X(t)* denotes the di fference between the prey's position and the position of the hawks in the present iteration. *J* is the strength of the prey's random escape.

When |*E*| < *0.5* and *r* ≥ *0.5*, as the prey has no energy to escape, a hard siege will be performed. Accordingly, the hawks hardly circle the prey to execute a surprise attack. Therefore,

$$X(t+1) = X\_{best}(t) - E[\Delta X(t)]\tag{36}$$

Di fferent stages of the HHO algorithm are represented in Figure 11. Moreover, detailed explanations and mathematical formulations of soft and hard siege manners are provided in [32,33].

**Figure 11.** The stages of the HHO algorithm [33].

#### *6.2. Atom Search Optimization (ASO)*

Atom search optimization is a novel physics-based method for global optimization problems [34,35]. ASO is a population-based technique, which imitates the motion of atoms subjected to interactions and other constraint forces [33,34]. The total forces of interaction that a ffect the *i*th atom in the *d*th dimension are described as follows:

$$F\_i^d(t) = \sum\_{j \in \text{Kbest}} rand\_j F\_{ij}^d(t) \tag{37}$$

where *randj* is a number randomly distributed in the range [0, 1] and *Kbest* is a group of *K* atoms having the best values fitness function. The value of K has to be decreased gradually as the iterations advance, to let ASO perform more exploitation in the latest iterations, as shown in the following expression:

$$K(t) = N - (N - 2) \times \sqrt{\frac{t}{T}}\tag{38}$$

where *N* denotes the total number of atoms forming an atomic structure, *t* denotes the present iteration, and *T* is the total number of iterations. The interaction force Fdij presents the gradient of the potential Lennard-Jones (L-J) and is described in the following formula [29]:

$$F\_{ij}^d = -\eta(t) \left[ 2\left(h\_{i\bar{j}}(t)\right)^{-13} - \left(h\_{i\bar{j}}(t)\right)^{-\mathcal{T}} \right] \frac{\overrightarrow{r}\_{ij}}{r\_{i\bar{j}}} \tag{39}$$

where η*(t)* presents the depth function, which controls the regions of attraction and repulsion; *hij(t)* is the ratio between the distance between two atoms *rij* and the scaled distance between them σ(t).

$$h\_{ij}(t) = \frac{r\_{ij}}{\sigma(t)}\tag{40}$$

$$\eta(t) = a \left( 1 - \frac{t - 1}{T} \right)^3 e^{-\frac{2\mu t}{T}} \tag{41}$$

where α denotes the depth weight. The function of the force of interaction under different values of η with respect to a variable scaled distance (*h* changed from 0.9 to 1) is presented in Figure 12.

**Figure 12.** The behavior of the interaction force function with the scaled distance (*h*) under different depth values [34].

Equation (42) presents the expression of the scaled distance between any two atoms *i* and *j*:

$$h\_{ij}(t) = \begin{cases} h\_{\text{min}} & \frac{r\_{ij}(t)}{\sigma(t)} < h\_{\text{min}} \\ \frac{r\_{ij}(t)}{\sigma(t)} & h\_{\text{min}} \le \frac{r\_{ij}(t)}{\sigma(t)} \le h\_{\text{max}} \\ \hline \end{cases} \tag{42}$$

where *hmin* and *hmax* describe the upper and lower limits of the distance (*h*), respectively, and are determined by the following formula:

$$\begin{cases} h\_{\text{min}} = \mathcal{g}\_0 + \mathcal{g}(t) \\ \qquad h\_{\text{max}} = u \end{cases} \tag{43}$$

$$g(t) = 0.1 \times \sin\left(\frac{\pi}{2} \times \frac{t}{T}\right) \tag{44}$$

$$\sigma(t) = \left\| \mathbf{x}\_{ij}(t) \right\| \frac{\sum\_{j \in \text{Kbest}} \mathbf{x}\_{ij}(t)}{K(t)} \Big\|\_{2} \tag{45}$$

where *g0* and *u* denote the lower and upper limits, respectively. *g(t)* denotes a drift factor to give the algorithm the ability to transform from the exploration to the exploitation stage.

Presuming that there is a covalent bond between each atom in ASO and the best one, the constraint force resulting from this combination be expressed as follows:

$$\mathbf{G}\_{i}^{d}(\mathbf{t}) = \boldsymbol{\lambda}(\mathbf{t}) \Big(\mathbf{x}\_{\text{best}}^{d}(\mathbf{t}) - \mathbf{x}\_{i}^{d}(\mathbf{t})\Big) \tag{46}$$

where

$$
\lambda(t) = \beta \times e^{-\frac{2\Lambda}{\tau}} \tag{47}
$$

where *xdbest*(*t*) denotes the location of the best atom in the *dth* dimension, λ*(t)* indicates the Lagrangian multiplier, and β denotes a multiplier weight. After defining the forces of interaction and L-J potential resultant constraint force, the acceleration, to which the *i*th atom is projected in the dimension *d* and iteration *t*, is determined as follows:

$$\begin{array}{lcl} a\_i^d(t) &= \frac{\overline{F}\_i^d(t)}{m\_i^d(t)} + \frac{G\_i^d(t)}{m\_i^d(t)} \\ &= \alpha \Big(1 - \frac{t-1}{T}\Big)^3 e^{-\frac{2\lambda t}{T}} \\ &\times \sum\_{j \in \text{Kbest}} \frac{\text{rand}\_j \Big[2\left(h\_{ij}(t)\right)^{-13} - \left(h\_{ij}(t)\right)^{-T}\Big]}{m\_i(t)} \\ &\cdot \frac{\left(\overline{x}\_j^d(t) - \overline{x}\_i^d(t)\right)}{\left\|\overline{x}\_j^d(t), \overline{x}\_j(t)\right\|\_2} + \beta e^{-\frac{2\mu}{T}} \frac{\left(\overline{x}\_{\text{best}}^d(t) - \overline{x}\_i^d(t)\right)}{m\_i(t)} \end{array} \tag{48}$$

where *mdi* (*t*) denotes the mass of the atom number *i* in the dimension *d* of iteration *t*, and is determined from the value of its fitness function as described below:

$$M\_i(t) = e^{-\frac{Fit\_j(t) - Fit\_{\text{heat}}(t)}{Fit\_{\text{mvert}}(t) - Fit\_{\text{heat}}(t)}}\tag{49}$$

$$m\_i(t) = \frac{M\_i(t)}{\sum\_{j=1}^{N} M\_j(t)}\tag{50}$$

where

$$Fit\_{\text{best}}(t) = \min\_{i \in \{1, 2, \dots, N\}} Fit\_i(t) \tag{51}$$

$$Fit\_{\text{worst}}(t) = \max\_{i \in \{1, 2, \dots, N\}} Fit\_i(t) \tag{52}$$

Finally, the velocity of the *i*th atom in the search space and the new position acquired at the next iteration (*<sup>t</sup>*+*1*) can be described by the following formulas:

$$v\_i^d(t+1) = rand\_i^d \cdot v\_i^d(t) + a\_i^d(t) \tag{53}$$

$$\mathbf{x}\_{i}^{d}(t+1) = \mathbf{x}\_{i}^{d}(t) + \upsilon\_{i}^{d}(t+1) \tag{54}$$

 For a deeper explanation of the atom search optimization technique, the procedure is discussed in detail in [34,35].

#### **7. Results and Discussion**

The simulation results have been obtained to validate the effectiveness of the proposed optimization techniques in determining the optimal values of the capacitance and the parameters of the DC voltage regulator. Moreover, to validate the accuracy of the those obtained from the results, these values have

been applied to the system under study under the integration of nonlinear loads and the condition of three-phase imbalance. The performance of the MMC STATCOM in the two study cases has been introduced in the simulation results.

#### *7.1. Nonlinear Load Case Study*

For the simulation purpose, a dedicated software program for the optimal estimation of the capacitor value, *<sup>K</sup>*p and *K*i, is developed in MATLAB/Simulink. The results of the proposed methods are compared with those of the conventional particle swarm optimization (PSO) [37–39]. The design parameters and the limits of the optimized variables are identical for the all optimization techniques proposed in this paper. The end value of iterations has been adjusted at 10 iterations and the search agents are proposed to be 10 agents. The convergence curves of the optimization algorithms are shown in Figure 13. The detailed results of the optimization process for ASO, HHO and PSO are summarized in Table 3. As shown from Figure 13, the HHO algorithm reached the optimum values of the capacitors and coefficients for the PI controllers within two iterations, and ASO reached the best value of the fitness function after four iterations. While the conventional PSO algorithm is the slowest among them, moreover, its finding of the best value of the fitness function (*F*) is achieved after eight iterations, and it is worse than that with the HHO and ASO. Both of the proposed algorithms were able to obtain convergen<sup>t</sup> results and satisfactory results to improve the voltage and reduce the THD of the current in the studied network, but HHO reached the optimum solution in less time.

**Figure 13.** The curves of ASO, HHO, and particle swarm optimization (PSO).


**Table 3.** Optimization parameters of the proposed optimization methods.

The optimal values for the capacitors of MMC and PI controller constants obtained using HHO and ASO have been applied to the network under study. The profile of the MMC DC voltage for both algorithms is shown in Figure 14, where it can be noticed that stability of the DC voltage is achieved in both cases.

**Figure 14.** The DC voltage of the MMC using HHO and ASO.

The THD of the current in the original state before optimization was 10.93% in the moment of the connecting MMC STATCOM, as shown in Figure 15, while the THD analysis of the grid current is shown in the Figure 16a, where the value of the THD is decreased to 0.60% with the optimal values of the capacitor converters and the optimal values of the PI controllers as determined using HHO. Moreover, the THD of the current with the optimal values of the capacitor converters and the optimal values of the PI controllers as determined using the ASO algorithm decreased to 1.21%, as shown in Figure 16b. Additionally, Figure 16c shows a comparison between the THD analysis of the current grid with di fferent methods (HHO, ASO and PSO).

**Figure 15.** Total harmonic distortion (THD) analysis of the grid current in the original case.

**Figure 16.** (**a**) The THD analysis of the current grid using HHO; (**b**) The THD analysis of the current grid using ASO; (**c**) Comparison between the THD analysis of the current grid with different methods.

The improvement of the voltage profile by compensating the reactive power is one of the main prospective results of MMC STATCOM. In order to examine the operation of the MMC STATCOM with the reactive power command, Qref is changed from (−5 MVAr); in this case the STATCOM operates in the inductive mode, to (+10 MVAr), where the STATCOM operates in the capacitive mode.

The grid side voltage and current of the MMC STATCOM are shown in Figure 17. As shown from Figure 17, when the set-point of the reference reactive power changed at 0.15 s, the current drawn from the grid steps from lagging the grid voltage by 90 degrees to leading the voltage by 90 degrees. The control system of the STATCOM introduces a rapid reaction for modifying the terminal voltage of the inverter in order to produce 10 MVAr of reactive power (capacitive mode). The load current and output voltage are represented in Figure 18.

**Figure 17.** The grid side voltage and current.

**Figure 18.** The load side voltage and current.

The active and reactive power with the connection of the MMC STATCOM is presented in Figure 19a,b respectively. From the figure, a high dynamic performance with acceptable raise time has been achieved for the reactive power variation. Moreover, the three-phase voltage of the converters in this case of study is shown in Figure 20.

(**a**) **Figure 19.** *Cont.*

**Figure 19.** The reactive and active power with connection to the MMC STATCOM. (**a**) Reactive power; (**b**) Active power.

**Figure 20.** Phase voltage converters.

#### *7.2. Unbalanced Load Case Study*

For more validation of the proposed configuration, the system has been tested under the case of unbalanced load. In this case of study, the load on phase C has been assumed to be changed at 1 s by 3 s, and the current at the grid side with imbalanced load is shown in Figure 21. The effectiveness of the proposed control scheme based on the optimization technique of HHO is examined. The voltages and currents at the load side are shown in Figure 22. Figure 23 shows the grid side voltages and currents, which indicate a balance in the grid currents and voltages. Moreover, Figure 24 shows the DC link voltage of the MMC for this case of study. As observed from the figures, in the case of imbalanced load, the proposed MMC STATCOM with the optimal design based on the optimization algorithms was able to maintain the voltage and currents as balanced as possible. Moreover, the maintenance of the voltage in the DC-link, constant without fluctuations, was also desired, and it has been achieved as observed in Figure 24.

**Figure 21.** The current at the grid side with an imbalanced load.

**Figure 22.** The magnitude voltage and current load with an imbalanced load.

**Figure 23.** The voltage and current at the grid side with an imbalanced load with MMC STATCOM.

**Figure 24.** The voltage DC link of the MMC with an imbalanced load.
