Electric Power Systems Research 84 (2012) 58–64Contents lists available at SciVerse ScienceDirectElectric Power SystemsResearchj o u r n a l h o m e p a g e :w w w.e l s e v i e r.c o m /l o c a t e /e p srDynamic equivalence by an optimal strategyJuan M.Ramirez a ,∗,Blanca V.Hernández a ,Rosa Elvira Correa b ,1a Centro de Investigación y de Estudios Avanzados del I.P.N.,Av del Bosque 1145.Col El Bajio,Zapopan,Jal.,45019,MexicobUniversidad Nacional de Colombia –Sede Medellín.Facultad de Minas.Carrera 80#65-223Bloque M8,114,Medellín,Colombiaa r t i c l ei n f oArticle history:Received 21May 2011Received in revised form 26September 2011Accepted 29September 2011Keywords:Power system dynamics Equivalent circuitPhasor measurement unit Stability analysisPower system stabilitya b s t r a c tDue to the curse of dimensionality,dynamic equivalence remains a computational tool that helps to analyze large amount of power systems’information.In this paper,a robust dynamic equivalence is proposed to reduce the computational burden and time consuming that the transient stability studies of large power systems represent.The technique is based on a multi-objective optimal formulation solved by a genetic algorithm.A simplification of the Mexican interconnected power system is tested.An index is used to assess the proximity between simulations carried out using the full and the reduced model.Likewise,it is assumed the use of information stemming from power measurements units (PMUs),which gives certainty to such information,and gives rise to better estimates.© 2011 Elsevier B.V. All rights reserved.1.IntroductionOne way to speed up the dynamic studies of currently inter-connected power systems without significant loss of accuracy is to reduce the size of the system model by means of dynamic equiva-lents.The dynamic equivalent is a simplified dynamic model used to replace an uninterested part,known as an external part,of a power system model.This replacement aims to reduce the dimension of the original model while the part of interest remains unchanged [1–6].The phrases “Internal system”(IS)and “external system”(ES)are used in this paper to describe the area in question,and the remaining regions,respectively.Boundary buses and tie lines can be defined in each IS or ES.It is usually intended to perform detailed studies in the IS.However,the ES is important to the extent where it affects IS analyses.The equivalent does not alter the transient behavior of the part of the system that is of concern and greatly reduces the dimen-sion of the network,reducing computational time and effort [4,7,8].The dynamic equivalent also can meet the accuracy in engineering,achieving effective,rapid and precise stability analysis and security controls for large-scale power system [4,8].However,the determi-∗Corresponding author.Tel.:+523337773600.E-mail addresses:jramirez@gdl.cinvestav.mx (J.M.Ramirez),bhernande@gdl.cinvestav.mx (B.V.Hernández),elvira.correa@ (R.E.Correa).1Tel.:+57314255140.nation of dynamic equivalents may also be a time consuming task,even if performed off-line.Moreover,several dynamic equivalents may be required to represent different operating conditions of the same system.Therefore,it is important to have computational tools that automate the procedure to evaluate the dynamic equivalent [7].Ordinarily,dynamic equivalents can be constructed follow-ing two distinct approaches:(i)reduction approach,and (ii)identification approach.The reduction approach is based on an elimination/aggregation of some components of the existing model [4,5,9].The two mostly found in the literature are known as modal reduction [6,10]and coherency based aggregation [2,11,12].The identification approach is based on either parametric or non-parametric identification [13,14].In this approach,the dynamic equivalent is determined from online measurements by adjusting an assumed model until its response matches with measurements.Concerning the capability of the model,the dynamic equivalent obtained from the reduction approach is considerably more reli-able and accurate than those set up by the identification approach,because it is determined from an exact model rather than an approximation based on measurements.However,the reduction-based equivalent requires a complete set of modeling data (e.g.model,parameters,and operating status)which is rarely avail-able in practice,in particular the generators’dynamic parameters [5,13,15,16].On the other hand,due to the lack of complete system data,and/or frequently variations of the parameters with time,the importance of estimation methods is revealed noticeably.Especially,on-line model correction aids for employing adaptive0378-7796/$–see front matter © 2011 Elsevier B.V. All rights reserved.doi:10.1016/j.epsr.2011.09.023J.M.Ramirez et al./Electric Power Systems Research84 (2012) 58–6459Fig.1.190-buses46-generators power system.controllers,power system stabilizers(PSS)or transient stability assessment.The capability of such methods has become serious rival of the old conventional methods(e.g.the coherency[11,12] and the modal[6,10]approaches).The equivalent estimation meth-ods have spread,because it can be estimated founded on data measured only on the boundary nodes between the study system and the external system.This way,without any need of informa-tion from the external system,estimation process tries to estimate a reduced order linear model,which is replaced for the external part.Evidently,estimation methods can be used,in presence of perfect data of the network as well to compute the equivalent by simulation and/or model order reduction[15].Sophisticated techniques have become interesting subject for researchers to solve identification problems since90s.For example, to obtain a dynamic equivalent of an external subsystem,an opti-mization problem has been solved by the Levenberg–Marquardt algorithm[17].Artificial neural networks(ANN)are the most prevalent method between these techniques because of its high inherent ability for modeling nonlinear systems,including power system dynamic equivalents[15,18–25].Power system real time control for security and stability has pro-moted the study of on-line dynamic equivalent technologies,which progresses in two directions.One is to improve the original off-line method.The mainstream approach is to obtain equivalent model based on typical operation modes and adjust equivalent parame-ters according to real time collected information[4].Distributed coordinate calculation based on real time information exchanging makes it possible to realize on-line equivalence of multi-area inter-connected power system in power market environment[4].Ourari et al.[26]developed the slow coherency theory based on the struc-ture preservation technique,and integrated dynamic equivalence into power system simulator Hypersim,verifying the feasibility of on-line computation from both computing time and accuracy[27].Prospects of phasor measurement technique based on global positioning system(GPS)applied in transient stability control of power system are introduced in Ref.[4].Using real data collected by phasor measurement unit(PMU),with the aid of GPS and high-speed communication network,online dynamic equivalent of interconnected power grid may be achieved[4].In this paper,the dynamic equivalence problem is formulated by two objective functions.An evolutionary optimization method based on genetic algorithms is used to solve the problem.2.PropositionThe main objective of this paper is the external system’s model order reduction of an electrical grid,preserving only the frontier nodes.That is,those nodes of the external system directly linked to nodes of the study system.At such frontier nodes,fictitious gen-erators are allocated.The external boundary is defined by the user. Basically,it is composed by a set of buses,which connect the exter-nal areas to the study system.There is not restriction about this set. Different operating conditions are taken into account.work reductionAfirst condition for an equivalence strategy is the steady state preservation on the reduced grid;this means basically a precise60J.M.Ramirez et al./Electric Power Systems Research 84 (2012) 58–64Fig.2.Proposed strategy’s flowchart.voltage calculation.In this paper,all nodes of the external system are eliminated,except the frontier nodes.By a load flow study,the complex power that should inject some fictitious generators at such nodes can be calculated.The nodal balance equation yields,j ∈Jp ij +Pg i +Pl i =0,∀i ∈I(1)where I is the set of frontier nodes;J is the set of nodes linked directly to the i th frontier node;p ij is the active power flowing from the i th to j th node;Pg i is the generation at the i th node;Pl i is the load at the i th node.Thus,the voltages for the reduced model become equal to those of the full one.For studies where unbalanced conditions are important,a similar procedure could be followed for the negative and zero equivalent sequences calculation.2.2.Studied systemThe power system shown in Fig.1depicts a reduced version of the Mexican interconnected power system.It encompasses 7regional systems,with a generation capacity of 54GW in 2004and an annual consumption level of 183.3TWh in 2005.The transmis-sion grid comprises a large 400/230kV system stretching from the southern border with Central America to its northern interconnec-tions with the US.The grids at the north and south of the country are long and sparsely interconnected transmission paths.The major load centers are concentrated on large metropolitan areas,mainly Mexico City in the central system,Guadalajara City in the western system,and Monterrey City in the northeastern system.The subsystem on the right of the dotted line is considered as the system under study.Thus,the subsystem on the left is the exter-nal one.There are five frontier nodes (86,140,142,148and 188)and six frontier lines (86–184,140–141,142–143,148–143(2)and 188–187).Thus,the equivalent electrical grid has five fictitious generators at nodes 86,140,142,148and 188.Transient stabil-ity models are employed for generators,equipped with a static excitation system;its formulation is described as follows,dıdt=ω−ω0(2)dωdt=1T j [Tm −Te −D (ω−ω0)](3)dE q dt=1T d 0 −E g −(x d −x d )i d +E fd(4)dE ddt =1T d 0−E d +(x q −xq )i q (5)dE fd dt=1T A−E fd +K A (V ref +V s −|V t |) (6)where ı(rad)and ω(rad/s)represent the rotor angular positionand angular velocity;E d (pu)and E q(pu)are the internal transient voltages of the synchronous generator;E fd (pu)is the excitationvoltage;i d (pu)and i q (pu)are the d -and q -axis currents;T d 0(s)and T q 0(s)are the d -and q -open-circuit transient time constants;x’d (pu)and x q(pu)are the d -and q -transient reactances;Tm (pu)and Te (pu)are the mechanical and electromagnetic nominal torque;Tj is the moment of inertia;D is the damping factor;K A and T A (s)are the system excitation gain and time constant;V ref is the voltage reference;V t is the terminal voltage;V s is the PSS’s output (if installed).The corresponding parameters are selected as typical [28].2.3.FormulationGiven some steady state operating point (#CASES )the following objective functions are defined,min f =[f 1f 2]f 1=#CASESop =1w 1opNg intk =1ωk ori (t )−ωkequiv (x,t )2(7)f 2=#CASESop =1w 2opNg intk =1Pe kori (t )−Pe k equiv (x,t )2(8)subject to:0=Ngen eqj =1H j −ni =1i ∈LH i(9)where ωk ori is the time behavior of the angular velocity of those generators in the original system,that will be preserved (Ng int),J.M.Ramirez et al./Electric Power Systems Research84 (2012) 58–6461Fig.3.Fitness assignment of NSGA-II in the two-objective space.Fig.4.Case1:from top to bottom(i)angular position37(referred to slack);(ii)angular speed28;(iii)electrical torque41,after a three-phase fault at bus172.62J.M.Ramirez et al./Electric Power Systems Research 84 (2012) 58–64Fig.5.Case 3:from top to bottom (i)angular position 40(referred to slack);(ii)angular speed 34;(iii)electrical torque 39,after a three-phase fault at bus 144.after a disturbance within the internal area;ωk equiv is the time behavior of the angular velocity of those generators in the equiv-alent system,after the same disturbance within the internal area;Pe k ori is the time behavior of the electrical power of those gen-erators of the original system within the internal area;Pe k equiv is the time behavior of the electrical power of those generators in the equivalent system within the internal area;H k is the k th genera-tor’s inertia;L is the set of generators that belong to the external system;Ngen eq is the number of equivalent generators [16,25].The set of voltages S ={V i ,V j ,...,V k |complex voltages stemming from PMUs }has been included in the solution.The main challenge in a multi-objective optimization environ-ment is to minimize the distance of the generated solutions to the Pareto set and to maximize the diversity of the developed Pareto set.A good Pareto set may be obtained by appropriate guiding of the search process through careful design of reproduction oper-ators and fitness assignment strategies.To obtain diversification special care has to be taken in the selection process.Special care is also to be taken to prevent non-dominated solutions from being lost.Elitism addresses the problem of losing good solutions dur-ing the optimization process.In this paper,the NSGA-II algorithm (Non dominated Sorting Genetic Algorithm-II)[29–31]is used to solve the formulation.The algorithm NSGA-II has demonstrated to exhibit a well performance;it is reliable and easy to handle.It uses elitism and a crowded comparison operator that keeps diver-sity without specifying any additional parameters.Pragmatically,it is also an efficient algorithm that has shown better results to solve optimization problems with multi-objective functions in a series of benchmark problems [31,32].There are some other meth-ods that may be used.For instance,it is possible to use at least two population-based non-Pareto evolutionary algorithms (EA)and two Pareto-based EAs:the Vector Evaluated Genetic Algorithm (VEGA)[36],an EA incorporating weighted-sum aggregation [33],the Niched Pareto Genetic Algorithm [34,35],and the Nondom-inated Sorting Genetic Algorithm (NSGA)[29–31];all but VEGA use fitness sharing to maintain a population distributed along the Pareto-optimal front.3.ResultsIn this case,the decision variables,x ,are eight parametersper each equivalent generator:{x d ,x d ,x q ,x q ,T d 0,T q 0,H,D }.In this paper,for five equivalent generators,there are 40parameters to be estimated.Likewise,in this case,a random change in the load of all buses gives rise to the transient behavior.A normal distribution with zero mean is utilized to generate the increment (decrement)in all buses.The variation is limited to a maximum of 50%.The disturbance lasts for 0.12s and then it is eliminated;the studied time is 2.0s.To attain more precise equivalence for severe operating conditions,this improvement could require load variations greater than 50%.However,this bound was used in all cases.Fig.2depicts a flowchart of the followed strategy to calculate an optimal solution.In this paper,three operating points are taken into account:(i)Case 1,the nominal case [37];(ii)Case 2,an increment of 40%in load and generation;(iii)Case 3,a decrement of 30%in load and gen-eration.To account for each operating condition into the objective functions,the same weighted factors have been utilized (w i =1/3),Eqs.(7)–(8).Table 1summarizes the estimated parameters for five equiv-alent generators,according to the two objective functions.TheJ.M.Ramirez et al./Electric Power Systems Research84 (2012) 58–6463 Table1Parameters of the equivalent under a maximum of50%in load variation.G equiv1G equiv2G equiv3G equiv4G equiv5f1f2f1f2f1f2f1f2f1f2x d0.1080.104 2.140 2.100 1.920 1.9100.1590.1790.3240.362x d0.1130.1130.8190.8690.3960.3960.08760.132 1.900 1.950T d010.7010.7039.6039.6018.8018.7011.5011.5011.7011.70x q0.7100.7500.5300.5190.2880.308 2.470 2.4600.7050.803x q0.3950.3850.8740.8050.9010.9260.9220.9530.8820.884T q0 4.950 4.82033.3033.40 4.280 4.10016.9016.9012.8012.90H 4.610 4.61022.2722.2747.2647.2669.2369.2333.9333.93D18.0318.40535.6536.6239.2239.141.9042.00705.1705.2 electromechanical modes associated to generators of the internalsystem are closely preserved.These generators arefictitious andbasically are useful to preserve some of the main interarea modesbetween the internal area and the external one[16].Thus,in order to avoid the identification of the equivalent gener-ators’parameters based on a specific disturbance,in this paper theuse of random changes in all the load buses is used.This will giverise to parameters valid for different fault locations.The allowedchange in the load(in this paper,50%)will result in a slight varia-tion of transient reactances.Further studies are required to assesssensitivities.Fig.3shows a typical Pareto front for this application.The NSGA-II runs on a Matlab platform and the convergence lastsfor3.25h for a population of200individuals and20generations.It is assumed that phasor measurement units(PMUs)areinstalled at specific buses(188,140,142,148,and86in the externalsystem,and141,143,145,and182in the internal system),whichbasically correspond to the frontier nodes.Likewise,it is assumedthat the precise voltages are known in these buses every time.Bythe inclusion of the PMUs,there is a noticeable improvement in thevoltages’information at the buses near them,due to the fact that itis assumed the PMUs’high precision.In this paper,the simulation results obtained by the full andthe reduced system are compared by a closeness measure,themean squared error(MSE).The goal of a signalfidelity measureis to compare two signals by providing a quantitative score thatdescribes the degree of similarity/fidelity or,conversely,the levelof error/distortion between ually,it is assumed that oneof the signals is a pristine original,while the other is distorted orcontaminated by errors[38].Suppose that z={z i|i=1,2,...,N}and y={y i|i=1,2,...,N}aretwofinite-length,discrete signals,where N is the number of signalsamples and z i and y i are the values of the i th samples in z and y,respectively.The MSE is defined by,MSE(z,y)=1NNi=1(z i−y i)2(10)Figs.4–5illustrate the transient behavior of some representa-tive signals after a three-phase fault at buses172(Fig.4)for theCase1;and144(Fig.5)for the Case3.Bus39is selected as theslack bus.Values in Table2show the corresponding MSE valuesfor the twelve generators of the internal system for each operat-ing case.Such values indicate a close relationship between the fulland reduced signals’behavior.In order to improve the equivalenceof a specific operating point,it is possible to weight it differentlythrough the factors w1–3,Eqs.(7)–(8).Table2shows the MSE’s val-ues when Case2has a higher weighting than Case1and Case3(w2=2/3,w1=w3=1/6).These values indicate that a closer agree-ment is attained between signals with the full and the reducedmodel.It is emphasized that the equivalence’s improvement couldrequire load variations greater than50%for off-nominal operatingconditions.4.ConclusionsUndoubtedly,the power system equivalents’calculationremains a useful strategy to handle the large amount of data,cal-culations,information and time,which represent the transientstability studies of modern power grids.The proposed approachis founded on a multi-objective formulation,solved by a geneticalgorithm,where the objective functions weight independentlyeach operating condition taken into account.The use of informationstemming from PMUs helps to improve the estimated equivalentgenerators’parameters.Results indicate that the strategy is ableto closely preserve the oscillating modes associated to the inter-nal system’s generators,under different operating conditions.Thatis due to the preservation of the machines’inertia.The use ofan index to measure the proximity between the signal’s behav-ior after a three-phase fault,indicates that good agreement isattained.In this paper,the same weighting factors have been usedto assess different operating conditions into the objective func-tions.However,depending on requirements,these factors can bemodified.The equivalence based on an optimal formulation assuresTable2MSE for Case2(Three-Phase Fault at Bus168).Angular position Angular speed Electrical powerw k=1/3w2=2/3,w1=w3=1/6w k=1/3w2=2/3,w1=w3=1/6w k=1/3w2=2/3,w1=w3=1/6 Gen39 2.13E−01 1.54E−01 6.36E−05 6.90E−059.40E−02 6.91E−02Gen28 1.79E−019.39E−02 3.43E−05 2.75E−05 4.09E−05 2.50E−05Gen29 1.12E−01 4.96E−02 4.84E−05 2.29E−059.25E−04 5.38E−04Gen30 1.01E−01 4.27E−02 3.31E−05 1.62E−059.75E−04 4.52E−04Gen32 2.80E−01 1.02E−01 4.34E−05 2.69E−05 4.26E−04 5.16E−04Gen33 2.97E−01 1.09E−01 4.81E−05 3.29E−059.85E−04 1.15E−04Gen34 1.70E−01 1.04E−01 4.03E−05 4.09E−059.19E−03 6.17E−03Gen37 1.68E−01 1.08E−01 5.23E−05 5.39E−05 1.12E−028.48E−03Gen38 1.53E−019.54E−02 4.47E−05 4.60E−05 1.29E−02 1.01E−02Gen40 1.91E−01 1.15E−01 5.27E−05 5.14E−05 2.26E−03 1.06E−03Gen41 1.93E−01 1.18E−01 5.52E−05 5.39E−05 2.50E−04 1.16E−04Gen42 1.68E−01 1.07E−01 4.17E−05 4.28E−05 5.23E−05 3.02E−0564J.M.Ramirez et al./Electric Power Systems Research84 (2012) 58–64proximity between the full and the reduced models.Closer prox-imity is reached if more stringent convergence’s parameters are defined,as well as additional objective functions,as line’s power flows,are included.References[1]P.Nagendra,S.H.nee Dey,S.Paul,An innovative technique to evaluate networkequivalent for voltage stability assessment in a widespread sub-grid system, Electr.Power Energy Syst.(33)(2011)737–744.[2]A.M.Miah,Study of a coherency-based simple dynamic equivalent for transientstability assessment,IET Gener.Transm.Distrib.5(4)(2011)405–416.[3]E.J.S.Pires de Souza,Stable equivalent models with minimum phase transferfunctions in the dynamic aggregation of voltage regulators,Electr.Power Syst.Res.81(2011)599–607.[4]Z.Y.Duan Yao,Z.Buhan,L.Junfang,W.Kai,Study of coherency-based dynamicequivalent method for power system with HVDC,in:Proceedings of the2010 International Conference on Electrical and Control Engineering,2010.[5]T.Singhavilai,O.Anaya-Lara,K.L.Lo,Identification of the dynamic equiva-lent of a power system,in:44th International Universities’Power Engineering Conference,2009.[6]J.M.Undrill,A.E.Turner,Construction of power system electromechanicalequivalents by modal analysis,IEEE Trans.PAS90(1971)2049–2059.[7]A.B.Almeida,R.R.Rui,J.G.C.da Silva,A software tool for the determinationof dynamic equivalents of power systems,in:Symposium-Bulk Power System Dynamics and Control–VIII(IREP),2010.[8]A.Akhavein,M.F.Firuzabad,R.Billinton,D.Farokhzad,Review of reductiontechniques in the determination of composite system adequacy equivalents, Electr.Power Syst.Res.80(2010)1385–1393.[9]U.D.Annakkage,N.C.Nair,A.M.Gole,V.Dinavahi,T.Noda,G.Hassan,A.Monti,Dynamic system equivalents:a survey of available techniques,in:IEEE PES General Meeting,2009.[10]B.Marinescu,B.Mallem,L.Rouco,Large-scale power system dynamic equiva-lents based on standard and border synchrony,IEEE Trans.Power Syst.25(4) (2010)1873–1882.[11]R.Podmore,A.Germond,Development of dynamic equivalents for transientstability studies,Final Report on EPRI Project RP763(1977).[12]E.J.S.Pires de Souza,Identification of coherent generators considering the elec-trical proximity for drastic dynamic equivalents,Electr.Power Syst.Res.78 (2008)1169–1174.[13]P.Ju,L.Q.Ni,F.Wu,Dynamic equivalents of power systems with onlinemeasurements.Part1:theory,IEE Proc.Gener.Transm.Distrib.151(2004) 175–178.[14]W.W.Price,D.N.Ewart,E.M.Gulachenski,R.F.Silva,Dynamic equivalents fromon-line measurements,IEEE Trans.PAS94(1975)1349–1357.[15]G.H.Shakouri,R.R.Hamid,Identification of a continuous time nonlinear statespace model for the external power system dynamic equivalent by neural net-works,Electr.Power Energy Syst.31(2009)334–344.[16]J.M.Ramirez,Obtaining dynamic equivalents through the minimizationof a lineflows function,Int.J.Electr.Power Energy Syst.21(1999) 365–373.[17]J.M.Ramirez,R.Garcia-Valle,An optimal power system model order reductiontechnique,Electr.Power Energy Syst.26(2004)493–500.[18]E.de Tuglie,L.Guida,F.Torelli,D.Lucarella,M.Pozzi,G.Vimercati,Identificationof dynamic voltage–current power system equivalents through artificial neural networks,Bulk power system dynamics and control–VI,2004.[19]A.M.Azmy,I.Erlich,Identification of dynamic equivalents for distributionpower networks using recurrent ANNs,in:IEEE Power System Conference and Exposition,2004,pp.348–353.[20]P.Sowa,A.M.Azmy,I.Erlich,Dynamic equivalents for calculation of powersystem restoration,in:APE’04Wisla,2004,pp.7–9.[21]A.H.M.A.Rahim,A.J.Al-Ramadhan,Dynamic equivalent of external power sys-tem and its parameter estimation through artificial neural networks,Electr.Power Energy Syst.24(2002)113–120.[22]A.M.Stankovic, A.T.Saric,osevic,Identification of nonparametricdynamic power system equivalents with artificial neural networks,IEEE Trans.Power Syst.18(4)(2003)1478–1486.[23]A.M.Stankovic,A.T.Saric,Transient power system analysis with measurement-based gray box and hybrid dynamic equivalents,IEEE Trans.Power Syst.19(1) (2004)455–462.[24]O.Yucra Lino,Robust recurrent neural network-based dynamic equivalencingin power system,in:IEEE Power System Conference and Exposition,2004,pp.1067–1068.[25]J.M.Ramirez,V.Benitez,Dynamic equivalents by RHONN,Electr.Power Com-pon.Syst.35(2007)377–391.[26]M.L.Ourari,L.A.Dessaint,V.Q.Do,Dynamic equivalent modeling of large powersystems using structure preservation technique,IEEE Trans.Power Syst.21(3) (2006)1284–1295.[27]M.L.Ourari,L.A.Dessaint,V.Q.Do,Integration of dynamic equivalents in hyper-sim power system simulator,in:IEEE PES General Meeting,2007,pp.1–6. [28]P.M.Anderson,A.A.Fouad,Power System Control and Stability.Appendix D,The Iowa State University Press,1977.[29]K.Deb,Evolutionary algorithms for multicriterion optimization in engineer-ing design,in:Evolutionary Algorithms in Engineering and Computer Science, 1999,pp.135–161(Chapter8).[30]N.Srinivas,K.Deb,Multiple objective optimization using nondominated sort-ing in genetic algorithms,put.2(3)(1994)221–248.[31]K.Deb,A.Pratap,S.Agarwal,T.Meyarivan,A fast and elitist multi-objectivegenetic algorithm:NSGA-II,IEEE put.6(2)(2002)182–197.[32]K.Deb,Multi-Objective Optimization using Evolutionary Algorithms,1st ed.,John Wiley&Sons(ASIA),Pte Ltd.,Singapore,2001.[33]P.Hajela,C.-Y.Lin,Genetic search strategies in multicriterion optimal design,Struct.Optim.4(1992)99–107.[34]Jeffrey Horn,N.Nafpliotis,Multiobjective optimization using the niched paretogenetic algorithm,IlliGAL Report93005,Illinois Genetic Algorithms Laboratory, University of Illinois,Urbana,Champaign,July,1993.[35]J.Horn,N.Nafpliotis,D.E.Goldberg,A niched pareto genetic algorithm formultiobjective optimization,in:Proceedings of the First IEEE Conference on Evolutionary Computation,IEEE World Congress on Computational Computa-tion,vol.1,82–87,Piscataway,NJ,IEEE Service Center,1994.[36]J.D.Schafer,Multiple objective optimization with vector evaluated geneticalgorithms,in:J.J.Grefenstette(Ed.),Proceedings of an International Confer-ence on Genetic Algorithms and their Applications,1985,pp.93–100.[37]A.R.Messina,J.M.Ramirez,J.M.Canedo,An investigation on the use of powersystem stabilizers for damping inter-area oscillations in longitudinal power systems,IEEE Trans.Power Syst.13(2)(1998)552–559.[38]W.Zhou,A.C.Bovik,Mean squared error:love it or leave it?A new look at signalfidelity measures,IEEE Signal Process.Mag.26(1)(2009)98–117.。
a r X i v :h e p -p h /9711228v 1 4 N o v 1997hep-ph/9711228October 1997O (α)QED Corrections to Polarized Elastic µe and Deep Inelastic lN ScatteringDima Bardin a,b,c ,Johannes Bl¨u mlein a ,Penka Christova a,d ,and Lida Kalinovskaya a,caDESY–Zeuthen,Platanenallee 6,D–15735Zeuthen,GermanybINFN,Sezione di Torino,Torino,ItalycJINR,ul.Joliot-Curie 6,RU–141980Dubna,RussiadBishop Konstantin Preslavsky University of Shoumen,9700Shoumen,BulgariaAbstractTwo computer codes relevant for the description of deep inelastic scattering offpolarized targets are discussed.The code µe la deals with radiative corrections to elastic µe scattering,one method applied for muon beam polarimetry.The code HECTOR allows to calculate both the radiative corrections for unpolarized and polarized deep inelastic scattering,including higher order QED corrections.1IntroductionThe exact knowledge of QED,QCD,and electroweak (EW)radiative corrections (RC)to the deep inelastic scattering (DIS)processes is necessary for a precise determination of the nucleon structure functions.The present and forthcoming high statistics measurements of polarized structure functions in the SLAC experiments,by HERMES,and later by COMPASS require the knowledge of the RC to the DIS polarized cross-sections at the percent level.Several codes based on different approaches for the calculation of the RC to DIS experiments,mainly for non-polarized DIS,were developped and thoroughly compared in the past,cf.[1].Later on the radiative corrections for a vast amount of experimentally relevant sets of kinematic variables were calculated [2],including also semi-inclusive situations as the RC’s in the case of tagged photons [3].Furthermore the radiative corrections to elastic µ-e scattering,a process to monitor (polarized)muon beams,were calculated [4].The corresponding codes are :•HECTOR 1.00,(1994-1995)[5],by the Dubna-Zeuthen Group.It calculates QED,QCD and EW corrections for variety of measuremets for unpolarized DIS.•µe la 1.00,(March 1996)[4],calculates O (α)QED correction for polarized µe elastic scattering.•HECTOR1.11,(1996)extends HECTOR1.00including the radiative corrections for polarized DIS[6],and for DIS with tagged photons[3].The beta-version of the code is available from http://www.ifh.de/.2The Programµe laMuon beams may be monitored using the processes ofµdecay andµe scattering in case of atomic targets.Both processes were used by the SMC experiment.Similar techniques will be used by the COMPASS experiment.For the cross section measurement the radiative corrections to these processes have to be known at high precision.For this purpose a renewed calculation of the radiative corrections toσ(µe→µe)was performed[4].The differential cross-section of polarized elasticµe scattering in the Born approximation reads,cf.[7],dσBORNm e Eµ (Y−y)2(1−P e Pµ) ,(1)where y=yµ=1−E′µ/Eµ=E′e/Eµ=y e,Y=(1+mµ/2/Eµ)−1=y max,mµ,m e–muon and electron masses,Eµ,E′µ,E′e the energies of the incoming and outgoing muon,and outgoing electron respectively,in the laboratory frame.Pµand P e denote the longitudinal polarizations of muon beam and electron target.At Born level yµand y e agree.However,both quantities are different under inclusion of radiative corrections due to bremsstrahlung.The correction factors may be rather different depending on which variables(yµor y e)are used.In the SMC analysis the yµ-distribution was used to measure the electron spin-flip asymmetry A expµe.Since previous calculations,[8,9],referred to y e,and only ref.[9]took polarizations into account,a new calculation was performed,including the complete O(α)QED correction for the yµ-distribution,longitudinal polarizations for both leptons,theµ-mass effects,and neglecting m e wherever possible.Furthermore the present calculation allows for cuts on the electron re-coil energy(35GeV),the energy balance(40GeV),and angular cuts for both outgoing leptons (1mrad).The default values are given in parentheses.Up to order O(α3),14Feynman graphs contribute to the cross-section forµ-e scattering, which may be subdivided into12=2×6pieces,which are separately gauge invariantdσQEDdyµ.(2) One may express(2)also asdσQEDdyµ+P e Pµdσpol kk=1−Born cross-section,k=b;2−RC for the muonic current:vertex+bremsstrahlung,k=µµ;3−amm contribution from muonic current,k=amm;4−RC for the electronic current:vertex+bremsstrahlung,k=ee;5−µe interference:two-photon exchange+muon-electron bremsstrahlung interference,k=µe;6−vacuum polarization correction,runningα,k=vp.The FORTRAN code for the scattering cross section(2)µe la was used in a recent analysis of the SMC collaboration.The RC,δA yµ,to the asymmetry A QEDµeshown infigures1and2is defined asδA yµ=A QEDµedσunpol.(4)The results may be summarized as follows.The O(α)QED RC to polarized elasticµe scattering were calculated for thefirst time using the variable yµ.A rather general FORTRAN codeµe la for this process was created allowing for the inclusion of kinematic cuts.Since under the conditions of the SMC experiment the corrections turn out to be small our calculation justifies their neglection. 3Program HECTOR3.1Different approaches to RC for DISThe radiative corrections to deep inelastic scattering are treated using two basic approaches. One possibility consists in generating events on the basis of matrix elements including the RC’s. This approach is suited for detector simulations,but requests a very hughe number of events to obtain the corrections at a high precision.Alternatively,semi-analytic codes allow a fast and very precise evaluation,even including a series of basic cuts andflexible adjustment to specific phase space requirements,which may be caused by the way kinematic variables are experimentally measured,cf.[2,5].Recently,a third approach,the so-called deterministic approach,was followed,cf.[10].It treats the RC’s completely exclusively combining features of fast computing with the possibility to apply any cuts.Some elements of this approach were used inµe la and in the branch of HECTOR1.11,in which DIS with tagged photons is calculated.Concerning the theoretical treatment three approaches are in use to calculate the radiative corrections:1)the model-independent approach(MI);2)the leading-log approximation(LLA); and3)an approach based on the quark-parton model(QPM)in evaluating the radiative correc-tions to the scattering cross-section.In the model-independent approach the QED corrections are only evaluated for the leptonic tensor.Strictly it applies only for neutral current processes.The hadronic tensor can be dealt with in its most general form on the Lorentz-level.Both lepton-hadron corrections as well as pure hadronic corrections are neglected.This is justified in a series of cases in which these corrections turn out to be very small.The leading logarithmic approximation is one of the semi-analytic treatments in which the different collinear singularities of O((αln(Q2/m2l))n)are evaluated and other corrections are neglected.The QPM-approach deals with the full set of diagrams on the quark level.Within this method,any corrections(lepton-hadron interference, EW)can be included.However,it has limited precision too,now due to use of QPM-model itself. Details on the realization of these approaches within the code HECTOR are given in ref.[5,11].3.2O (α)QED Corrections for Polarized Deep Inelastic ScatteringTo introduce basic notation,we show the Born diagramr rr r j r r r r l ∓( k 1,m )l ∓( k 2,m )X ( p ′,M h )p ( p ,M )γ,Z ¨¨¨¨B ¨¨¨¨£¢ ¡£¢ ¡£¢ ¡£¢ ¡£¢ ¡£¢ ¡£¢ ¡£¢ ¡£¢ ¡£¢ ¡z r r r r r r r r r r r r r rr ¨¨¨¨B ¨¨¨¨r r r r j r r r r and the Born cross-section,which is presented as the product of the leptonic and hadronic tensordσBorn =2πα2p.k 1,x =Q 2q 2F 1(x,Q 2)+p µ p ν2p.qF 3(x,Q 2)+ie µνλσq λs σ(p.q )2G 2(x,Q 2)+p µ s ν+ s µ p νp.q1(p.q )2G 4(x,Q 2)+−g µν+q µq νp.qG 5(x,Q 2),(8)wherep µ=p µ−p.qq 2q µ,and s is the four vector of nucleon polarization,which is given by s =λp M (0, n )in the nucleonrest frame.The combined structure functions in eq.(8)F1,2(x,Q2)=Q2e Fγγ1,2(x,Q2)+2|Q e|(v l−p eλl a l)χ(Q2)FγZ1,2(x,Q2)+ v2l+a2l−2p eλl v l a l χ2(Q2)F ZZ1,2(x,Q2),F3(x,Q2)=2|Q e|(p e a l−λl v l)χ(Q2)FγZ3(x,Q2),+ 2p e v l a l−λl v2l+a2l χ2(Q2)F ZZ3(x,Q2),G1,2(x,Q2)=−Q2eλl gγγ1,2(x,Q2)+2|Q e|(p e a l−λl v l)χ(Q2)gγZ1,2(x,Q2),+ 2p e v l a l−λl v2l+a2l χ2(Q2)g ZZ1,2(x,Q2),G3,4,5(x,Q2)=2|Q e|(v l−p eλl a l)χ(Q2)gγZ3,4,5(x,Q2),+ v2l+a2l−2p eλl v l a l χ2(Q2)g ZZ3,4,5(x,Q2),(9) are expressed via the hadronic structure functions,the Z-boson-lepton couplings v l,a l,and the ratio of the propagators for the photon and Z-bosonχ(Q2)=Gµ2M2ZQ2+M2Z.(10)Furthermore we use the parameter p e for which p e=1for a scattered lepton and p e=−1for a scattered antilepton.The hadronic structure functions can be expressed in terms of parton densities accounting for the twist-2contributions only,see[12].Here,a series of relations between the different structure functions are used in leading order QCD.The DIS cross-section on the Born-leveld2σBorndxdy +d2σpol Borndxdy =2πα2S ,S U3(y,Q2)=x 1−(1−y)2 ,(13) and the polarized partdσpol BornQ4λp N f p S5i=1S p gi(x,y)G i(x,Q2).(14)Here,S p gi(x,y)are functions,similar to(13),and may be found in[6].Furthermore we used the abbrevationsf L=1, n L=λp N k 12πSy 1−y−M2xy2π1−yThe O(α)DIS cross-section readsd2σQED,1πδVRd2σBorndx l dy l=d2σunpolQED,1dx l dy l.(16)All partial cross-sections have a form similar to the Born cross-section and are expressed in terms of kinematic functions and combinations of structure functions.In the O(α)approximation the measured cross-section,σrad,is define asd2σraddx l dy l +d2σQED,1dx l dy l+d2σpol radd2σBorn−1.(18)The radiative corrections calculated for leptonic variables grow towards high y and smaller values of x.Thefigures compare the results obtained in LLA,accounting for initial(i)andfinal state (f)radiation,as well as the Compton contribution(c2)with the result of the complete calculation of the leptonic corrections.In most of the phase space the LLA correction provides an excellent description,except of extreme kinematic ranges.A comparison of the radiative corrections for polarized deep inelastic scattering between the codes HECTOR and POLRAD[17]was carried out.It had to be performed under simplified conditions due to the restrictions of POLRAD.Corresponding results may be found in[11,13,14].3.3ConclusionsFor the evaluation of the QED radiative corrections to deep inelastic scattering of polarized targets two codes HECTOR and POLRAD exist.The code HECTOR allows a completely general study of the radiative corrections in the model independent approach in O(α)for neutral current reac-tions including Z-boson exchange.Furthermore,the LLA corrections are available in1st and2nd order,including soft-photon resummation and for charged current reactions.POLRAD contains a branch which may be used for some semi-inclusive DIS processes.The initial state radia-tive corrections(to2nd order in LLA+soft photon exponentiation)to these(and many more processes)can be calculated in detail with the code HECTOR,if the corresponding user-supplied routine USRBRN is used together with this package.This applies both for neutral and charged current processes as well as a large variety of different measurements of kinematic variables. Aside the leptonic corrections,which were studied in detail already,further investigations may concern QED corrections to the hadronic tensor as well as the interference terms. References[1]Proceedings of the Workshop on Physics at HERA,1991Hamburg(DESY,Hamburg,1992),W.Buchm¨u ller and G.Ingelman(eds.).[2]J.Bl¨u mlein,Z.Phys.C65(1995)293.[3]D.Bardin,L.Kalinovskaya and T.Riemann,DESY96–213,Z.Phys.C in print.[4]D.Bardin and L.Kalinovskaya,µe la,version1.00,March1996.The source code is availablefrom http://www.ifh.de/~bardin.[5]A.Arbuzov,D.Bardin,J.Bl¨u mlein,L.Kalinovskaya and T.Riemann,Comput.Phys.Commun.94(1996)128,hep-ph/9510410[6]D.Bardin,J.Bl¨u mlein,P.Christova and L.Kalinovskaya,DESY96–189,hep-ph/9612435,Nucl.Phys.B in print.[7]SMC collaboration,D.Adams et al.,Phys.Lett.B396(1997)338;Phys.Rev.D56(1997)5330,and references therein.[8]A.I.Nikischov,Sov.J.Exp.Theor.Phys.Lett.9(1960)757;P.van Nieuwenhuizen,Nucl.Phys.B28(1971)429;D.Bardin and N.Shumeiko,Nucl.Phys.B127(1977)242.[9]T.V.Kukhto,N.M.Shumeiko and S.I.Timoshin,J.Phys.G13(1987)725.[10]G.Passarino,mun.97(1996)261.[11]D.Bardin,J.Bl¨u mlein,P.Christova,L.Kalinovskaya,and T.Riemann,Acta Phys.PolonicaB28(1997)511.[12]J.Bl¨u mlein and N.Kochelev,Phys.Lett.B381(1996)296;Nucl.Phys.B498(1997)285.[13]D.Bardin,J.Bl¨u mlein,P.Christova and L.Kalinovskaya,Preprint DESY96–198,hep-ph/9609399,in:Proceedings of the Workshop‘Future Physics at HERA’,G.Ingelman,A.De Roeck,R.Klanner(eds.),Vol.1,p.13;hep-ph/9609399.[14]D.Bardin,Contribution to the Proceedings of the International Conference on High EnergyPhysics,Warsaw,August1996.[15]M.Gl¨u ck,E.Reya,M.Stratmann and W.Vogelsang,Phys.Rev.D53(1996)4775.[16]S.Wandzura and F.Wilczek,Phys.Lett.B72(1977)195.[17]I.Akushevich,A.Il’ichev,N.Shumeiko,A.Soroko and A.Tolkachev,hep-ph/9706516.-20-18-16-14-12-10-8-6-4-200.10.20.30.40.50.60.70.80.91elaFigure 1:The QED radiative corrections to asymmetry without experimental cuts.-1-0.8-0.6-0.4-0.200.20.40.60.810.10.20.30.40.50.60.70.80.91elaFigure 2:The QED radiative corrections to asymmetry with experimental cuts.-50-40-30-20-100102030405000.10.20.30.40.50.60.70.80.91HectorFigure 3:A comparison of complete and LLA RC’s in the kinematic regime of HERMES for neutral current longitudinally polarized DIS in leptonic variables.The polarized parton densities [15]are used.The structure function g 2is calculated using the Wandzura–Wilczek relation.c 2stands for the Compton contribution,see [6]for details.-20-100102030405000.10.20.30.40.50.60.70.80.91HectorFigure 4:The same as in fig.3,but for energies in the range of the SMC-experiment.-20-10010203040500.10.20.30.40.50.60.70.80.91HectorFigure 5:The same as in fig.4for x =10−3.-200-150-100-5005010015020000.10.20.30.40.50.60.70.80.91HectorFigure 6:A comparison of complete and LLA RC’s at HERA collider kinematic regime for neutral current deep inelastic scattering offa longitudinally polarized target measuring the kinematic variables at the leptonic vertex.。
