Comparative study of SBOAs on the tuning procedure of the designed SMPI controller for MIMO VSP/HVDC interconnected model

This paper presents a new application of advanced SMPI controller for a newly introduced interconnected dynamic system with VSP based HVDC links for frequency control problem. This work presents an outgrowth of analysis about the Swarm – Based Optimization Algorithms (SBOAs) in the tuning process of Multivariable Proportional – Integral (MPI) controller, which is called swarm – based MPI (SMPI). PSO, GOA and GWO algorithms are used for tuning process of the designed SMPI. The VSP based HVDC model is added for mitigation of system frequency dynamics with emulating virtual inertia. The proposed SMPI controller are designed for enhancing the dynamic performance of this system ’ s states during contingencies and they are compared with the conventional designed MPI controller. Deviation characteristics of the step function in MIMO transfer function of the VSP based HVDC model is considered as the common performance index in the proposed algorithms. On the other hand, the role of the proposed SMPI controller is to stabilize such interconnected system while minimizing the associated cost function. Matlab simulations next to the performed Nyquist ’ s stability anslysis demonstrate how the tuned SMPI can remarkably improve the frequency deviations and the damping of the inter-area oscillations excited during a fault. This enhancement is more obvious especially when SMPI controller for a VSP based virtual inertia emulation is tuned using GWO method.


Introduction
The rapid evolution and expansion of the existing AC network into an interconnected low-inertia scenario with a high share of renewablebased generation units, is leading to a system with high complexities, which can impact the dynamics of the overall network and particularly the frequency stability of AC/DC interconnected modern system [1].AC/DC system mainly consists of hybridization of the common Alternative Current (AC) based network and the Direct Current (DC) based components.Nowadays, thanks to significant advancement in power electronic technologies, it became possible to use high-voltage DC (HVDC) transmission systems and power electronic-based components in a very effective manner.HVDC systems overcome the main drawbacks of HVAC transmission, especially for a very long distance application.The frequency control mechanism for such a multi-area modern system, during load and resource variation, is very important.It can consist of primary and secondary frequency control.Primary control is referring to the inertial response of the governors against the frequency mismatch.Its main role is to provide enough inertia for maintaining the stability of the system especially within the next couple of seconds after the frequency drop.Secondary frequency control or AGC is a higherlevel control mechanism that can facilitate different tasks such as frequency restoration, tie-line power flow control between authority areas and economic dispatch of generation units [2]- [3].The interconnection of multi areas can be supported by AC lines of HVDC tie-lines which can enable the transfer of the scheduled power between different zones.It can also give more support and back up to the neighboring area in abnormal situations.Considering the operation limitations of AC lines especially for long-distance, the DC transmission is receiving increasing attention over the last recent years.
HVDC links are known as one of the main utilization of power electronics in modern power systems, which could convey several advantages such as fast and bidirectional control on active and reactive power, power oscillation damping and frequency stability support [4,5].With proper coordination, they can improve the reliability of complex systems and can act as a kind of firewall against cascading disturbances.Due to these advantages, in many projects, the DC links or hybrid AC/DC interconnections are becoming the preferred option [6][7].
Recently, several attempts have been made in the control of grid connected power converted for enhancing the dynamics of low-inertia AC/DC systems by providing virtual inertia to the grid [8][9].
The virtual inertia is emulated by using advanced control of power converters considering their short term energy storage, [10] which gives rise to the possibility of having a higher amount of distributed generation systems, connected to the grid through power converters, without hindering the system stability.However, the application of this kind of converters are not restricted to low power applications or generation parts.For instance, it can be a beneficial option for HVDC power converters.As reported in [11], a novel technique, named VSP based HVDC system, is introduced for inertia emulation capability.Up to now, only the basic model and its mathematical model are reported while no further improvements in tuning and its optimal operation have been reported yet.The AC/DC interconnected system with VSP capabilities is a complex MIMO system that suffers from the lack of advanced control methodology for better operation of the VSP based HVDC system.This paper is trying to fill this gap by proposing a heuristic-based MPI controller for such a system.
Recently, different control methods are presented to achieve accurate and suitable damping effects, such as sliding mode control (SMC) [12][13][14].SMC is basically a roust control technique that is used to control systems with uncertainties.The biggest advantage an SMC possess over PID is the finite-time stability of the closed-loop system (may not be true for discrete time systems).Another important point to remember is that the SMC can reject the disturbances (disturbance decoupling, when properly designed), whereas PID can reduce or remove the effect of uncertainties.In the most common sliding mode control schemes the main limitation is that the performance of the controller relies on the sliding surface.If the sliding surface of SMC is not designed properly, it may lead to unacceptable performance while the applications in real practice is not easy.
Among the controller variety, Proportional-integral (PI) becomes the controller that is most common and the most applied in a physical system [15].The reason is that it has a characteristic that offers simplicity, clear functionality, and ease of use [16].However, Ho et al. [17] reported that only one-fifth of PI control loops are in good condition.The others are not, where 30% of PI controllers are not able to perform well due to the lack of tuning parameters, 30% due to the installation of a controller system operating manual, and 20% due to the use of default controller parameters.
In recent years, many researchers have paid attention to the MPI controllers design for various systems such as in Industrial Scale -Polymerization Reactor [18], Coupled Pilot Plant Distillation Column [19], Narmada Main Canal [20], Quadruple-Tank Process [21], Boiler-Turbine Unit [22], and Wood-Berry Distillation Column [23].Research by Kumar et al. [18] had proposed a synthesis based on the approximation of the relative gain array (RGA) concept to the multivariable process.
The method was further improved by relative normalize gain array concept (RNGA).The controller based on RNGA concept provides better performance than the RGA concept.Both concepts use the nonstandard PI controller which requires Maclaurin series expansion [24].Multivariable PI or PID controller utilization, by attention to all its advantages, has been reported, specifically on industrial applications, as a method in several scientific societies.For achieving a safe, robust and flexible strategy for optimal coefficients arrangement a metaheuristic based algorithms, has been discussed frequently [25][26][27][28][29][30][31].It should be mentioned that the utilized cost functions in these works are focused on the different index such as ITAE and ISE.Some of the main limitations in these mentioned papers are the need for high availability and direct access to the state variables while the distortion of the measurable state variables with the noise is an additional burden in actual implementation.
But in this paper, by attention to the determined features of the step response and its accessibility in every iteration of the simulation, a fitness index is used according to the extracted features of the step function.Also by considering the variety of intelligent optimization algorithms in the tuning process, it is necessary to review all methods as comparative researches.The utilization of swarm-based optimization algorithms, specifically GOA and GWO, will increase the convergence performance in comparison with the applied typical MPI tuning [32][33][34][35][36][37].MPI controller has few parameters to adjust and is relatedly simple to implement.It's robust and can be used in parallel computation.The advantage and necessity of MPI controller are on its feasibility and easy to be implemented.The MPI gains can be designed based upon the MIMO system parameters if they can be achieved or estimated precisely.Moreover, the MPI gain can be designed just based on the system tracking error and treats the system to be "black-box" if the MIMO system parameters are unknown.However, the MPI controller generally has to balance all matrix-gains impact on the whole system and may compromise the transient response.Thus, a complementary action from heuristic based algorithm for such a controller can advance the performance of the controller.This is another important topic that going to be answered, for interconnected system with VSP, in this paper.
The application of swarm based algorithm for better tuning of MPI approach is a very interesting idea which can further enhance the performance of the controller.The first and the most basic presented definition for "swarm" is based on a similar and simple set of agents that are interconnected locally among themselves.So by considering the above explanation, swarmbased optimization algorithms (SBOAs) have been introduced as a subset of natureinspired, crowdbased algorithms which are able to generate low cost, fast and robust solutions for an extensive set of the multiplex problems.Simultaneous to introduce swarmbased computational models, there has been a mutation in the number of scientific activities which presented the successful utilization of SBOAs in a various set of optimization tasks and complex problems.SBOA principles have been successfully achieved in a variety of problems consist of function and structural optimization, finding optimal routes, scheduling and image and data analysis.Swarmbased Computational models have been further used in a wide-range of diverse domains, involving machine learning, bioinformatics and medical informatics, dynamical systems and operations research.Until now, several AI algorithms are performed for achieving to an optimal tuning coefficients method, such as Artificial Bee Colony [38] and Ant Colony algorithm [39][40][41]  controllers.These AI algorithms are created based on a simple structure and therefore they cover some of the optimization problems which are performed for uncomplicated systems.
In this work, a new application of heuristic-based MPI controller design for AC/DC interconnected system with VSP capabilities is proposed.This control approach is based on the evolutionary-based swarm intelligence (SI) optimization algorithms.Particle Swarm Optimization (PSO), Grasshoppers Algorithm (GOA) and Grey Wolf Optimization (GWO) are the main SBOAs, which are used for the tuning process of MPI controller.The focus of this study is on a new design of advanced control for VSP concept through the converter stations of HVDC link, which is achieved by means of an SMPI controller in a multi-machine AGC system.
Another importance of the proposed approach is on the proposition of a unique performance index as part of the optimization process.The performance index is defined as a function consists of a linear combination of step response features such as stability index (SI) and settling time (ST).Finding the best MPI coefficients, as the main target among this process, is realized using the minimization of the proposed performance index.
The proposed advanced control approach is validated through the Matlab simulations of an interconnected AGC dynamic model with HVDC with the VSP capabilities.The proposed method is based on the minimization of a performance index, which is designed by attention to the transfer function characteristics matrix in an MIMO system.The MPI controller design method is performed on the systems to enhance the characteristics performance model versus the conventional PI controller method.
The rest of the paper is organized as the following.In the next section, Section 2, the mathematical formulation of the studied two-area interconnected AGC dynamic model, consists of a detailed model with DC link with VSP capabilities, with an MIMO transfer function characteristics are discussed.Then, in the next Section 3, a typical MPI controller design considering the tuning of its coefficients are presented.In Section 4 swarm-based optimization algorithms (SBOA) considering PSO, GOA and GWO insertions are explained.The simulation results for VSP-HVDC and AC/HVDC model are presented in Section 5.In this section, the robustness of the controller against parameter changes with Nyquist analysis is also presented and compared.Finally, the conclusions and controller remarks are presented in Sections 6.

VSP-based AC/HVDC interconnected system
In this section, the concept of VSP for inertia emulation through the converter stations of the HVDC links in a two-area system will be presented.

Dynamic model of the VSP facilitated system
We deploy a high-level control architecture to simulate and analyze the frequency behavior of the AC/HVDC two-area system equipped with VSP functionalities.In general, the VSP strategy is based on the control of voltage source converters through an active power synchronization loop and a virtual admittance [11].
After an analysis of the electromechanical model of the VSP, the dynamic relationship between input and output power of each converter station could be a second-order transfer function.Then a characteristic equation of the VSP dynamic model can be derived for the description of this second-order system, namely, where ζ is the damping factor,ω n is the natural frequency.The input signal, u, is the reference DC power which is related to other states of the global multi-area system.The output signal consists of two states variables: where Δx p,vsp represents the injected DC power and Δx f,vsp is the derivative term of this power for each VSP, respectively.Based on classic control concepts, this second-order system could be represented by a set of two linear state equations: where Δu denotes the input signal of the DC power reference (ΔP DC, ref ) that comes from a higher-level control.Here Δu is generated by other control signals.For a two-area interconnected system with Area i and Area k, let us assume that a VSP based HVDC link is implemented.As the HVDC line is located between these two areas, the frequency deviation of each area is the most suitable control signal.In the case of parallel AC/DC lines, the AC tie-line power deviation could be used as another control signal in order to achieve suitable coordination between these two lines.Thus we can have: where K fi,vsp , K fk,vsp and K AC are the proportional coefficients for the frequency deviations Δω i , Δω k and AC power flow deviations ΔP tieAC,ik of these two areas, respectively.This state space presentation in (3) will be part of the global model for the frequency dynamics of the two-area AC/DC system with VSPbased inertia emulation.Fig. 1 depicts the high-level control structure of the two-area system with a VSP-based AC/DC transmission.For the frequency regulation, typically a low-order linearized model could be used for modeling the load-generation dynamic behavior, namely the AGC model.Thus to study the small-signal modeling and stability of the high-level control design, we can describe the dynamic equations of these two AC/DC interconnected areas considering the dynamics of the VSP-based HVDC link in the Laplace domain as follows: ) ) 1 + sT pk (6) where ΔP mi and ΔP mk denote total active power from all generation units (GENs) within Area i and k, respectively.The part of AGC control calculates the Area Control Error (ACE) of Area i.The ACE signal now contains the frequency deviations of that area and both AC/DC power flow deviations, and acts as the input for an integral control action, ) where for Area i, β i is the frequency bias.In ( 5) and ( 6), the parts Δx p,vspi and Δx p,vspk denote the total output of the VSP based DC transmitted power.Thus for this one VSP station on the HVDC line in the two-area system, the following equations can be further derived from (3), where and

Mathematical state-space VSP-HVDC model
As it is mentioned in Section 2.1, we can have a global linearized mathematical representation of the two-area AC/DC interconnected power system with VSP-based inertia emulation strategy: All of the required coefficients and constant values are available in [42][43].

MIMO transfer function characteristics
As it was mentioned in this work, simulation results and conclusions will be presented as a MIMO system.Therefore, it is necessary to review this type of systems, specifically as a two-input two-output model.For achieving a common understanding about MIMO system, assume the following block diagram shown in Fig. 2. As it is observed, this figure represents a MIMO system as a two-input two-output model.
By attention to above block diagram, it can obtain: So ( 19) can be rewrite as a matrix form: By considering to the applied controllers in this MIMO model, it can obtain: By replacing (20) in ( 21): Also: Therefore, it can be obtained as follows: and: By replacing and simplifications ( 19) -( 23), finally we can achieve to the MIMO transfer function for a twoinputs twooutputs model: where T ij (i,j = 1,2) are the main elements of transfer function matrix which are calculated individually as follow: where And for MPI controller: Finally, by considering the described equations, MIMO transfer function for VSP/HVDC model can be driven as follows: It is possible to calculate MIMO transfer function for VSP/HVDC model by attention to all constant values in the steady state model (15) and calculation matrixes A vsp , B vsp and C vsp .The obtained transfer functions are presented in Appendix A.

Typical MPI controller design
The typical tuning process of an MPI controller for a linear multivariable stable plant is introduced in [44][45].Therefore, the obtained MPI controller coefficients from this method, (K P and K I ), can be determined by notification to the tuning MP controller coefficient (K P ).

Tuning of MP coefficient
As it was mentioned, for finding K P and K I in MPI controller based on the typical method, it is necessary to determine K P as the main coefficient in MP controller.As Fig. 3, the MIMO system in the existence of a general controller, as G C , can be illustrated.
By attention to Fig. 3 and considering MP controller for G C (s), it can be seen: Then by considering steady state model for the main plan, it is represented by: By assumption initial conditions: Therefore, we have: By replacement to the state equation: By attention to Fourier series definition: By considering t = 0: where K P matrix is given by: where (CB) + is the left inverse matrix for CB and it is given by: It is important to know that in the square matrix: Finally, by attention to ( 41) -(43) and selection k 1 ,…,k m , K P will be achieved in (41).

Tuning of MI coefficient
In this section the control signal, considering the assumed MI controller in Fig. 3 and the equation of G C (s) = K I /s, will be achieved: Now by definition z(t) = ∫ t 0 e(τ)dτ, the closed-loop system equations will be arranged as: So the stable of closed-loop system will be guaranteed by: And finally, when the state matrix is square then the coefficient K I is achieved by: where, Cnsidering Section 3.A and 3.B, also ( 40) -( 44), the conventional MPI (MP + MI) coefficients are available for achieving.

Swarm -Based optimization algorithms (SBOAs)
As explained in the introduction, SBOAs can be called as an individual branch of artificial intelligence (AI) which are applied to implement the communal behavior of social swarms in nature.The social interconnections between swarm individuals can be either direct or indirect [46][47][48].The waggle dance of honey bees can be introduced as an example of direct interconnection which is through visual or audio contact.When one individual can change the situations in its around environment, actually indirect interconnection happens.Also, the other individuals answer to the new environment, such as the pheromone trails of ants that they settle on their way to search for food sources.
In every swarmbased algorithm, it is actually important to address the exploration and exploitation of a search space.The process of encountering new regions of a search space is called the exploration process.Also the process of encountering those regions of a search space within the neighborhood of previously visited points, will be introduced as exploitation process.For achieving to the specific targets and have a successful search algorithm, it is necessary to design a suitable ratio between exploration and exploitation.Fig. 4 can describe these two management concepts by attention to the nature of algorithms.As it is mentioned, some algorithms, such as Steepest Descend, Gradient and Newton Raphson, are located so closed to the vertical axis (exploitation parameter).These algorithms are so depended on the start point.In the other hand, if they have been tuned from suitable start point then the optimization process will be finalized to the best answer.Therefore, such algorithms direct the process greedily.But on the contrary, those algorithms which are so close to the horizontal axis (exploration parameter), will be based on random searches.In the other hand, every algorithm which emphasizes exactly to the exploration parameter, will be impractical in the optimization process.
By considering the above descriptions, it is obvious that a SBOA can be achieved to find the best answer when it uses a correct combination between these two factors.
SBOAs are known as a subset from a bigger set which are called meta -heuristic methods.A meta -heuristic is a high-multiplex problem that presents a set of strategies to expand heuristic optimization algorithms.
Meta -heuristics have been demonstrated by the scientific community to be a viable, and often superior, alternative to more traditional (exact) methods of mixed integer optimization such as dynamic programming.By attention to the recently obtained experiences in complicated or large problem, meta -heuristics are capable to give a suitable trade -off between solution quality and computing time.Fig. 5 presents a total scheme from the metaheuristic algorithms classification based on the nature of applied methods.
As it was described, SBOAs will be classified as one of the most important subsets.By attention to this scheme, Particle Swarm Optimization (PSO), Grasshoppers Algorithm (GOA) and Grey Wolf Optimization (GWO) are the main SBOAs, which are performed in this work.

Particle swarm optimization (PSO)
This methodology, as the first and the most developed SBOAs, is a meta -heuristic approach inspired by the individual and social behavior of the flocking birds [49].In PSO algorithm, particles are known as the population of candidate solutions in optimization problem.Each one of particles is determined by means of two important components which are position and velocity vectors.In this method, every particle tracks its direction towards the total optimum based on Newtonian mechanics.The particle components are defined for the i th particle and in a Ddimensional space: Also, the upgrade rules in this algorithm, according to the components are illustrated below equations respectively: x where × and v are the main components in PSO algorithm.w is applied inertia to achieve the convergence speed and balance between exploitation and exploration.c 1 and c 2 are two positive constants and r 1 and r 2 are random numbers obtained from a uniform random distribution function in the interval [0,1].x pbi is the best position that i th particle have seen and x g is the position of the best particle seen so far among the upgrade process.

Grasshopper optimization algorithm (GOA)
Grasshoppers are classified as a set of insects which generally harm to the crop production as they consider as a kind of agriculture pest.These kinds of insects are usually seen individually in out around environment but mostly they fly as a large swarm in nature.The swarm of grasshopper has one distinctive feature which is its swarming behavior in both the nymph as well as adulthood in grasshopper.The nymph grasshopper flies like rolling cylinders in millions of numbers.Grasshoppers destroy all the vegetation which exists on their path during their motion.After nymph -hood, they become adult and gradually create a swarm in air and then they flight over a very large distance.
In the swarm larvalhood, the swarm of grasshopper generally has a very slow moving.The basic feature in the larvalhood is the small step of the grasshopper.Unlike this duration, the most important feature of swarm in adulthood is long range and sudden motion of swarm.The grasshopper swarm mostly create for the finding of food source.This searching process is another important characteristic for the swarm of grasshopper.As it was explained in the introduction paragraph, exploration and exploitation are two important determined factors in SBOA.Both tendencies are enhanced by grasshoppers naturally in which they fly abruptly as well as locally in small areas.
A scientific model is performed of this behavior of grasshopper swarm to design a nature-inspired optimization algorithm.This mathematical model utilizes to simulate the behavior of the swarm grasshoppers is illustrated as follows [50]: where X i denotes the position of the i th grasshopper, S i is the social interconnection, G i is the gravity force of the i th grasshopper and A i specifies the wind advection.For achieving to the random behavior, (54) can be written as: where r 1 , r 2 and r 3 are random numbers between [0,1].
where d ij is defined as the distance between the i th and j th grasshopper, calculate as: s is a function to explain the strength of social forces, is shown in (58) and d ij is a unit vector from the the i th to the j th grasshopper which can be shown as The s function, which denotes the social forces, is calculated as follows: where f is introduced as the intensity of absorption and is the attractive length scale.The function s indicates the effects on the social interconnection (absorption and repulsion) of grasshopper.The distances are determined from as 0 to 15.The suitable distance for the interlude of repulsion is considered between (0, 2.079).Also the comfortable distance of a grasshopper is 2.079 units from other grasshopper, as there is neither absorption nor repulsion for a grasshopper when it is away from other grasshopper by 2.079 units.This is also called as comfortable zone.Equation (59) also describes a kind of social behaviors which exist among intelligent grasshoppers.This factor changes by means of two parameters, l and f in equation ( 59).The consequences of l and f can be seen on function s, in (59), by considering to vary these two parameters.Therefore, the comfort zone can be changed by tuning of l and f.It is  important that the absorption or repulsion regions are very small for some values (l = 1.0 or f = 1.0 for instance).From all above values, l = 1.5 and f = 0.5 [50].
As it was indicated, totally there are three forces according to location for every grasshopper: absorption, repulsion, and neutral.By considering to these forces, the around space is divided into before the comfort zone, inside the comfort zone, and after the comfort zone.To show the interconnection between grasshoppers by attention to comfort area, Fig. 6 shows a conceptual model.
Also in facilitated structure, the social interconnection can be as the motivating force [51].This fact is observed particularly in some earlier locust swarming models.By attention to the previous explanations, two grasshoppers can divide their around space into comfort zone, absorption and repulsion region by means of function s.There is a zero cross when the distances are greater than 10 the value of function.Also there is a specific limitation in strong forces for the large distances between grasshoppers.In this situation, the force component cannot be utilized using this function.Then it is necessary to keep the distance of grasshoppers and actually map it in the interval.The G component in (54) is calculated as follow: where g is the gravitational constant and e → g shows a unity vector towards the center of the earth.Also for calculation of the A component in (54) as follows: where u is a constant drift and e → w is a unity vector in the direction of wind.It is important that the maximum correlation happens with the direction of wind for the motion of the nymph grasshoppers when they don't have wings.According to the applied substitutes in (54), this equation can be rewrite as: where N is the number of grasshoppers.It is important to keep the position of nymph grasshoppers above the threshold value when they land on the ground.As it is cleared, the simulation process installation is not achieved by attention to (62) because of inaccessibility to the suitable exploring and exploitation the search space around the solution.Therefore, this equation can be utilized just in the free space for the swarm grasshoppers.The interconnection between grasshoppers in swarm can be calculated by means of (62).
The obtained experiences indicate that the speed of the grasshopper is so high for accessibility to the comfort zone and then the convergence process is not enhanced.For this reason, usually the mathematical model can't be installed in the optimization problem directly.To overcome this problem, a facilitated version for (62) is presented as follows: In (63), the upper and lower band are defined by ub d and lb d in Ddimension, respectively.The value of the Ddimension for the best obtained solution so far, is introduced by Td .Also c is a parameter that should be decreased because of shrinking three zones (comfort, repulsion and absorption).
As it is seemed in (63), the upgrade process of the position of grasshopper will be achieved according to the current position, the position of target and the all other grasshoppers position.It can be concluded that the situation of all grasshoppers is noticed in finding the location of search agents around the target.
The parameter c which is in the inner part of (63), actually is utilized as a setting factor for tuning of repulsion/absorption forces between grasshoppers.As it has seemed, the inner c is corresponding to the number of iterations.Also the parameter c which is in the outer part of (63), effectively is useful as an important factor for tuning of the search coverage around the target when as the iteration is increasing.

E. Rakhshani et al.
Therefore, it is cleared that the first part in equation ( 63), which is related to summation term, will be focused on the position of other grasshoppers and finally performs the grasshopper interconnection in nature.
The other main part in (63), Td , points to the variation rate of the tendency for grasshoppers to fly towards the food source.Also as it was explained in above, the parameter c can simulate the decreasing process of grasshopper which is proceeded towards the food source.All existent parts in (63), Td and the summation part, are multiplied by the random values which create more random behavior for grasshoppers.
The obtained formulation in (63), by using the explained two terms can explore and exploit the search space.The required structure in the level tuning of exploration and exploitation for search agents is provided as follows: As it is observed in (64), the maximum and minimum values for c are introduced by c max and c min respectively.Also the current of iteration is shown using l and L is the biggest number of iterations.By attention to (64), it is indicated that the number of iteration and the parameter c are related together in order to make balance between the exploration and exploitation process.
By considering to above explanations, it is obvious that the proposed model can enhance gradually grasshopper to fly towards a task.In grasshopper optimization algorithm, the grasshopper with best objective value is assumed as the fittest grasshopper during optimization process.This rule can cause to keep the best solution in every iteration for GOA.

Grey wolf optimization algorithm (GWO)
Comprehensive analysis about the taken processes from nature, finally focused the scientists and researchers to the swarmbased life of grey wolf.These studies lead to appear a well -regarded optimization algorithm that is called GWO.According to the achieved researches, the social behavior of grey wolves has a hierarchical structure [52].The main part of swarmbased algorithm in GWO imitates by the social interconnection of grey wolves in nature.Almost a crowd of wolves is classified into specified types of members according to the command level including α, β, δ and ɷ.The most presiding wolf is α and this ability will be reduced from α to ɷ.The hunt process in wolves usually performs as a crowd.In the other hand, wolves cooperate in a smart way to hunt a prey (victim).At first, the hunt process is started by pursuing the prey by Grey wolves as a team.Then they attempt to surround it circularly.This strategy will be continued in order to increase the chance of hunt.A conceptual diagram from this process is indicated in Fig. 7.
Above diagram insists that various positions can be enhanced by attention to the location of the grey wolf (x , y) and the situations of the prey (x* , y*).This idea can be extended for more than the 2th dimension.It can be possible to create a mathematical model as follows: In ( 65) and (66), X → (iter + 1) emphasizes on the new position of the wolf.Also X → p (iter), by subscript p, indicates to the current location of prey (victim).

According to above equations, For calculation of A
→ and C → , as specified coefficients, the equations are presented by (67) and (68).It is obvious that the distance which relates to the location of the prey (victim), will be shown by D → .
where a → is a tuning parameter that reduces linearly from 2 to 0 during the process.Also two random parameters, r → 1 and r → 2 , are produced from the interval 0 to 1.
In the optimization process, the earliest three best solutions are selected and kept.Then the update process of their locations will be continued based on the best search agents.Finding a mathematical model for the hierarchical behavior formulation of wolves in enhancement of GWO algorithm, will be realized by classification their population into: α, β, δ and ɷ.Fig. 8 indicates the upgrade process of positions for α, β and δ, respectively, in a two-dimensional space.
The other nominee for solutions are focused on ɷ.Actually in GWO, the optimization process is managed by α, β, δ and ɷ.So by keeping the first three best obtained solutions so far and enforcing the other search agents (including the ɷ) to upgrade their positions based on the location of the best search agent.Also in this process, the main below equations are presented.
In fact, this swarmbased algorithm will be completed using ( 69) -(75) and therefore the upgrade process of position is achieved based on to α, β and δ in a dimensional search space.In (67) and (68), A → and C → are enforced this process to explore and exploit the search space.
Reducing A → will be lead that half of the iterations is devoted to exploration ( A → ≥ 1) and the other half is dedicated to exploitation

Simulation results
This section is presented with two main parts.The first part is focused on preparation and cost function determination.The last part will focus on the desing of the controller, test and simulation.To evaluate the performance of SBOAs in finding the best coefficients of multivariable PI (MPI) for VSP/HVDC model, at first, it is necessary to propose a new cost function according to desired MIMO characteristics.The applied minimization process will be achieved using SBOAs for the proposed cost function and it will be compared with the typical tuned SMPI based on Section.3. Also for achieving a comprehensive comparison between different models, the model with AC/HVDC system is also presented in an individual section and the proposed SMPI controller is performed and discussed.Simulation time range is equal to 50 s.

Proposed cost function
In this section, it is necessary to select a specified approach among the optimization process.By attention to the determined features of step response and its accessibility in every iteration of the simulation, a fitness index is used according to the extracted features of the step function.So the calculation process of the cost function is performed according to the obtained MIMO transfer function in this part.
It should be noted that the proposed fitness index actually is a function, including summation of the selected step response features, which they are formulated as a linear combination by weighted coefficients.The utilized features, which are mentioned in this work, are  The movement manner of closedloop poles as an effective parameter should be mentioned in simulations.As it is cleared, the pole position can be affected by the system's stability.In this study, the most attention should be given to adding an equal index for managing stability during the optimization process.To determine this component, at first it is important to find the real value of poles.After that, we can select the most positive value between the obtained reals while determining their maximum.It is noted that when these poles locate on the imaginary axes or in the right hand of s (unstable poles),  the calculated parameter will be zero or bigger than zero.To overcome this problem, we can write: By attention to (76), the maximum real part of closed loop poles will be negative.Finally, the SI parameter can be defined as follows.It is obvious that the lower SI the better: In (76), p 1 , p 2 ,…,p i are the obtained closed-loop poles.

Settling Time Index (ST):
One of the most specified features for step response in a system is the minimum time that the amplitude gets to correct error of 0.05.On the other hand, according to (20), VSP/HVDC model, as the basic system, has MIMO transfer function.Then by implementation of step unit function as the input signal, the obtained step response will be calculated as a 2 × 2 matrix.Therefore, in this work we will confront four ST parameters which every one of them belongs to one of the element matrices of the step response.On the other hand, it is clear that for MIMO step response: Eventually, by considering the worst possible situations, and the obtained elements of the matrix (78), this index can be determined by:

Rise Time Index (RT):
The other effective feature of step response, can be determined as the time it takes for the response to rise from 10% to 90% of the steady-state response.Similarly, there is a matrix for the determined RT as follows: and finally this index will be calculated as: By notification to (77), ( 79) and (81), the cost function can be defined as a linear combination of above indexes: where k 1 , k 2 and k 3 are constant and real positive digits which are selected by attention to the target design.By definition of their values, we can change the role of every above index in the minimization of a cost function.It is notified that in this case, all k i are selected equal to 1.
As it is indicated in (82), the fitness index is produced by three homogenous terms with the same dimension.This important point is applied by selecting the manner of three indexes, including SI, ST and RT.
Therefore, SBOAs will be completed by implementing this cost function during the optimization procedure.Fig. 9 illustrates a conceptual process from the optimization method for the proposed cost function.This task is achieved using swarmbased strategies.
As shown in Fig. 9, the proposed optimization algorithm can be performed as follows: a) Set up Initial parameters for swarmbased optimization methods, consist of: • Number of decision variables (nVar = 4).
• Population size (nPop = 80).b) Selection of decision variables randomly and replacing in MPI structure as follows:    (82) in (e), before achieving to the optimal value, is fixed on a specified value and then this situation is continuously repeated more than 15 iterations.• Stability limitations: It is important to stay all the obtained poles in the closedloop transfer function on the left hand of S pane in all iterations.
• Possibility of realization: The industrial implementation of MPI controller will be related to the obtained actuator signal as U according to ( 22) and ( 34).
• Final Decision about above limitations: If every set of the selected Ki is lead to one of above limitations, that set will be deleted from the collection of solutions.g) New Iteration = past Iteration + 1 h) If (new iteration < max iteration) then update parameters and go back (b) i) Show the last value for Z and MPI coefficients as follow: j) End."

Preparation process
By attention to Section 2, input signals for MIMO VSP/HVDC model are considered as (17) and Fig. 10 shows the variation of the input variables.
Namely, it is assumed that contingencies happened as step load changes in the first area by rising to 0.03p.u. at t = 2 sec, and then the load value deviation becomes 0.015p.u. at t = 30 sec.Then for simulation of these step changes, r(t) can be defined as input signals: For high availability to all parameters and the obtained separated results in simulations, it is better to program as an extendable structure.As it was described in PSO algorithm, every particle can carry out structures, consisting of the position, objective value (cost value), velocity and the best structure.In the same way, the best structure will be included position and objective value (cost value) part too.This separated structure will be achieved on others SBOAs, such as GOA and GWO.
Table 1 shows common tuning parameters for all SBOAs.Additional parameters for different numbers of iteration with different algorithms are presented in Table 2 as well.
In this scenario, all parameters have to be designed and installed in a way that the gains of K i (i = 1,2,3,4) in the below MPI controller can be calculated:

First Scenario: Simulation of the designed SMPI for VSP/HVDC model
In this section, a VSP/HVDC model, as it was defined in Section 2, will be considered as the base test system to evaluate the achievement SMPI controller compared to a typical setting method.It is mentioned that two load step changes are happening at 3 sec and 30 s as shown in Fig. 10.The performance of the developed model is simulated by developing both MATLAB scripts and SIMULINK R2016b® environment which is a programming platform for numerical computation and data.For implementing SMPI with different SBOAs, PSO, GOA and GWO are designed and installed to calculate the gains of K i (i = 1,2,3,4).
According to Fig. 9 and Section 5.1, basically, a MIMO matrix will be determined as the step response.Then the obtained features, including ST, SI and RT, are extracted and entered in the optimization algorithm.
Fig. 11 shows the variation of cost function value versus the number of iterations.In this figure, the cost function of deviations is illustrated based on SBOAs.According to Fig. 11, GWO scenario presents faster convergence by less iteration.Table .2 can provide the possibility of comparison between different swarmbased intelligence methods and their convergence speed in the same situations.As it is observed, this comparison is applied in various population size.By considering the recorded values, the GWO and GOA characteristics are so closed and more suitable than PSO.Also the obtained PSO features are slower and farther than other two SBO algorithms.By comparison between the final cost function, it is obvious that the optimization process in GWO and GOA has achieved more accurate than PSO.So the obtained results in Table .2, can satisfy the cost function deviation in Fig. 11.
For having a better comparison, the tuning process of tuning SMPI coefficients is applied with four SBOAs, consist of a typical tuned method, tuned using PSO, GOA and GWO strategy.For the first and the second state variables in VSP/HVDC model, Fig. 12 (a) and (b) exhibit the frequency deviations in area 1 and 2 (Δɷ 1 , Δɷ 2 ) using typical tuned MPI controller method and the proposed SBOAs.
Based on the simulation results, it is clear that the conventional tuning method, presented in Section 3, has the least effect on output signal.Just the reverse of usual tuning method, Fig. 12 shows that the applied swarm-based tuning method can represent a more accurate performance from the old setting strategy.This fact is more cleared during the optimization process by means of GWO and GOA.
Fig. 13 discloses tieline AC power between areas (ΔP tieAC12 ) deviations in VSP -HVDC model.It is absolutely clear that ΔP tieAC12 is controlled by the acceptable damping using GWO versus other SBOAs basis on the simulation results.
Also in this case, the conventional MPI controller has more suitable performance than PSO.Convinces the precise performance of estimation stage too.
Fig. 14 (a) and (b) indicate the injected DC power deviation and its derivative in VSP 1 for VSP/HVDC model (Δx p,vsp1 and Δx f,vsp1 respectively).In both (a) and (b), the SMPI controller design process has more suitable performance using GWO.Also the typical MPI tuning method has less direct effect on system without any controller.
It is clear that GOA has more closed results to GWO method.This conclusion is seen in all simulation results.Fig. 15 (a) and (b) indicate the injected DC power deviation and its derivative in VSP 2 for VSP/HVDC model (Δx p,vsp2 and Δx f,vsp2 respectively).In both, the SMPI controller design process has more suitable performance using GWO.The obtained values for the SMPI controller gains (K P and K I matrixes) with different SBOAs, consist of PSO, GOA and GWO, are given in Table .3. Also the final performance index values for SBOAs are described in Table .3.

Sensitivity analysis for VSP/HVDC case
As it is obvious, the stability condition for a closedloop structure is described by this way that if and only if the trajectory of the Nyquist diagram of G(sjɷ)H (jɷ) from -∞ < ɷ < ∞ encircles the point (-1,0) in a counter-clockwise direction as much times as G(s)H (s) has unstable poles.
Above explanation can be extended by the generalized Nyquist's stability criterion for MIMO systems, which is a modification from the Nyquist's stability criterion for SISO systems [53][54].
By attention to Section II and (20), the generalized Nyquist's stability criterion for SBOAs, including PSO, GOA and GWO can be achieved by the below figures.From the obtained graphs in Fig. 16, it is observed that GWO and GOA methods are more stable than PSO tuning method.In this figure, the conventional MPI controller design is shown by the 'No Control' legend.

Robustness analysis for VSP/HVDC case
In this section, it is possible to analysis the behavior of the tuned SMPI controller using GWO method, as the most optimal SBOAs on the VSP/HVDC model, against parameter changes.By considering a common setup parameter for both VSPs and a direct relation between selected parameter and tieline AC power deviation among area 1 and 2, we can choose k acVPS as the most suitable parameter.In the other hand, we will able to study robustness analysis on VSP/HVDC model, as the structured uncertainty.Therefore, we present the reaction of Δɷ 1 , Δɷ 2 and ΔP tieAC12 after applying this robust evaluation strategy.As it seems, by decreasing k acVSP from the standard value, frequency deviations in both areas start to decrease step by step.This change is shown in Fig. 17 .Therefore, as it is cleared, despite the implementation of k acVSP variations in VSP/HVDC model, the proposed SMPI controller stay in the stability range.

Second Scenario: Design of SMPI on AC/HVDC model
As it was described in the previous section, VSP/HVDC model, will be obtained by means of utilization a higher level control technique for emulation using VSP.In the other hand, this control structure performs by programming the electrical performance of a synchronous generator (SG) in a HVDC line which is responsible for controlling the VSCconverter.Therefore, every tuning on this model, such as tuning SMPI controller basis on this paper, as separately can impress the performance of HVDC model.So the second scenario will be centralized on the SMPI controller effects for HVDC model.
By attention to this model, the converter stations of the HVDC links in a multi -area AC/DC interconnected system, as a dynamic model of the facilitated MIMO system by means of following equations: Finally, by considering described equations, MIMO transfer function for HVDC model can be driven as follows: The MIMO transfer function in (92) for the AC/HVDC model shown in equation (93) considering its state space matrices is presented in Appendix A.
As it was mentioned for VSP/HVDC model, by considering to the same initial parameters for swarmbased optimization algorithms and the same load step change as input deviations, simulation results for HVDC model will be presented as follow: By attention to Table .4 and Fig. 18, it is observed that the cost function value for the tuning process using PSO method has performed so closed values to the output results using GWO method.Also in HVDC model, for the initial iterations, GOA method performs slowly but in the other iterations this method can compensate its backwardness suitably and tends to the same permanent value.
By attention to the obtained values for PI coefficients, (K i , i = 1,2,3,4), it seems that tuned SMPI controller in HVDC model tends to behave similar to a MP controller.
By considering Figs.19 and 20, it is clear that the typical tuned MPI controller is slower than SMPI controllers (using PSO, GOA and GWO).Also it seems that the obtained results from the tuned SMPI controllers, using PSO and GWO, are more suitable than the proposed controller tuning by GOA.
The obtained simulation results from HVDC model can give the comparison opportunity to VSP/HVDC model.This comparison can perform by attention to the following tables for every swarmbased optimization algorithm.As it is observed, frequency deviation in both twoareas and tie line AC power deviation between areas, as the same state variables, are selected from both models and then they are compared in 10 sample time, with 5 s for every range time.These tables can represent the comparison ability for finding the leveling effect of utilization swarmbased optimization algorithms in the tuning process of MPI controller.
By attention to Table .5, 6 and 7, it is observed that the obtained simulation results in the SMPI tuning process from VSP/HVDC model are closer to zero value and then all swarmbased optimization algorithms, PSO and GOA and GWO, have more suitable values for VSP/HVDC model.

Sensitivity analysis for AC/HVDC model
As it was illustrated in the first scenario.In this section and about HVDC model, we can utilize the generalized Nyquist's stability criterion for the second scenario's model.
Therefore, the sensitivity analysis frequency for SBOAs, including PSO, will be obtained by the below figures.Fig. 21 presents Nyquist stability for all four transfer functions elements in G HVDC Total from (93).From this figure, it is obvious that the tuning methods, including PSO and GWO, perform so closed and they can create a suitable performance for SMPI controller.Absolutely, the GOA method creates more sensitive Nyquist diagram which is related to be weaker its performance in comparison with two other SBOAs.

Conclusion
A new high-level control application for a heuristic-based multivariable PI (MPI) controller design for an interconnected dynamic AGC model with VSP is presented and discussed.
Swarm-based optimization algorithms (SBOAs), as the most applied, naturalbased and flexible metaheuristic methods, is utilized to achieve the optimized solution where exact optimization techniques fail to obtain the satisfying results.In this work, it is tried to make a new challenge for comparative analysis SBOAs and their performance to adjust a Multivariable PI (MPI) controller, which is called Swarm -Based Multivariable PI (SMPI) controller.By considering more experienced background in the optimization problems, recently applied researches and well-known structure in scientific societies, particle swarm optimization (PSO), grasshopper algorithm (GOA) and also grey wolf optimization (GWO) are mentioned, among all SBOAs, to implement in this tuning MPI process.These selected SBOAs are setup and tested on a multiinput multi output (MIMO) AGC model with high voltage direct current (HVDC) with the virtual synchronous power (VSP) strategy, which is called VSP/HVDC model.The proposed method is based on the minimization of a performance index, which is designed by considering the transfer function characteristics matrix in a MIMO system.The SMPI controller design method is performed for the systems to improve characteristics performance of the model versus the conventional MPI controller method.For achieving to the best comparison between methods: 1-A conventional tuning method is selected to design MPI gains, consist of K P and K I matrix, in VSP/HVDC model.2-A performance index is selected based on the step response information, especially important characteristics such as settling time (ST), stability index (SI) and rise time (RT) are considered to minimize using SBOAs.3-Set-up parameters for SBOAs, consist of PSO, GOA and GWO, are tried to be close together for more accurate comparison between methods.4-The obtained results exhibit that the proposed SBOAs, especially GWO and GOA, have more suitable performance than the conventional tuning method.
It should be noted that in the presented study PSO has been selected the base swarm-based optimization algorithm for comparing with other methods.GOA is selected at one of the newest intelligent algorithm among the other.Respect with other existing methods, application of other swarm based algorithms, such as monarch butterfly optimization (MBO), earthworm optimization algorithm (EWA), elephant herding optimization (EHO), can be considered as part of our future research.
Performed studies show that the proposed controller has few parameters to adjust and is relatedly simple to implement.Its robust and can be used in parallel computation.The speed of its convergence for the studies cases are fast with short computational time.But despite of the advantages, it is important to take into account the importance of initial designed parameters.Improper values for the initial parameters might lead to premature convergence reaching to a local minimum especially with very big system with high complexities.The industrial implementation of MPI controller will be related to the obtained actuator signal as U, shown in equation (22), and during the design process it is important to stay all the obtained poles in the closedloop transfer function on the left hand of S-plan in all iterations.

Fig. 1 .
Fig. 1.Basic frame of AGC in multi-area systems with a VSP based AC/DC transmission.

Fig. 4 .
Fig. 4. Search space management in algorithms based on exploration and exploitation.

Fig. 7 .Fig. 8 .
Fig. 7.The position of a wolf can be defined around a prey with considering a 2D mesh around the prey.

Fig. 13 .
Fig. 13.Variation of tie line AC power between areas.
2 and the vector C → also modified exploration when C → > 1.Also for C → < 1, exploitation will be highlighted.According to α, β and δ in a dimensional search space.In these formulas, A → and C → are obliging the GWO algorithm to explore and exploit the search space.With decreasing A → , half of the iterations is devoted to exploration ( A → ≥ 1) and the other half is dedicated to exploitation ( A → < 1).The range of C → is 0 ≤ C → ≤ 2 and the vector C → also improves exploration when C → > 1. Exploitation is emphasized when C → < 1.
Some of the main assumptions which are utilized in the simulation process of the proposed swarmbased optimization MPI (SMPI) controller are as follow: a) To set up initial parameters for swarmbased optimization methods: Number of decision variables (nVar = 4).Lower and upper bond of decision variables (VarMin = -1 and Var-Max = 1).Maximum number of iterations (MaxIt = 300).Population size (nPop = 80).b) To set up VSP/HVDC model: A VSP , B VSP and C VSP as the constant state, input and output matrices.Implementation of input reference signal at r(t).

E
.Rakhshani et al.
Calculate the closed loop transfer function by the obtained decision variables from (b) and equations (23) -(32).d) Extract the proposed cost function terms, for achieving to SI, ST and RT, by means of the obtained result from (c).

Fig. 20 .
Fig. 20.Variation of Tie line AC power between areas (p.u.MW) for HVDC model.
the final cost function by replacing the obtained results from (d) into equation (82).f) Consider the most important limitations in this process, as it follows: • Proper Local Minimizations: If the obtained cost function, according to

MPI Controller Multivariable PI Controller PI Controller Proportional -Integral Controller PSO Particle Swarm Optimization SBOAs Swarm -Based Optimization Algorithms SI Stability Index SMPI Controller Swarm -Based Multivariable Proportional -
which are applied in the tuning process of MPI

Table 1
Setting The Same Parameters of SBOAs for VSP/HVDC Model.

Table 2
Comparative transient performance for VSP/HVDC system, tuned by SMPI.

Table 3
The best performance index for SBOAs.

Table 4
Comparative transient performance for HVDC model, tuned by SMPI.

Table 5
Comparison Analysis Between State Variables in The Tuning SMPI Process Using PSO.

Table 6
Comparison Analysis Between State Variables in The Tuning SMPI Process Using GOA.

Table 7
Comparison Analysis between State Variables in The Tuning SMPI Process Using GWO.