Optimal Allocation of STATCOM to Enhance Transient Stability Using Imperialist Competitive Algorithm

1  Introduction

Generally, the concept of FACTS is utilizing the power electronics-based devices on the high-voltage side of the grid, which can be electrically controlled. In this regard, some specific methods are proposed for developing and using the power electronic-based devices aiming for controlling power distribution on the high-voltage side of the grid. The FACTS devices have two main goals as follows [1]:

•   Enhancing the power transmission ability of transmission systems

•   Passing power through desired paths

The first goal means that if the system stability can be maintained during the error occurrence, the transmission power can be increased until approaching the thermal limit by applying the real-time power control. This does not mean that the lines always operate at their thermal limit (because the power losses will be high and unacceptable); however, this can be applied in emergency conditions. Using FACTS controllers instead of considering a large stability margin will significantly increase the power transmission in normal condition. One of the limitations of exploiting the transmission system is that increasing the transmission power of a line in a system will decrease its stability [2].

If the transmission power gradually increases, the system suddenly becomes unstable at a specific level of transmission power. This level is the maximum power that can be transmitted in a steady state. This bound is not an irreversible value stabilized by designing a synchronous motor and transmission line; this can change by different factors. The main factors affecting the bound are synchronous motor excitation, number and connection of lines, types of synchronous motors connected to the grid and compensating devices. Compensators enable the system to pass more power through the lines in the stability limit by injecting the reactive power; as a result, the loadability of the power system is increased. The second goal defines that if the line current is controllable, for example, by changing effective impedance of the line, the power can be transmitted through arbitrary power transmission corridors and restricts the power transmission. Besides, this goal implicitly expresses that a rapid change of power transmission path should be possible in the case of event occurrence so that the power transmission can be done desirably in the whole transmission system. Fulfilling the two above-mentioned goals demand the improvement of controllers and high-power compensators [3–6].

In recent years, voltage collapse has become one of the main concerns in the field of power grids. There is a direct relationship between several large turnouts in countries like Japan, Italy, and US and the voltage collapse. The voltage instability occurs whether a high load is fed into the grid or a local reactive power shortage occurs. Several attempts have been made to overcome this challenge and ensure voltage stability in the grid, and eventually, the FACTS devices were invented. It was found that these devices are useful tools for enhancing the voltage static stability margin in the power grid and can operate near the steady-state stability limit. Another advantage of the FACTS devices is increasing the capacity of transmission lines. Of course, exploiting the advantages of FACT devices is highly dependent on their placement. The optimal placement of FACTS devices has attracted the attention of many researchers in the field of power systems, and the researchers have proposed different solutions [7–10].

There are different types of mathematical methods for solving optimal placement such as the modal analysis method. The modal analysis or eigenvalues is used to identify the busses near the voltage collapse point or recognize the point at which the power feeding is more useful [11]. In [12], the authors have used a combination of dynamic and static methods to optimize the compensators’ location, type, and size. This method is based on the continuous and optimal load distribution and analyzes the eigenvalues. The results show that the proposed method was useful in coordinating the performance of FACTS performance to increase the power system dynamic. Besides, in [13], authors have investigated the location and number of compensators in the power market using the mixed integer nonlinear program (MINLP) method. It should be noted that this method has some problems in detecting optimal points when there are some local extremums. Despite the methods mentioned above, optimization methods have also been used for solving this problem. The basic of these methods is using optimization algorithms. The genetic and particle swarm algorithms fall in this type of method [14–16].

In [17], the ant colony and genetic algorithms are used to optimize the location and capacity of STATCOM for enhancing voltage stability. The voltage stability is a critical problem in power systems, especially in the load peak. Therefore, the authors have used the STATCOM with equivalent-current-inject (ECI) to compensate for the reactive power and have optimized it using the above-mentioned algorithms. Finally, they evaluate IEEE 30-bus test, which shows enhanced voltage stability. In [3], the particle swarm optimization (PSO) method is used to determine the optimal location and controller parameters of STATCOM. First, a systematic procedure has been provided to determine the optimal location of STATCOM in order to enhance the transient stability after severe disturbances. The proposed algorithm is applied in a two-area system that uses the PSO for finding the optimal location of STATCOM. In addition, the STATCOM-based damping controller is presented, and the controller parameter along with PSSs was determined synchronously by using the PSO. The nonlinear simulation results reveal that the optimal STATCOM has extended the cleaning time of the critical error, which as a result, enhances the transient stability of the system. Besides, a coordinated design of STATCOM-based controllers and power system stabilizer (PSS) which use this method can significantly enhance the system stability. This is done by rapid damping of the power system oscillations under small and severe disturbance conditions. Eventually, the coordinated design problem has been extended to the four-machine two-area power system.

In [4], the advanced PSO algorithm has been used to determine the desired size and appropriate assignment of STATCOM in a power system. This algorithm has successfully found the optimal settings for current infrastructures as well as the optimal locations, size, controller settings, and STATCOM settings. The effectiveness of this algorithm has been tested in an IEEE 30-bus test standard system. The results show that the proposed algorithm is able to find the best solution with high statistical significance and high-level of convergence. In [5], a time-varying acceleration coefficients particle swarm optimization (TVAC-PSO) method has been proposed to find the optimal location and size of STATCOM to increase the voltage safety in probable line outage conditions in the power systems. The most severe line outage conditions were considered, then the PSO-TVAC optimization algorithm was used, and the efficiency of the method was evaluated in IEEE 30-bus test standard. This method can be performed for large-size power systems. More target functions can be considered instead of the two target functions of this study.

In [6], a novel method has been proposed for detecting the optimal location of STATCOM and static var compensator (SVC). The power system which operates in a stability boundary lost its stability by occurring disturbance. Given the probable outage, knowing the instability voltage in the assignment of FACTS devices is significant. This article defines the optimal location for installing SVC and STATCOM by considering the power system possibilities and their effects on the constant voltage. The stability voltage of the Khuzestan’s power grid, as a case study, and the optimal location for installing have been studied. By analyzing the location of the weakest bus in constant conditions, it was found that the STATCOM has a higher voltage stability margin than the SVC.

In [7], the quantitative evaluation and appropriate installation of SVC have been addressed to enhance the transient stability in a post. The evaluation procedure includes recognizing areas with low voltage safety and key errors, creating unsafety voltage in power system based on the quantitative analysis and installing SVC, which brings stability and safety for a power system.

In [8], the unified power flow controller (UPFC) locating has been performed to increase the voltage safety and reduce the congestion using the PSO algorithm. In addition to the optimal location, the optimization of size and number of UPFC has been performed. This method recognizes the compact transmission line and weak BUSs as the candidate places for UPFC installation based on the maximum probable intensity index and the L index. The optimal location of UPFC for reducing the congestion and enhancing the safety of voltage is modeled as an optimization problem and has been solved by the PSO algorithm. This algorithm is an effective and practical method for selecting and assigning UPFC in electrical power systems.

Damping the oscillations of the power system between the interconnected areas is vital for the safe operation of the system. In addition to PSSs, the FACTSs are also used for improving system stability. Specifically, utilization of traditional PSS in multi-machine systems cannot lonely provide enough damping for inter-area oscillations. In such cases, the flexible AC transmission system (FACTS) with power oscillation damping (POD) have been considered as a more effective solution. In addition, given that the electricity market has been deregulated in recent years, the traditional concepts and procedures used in the power systems have been changed. Therefore, the available power systems should be used in a more efficient way so that their capacities increase by installing FACTS devices. These devices play an increasingly significant role in the operation and control of power systems. Nevertheless, the non-coordinated local control of PSSs and FACTS devices can lead to unstabilize interactions. Various works have been done regarding the coordination between FACTS POD and PSSs controller to enhance overall system performance. Some of these methods are based on complicated nonlinear simulations, while the others are based on the linear power system model.

In [9], an optimization-based algorithm has been proposed for simultaneous coordination between some controllers. This algorithm uses a hybrid second-order programming method to enhance the system’s overall performance. The interaction between these controllers will enhance by minimizing the target function in which the effects of PSSs and FACTS POD controllers have not been considered. Therefore, the overall performance of the system will be optimized.

In [10], the investigation of three heuristic methods (GA, TS, SA) in the optimization of FACTS devices in a power system has been addressed. The optimization has been done on three factors, including location, size, and type of FACTS devices. Five FACTS controllers have been used for modeling the steady-state. The proposed method was performed in the IEEE 118-bus test. The results show that these three convergence algorithms have similar optimal solutions. However, the TS and GA algorithms are converged faster than the SA. The simulation shows that the FACTS devices can enhance the system stability margin, while some limitations may be observed.

In [16], the nondominated sorting genetic algorithm (NSGA) has been proposed to optimize the location and size of SVC to minimize power loss and enhance the voltage profile. The desired placement of SVC can considerably increase the voltage safety in a power system. The effectiveness of the proposed method has been tested in the IEEE 30-bus test. It was found that this algorithm could be implementable in the power system.

In this paper, the ICA algorithm has been used to determine the optimal location and size of STATCOM by considering the transient stability of the power system. In Section 2, the STATCOM is introduced. The introduction of the ICA algorithm and the problem analysis procedure are addressed in Section 3. In Section 4, the modeling of the problem is presented. In Section 5, the discussion of the results of simulations performed in MATLAB software and the plots’ analysis are provided. Finally, the conclusions and future works are provided in Section 6.

2  Static Synchronous Compensator (STATCOM)

The basis of STATCOM performance is the same as the synchronous condenser. Since the power electronic devices are used in building such devices, they are called static compensators. The converters used in this compensator provide the required reactive power locally (at the grid-STATCOM connection location), and its output can be continuously regulated. For this reason, in some scenarios in which the power grid voltage is widely changed (when a disturbance occurs or after error correction), this compensator is used. Fig. 1 shows a schematic of STATCOM and its V-I characteristic. Production or absorption of reactive power by the voltage source converter (VSC) is done by regulating the Vref voltage.

Figure 1: STATCOM and its V-I characteristic

The most important applications of STATCOM include:

•   Voltage dynamic control

•   Transient stability improvement

•   Eliminating power oscillations in power grid

•   Controlling both reactive and real powers

3  Imperialist Competitive Algorithm (ICA)

In optimization, the goal is to find an optimal answer in terms of the problem variables. Here, an array of problem variables should be created, optimized by the target. In the ICA, this array is called a country. In an optimization problem with N variables, a country has N features, each being equal to a variable. These features determine the power of that country (the corresponding cost function). In the optimization problem, this algorithm seeks for the best country. In fact, finding this country is equal to finding the best parameters of the problem, which produces the minim cost function.

To start the algorithm, an N country should be created. These countries are the possible answers to the problem and are equal to the chromosome in the genetic algorithm and particle in the particle swam optimization.

The Nimp of best members of this population, i.e., countries with a minimum cost function, are selected as imperialists. The remaining countries form the colonies, each belonging to an empire. First, some colonies, based on the power of the imperialist, are given to each imperialist. Therefore, the total countries are divided into two groups: the majority group called a colony, and the minority group called empire or colonialist. An emperor, along with its colonies, is called an empire. As stated earlier, this algorithm seeks to find the best country, and to this end, some policies (operators) are defined. The main policies (operators) of this algorithm are:

•   Assimilation (attraction)

•   Translational motion or revolution

•   Colonial competition

•   Eliminating the weakest empire

The two first policies are done in every empire, while the third one is done between empires. Finally, the fourth policy is applied to the weakest empire.

Under the assimilation policy, every colony moves toward its colonialist to obtain the empire features. Under the revolution policy, some of the colonies of every empire are located in a new location which is adopted randomly. Under the colonial competition, the weakest empire loses one of its colonies, and as a result, its power will degrade. Finally, under the elimination policy, the empire with no colony will be eliminated and itself become a colony of other empires.

The power of an empire depends on the power of the empire and all its colonies. Finally, the power of every empire is equal to the power of empire plus to a percentage of its average colonial power. The pseudocode of the ICA flowchart shown in Fig. 2 is as follows:

Figure 2: Flowchart of ICA

•   Select some random points on the function and forms the initial empires.

•   Moves the colonies toward the imperialist countries (assimilation policy).

•   A small number of colonies of each empire are moved to a new position that is selected as randomly (revolution policy).

•   If there is a colony in an empire who has a cost less that the imperialist, switch the location of colony and imperialist.

•   Calculate the cost (power) of whole empire (by considering the cost of imperialists and their colonies).

•   Select one colony from the weakest empire, and give it to the empire with maximum possibility for taking it up.

•   Eliminate the weakest empire who have no colony.

•   If two empires are getting very close to each other, merge them and forms a single empire.

•   If only a single empire is left, stop, otherwise move back to 2.

By having cost (power) of all imperialists, the normalized power can be considered as follows:

Cj=fcost(imp,j)−(fcost(imp,i))(1)

where fcost(imp,j) is the power of nth imperialist, and Cj is the normalized power of jth imperialist. The relative normalized power is calculated by Eq. (2) so that one can divide the colonies between imperialists.

Normalized Powerj=Cj∑i=1NimpCi(2)

where Nimp is the number of empire. From another perspective, the normalized power of an imperialist is the ratio of colonies that are administered by the imperialist. So, the initial number of colonies of an imperialist is expressed as follow:

NCj=Round(|Normalized Powerj|.Ncol)(3)

where NCj is the initial number of colonies of an empire, and Ncol is total number of colonial countries. Besides, Round is a function that gives the nearest integer to a decimal number.

The power of jth empire, TCj, is obtained as follows:

TCj=fcost(imp,j)+ξ.∑i=1NCjfcost(col,i)NCj(4)

where is a value between 0 and 1. The normalized power of jth empire, NTCj, is expressed by:

NTCj=TCj−(TCi)(5)

By using the power of jth empire, the possibility of conquering a colony by the jth empire can be obtained by following equation.

Pj=|NTCj∑i=1NimpNTCi|(6)

Pj is used in empires’ competition policy to conquer a colony of the weakest empire. By having the possibility of conquering each empire, a mechanism like the roulette wheel in the genetic algorithm is needed to give the colony to an empire based on the power of empires.

The assimilation policy (see Fig. 3) is defined as Eq. (7) in the ICA, which leads to move the colonies toward emperor:

{x}new={x}old+U(0,β×d)×{V1}(7)

where {x}new and {x}old are the new and old location of each colony. Parameter β is a parameter higher than one denoting the relocation step. Parameter d represents the distance between a colony and its empire. U(0,β×d) is a random number with a uniform distribution between 0 and β×d. Finally, {V1} is a unit vector with same direction with the path between colony and its corresponding empire.

Hint: in the modified method, the obtained value of U every time is minus from one (1-U), and for these points, the algorithm also performs the seeking and calculating.

Figure 3: Assimilation policy

The algorithm is repeated until the convergence condition is fulfilled and/or the total number of iterations is reached. After some time, all empires fall, and only a single one remains, and the other countries will be under the control of the remained empire. In this new ideal world, all colonies are administered by a single empire, and the location and cost of colonies are equal to the location and cost of the imperialist country. In this new world, there is no difference between the colonies and between the colonies and the imperialist country. In other words, all countries at the same time are colonies and colonialists. In such situations, the imperialistic competition is over, and the algorithm is stopped [18–26].

4  Problem Formulation

The power transmission capacity and the system’s transient stability can be enhanced by changing the location of STATCOM from the middle point toward the sending point. For certain initial operational conditions, there is an optimal location for STATCOM at which the maximum power of sending point related to Section 1 of the transmission line is equal to the maximum power of receiving point related to Section 2 of the transmission line. In view of transient stability, the maximum power of sending point is considered a power of productive terminals that can be received immediately after an error affecting the transient stability. So, the target function for preserving the transient stability can be defined as follows:

J=|MAX(Δδ1−Δδ2)|(8)

where Δδ1, and Δδ2 are equal to the angular deviation of the rotor after the generators’ disturbance in areas 1 and 2, respectively. Besides, |MAX(Δδ1−Δδ2)| is the absolute value related to the maximum difference between the rotor’s angular deviations for both areas.

If |MAX(Δδ1−Δδ2)| is less than 180 degree, the system is considered as stable. The time-domain simulation of the nonlinear system model is used to calculate the target function. The limitations of the problems include the bounds of the desired location. Therefore, the design problem can be formulated in the form of an optimization problem under the following limitations.

Lmin≤L≤Lmax(9)

where L is a part of length from the sending point to the STATCOM location.

The following steps should be followed for seeking the optimal location for STATCOM in order to enhance the transient stability.

Step 1: firstly, a fault-clearing time (TFC) is selected high enough so that the system be unstable in all STATCOM locations.

Step 2: using the ICA to minimize the target function J

Step 3: examining the system stability

Step 4: if the system is unstable, the TFC value is reduced and Step 2 is repeated. Stop if the system becomes stable.

At TFC = TFCF, system is stable if the STATCOM is located in the optimal location obtained from the above-mentioned method, and for TFC > TFCF, the system is unstable in all locations.

The PSS block related to Sim Power System (SPS) tool is used to add damping to rotor oscillations related to the machine by controlling its stimulation. The output signal related to PSS is used as the additional input (VS) for the stimulation system block. The PSS input signal can be either machine speed deviation or the accelerator power. The PSS model uses a low pass filter, an overall gain, washout high pass filter, a phase compensator system, and an output limiter (see Fig. 4). The overall gain, KP, determines the damping produced by stabilizing. The washout filter eliminates the low frequencies that exist in the input signal and enables PSS to respond only to the input changes. The compensator system is used to compensate for the phase delay between phase stimulation and the electrical torque of the synchronous machine [27–32]. It should be noted that the PSS and STATCOM-based controllers have been designed to minimize the power system oscillations after a huge oscillation to enhance the power system stability (Fig. 4).

Figure 4: Structure of power system stabilizer

5  Simulation and Results Analysis

The MATLAB/SIMULINK model has been designed using a set of SPS blocks to determine the optimal location and size of STATCOM. The system under study comprises two hydraulic generator units, one with power of 1000 MVA in area 1 and the other with the power of 5000 MWA in area 2. A six-order model (shown in Fig. 5) has represented the generators. Both generators are equipped with a hydraulic turbine and governor (HTG) and a stimulation system. The HGT represents a nonlinear hydraulic turbine model, and a proportional–integral–derivative controller (PID) system, and servomotor. The stimulation system is composed of a DC voltage regulator without an actuator saturation function. Generators with output voltages of 13.8 kV are connected by a 700 km three-phase transmission line with zero sequence impedance 0.01755Ω/km, inductance line 0.873 mH/Km, capacitor line 13.33 nf/km through three-phase transformators. The output voltage of transformator is 500 kV. There is a three-phase error with error resistance of 0.001 Ω in area 1 and a resistance load with active power 5000 MW and voltage 500 KV in area 2. The main system is returned to its initial condition after cleaning the error. The initial power outputs related to the selected generators are P1 = 0.95 pu, and P2 = 0.80 pu. The STATCOM used for this model is a phasor model with nominal power of 200 MVA. The total frequency of system is 60 HZ.

Figure 5: Two-motor system designed with MATLAB simulink

First, the STATCOM is eliminated from the system, then a three-phase error is applied to the system at t = 1 s, and the error is cleaned at t = 1.1 s. It can be seen that the system represents stable behavior after cleaning the error. As can be seen from Figs. 6–11, after the occurrence of an error, the rotor angular difference, rotor speed, voltage of terminals 1 & 2, and buss’s voltages move toward instability, and as a result, the system become unstable. The wave shape of voltages and power after error cleanness will be as undamped oscillations, denoting the system’s instability (Figs. 6 to 11).

Figure 6: Variations of rotor angular difference in the absence of STATCOM

Figure 7: Variations of rotor speed (W1,W2) in the absence of STATCOM

Figure 8: Variations of voltage of terminals 1 and 2 in the absence of STATCOM

Figure 9: Variations of rotor speed difference in the absence of STATCOM

Figure 10: Variations of positive sequence voltage B1, B2, and B3 in the absence of STATCOM

Figure 11: Variations of line power in the absence of STATCOM

In the next step, the STATCOM is embedded in the system, and the cleaning time will be longer. It can be seen that, for TFC > TFCF, the system is unstable in all locations. As can be observed from Figs. 12 to 17, after the occurrence of error and by considering the cleaning time of 1.15 s, the rotor angular difference, rotor speed, line power, voltage of terminals 1 & 2, and the bus’s voltage are unstable. In this case, if STATCOM is placed anywhere in the system, the system will not be stable.

Figure 12: Variations of rotor angular difference in the presence of STATCOM for TFC > TFCF

Figure 13: Variations of rotor speed (W1,W2) in the presence of STATCOM for TFC > TFCF

Figure 14: Variations of voltage of terminals 1 and 2 in the presence of STATCOM for TFC > TFCF

Figure 15: Variations of rotor speed difference in the presence of STATCOM for TFC > TFCF

Figure 16: Variations of positive sequence voltage B1, B2, and B3 in the presence of STATCOM for TFC > TFCF

Figure 17: Variations of line power in the presence of STATCOM for TFC > TFCF

In the next step, the STATCOM is located in the middle point of the transmission line, and the cleaning time is reduced. The following values are used:

L = 700 M    L1 = 350 M    L2 = 350 M

SF = 200 MVA  C = 350 μF    Tclear = 1.1 s

It can be observed that the waves are damped in an oscillation form, and the system will be stable (Figs. 18 to 23). According to the target function of the system transient stability, a system is stable if:

Figure 18: Variations of rotor angular difference in the presence of STATCOM without optimization

Figure 19: Variations of rotor speed (W1,W2) in the presence of STATCOM without optimization

Figure 20: Variations of voltage of terminals 1 and 2 in the presence of STATCOM withoutoptimization

Figure 21: Variations of rotor speed difference in the presence of STATCOM without optimization

Figure 22: Variations of positive sequence voltage B1, B2, and B3 in the presence of STATCOM without optimization

Figure 23: Variations of line power in the presence of STATCOM without optimization

|MAX(Δδ1−Δδ2)|<180

It can be observed from Fig. 18 that MaxΔδ=109.2 deg, and the setting time is 4.48 s. In Fig. 19, MAXW1=1.011,MAXW2=1.001 pu. In Fig. 21, MAX(W1−W2)=0.012 pu. Finally, in Fig. 23, the maximum line power is 1405 MW, and the setting time is 2.82 s.

By observing the wave shapes and creating the error at t = 1 s, some oscillations and noises have occurred in the system, resulting in instability in the grid. However, these oscillations have been damped by the presence of STATCOM, and the grid reaches stability. As can be seen from the above figures, at the moment of error occurrence, a severe increase in generators’ speed is occurred, with or without STATCOM. Besides, the angular difference of generator rotors is increased, but this increase is not applied suddenly and as impact because of generators’ inertia; however, the buss voltage and generator terminals experience a severe step drop, and some electromagnetic transients are created. It should be noted that in mechanical systems, these transients are significantly lower due to the inertia of such systems.

In the next step, the optimal location and size of STATCOM were obtained by using ICA and the steps provided in Section 4. The procedures of the ICA toolbox have been used to optimize the target function. Some certain numbers of parameters are needed to implement the ICA technique. Table 1 shows the parameters used in the ICA algorithm.

After implementing the ICA optimization procedure, the optimal size of STATCOM and longest critical error cleaning time are shown in Table 2. Figs. 24 to 29 represent the simulation results considering the parameters optimized by the ICA method.

Figure 24: Variations of rotor angular difference considering for the parameters optimized by ICA

Figure 25: Variations of rotor speed (W1,W2) for the parameters optimized by ICA

Figure 26: Variations of voltage of terminals 1 and 2 for the parameters optimized by ICA

Figure 27: Variations of rotor speed difference for the parameters optimized by ICA

Figure 28: Variations of positive sequence voltage B1, B2, and B3 for the parameters optimized by ICA

Figure 29: Variations of line power the parameters optimized by ICA

By optimizing the problem with the ICA method, the responses of the grid parameters become stable, which is obvious in the above figures. In this case, it is observed that the stability criteria are improved compare to the case without optimization. In other words. By comparing the rotor angular of the scenario in which the STATCOM is traditionally located on the line, it can be observed that the maximum rotor angular different, i.e., overshoot factor, is considerably reduced in ICA method than the traditional method. Besides, the setting time is also reduced for ICA than the traditional ones. These two factors represent the higher stability of the system in the optimization scenario than the traditional ones.

6  Conclusions

In this paper, the ICA technique has been used to determine the optimal location and size of STATCOM. Firstly, a systematic procedure has been used to determine the optimal location of STATCOM for enhancing the transient stability after a severe disturbance. Then, the ICA algorithm is used to find the optimal location and size of STATCOM in a two-area system. The nonlinear simulation results show that the STATCOM at the optimal location results in an increased error cleaning time, and as a result, an enhanced system’s transient stability.

Funding Statement: This work was supported in part by an International Research Partnership “Electrical Engineering—Thai French Research Center (EE-TFRC)” under the project framework of the Lorraine Université d’Excellence (LUE) in cooperation between Université de Lorraine and King Mongkut’s University of Technology North Bangkok and in part by the National Research Council of Thailand (NRCT) under Senior Research Scholar Program under Grant No. N42A640328.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.