Comparison of two metaheuristics with mathematical programming methods for the solution of OPF

Comparison of two metaheuristics with mathematical programming methods for the solution of OPFP.N.Biskas,N.P.Ziogos,A.Tellidou,C.E.Zoumas,A.G.Bakirtzis and V.PetridisAbstract:Different optimisation methods developed for the solution of the nonlinear OPF problem with both continuous and discrete variables are compared.Two mathematical programming methods are compared with two metaheuristics,an enhanced genetic algorithm and a particle swarm optimisation implementation.Test results from the application of the methods to several IEEE systems are presented and eful conclusions are drawn concerning the execution times and the ‘optimum’costs provided by all four tested methods.List of symbols P G unit active power output vector Q G unit reactive power output vector V bus voltage magnitude vector t transformer tap settings vector b SH bus shunt admittance vector u transformer phase shifts vector hbus voltage phase angle vector1IntroductionSince its introduction in the 1960s optimal power flow (OPF)has been the subject of intensive research due to its ability of integrating the economic and security aspects of a power system into one mathematical formulation.In its most general formulation,OPF is a nonlinear,nonconvex,large-scale static optimisation problem,with both contin-uous and discrete control variables.Numerous optimisation methods have been proposed for the solution of the non-linear OPF problem,such as mathematical programming methods,stochastic global optimisation approaches and metaheuristics.The mathematical programming methods can be categorised as gradient techniques such as the generalised reduced gradient technique,successive quadratic programming (SQP),Lagrange Newton approaches,suc-cessive linear programming (SLP)and interior-point (IP)methods.Metaheuristics include simulated annealing type methods,genetic algorithms and evolutionary programming methods.The gradient techniques [1,2]were the first approaches proposed to solve OPF problems.These approaches exhibit slow convergence,making small steps near the optimum.The SQP approaches [3,4]use the second-order derivatives to improve the rate of convergence of thegradient approaches.Their modelling is based on the quasi-Newton process,in which the approximation of the hessian matrix of the Lagrangian function is used to overcome the difficulties encountered in QP problems.However,as in all quasi-Newton methods,the reduced hessian built iteratively is dense which may make these methods too slow as the number of control variables increases.Instead of solving the original problem,Newton methods [5,6]solve the problem’s Karush–Kuhn–Tucker (KKT)optimality conditions.In [5]the Newton method (for unconstrained optimisation)is combined with a Lagrange multiplier method (for optimisation with equalities),while the inequality constraints are added as quadratic penalty terms to the problem objective,multiplied by appropriate penalty multipliers.The major difficulty of Newton methods turned out to be the efficient identification of binding inequalities,an issue later studied in [6].The Newton method is favoured for its quadratic convergence properties,but its solution process remains time-consuming.The SLP approaches [7,8]are based on the linearisation of OPF constraints and the objective function.The central modelling issue is the linear representation of network equations.An incremental linear model is adopted,and the simplex algorithm (or a variant of it)is used.Due to the fact that LP approaches provide solutions jumping from one vertex to another,during the solution process a successive refinement of cost curves into smaller segments is needed to overcome the discontinuities.The execution time required for the solution obtained by SLP methods increases substantially with the system size and the number of constraints.IP methods have been attracting interests in the research community since their introduction and several IP variants have been developed during the last two decades for the solution of the nonlinear OPF problem [9–11].The theoretical foundation behind IP methods consists of three basic building blocks:Newton’s method for solving nonlinear equations (for unconstrained optimisation),Lagrange’s method for optimisation with equalities,and a barrier method for optimisation with inequalities.Among the many variants,the predictor–corrector primal–dual IP method has been computationally the most efficient for NLP problems.IP methods initially convert the inequality constraints to equalities by the introduction of nonnegative slack variables.A logarithmic barrier function of the slackE-mail:pbiskas@egnatia.ee.auth.grThe authors are with the Department of Electrical &Computer Engineering,Aristotle University of Thessaloniki,54124Thessaloniki,Greece r IEE,2006IEE Proceedings online no.20050047doi:10.1049/ip-gtd:20050047Paper first received 7th February 2005and in final revised form 12th September 2005variables is then added to the objective function,multiplied by a barrier parameter.The barrier parameter is analogous to the problem duality gap,therefore it is expressed as a function that is gradually reduced to zero during the solution process.The KKT conditions of the modified problem are formulated as a set of nonlinear equations that can be solved by the Newton’s method through an iterative process.One attractive feature of IP methods is their ability to handle nonlinear inequalities without requiring an active set identification,as in Newton methods.Additionally,a strictly feasible starting point is not mandatory;only the strict positive conditions on a subset of variables must be satisfied at the initial point and subsequent iterations. OPF programs based on mathematical programming methods have been designed for purely continuous-variable OPF.However,OPF is a mixed-integer NLP problem with discrete control variables,such as transformer tap settings and phase shifts and bus shunt admittances.The presence of discrete control variables and the inclusion of the nonlinear power-flow equality constraints make OPF a nonconvex problem.To avoid the prohibitive computational require-ments of mixed-integer programming,discrete control variables are initially treated as continuous,and postproces-sing discretisation logic is subsequently applied[12,13]. Whereas the effects of discretisation on load tap-changing transformers are small and usually negligible[14],the rounding of switchable shunt devices may lead to voltage infeasibility,especially when the discrete VAr steps are large,and requires special logic[13].In recent attempts to overcome the limitations of the mathematical programming methods,metaheuristic methods have been developed that treat the discrete control variables in their exact form and take the feasibility of OPF into account.Metaheuristics include SA-type methods[15], GAs and evolutionary computation techniques,such as the recently presented method called particle swarm optimisa-tion(PSO).Metaheuristics are methods that are based on processes observed in physics,biology or sociology.They are stochastic search techniques that use the mechanics of evolution to produce optimal solutions to a given problem. They are classified as heuristics,since,with the exception of SA,no formal proof of their convergence to the global optimum exists.They are particularly well-suited to nonmonotonic solution surfaces,where many local optima may exist.One basic disadvantage of these methods is that they do not compute nodal prices,which are essential for pricing energy under the deregulated environment of the electric energy sector.GA approaches[16,17]overcome the limitations of the mathematical programming approaches in the modelling of nonconvex cost functions,discrete control variables and prohibited unit operating zones,showing superb modelling flexibility.The main disadvantage of GAs is that they are stochastic algorithms and the solution they provide to the OPF problem is not guaranteed to be the optimum.The PSO technique is based on social behaviour of natural creatures,such as birds orfish.It combines social psychology principles in socio-cognition human agents and evolutionary computations.The literature on PSO approaches for the solution of the OPF problem is limited [18–20].The PSO method has the modellingflexibility of GAs and is generally simple in concept,easy to implement and computationally efficient.PSO has theflexibility to control the balance between the global and local exploration of the search space,which enhances its search capabilities avoiding premature convergence.Unlike mathematical programming methods(and like GAs)the solution quality of PSO does not rely on the starting point in the search space.However,it inherits the disadvantage of GAs in not guaranteeing the global optimum.In this paper a comparison of the solutions provided by mathematical programming methods with the corres-ponding solutions of metaheuristics is presented.For the case studies of mathematical programming methods,two commercial NLP packages,GAMS 2.50(MINOS)[21] and LINDO API2.0(CONOPT)[22],were employed. Metaheuristics comprise an enhanced genetic algorithm[17] and a PSO implementation.Several test systems,created by replicating the IEEE one-area RTS96and the IEEE 118-bus system,are used for the comparison.The execution times and the optimum costs provided by all four tested methods are reported.Our goal is to provide results on the scalability of the two metaheuristics,by stressing both to their limits.2OPF problem formulationThe OPF problem can be formulated as a mathematical optimisation problem as follows:min fðzÞð1Þsubject to gðzÞ¼0ð2ÞhðzÞ0ð3Þz min z z maxð4Þwhere z¼b P TGQ T G V T h T t T u T b T SH c T is the vector of unknown variables,function f in(1)represents the power system’s operation optimisation goal,equality constraints (2)are the nonlinear power-flow equations,inequality constraints(3)are nonlinear functional operating con-straints that correspond to branchflow limits(MVA),and constraints(4)define the feasibility region of the problem variables,such as:unit active and reactive power output limitsbus voltage magnitude limitsreference bus voltage phase angle limits(set to0.0)transformer tap setting and phase shift limits(discrete values)andbus shunt admittance limits(continuous or discrete control).In our case,the objective function f(z)is the power system’s operating cost.Bus shunt admittances,transformer tap settings and phase shifts are considered here as discrete variables.The handling of discrete variables is described in the following Section.3Mathematical programming methods3.1Solvers3.1.1MINOS under GAMS:For problems with nonlinear objective function and nonlinear constraints, a projected augmented Lagrangian algorithm is used[23]. MINOS solves a sequence of subproblems,in which the objective is an augmented Lagrangian function,and the constraints are the true linear constraints plus linearisations of the nonlinear constraints.Experimental runs using SNOPT and CONOPT under GAMS have shown that these solvers provide similar results with MINOS,concerning execution time and number of solver iterations required for convergence to the optimal solution.Therefore only MINOS(under GAMS)solutions are presented in this paper.3.1.2CONOPT under LINDO API :LINDOAPI [22]uses the CONOPT solver for NLP problems.CONOPT is a feasible-path solver based on the generalised reduced gradient method.CONOPT is based on the usual NLP model in which all variables are continuous and all constraints are smooth with smooth first derivatives.In addition,the Jacobian matrix is stored and processed as a sparse matrix.CONOPT attempts to find a local opti-mum satisfying the usual Karush–Kuhn–Tucker optimality conditions.3.2Implementation details3.2.1Constraint matrix representation :Inour LINDO API implementation the sparse matrix representation of the constraint matrix of the problem is provided by the user to the solver,whereas in GAMS it is built by the modelling interface.Therefore the LINDO API implementation avoids the additional modelling time of GAMS,as discussed in Section 5.3.3.2.2Jacobian matrix computation :In GAMSthe Jacobian matrix is computed by the solver (using symbolic differentiation),whereas in the LINDO API implementation it is provided by the user to the solver.Experimental runs have shown that the speed-up attained by the computation of the Jacobian matrix in our LINDO API implementation reaches 300%in the total execution time.3.2.3Discrete variables handling :The pro-blem discrete variables are handled in both implementations as follows.The continuous OPF problem is initially solved (all variables are considered as being continuous).A simple postprocessing discretisation logic is subsequently applied,in which the values of bus shunt admittances,transformer tap settings and phase shifts are rounded to their closest discrete value.The OPF problem is solved again with the values of the variables fixed to the previously computed discrete values.The solution provided by this process is suboptimal but very close to the solution provided by a MINLP solver,as shown in Section 5.2.4Metaheuristics4.1OPF problem formulationA slightly different OPF formulation is used for the metaheuristic methods by separating the problem decision variables z to state variables x and control variables u as follows:min f ðx ;u Þð5Þsubject to g ðx ;u Þ¼0ð6Þh ðx ;u Þ 0ð7Þu 2Uð8Þwherex ¼h_TV T Ljk Tsystem state vector ð9Þu ¼P_TGV T Gt T u Tb T SHjk Tsystem control vectorð10Þand V L is the load (PQ )bus voltage magnitude vector,while V G is the generation (PV )bus voltage magnitude vector.A ‘hat’above vectors h and P G denotes that the entry corresponding to the slack bus is missing.The equality constraints (6)are the nonlinear power-flow equations.The inequality constraints (7)are the functional operating constraints,such as branch flow limits (MVA,MW or A),load bus voltage magnitude limits,generator reactive capabilities and slack bus active power output limits.Constraints (8)define the feasibility region of the problem control variables,such as unit active power output limits,generation bus voltage magnitude limits,transformer tap setting and phase-shift limits and bus shunt admittance limits.4.2Enhanced genetic algorithmThere are three major differences between GAs and mathematical programming methods.First,GAs operate on the encoded string of the problem parameters rather than the actual parameters of the problem.Each string can be thought of as a chromosome that completely describes one candidate solution to the problem.Secondly,GAs use a population of points rather than a single point in their search.This allows the GA to explore several areas of the search space simultaneously,reducing the probability of finding local optima.Thirdly,GAs do not require any prior knowledge,space limitations or special properties of the function to be optimised such as smoothness,convexity,unimodality or existence of derivatives.They only require the evaluation of the so-called fitness function (FF )to assign a quality value to every solution produced.In this Section the enhanced genetic algorithm (EGA)of [17]is implemented.Each chromosome is formed using 12bits for each unit active power output,ten bits for each bus voltage magnitude,five bits for each transformer tap setting and two or three bits for each bus shunt admittance.The EGA population is consisted of 100chromosomes.All the members of the population are initially generated at random,except from one,which is initialised using the corresponding genotypes of the upper bound of the control variables (a new feature here as compared with [17]).Simple genetic evolution takes place by means of three basic genetic operators:Parent Selection :The selection rule used in the EGA is a simple roulette-wheel selection.Crossover :Uniform crossover has been used in the EGA,and is applied with initial probability of 0.9.Mutation :The mutation operator is applied to every bit of one chromosome with initial probability of 0.001.A set of advanced and problem-specific genetic operators [17]has been added to the simple genetic evolution to increase its convergence speed and improve the quality of solutions.The advanced features included in our EGA implementation are fitness scaling (a linear transformation here),elitism,multiparameter hill-climbing (number of parameters is set to three)and elite self-fertilisation.The problem-specific operators,which are gene swap operator,gene cross-swap operator,gene copy operator,gene inverse operator and gene max–min operator,are applied to all chromosomes with a probability of 0.3.The OPF functional operating constraints (7)are included in the GA solution by augmenting the GA FF by appropriate penalty terms for each violated functional constraint.Constraints on the control variables (8)are auto-matically satisfied by the selected GA encoding/decodingscheme.Therefore the GA FF is formed as follows:FF ¼AP N G i ¼1F i P Gi ðÞþB P N C i ¼1w j Pen jð11ÞPen j ¼j h j x ;u ðÞj ÁH h j x ;u ðÞÀÁð12Þwhere FF fitness functionA ,B positive constant numbersF i (P Gi )fuel cost function of unit i ,in our case it is a quadratic functionw j weighting factor of functional operating constraint jPen j penalty function for functional operating constraint jh j (x ,u )violation of j th functional operating constraint,if positiveH (Á)the Heaviside (unit step)function N G number of unitsN Cnumber of functional operating constraintsThe weighting factors w j of the penalties that correspond tothe violation of the inequality constraints (7)are adjusted appropriately in each case study to avoid violations of constraints (7)in the final solution.Given a candidate solution to the problem,represented by a chromosome,the FF is computed as follows:(i)Decode the chromosome to determine the values of the control variables u .The computed control vector satisfies,by design,constraints (8).(ii)Solve the power-flow equations (6)to compute the state vector x .(iii)Determine the violated functional constraints (7)and compute associated penalty functions (12).(iv)Compute the FF using (11).In step (ii)a simple fast decoupled load flow (FDLF)is used with no PV–PQ bus-type switching since generator reactive capabilities are incorporated in the functional operating constraints,and no local control adjustments,such as tap and switchable shunts,since the settings of these controls are determined by the GA.Therefore only a few load flow iterations are required for convergence.The FDLF B 0and B 00matrices are formed and factorised only once in the beginning;the effect of the changes of shunt admittances on the B 00matrix is neglected.If the OPF is used for the calculation of maximum power transfers or maximum system loadability limits,where the active and reactive power are strongly coupled near the optimum and the decoupled load flow may experience convergence problems,a Newton load-flow solution (with no PV–PQ bus-type switching and no local control adjustments)can be used at the expense of increased computational requirements.A predefined number of generations is used as a termination criterion for the EGA.4.3Particle swarm optimisationIn this Section the standard PSO (SPSO)technique is initially presented.A modification of the SPSO (known as dissipative PSO,DPSO)is subsequently analysed,which has been proved to perform better than the SPSO technique [24].4.3.1Standard PSO :The PSO technique employsan iterative process which uses a population of particles to model the social behaviour of individuals of a swarm.Each particle is a vector containing the control vari-ables,suggesting a possible solution to the OPF problem.The location of the i th particle is represented asu i ¼u i 1;u i 2;...;u iD ðÞ,where u id 2b u min d ;u maxd c ,d 2½1;D ;u min d and u maxdare the lower and upper bounds of variable d ,respectively,and D is the number of control variables.The particles are moving in a D -dimensional space.Each particle attempts to minimise the following objective function:Obj ¼X N G i ¼1F i P Gi ðÞþBX N C i ¼1w j Pen j ð13Þwhere Pen j is given by (12).Each particle updates itsposition in each iteration t according to the following equation:u t þ1id¼u t id þv t þ1id ð14Þwhere v id is the velocity of the i th particle in the d th dimension representing the change in position between two iterations.In this updating process each individual takes into consideration two kind of information;his own experience,represented by his current position and its best position so far (pbest ),and other individuals’experience,represented by the best position visited so far by the entire swarm (gbest ).Three terms are used for the modelling of the particle’s knowledge,the summation of which results in particle’s velocity vector:its current velocities,its distance from pbest and its distance from gbest.Hence,the particle’s velocity can be expressed as follows:v t þ1id¼w Áv t id þc 1Árand 1Ápbest tid Àu t id ÀÁþc 2Árand 2Ágbest td Àu t id ÀÁð15Þwhere w inertia weight,which controls the impact of the previous velocity on the current one c iacceleration coefficientsrand i random numbers in the range [0,1]pbest t id pbest of particle i in dimension d until iteration t gbest t d gbest of the swarm in dimension d until iteration t u t idcurrent position of the i th particle in the d th dimensionNote that gbest is not updated in the end of each iteration but it is updated asynchronously [25];after the evaluation of the i th particle’s objective function,i th particle’s pbest is computed and gbest is updated and provided to (i +1)th particle to update its velocity using (15).4.3.2Dissipative PSO :The standard PSO exhibitssatisfactory performance in the first iterations and a nearly stagnated evolution in the next iterations of the algorithm.Therefore standard PSO does not possess the ability to improve the quality of the solutions,as the number of iterations increases [24].This is more evident in large,real-life optimisation problems.To prevent the system from being in a stationary state an opening dissipative system is constructed by introducing negative entropy through additional chaos in velocity and location of the particles [24].The DPSO exhibits better performance as compared to SPSO,due to its ability to achieve sustainable development,especially when the dimensionality of the problem increases.The condition that introduces chaos in velocity isIf ðrand i o c v Þthen v t þ1id ¼rand i Áv max d ;for i ¼1;2ð16Þwhereas in location it isIf ðrand i o c ‘Þthen u t þ1id ¼Random u min d ;u maxdÀÁð17Þwhere c v ,c ‘,are the chaotic factors in the range [0,1]for velocity and location respectively,‘Random’is a randomnumber between u min d u max d ,and v maxdis the maximum velocity in the d th dimension.Conditions (16)and (17)are checked and executed after (15)and (14),respectively,of the SPSO.In contrast to [24],in which only one of the conditions (16),(17)is applied in each case study,in this paper both conditions are applied in all cases,since several empirical tests have indicated better performance of this scheme in the OPF problem.4.3.3Discretisation logic :The discrete controlvariables are handled by using discrete numbers instead of continuous to express the current position and velocity in (14),(17)and (15),(16),respectively [19].The values of the velocity and position for these control variables are rounded to their closest discrete value.4.3.4Objective function computation :Thecomputation of each particle’s objective function is similar to the process followed by the EGA (Section 4.2).The objective function is computed as follows:(i)Acquire the values of the control variables u from theparticle.(ii)Solve the power-flow equations (6),using the FDLFdiscussed in Section 4.2,to compute the state vector x .(iii)Determine the violated functional constraints (7)andcompute associated penalty functions (12).(iv)Compute the objective function using (13).4.3.5Parameter selection :The values of theparameters have been selected after numerous tests andtaking into consideration the literature on PSO [18–20]and several empirical studies [26].The DPSO population consists of 100particles,all of which are randomly initialised,except from one,in which all control variables are initialised to their upper limit,as in the EGA.The values of both acceleration coefficients in (15),c 1,c 2,are taken equal to 2.The inertia weight w ,which is used to regulate the balance between global and local exploration,takes an iteration dependent value in the standard PSO process [18,20].However,in DPSO it has been shown [24]that a fixed value of 0.4provides better performance to the evolution process.The values of the chaotic coefficients,c v in (16)and c ‘in (17),are taken equal to 0.002and 0.001,respectively [24].The maximum velocity in the d th dimension in (16)is taken asv max d¼u max d Àu min dÀÁNð18Þwhere N ¼10[18].4.3.6Termination criterion :A predefined num-ber of iterations is used as a termination criterion for the DPSO algorithm.The value of the termination criterion is increased in accordance with the size of the problem search space.5Test results5.1Test systemsSeveral test cases are considered to evaluate the presented OPF methods.Several systems with increasing number of areas,derived by replicating the IEEE one-area RTS96[27]and the IEEE 118-bus system [28],are tested by all methods.The unit-cost data for the IEEE RTS96are derived from the heat-rate data provided in [27](fitted by quadratic functions,shown in columns 6–8of Table 1)and the fuelTable 1:Unit data for IEEE 118-bus systemUnit groupFuel/unit typeFuel cost [h /Gcal]P min G [MW]P maxG[MW]Heat-rate coefficients Generator bus numbersa [Gcal/h]b [Gcal/MWh]c [Gcal/MW 2h]U12oil/steam 20.874 2.4012.00 4.2387432.18752110.03987821,4,6,8,15,18,19,24,27,32,34,36,40,42,55,56,62,70,72,73,74,76,77,85,90,91,92,99,104,105,107,110,112,113,116U20oil/CT 46.04615.8020.0054.226666À1.81488740.137875531,46,87U50hydro 0.00015.0050.000.0000000.00000000.000000054,103,111U76coal/steam 6.66015.2076.0033.798806 1.98198230.007796912U100oil/steam 20.87425.00100.0032.697565 1.89244290.0030066–U155coal/steam 6.66054.25155.0046.463859 1.90170100.001402359,61U197oil/steam 20.87468.95197.0047.783680 1.92958620.001253925,49,100U350coal/steam 6.660140.00350.0081.532894 1.87312290.000822326,65,66U400nuclear 6.549100.00400.00102.538242 2.16416890.000247810,80U550coal/steam6.660165.00550.00337.4085001.75234650.001115469,89cost data listed in the third column of Table1.The value of water is zero assuming excessive inflows.The lower and upper bounds of the bus voltage magnitudes are0.9and 1.1,respectively.The transformer tap settings take17 discrete values in the range[0.9,1.1]with a step equal to 0.0125and the bus shunt admittances take four discrete values in the range[À150,0]MVAr with a step equal to50MVAr.Unit data for the IEEE118-bus system(not provided in [28])are given in Table1.All but the last unit types of the IEEE118-bus system are identical with the ones of the IEEE RTS96.One new unit type,U550,is introduced. The last column of Table1shows the buses of the IEEE-118system at which units of each unit type are connected. The lower and upper bounds of the bus voltage magnitudes are0.9and1.1,respectively.The transformer tap settings take17discrete values in the range[0.9,1.1]with a step equal to0.0125and the bus shunt admittances takefive discrete values in the range[À40,40]MVAr with a step equal to20MVAr.The tie-line data of IEEE-RTS96and IEEE118-bus systems are shown in Table2.Thefirst digit in each bus number denotes area number,while the remaining digits denote the bus number inside each area. The characteristics of all test systems are shown in Table3.All cases were run on a1.53GHz AMD Athlon/ 256MB RAM PC.5.2Discretisation logic checkAn experiment is performed to determine the difference in the solution attained by the postprocessing discretisation logic discussed in Section3.2.3(suboptimal solution),with the solution provided by a MINLP solver.For comparison the LINDO API MINLP solver[22]is used,which employs the branch-and-bound method.The tests are performed using the IEEE3-area RTS96(case3).The discretisation logic used in this paper provides a feasible solution with system cost equal to212805.47h/h running for12seconds, while the MINLP solver provides a system cost equal to 212805.17h/h running for30.5min.The difference in system cost is negligible,whereas the run-times are vastly different.5.3ResultsTable4presents the results from the application of the mathematical programming methods to the test systems. In columns2–4the size of the OPF problem(1)–(4)is presented,including the number of unknown variables,the number of constraints and the number of nonzero elements of the sparse constraint matrix of the problem.Thefifth column reports the operation cost for each test system, which is the same for both mathematical programming methods.In columns6–9the execution times of the two implementations are presented,whereas columns10–13 report the number of solver iterations for the continuous relaxation and the discretisation step for both mathematical programming methods.The execution times reported in the sixth and eighth columns for the two methods include the postprocessing discretisation step(Section3.2.3).The modelling time reported in the seventh column includes the time to build the sparse constraint matrix of the problem(to be processed by the solver)and the time to compute the Jacobian matrix,for both the continuous relaxation and the discretisation process.The great difference in the run-times of these methods(sixth and ninth columns)is attributed to the LINDO API program structure,which provides all necessary information(Jaco-bian matrix,sparse matrix representation of the constraint matrix)to the solver ready to be processed,thus avoiding the additional modelling time of GAMS.Table5summarises the results from the application of the metaheuristic methods to the test systems.Twenty runs have been performed for each case examined.The second column reports the number of control variables,which are the same for both metaheuristic methods.The third column reports the average(over the20runs)number of iterations (NI)to arrive at a good-quality OPF solution with the DPSO.A good-quality OPF solution is one with objective function within0.1%of the objective obtained after allowing the DPSO to evolve for the maximum number of iterations in each case,which is shown in the fourthTable2:Test systems tie-line dataTie-lines r[p.u.]x[p.u.]b[p.u.]Rating[MVA]IEEE RTS96107–2030.0420.1610.044175 123–217,103–4210.0100.0740.155500113–215,202–4130.0100.0750.158500101–508,223–318,402–523,523–6180.0130.1040.218500121–325,221–6250.0120.0970.203500IEEE118-bus systems1090–20400.026100.07030.01844500 1105–20560.039060.18130.04610500 1027–3042,1040–3085,1085–4105,2027–5105,1012–6090,4012–61050.053000.18300.047205002085–3090, 3105–4085, 3027–4040, 4090–5085, 1042–60850.010500.02880.00760500Table3:Test systems characteristicsCaseno.Test system Buses Units Lines Trans-formersShunts1IEEE1-area RTS96243338512IEEE2-area RTS964866791023IEEE3-area RTS9673991201634IEEE4-area RTS96971321602145IEEE5-area RTS961211652002656IEEE6-area RTS961461982413267IEEE1Â118-bus118541869148IEEE2Â118-bus23610837418289IEEE3Â118-bus354162563274210IEEE4Â118-bus472216752365611IEEE5Â118-bus590270940457012IEEE6Â118-bus70832411295484。